# Diffuse optics for glaciology

Markus Allgaier<EMAIL_ADDRESS>Department of Physics and Oregon Center for
Optical, Molecular, and Quantum Science, University of Oregon, Eugene, Oregon
97403, USA  Brian J. Smith  Department of Physics and Oregon Center for
Optical, Molecular, and Quantum Science, University of Oregon, Eugene, Oregon
97403, USA

###### Abstract

Optical probing of glaciers has the potential for tremendous impact on
environmental science. However, glacier ice is turbid, which prohibits the use
of most established optical measurements for determining a glacier’s interior
structure. Here, we propose a method for determining the depth, scattering and
absorption length based upon diffuse propagation of short optical pulses. Our
model allows us to extract several characteristics of the glacier. Performing
Monte Carlo simulations implementing Mie scattering and mixed boundary
conditions, we show that the proposed approach should be feasible with current
technology. The results suggest that optical properties and geometry of the
glacier can be extracted from realistic measurements, which could be
implemented with low cost and small footprint.

##  I Introduction

The optical properties of glacial ice encode information on paleoclimate [1,
2], ice crystal orientation and flow [3, 4, 5], and hold the key to precise
modeling of radiative transfer [6, 7, 8]. Optical remote sensing techniques
are commonly employed to monitor snow cover [9]. While lidar cannot penetrate
more than a few centimeters into snow [9], exact knowledge of optical
properties of snow and ice is nevertheless required to achieve high precision
[8]. The transport of light in glacier ice is dominated by scattering from
bubbles, grain boundaries and dust [10], while pure ice absorption is small in
comparison. In deep glacier ice, diffuse optics have been employed to
successfully measure scattering and absorption coefficients independently of
each other from the shape of laser pulses transmitted through a large volume
[11, 12, 13]. Understanding of diffuse optics is tantamount to perform such
measurements. Penetrating, non-invasive optical measurements on the surface
would simplify data collection and allow the usage of low-cost, portable off-
the-shelf components. In addition, similar techniques are known in medical
imaging as diffuse optical imaging [14], where appropriate boundary conditions
for diffuse optics have long been known [15]. The theory of diffuse transport
of light with boundaries holds the key to measuring not only local optical
properties, but also to extracting information about the volume geometry, such
as the thickness of a glacier. Ice thickness is typically measured using
radar-sounding[16]. Such techniques are successfully used on cold glaciers,
but only low-frequency systems, which are typically not commercially
available, are of use for temperate glaciers [16, 17]. The use of commercially
available optical components could provide similar performance at a much
reduced system size and cost.

While many optical properties of snow and ice are well studied [7, 6, 18, 19,
20, 8], precise data on the scattering and absorption properties of near-
surface ice is scarce. Scattering coefficients of glacier ice has been
measured in deep ice ($>1000$ m) at the South Pole [11, 12]. By extrapolating
data collected in the arctic region one can conclude that scattering near the
surface is too strong for penetrating measurements [10]. However, recent
measurements in the Greenland ablation zone revealed that this is not the case
in warmer regions, with scattering lengths on the order of $1$ m [8]. Recently
extracted ice cores from the Rhône glaciers also show a low density of bubbles
of small size [5], corroborating the order of magnitude of scattering observed
in Greenland. This prompts a re-evaluation of the feasibility of a penetrating
optical measurement technique.

In this paper we present a simple model for diffusive transport of light in a
shallow glacier and outline the very basic principles needed to describe
diffuse optical measurements of glacier ice. We describe the glacier as a
semi-infinite slab, where transport of light is governed by both absorption
and scattering. We evaluate appropriate boundary conditions for both an
idealized, totally reflective bottom termination and Fresnel-reflection at the
top ice-air interface, which is known in the context of diffuse optical
medical imaging [15]. With these boundary conditions, we show that information
about the geometry becomes easily accessible by constructing solutions to the
diffusion equation using the well-known method of images [21]. This serves as
an extension of the infinite-volume random walk model developed for AMANDA
[22]. In the second half of the paper we employ Monte Carlo simulations of a
random walk to verify the validity of the model and give an estimate of
measurement performance and limitations under realistic conditions. The
results show that scattering and absorption length are accessible through both
time-resolved measurements and simple intensity measurements on the surface,
while information about thickness and geometry in general are accessible with
time-resolved measurements only. We discuss challenges arising for a possible
experimental implementation, such as laser sources, detectors, realistic
glacier ice structure and bedrock reflectivity. While more research is needed
to understand all these details as well as test the robustness of the
presented model, these results do outline the potential of diffuse optical
measurements for glaciology.

##  II Theory

###  II.1 Optical properties of glacier ice

The optical properties of glacier ice are governed by both the properties of
pure ice, as well as the effect of trapped air bubbles and impurities. For
example, the absorption coefficient $\beta$ of pure ice, which attenuates
transmitted light after length $d$ by the factor $\mathrm{exp}(-\beta d)$, is
small in the visible spectrum, measurements give values between
$10^{-1}\mathrm{m}^{-1}$ in a laboratory setting to $10^{-3}\mathrm{m}^{-1}$
in the field [7]. While measuring the exact value is a challenge, it is
evident that pure ice is optically very clear, and absorption at wavelength
above 600nm seems unaffected by impurities [7, 8]. Transport of light is
instead strongly governed by scattering at bubbles, impurities, crystal
boundaries, and various other sources, with bubbles giving the dominant
contribution [10]. If the transmission through snow or ice is affected by both
scattering and absorption, transmitted intensity $I(d)$ at depth $d$ can be
described by the attenuation coefficient $\alpha=\sigma+\beta$:

$\frac{I(d)}{I_{0}}=e^{-\alpha d}=e^{-(\sigma+\beta)d},$ (1)

where $I_{0}$ is the input light intensity and $\sigma$ is the isotropic
scattering coefficient, which is much larger in snow than it is in ice. After
a distance of $1/\alpha$ in the ice, the initial intensity will have dropped
to $1/e$. The size and character of $\sigma$ varies strongly throughout the
literature. Attenuation, scattering and absorption coefficients have recently
been measured in the Greenland ablation zone [8]. Without snow cover in the
ablation zone of a glacier, the bare ice is accessible to optical measurements
and appears to be mostly uniform at depths of more than a few centimeters. The
scattering and absorption coefficients were found to be of order
$\sigma\approx 1\mathrm{m}^{-1}$ and $\beta\approx
10^{-2}...10^{-1}\mathrm{m}^{-1}$, respectively, corresponding to scattering
and absorption length (mean free path) of $1m$ and $10...100m$. These findings
suggest that diffusive transport of light can be quite efficient, especially
in the spectral region around 400 nm. This implies that a small portion of
light injected at the surface reaches deep into the volume, and some signature
of the volume depth may be encoded in light detected at the surface.

From the existing literature it is evident that absorption in the visible
spectrum depends strongly on the content of light absorbing impurities, while
the physical scattering coefficient is largely independent of wavelength.
Independent measurement of $\sigma$ and $\beta$ without the use of reference
values, i.e., not relying on the attenuation coefficient from Eq. 1, is
challenging. To the best of the authors’ knowledge this has only been achieved
in the framework of the random-flight model employed at AMANDA [22, 12], where
the distribution of photon arrival times between a source and detector
embedded deep in the ice was measured.

Exact knowledge of $\sigma$ and $\beta$ are tantamount to establish a useful
and realistic model of the transport of light in a glacier. If these
parameters can be extracted from measurement, such a model can in turn be used
to obtain information about the geometry of a glacier, namely its thickness.
These questions are similar to diffuse optical imaging in the context of life
science [14].

###  II.2 Random flights: From multi-scattering to diffusion

To derive a simple model of diffuse light transport in a glacier, we model the
glacier as a semi-infinite slab. The system is symmetric around the $z$-axis.
The volume is uniformly turbid, meaning the scattering coefficient is larger
than the absorption coefficient, and both, as well as the index of refraction
$n$, are assumed to be constant. The volume has a reflective boundary at the
bottom ($z=-Z$), and a semi-reflective boundary at the top ($z=0$). An
instantaneous (infinitely short pulsed) source of light is injected at
$z=\rho=0$, a detector is placed at the surface to collect light at the point
of interest $z=0,\rho$, where we are using cylindrical coordinates. We assume
that light occupies a narrow spectral range and spectral dependencies can
hence be neglected. The geometry is depicted in Figure 1.

Figure 1: Sketch of transport of light in glacier ice of thickness $Z$. Light
gets injected downward from the surface and eventually gets scattered into a
random direction by bubbles in the ice. Some light makes it back to the
surface and gets detected by a detector separated from the source by the
distance $\rho$.

After injection, photons travelling on average for
$l_{\mathrm{sca}}=1/\sigma$, where they are scattered into a random direction,
as prescribed by the probability density $p(\theta)$ or scattering angle
$\theta$. The photon will travel a random distance $r$ before it is scattered
again. Each photon is scattered many times before reaching a boundary or being
absorbed, performing a ”random flight”, while on average covering the distance
$l_{\mathrm{sca}}$ with every scattering step, but in a random direction.
Finding the photon after a certain number of scattering steps becomes a
statistical problem instead of an optical one that can be described simply by
propagation. The probability density $P(R)$ to find a photon at distance $R$
from the injection point after $N$ steps is described by the distribution [23]

$P(R)=\frac{1}{\left(\frac{2}{3}\pi N\langle
r^{2}\rangle\right)^{3/2}}\mathrm{exp}\left(-3R^{2}/2N\langle
r^{2}\rangle\right),$ (2)

where $\langle r^{2}\rangle$ described the mean square displacement in the
random flight. Here, we first establish the connection between the governing
length scale of the random flight, $\langle r^{2}\rangle$, and the transport
parameter of diffusion, the diffusion constant $D$. Using this parameter, we
rewrite the above probability density into the well-known solution to the
diffusion equation in an infinite volume. Again, it represents a probability
density to find a photon at distance $R$ from the source at time $t$ after
being emitted:

$P(R,t)=\frac{1}{\left(4\pi
Dt\right)^{3/2}}\mathrm{exp}\left(-3R^{2}/4Dt\right),$ (3)

with

$D=\frac{N\langle r^{2}\rangle}{6t}=\frac{c\langle
r^{2}\rangle}{6l_{\mathrm{sca}}},$ (4)

where we introduced the stepping rate $N/t$. The distance traveled is
$Nl_{\mathrm{sca}}=ct$ by design, with the speed of light in the medium, $c$,
which removes the apparent time-dependence of the diffusion constant. The
diffusion constant is independent of the actual distribution of scattering
lengths and depends only on a mean transport parameter $\langle r^{2}\rangle$.

Askebjer et al.[22] uses an additional dampening term to account for
absorption. This term introduces a weight to components at long times. With
the absorption term, we modify the previous density function. Since the
probability to find the photon at large times is exponentially decreasing, the
following density function is formally not a probability density function,
because it does not integrate to 1, but it does still have the same form and
information about where and when to find the photon. It reads

$P(R,t)=\frac{1}{\left(4\pi
Dt\right)^{3/2}}\mathrm{exp}\left(-\frac{R^{2}}{4Dt}-\beta t\right).$ (5)

This is a solution to the photon diffusion equation. Transitioning to
macroscopic quantities, from a probability density $P(R,t)$ for the location
of a single particle to a fluence rate $\phi(R,t)$, describing the isotropic
fluence rate in units of [${\rm{Wm}}^{-2}$], originating from an isotropic
point source of unit energy. The diffusion-kernel with absorption is a well-
known solution to the isotropic, infinite-volume case of the photon diffusion
equation

$\frac{\partial\phi(R,t)}{\partial
t}=D\nabla^{2}\partial\phi(R,t)-c\beta\phi(R,t)+cS(R,t),$ (6)

where $S(R,t)$ is the source term.

###  II.3 The role of anisotropy

Before we move on, it is worthwhile to inspect the diffusion constant $D$ and
the physical interpretation of the mean distance and mean square displacement.
For a ”truly random” flight, the scattering angle $\theta$ is uniformly
distributed. It is straightforward that mean distance and mean square
displacement are then

$\langle r\rangle=\frac{1}{N}\sum_{i=1}^{N}r_{i}=0,\quad\langle
r^{2}\rangle=\frac{1}{N}\left(\sum_{i=1}^{N}r_{i}\right)^{2}=l_{\mathrm{sca}}^{2}.$
(7)

This implies that the distribution is always centered around the source point,
since the average distance traveled is 0, while the mean square displacement,
which is a measure of the extent of the distribution, is the scattering length
squared. Scattering by bubbles in particular is not isotropic, but forward-
peaked, with the average cosine of the scattering angle being
$\langle\cos\theta\rangle=0.75$ [22]. This is commonly known as the anisotropy
factor first introduced in meteorology [24] to describe attenuation in terms
of an effective scattering constant

$\sigma_{\mathrm{eff}}=\sigma(1-\langle\cos\theta\rangle).$ (8)

This effective scattering coefficient can successfully describe attenuation in
the sense of Eq. 1, even when the distribution of scattering angles is not
uniform, i.e., when scattering is anisotropic. While attenuation may behave
the same for an effective scattering constant compared to an actual isotropic
scattering constant, mean distance and mean square displacement differ[22]:

$\langle r\rangle=l_{\mathrm{sca}},\quad\langle
r^{2}\rangle=2l_{\mathrm{sca}}^{2}$ (9)

This yields a factor of 2 difference for the diffusion constant between truly
isotropic scattering and effective isotropic scattering. Using the mean square
displacement for forward-peaked scattering and substituting into Eq. 4 along
with the effective scattering length we find that the diffusion constant for
the case of scattering on bubbles or other forward-peaked scattering processes
is

$D=\frac{c\langle
r^{2}\rangle}{6l_{\mathrm{eff}}}=\frac{2cl_{\mathrm{eff}}^{2}}{6l_{\mathrm{eff}}}=\frac{cl_{\mathrm{eff}}}{3}.$
(10)

###  II.4 Boundary condition

The volume we are employing to model the glacier has two boundaries. The top
boundary is semi-reflective, the bottom boundary is reflective. As mentioned
before, assuming unit reflectivity at the base is an oversimplification, but
it allows for a straightforward, intuitive solution. This bottom boundary can
easily be modeled using the method of images, which is well known in the
context of electromagnetism and conduction of heat, which is also described by
diffusion [21]. The method is based on the assumption that any sum of
solutions $\phi_{1},\phi_{2}$ to the diffusion equation is itself a solution
to the diffusion equation

$\frac{\partial\left(\phi_{1}(R,t)+\phi_{2}(R,t)\right)}{\partial
t}=D\nabla^{2}\left(\phi_{1}(R,t)+\phi_{2}(R,t)\right)-c\beta\left(\phi_{1}(R,t)+\phi_{2}(R,t)\right)+cS(R,t).$
(11)

This follows from the linearity of the diffusion differential equation. The
method of images allows the construction of a solution satisfying various
boundary conditions from a source term, and several mirror sources. For a
fully reflective barrier at $z=-Z$, the mirror source term is identical to the
source, but mirrored around the boundary, giving the solution with two source
terms, where one is mirrored around the boundary and hence placed at $z=-2Z$:

$\phi(R,t)=\frac{1}{\left(4\pi
Dt\right)^{3/2}}\mathrm{exp}\left(-\frac{\rho^{2}}{4Dt}-\beta
t\right)\left(\mathrm{exp}\left(-\frac{z^{2}}{4Dt}\right)+\mathrm{exp}\left(-\frac{(z+2Z)^{2}}{4Dt}\right)\right).$
(12)

We simplify this expression by defining a point-source term as

$\phi_{+}(z)=\frac{1}{\left(4\pi
Dt\right)^{3/2}}\mathrm{exp}\left(-\frac{\rho^{2}}{4Dt}-\beta
t\right)\mathrm{exp}\left(-\frac{z^{2}}{4Dt}\right),$ (13)

which allows us to rewrite the above equation

$\phi(R,t)=\phi_{+}(z)+\phi_{+}(z+2Z).$ (14)

Here, it is assumed that the source is isotropic, which isn’t the case for
injection of a laser beam from the surface. We can approximate isotropic
behavior by replacing the source (and mirror source) by effective sources,
shifted by one effective scattering length into the slab:

$\phi(R,t)=\phi_{+}(z+l_{\mathrm{eff}})+\phi_{+}(z+2Z-l_{\mathrm{eff}}).$ (15)

This equation would be sufficient if all light escaped the top boundary, in
which case it completely describes the time-dependent power density
distribution in the slab in the presence of a reflective boundary.
Unfortunately, behavior at the top boundary isn’t that simple. To introduce a
mixed boundary condition that describes the Fresnel reflections at the
boundary between two media with different refractive index, ice and air, we
follow the derivation by Haskell et al.[15]. To find a suitable expression, we
need to introduce more radiometric quantities, namely the radiance $L$ and
flux $\vec{j|}$. Radiance is a directional power density per unit area and
therefore describes the directional transport of optical power. In a diffuse
medium, the total radiance is composed of the isotropic fluence rate and a
small directional flux:

$L(R,\vec{s},t)=\frac{\phi(R,t)}{4\pi}+\frac{3}{4\pi}\vec{j}(R,t)\cdot\vec{s}.$
(16)

The general diffusion approximation which allows us to describe the transport
of light with such equations in the first place dictates that the radiance is
dominated by the fluence rate, and the flux component is small. All outward
radiation that is repelled by Fresnel reflection can be replaced by an inbound
irradiance, analog to the method of images:

$L_{in}=RL_{out}=\iint_{\vec{s}\cdot\vec{n}>0}R(\vec{s})L(\vec{s})\vec{s}\cdot\vec{n}d\Omega,$
(17)

with

$R(\theta)=\frac{1}{2}\left(\frac{n\cos\theta^{\prime}-\cos\theta}{n\cos\theta^{\prime}+\cos\theta}\right)^{2}+\frac{1}{2}\left(\frac{n\cos\theta-\cos\theta^{\prime}}{n\cos\theta+\cos\theta^{\prime}}\right)^{2},\quad\theta<\theta_{c},$
(18)

where $\theta_{c}$ is the critical angle above which total internal reflection
occurs, $\theta^{\prime}$ is the external refracted angle. Note that we are
averaging the Fresnel reflection for both linear polarizations, assuming that
polarization is random after many scattering steps.

Left and right side of the reflected radiance equation evaluate as follows:

$\frac{\phi(R,t)}{4}+\frac{j_{z}}{2}=R_{\phi}\frac{\phi(R,t)}{4}-R_{j}\frac{j_{z}}{2},$
(19)

where the integrals for fluence rate and flux are

$R_{\phi}=\int_{0}^{\theta_{c}}2\sin\theta\cos\theta R(\theta)d\theta,$ (20)
$R_{j}=\int_{0}^{\theta_{c}}3\sin\theta\cos^{2}\theta R(\theta)d\theta.$ (21)

With these definitions, we can write

$j_{z}=-\frac{1}{2}\frac{1-R_{\phi}}{1+R_{j}}\phi.$ (22)

Using Fick’s law, which links fluence rate and flux in the diffusive regime,
the boundary condition can finally be identified as

$-\frac{D}{c}\nabla\phi(R,t)=j=-\frac{1}{2}\frac{1-R_{\phi}}{1+R_{j}}\phi,$
(23)

with

$\phi(R,t)=\frac{2l_{\mathrm{eff}}}{3}\frac{1+R_{j}}{1-R_{\phi}}\frac{\partial\phi}{\partial
z}\equiv h\frac{\partial\phi}{\partial z}$ (24)

at the top of the volume (z=0). Ignoring the bottom boundary entirely for the
moment, the mirror source’s fluence rate $\phi_{m}$ can be calculated using
the boundary condition and infinite-volume solution [21]:

$\displaystyle\frac{\partial\phi_{m}}{\partial
z}-\frac{\phi_{m}}{h}=-\left(\frac{\partial}{\partial
z}-\frac{1}{h}\right)\frac{1}{\left(4\pi
Dt\right)^{3/2}}\mathrm{exp}\left(-\frac{\rho^{2}}{4Dt}-\beta
t\right)\mathrm{exp}\left(-\frac{(z+l_{\mathrm{eff}})^{2}}{4Dt}\right)$ (25)

at $z=0$, at which point we may use
$-(z+l_{\mathrm{eff}})=z-l_{\mathrm{eff}}$. Multiplying both sides by
$\mathrm{exp}(-z/h)$ gives:

$\displaystyle\frac{\partial}{\partial z}\left(e^{-z/h}\phi_{m}\right)=$
$\displaystyle\frac{\partial}{\partial z}\left(e^{-z/h}\frac{1}{\left(4\pi
Dt\right)^{3/2}}\mathrm{exp}\left(-\frac{\rho^{2}}{4Dt}-\beta
t\right)\mathrm{exp}\left(-\frac{(z-l_{\mathrm{eff}})^{2}}{4Dt}\right)\right)$
$\displaystyle\quad-\frac{2}{h}e^{-z/h}\frac{1}{\left(4\pi
Dt\right)^{3/2}}\mathrm{exp}\left(-\frac{\rho^{2}}{4Dt}-\beta
t\right)\mathrm{exp}\left(-\frac{(z-l_{\mathrm{eff}})^{2}}{4Dt}\right).$ (26)

This expression can be integrated, moving the new dampening terms
$\mathrm{exp}(-z/h)$ to the right side:

$\displaystyle\phi_{m}(R,t)$
$\displaystyle=\phi_{+}(z-l_{\mathrm{eff}})-\frac{2}{h}\int_{z}^{\infty}e^{-(z-z^{\prime})/h}\phi_{+}(z^{\prime}-l_{\mathrm{eff}})dz^{\prime}$
$\displaystyle=\phi_{+}(z-l_{\mathrm{eff}})-\frac{2}{h}\int_{0}^{\infty}e^{-l/h}\phi_{+}(z-l-
l_{\mathrm{eff}})dl,$ (27)

which consists of a mirror source at $z=-l_{\mathrm{eff}}$, and an
exponentially dampened line of sinks extending from $z=-l_{\mathrm{eff}}$ to
$-\infty$.

So far, we have extended the infinite-volume description to non-isotropic
injection and treated the two interfaces separately with different boundary
conditions. For a general solution, both need to be combined. Since both
interfaces are at least partially reflective, the solution will have to
consist of an infinite sum of sources, sinks and their mirror images to
account for the repeated reflections between both. We shall restrict the
treatment of this problem to the surface, $z=0$, and calculate the ratio
between two consecutive round trips through the volume at the surface. Using a
single source at $z=0$, the ratio between the signal at $\rho,z=0$ and the
signal from a mirror source at $(\rho,z=-2Z)$ is then

$\frac{\phi_{m}}{\phi}=\mathrm{exp}\left(-\frac{Z^{2}}{Dt}\right)=\mathrm{exp}\left(-\frac{3Z^{2}}{cl_{\mathrm{eff}}t}\right).$
(28)

Assuming typical diffuse optical transit times $ct\approx 1$ m and a
sufficiently large volume to justify the diffuse treatment,
$Z/l_{\mathrm{eff}}\gg 1$, this ratio will be small. In addition, for the
signal at the surface, only full integer reflections will contribute, and it
is reasonable to only consider the first roundtrip reflection. Hence, we can
write the full solution as the effective source distribution accounting for
the partial boundary condition at the surface, plus a mirror distribution
centered around $z=-2Z$, consisting of 6 terms in total:

$\displaystyle\phi(\rho,t,z=0)=$ $\displaystyle\
\phi_{+}(l_{\mathrm{eff}})+\phi_{+}(-l_{\mathrm{eff}})+\phi_{+}(l_{\mathrm{eff}}+2Z)+\phi_{+}(-l_{\mathrm{eff}}+2Z)$
$\displaystyle-\frac{2}{h}\int_{0}^{\infty}e^{-l/h}\phi_{+}(-l-l_{\mathrm{eff}})dl-\frac{2}{h}\int_{0}^{\infty}e^{-l/h}\phi_{+}(l+l_{\mathrm{eff}}+2Z)dl,$
(29)

once again using the compact definition of a source term $\phi_{+}$. If higher
precision is needed, the bottom mirror source distribution can be mirrored
once more by applying the same method that was used to obtain the top source
distribution, and so forth.

While constructing a mirror source image to account for a totally reflective
boundary is both simple and instructive, it may not be applicable in many
real-world applications. In the presence of loss at the lower boundary,
introducing a mixed boundary may be more appropriate. If optical properties of
the interface are known, this can be accomplished by following the derivation
shown here but may prove challenging in the presence of a lossy medium.
Employing a so-called extrapolated boundary [15], which involves shifting the
mirror source image further away from the physical boundary, may prove
valuable to approximate a lossy interface sufficiently. We discuss realistic
values for the interface’s reflectivity in section IV.4.

###  II.5 Instantaneous source and time-resolved measurement

The above formula clearly shows the signature of three properties of the
volume: diffusion constant, absorption constant, and thickness. Assuming the
pure ice refractive index is known, the diffusion constant is directly related
to only the scattering coefficient. The distribution $\phi(\rho,t)$ on the
surface can in principle be measured using photon counting detectors, which
offer time-resolved detection of even dim light with resolution of the order
of $\sim 100ps$. Since the diffusion approximation is only valid for many
scattering steps, and the effective scattering length is of the order of
$\sim$1m, it can be assumed that typical arrival times are at least of the
order of many nanoseconds. Using these assumptions and a refractive index of
$1.31$ [19], we can calculate the arrival time distributions at a hypothetical
detector for different offsets of the detector from the light source, and for
different thicknesses, shown in Figure 2.

Figure 2: Time-dependent fluence rate for an instantaneous point-source and a
detector separated by $\rho$. (a) Different detector offsets, showing the
fluence rate distribution for an infinitely deep slab (black lines) and with a
bottom boundary 10m deep (blue line). (b-d) Distribution for different
thicknesses, shown for detector separation of $\rho$=10m, 30m, 50m,
respectively.

Figure 2a shows the fluence rate as a function of time at four different
detector positions for an infinitely deep volume. The blue curves are the
corresponding fluence rates for a 10m-deep volume. The curves are all
normalized to unit area for the infinite volume. It can be seen that the
presence of a shallow boundary not only increases the fluence rate at the
surface, but it also shifts the distribution to longer arrival times. Assuming
a priori knowledge of optical properties, this is a clear signature of volume
thickness that could potentially be used to invert a time-resolved
measurement. With increasing distance between source and detector, this
signature apparently grows larger as well. This can be understood in the
picture of source and mirror source: Moving the detector further away from the
source exponentially dampens the fluence rate from the source but does not
proportionally change the separation from the mirror source at $z=-2Z$, hence
increasing the portion of the signal that bears the signature of volume
thickness. In Figure 2b-d we show the fluence rate distribution for different
volume thickness in three different detector offsets, 10m, 30m and 50m
respectively. Again, at large detector separation the difference between
different volume thicknesses becomes more pronounced, but takes the shape of
an almost symmetrical offset. With the parameters we assumed, it seems likely
that thicknesses around 25m could be measured in principle. For larger
thickness, a significant change in the shape of the curves is only present at
very large distance from the source, where any remaining signal would be weak.
It should be noted that normalizing the curves in Figure 2a individually is
not representative of an actual measurable signal, since absorption dominates
the model behavior at larger source-detector separations, and actual
measurable signal per injected pulse decreases. The magnitude of this decrease
is better quantified in terms of time-integrated fluence rate.

In addition to the signature of volume size there is a clear indication that
the effective scattering and absorption coefficients can be extracted
independently from a measured arrival time distribution, as was done for the
infinite-volume situation [10, 12]. Close examination of equation 29 reveals
that all terms present in the bounded volume model are of similar nature to
the infinite volume solution from equation 5. Diffusive terms
$\sim\mathrm{exp}(-1/Dt)$ dampen short arrival times and bear the signature of
the diffusion constant and hence the scattering coefficient. Absorptive terms
$\sim\mathrm{exp}(-\beta t)$ dampen dominantly long arrival times. We can
therefore conclude that in a similar fashion to the deep measurements [10,
12], scattering and absorption coefficient can be independently extracted from
measurements at the surface, assuming that the ice thickness is either known,
or detector separation is small, in which case thickness has a minor influence
on the detector signal.

Another point worth noting is that the modification of the infinite solution
achieved by introducing the mixed boundary condition at the surface is most
prominent at small detector separation, mostly influencing the peak, and
becomes negligible otherwise.

###  II.6 Stationary sources and stationary measurement

While time-resolved measurements have the potential of proving abundant
information, they are not as simple as time-integrated intensity measurements.
A stationary measurement, i.e., a continuous source and non-time-resolving
detector that only measures intensity, would be much more desirable for a
portable, low-cost measurement apparatus. Using non-time-resolved detection
would yield the time-integrated, stationary fluence rate, or short fluence
$\phi_{S}(\rho)$ at the surface:

$\phi_{S}(\rho)=\int_{0}^{\infty}\phi(\rho,t,z=0)dt.$ (30)

Employing a stationary source gives an expression that can be shown to be
identical with the stationary fluence. For this, the fluence rate at time $t$
is the instantaneous fluence rate integrated over all past times from
$-\infty$:

$\int_{-\infty}^{t}\phi(\rho,t-t^{\prime},z=0)dt^{\prime}=-\int_{\infty}^{0}\phi(\rho,t,z=0)dt=\int_{0}^{\infty}\phi(\rho,t,z=0)dt=\phi_{S}(\rho).$
(31)

We find that all time-dependencies in the fluence rate (equation 29) are of
the form

$C\int_{0}^{\infty}dt\frac{1}{t^{3/2}}\mathrm{exp}\left(-\left(\frac{a}{t}+bt\right)\right)=C\sqrt{\frac{\pi}{a}}\mathrm{exp}\left(-2\sqrt{ab}\right),$
(32)

where $C,a,b$ are constant. Calculating the fluence from equations 32 for one
term of the above form yields a stationary source term:

$\phi_{S+}(\rho)=\frac{1}{4\pi
D}\frac{1}{\sqrt{\rho^{2}}}\mathrm{exp}\left(-\sqrt{\frac{3\rho^{2}}{l_{\mathrm{eff}}l_{\mathrm{abs}}}}\right).$
(33)

Applying this recipe to equation 29 can be done term-by-term, yielding the
complete fluence:

$\displaystyle\phi_{S}(\rho)=$ $\displaystyle
2\phi_{S+}(\rho^{2}+l_{\mathrm{eff}}^{2})+\phi_{S+}(\rho^{2}+(l_{\mathrm{eff}}+2Z)^{2})+\phi_{S+}(\rho^{2}+(-l_{\mathrm{eff}}+2Z)^{2})$
$\displaystyle+\frac{2}{h}\int_{0}^{\infty}e^{-\frac{l}{h}}\phi_{S+}(\rho^{2}+(-l-l_{\mathrm{eff}})^{2})dl+\frac{2}{h}\int_{0}^{\infty}e^{-\frac{l}{h}}\phi_{S+}(\rho^{2}+(l+l_{\mathrm{eff}}+2Z)^{2})dl.$
(34)

The result consists of the beforementioned 6 terms: 2 sources, 1 sink, and
their mirror images, whereby the first two source terms are identical at
$z=0$. In terms of a injected energy, the time-integrated fluence rate has the
units of $\mathrm{[Jm^{-2}]}$. Integrating over a hypothetical detector
aperture would yield the detected energy. Each term displays an exponential
decay with separation $\rho$ from the source, with each term again bearing
different offsets to this value due to different effective source position.
The exponentials fall off as
$\mathrm{exp}(-\sqrt{3\rho/l_{\mathrm{eff}}l_{\mathrm{abs}}})$. Hence, a fit
to data would reveal the product of scattering and absorption coefficients.

Figure 3: Time-integrated fluence rate for an instantaneous point-source and a
detector separated by $\rho$, compared for (a-b) different scaling parameters
$l_{\mathrm{eff}}l_{\mathrm{abs}}$ on linear and log-scale, respectively, (c)
infinite and 10m-deep volume, and (d) for constant scaling constant
$l_{\mathrm{eff}}l_{\mathrm{abs}}$ but varying scattering length
$l_{\mathrm{eff}}$.

The integrated fluence rate for different system parameters is shown in Figure
3. Panel a shows the fluence rate for constant scattering length
$l_{\mathrm{eff}}=1$ m, but different absorption lengths, resulting in a
different scaling parameter $l_{\mathrm{eff}}l_{\mathrm{abs}}$. Panel b shows
the same data on a log-scale. The difference is most apparent in the tails of
the distribution, where the curves cover about one order of magnitude for a
realistic range of parameters, displaying different slopes on a log-scale.

Figure 3c shows the same curve for two different geometries, the infinitely
deep volume, and a shallow one with a depth of 10m. It is clear that the
difference is minor and only shows up in the tail of the distribution. This
hints that the integrated fluence rate is not an ideal measure of volume
geometry, especially compared to the time-dependent fluence rate. This
highlights that stationary measurements are more useful for extracting optical
properties rather than measuring thickness - especially so since the
measurement is largely independent of thickness.

Figure 3d shows the distribution for different scattering lengths, but
constant overall scaling constant $l_{\mathrm{eff}}l_{\mathrm{abs}}$. In
comparison to panel b the slope in the tail (corresponding to the exponential
term) is retained, while the behavior at small detector separation is
different for each set of parameters. A difference is only visible for small
detector separation, where the diffusion approximation may not hold. This
implies that independent measurement of scattering and absorption is not
possible in this configuration.

Lastly, the integrated fluence rate can answer question of the absolute
measurable signal level in a time-resolved measurement, and its fall-off with
detector position $\rho$. According to equation 31, the integrated fluence
rate is identical to the integral over the entire arrival time distribution
shown in Figure 2. Therefore, the integrated fluence rate is a measure for the
total measured signal at position $\rho$ in a time-resolved measurement,
revealing the overall system loss.

##  III Monte Carlo Simulations

The model developed so far is based on several assumptions and approximations.
We’re assuming that the physical Mie scattering process can be approximated
with an effective isotropic process, and that the diffusion approximation
holds, specifically at the boundaries. The bottom boundary is modeled as a
perfect mirror, which surely is far from reality. To test if these assumptions
are reasonable, we are comparing the model to Monte Carlo simulations
implementing a more physical process. However, these simulations have their
own limitations, especially since a computationally efficient physical
scattering phase function has to be chosen. Agreement with a realistic,
physical scattering process is only guaranteed if the phase function is well
known.

The simulations calculate the trajectory of an individual photon. Scattering
is implemented using a Henyey-Greenstein distribution [25] around the pitch-
axis of the photon flight vector, and a uniform distribution around the roll
axis, which is the vector itself, using the Rodrigues’ rotation formula to
calculate. The asymmetry factor for the employed Henyey-Greenstein
distribution is $g=\langle\cos\theta\rangle=0.75$. The two random angles
$\theta,\gamma$ in each scattering step are obtained from uniformly
distributed random numbers $\xi_{1},\xi_{2}\in\\{0,1\\}$ using the inverse
cumulative Henyey-Greenstein density function

$\displaystyle\theta=\mathrm{cos}^{-1}\left[\frac{1}{2g}\left(1+g^{2}-\left(\frac{1-g^{2}}{1+g-2g\xi_{1}}\right)^{2}\right)\right],$
(35) $\displaystyle\gamma=2\pi\xi_{2}.$ (36)

Since this approximation to the true Mie scattering distribution only accounts
for the forward component of scattering, artifacts around the injection point
have to be expected. While the bottom boundary is assumed to be lossless and
cosine-distributed

$\displaystyle\theta=\mathrm{sin}^{-1}\left(\xi_{1}\right),$ (37)
$\displaystyle\gamma=2\pi\xi_{2}.$ (38)

Given that the bedrock of a glacier will have some degree of microscopic and
macroscopic roughness, this approximation seems appropriate, and is in line
with experimentally observed emission profiles [26]. In the presence of a
water or air layer below the ice, a (partial) specular reflection may be more
appropriate. For cases $l_{eff}\ll Z$, the distribution will not matter, since
randomization of the distribution will occur regardless of the scattering
distribution at the interface. For all cases under investigation in the
following chapter we observe no significant differences in the simulated
distributions for reflective, cosine or uniformly distributed angles at the
lower interface.

A Fresnel reflection is implemented at the top boundary, assuming unpolarized
light and summing the contribution of both linear polarizations, assuming a
refractive index of ice of 1.31 [19]. For this refractive index, total
internal reflection takes place for trajectories that intersect the surface at
an angle larger than $\theta_{C}=50^{\circ}$. Each step implements absorption
in the form of a survival probability of
$\mathrm{exp}(-l_{\mathrm{sca}}/l_{\mathrm{abs}})$ and trajectories are
terminated either by absorption or return to the surface. Figure 4 shows 2000
trajectories for a deep and a shallow volume and denotes injection point,
effective source position and boundaries. The overlayed trajectories display
the usual pear-shaped distribution.

Figure 4: Simulated random flights for 2000 photons with constant scattering
length $l_{\mathrm{eff}}=1$m and absorption length $l_{\mathrm{abs}}=100$m.
The bottom boundary implements a total, diffuse (cosine-distributed)
reflection, the top boundary implements Fresnel reflection. Light is injected
straight down, getting scattered after one physical scattering length, whereas
the model’s effective source would be one effective scattering length down.
(a) Shallow volume, $Z=10$m, (b) Deep volume, $Z=50$m.

For each photon trajectory the endpoint and time-of-flight are accessible. A
subset of trajectories ending in a defined detector area of
0.01$\mathrm{m^{2}}$ \- achievable with off-the-shelf optics - is evaluated in
histograms to simulate the outcome of a time-of-flight measurement. Surface
intensity measurement can be simulated easily by counting all trajectories
ending within detector areas along a line on the surface.

##  IV Feasibility of measurements

The simulations can be compared to the theoretical model. The simulations can
easily predict photon count rates and measurement duration for a given source
brightness. For a successful measurement of any parameters, it should be
possible to extract these parameters by fitting the theoretical curves to the
simulated data.

###  IV.1 Intensity measurements

The easiest measurement to implement is an intensity measurement as a function
of detector position. Here, photon counting detectors such as photo multiplier
tubes or avalanche photo diodes offer good efficiency, temporal resolution,
and more importantly a dynamic range over many orders of magnitude. This helps
when extracting meaningful data from measurements. To ascertain how well this
technique would fare given the complicated nature of equation 34, we performed
simulations for different combinations of parameters $l_{\mathrm{eff}}$ and
$l_{\mathrm{abs}}$, corresponding to a scaling parameter
$l_{\mathrm{eff}}l_{\mathrm{abs}}$ between 20$\mathrm{m^{2}}$ and
250$\mathrm{m^{2}}$, covering a wide range of expected values. The results are
shown in Figure 5.

Figure 5: Simulated time-integrated detector counts per injected photon on the
surface for 100cm2 detector area with $10^{8}$ simulated photon trajectories.
(a) Reconstructed vs simulated scaling constant (blue curve) over the ideal
ratio of 1 (black curve). Points shown contain combination of scattering and
absorption lengths in the range of $l_{\mathrm{eff}}=0.5...1.5m$ and
$l_{\mathrm{abs}}=20...250m$. Data is averaged in bins of width $20m^{2}$ with
$R^{2}=0.99$ ($R^{2}=0.93$ for the raw data). (b) Simulation with scaling
constant $l_{\mathrm{eff}}l_{\mathrm{abs}}=130\mathrm{m^{2}}$ and $10^{8}$
photon trajectories. The solid line is a fit of the mixed-boundary diffusion
model to the simulated histogram, excluding the center portion. Reconstructed
scaling constant is $l_{\mathrm{eff}}l_{\mathrm{abs}}=134(\pm
7)\mathrm{m^{2}}$. The inset shows the two-dimensional distribution on the
surface with the dashed line indicating where the cut was taken.

For fitting the theoretical model, we scale it by detector area. Simulated
data points $y$ are inversely weighted as $1/y$ to account for the large
dynamic range of 4-5 orders of magnitude. As a result, the scaling parameter
used for the simulation can be extracted. Accuracy is in the range of 10% for
smaller values. Faithful fitting seems increasingly difficult at larger values
for $l_{\mathrm{eff}}l_{\mathrm{abs}}$. While scattering and absorption length
are used as independent fitting parameters as prescribed by equation 34, only
their product yields reasonable agreement with the expected value. We believe
that this is caused by the invalidity of the diffusion approximation in
proximity to the detector. In Figure 5b we show one example of simulated
measurement data for a scaling constant of
$l_{\mathrm{eff}}l_{\mathrm{abs}}=130\mathrm{m^{2}}$. The fit retrieves a
value of $l_{\mathrm{eff}}l_{\mathrm{abs}}=134(\pm 7)\mathrm{m^{2}}$. In the
lowest counting detector positions, detection probability per injected photon
is of the order of $10^{-9}$. Assuming injection rates of
$10^{15}\mathrm{s}^{-1}$, which corresponds to roughly $1$ mW of optical
power, we can expect count rates of the order of $10^{6}s^{-1}$. This is
equivalent of about $1$ pW of optical power. Accounting for limited detector
and collection efficiency, these estimates leave ample room to clear ambient
nighttime background using photon counting detectors or highly efficient photo
diodes.

Note that product of absorption and scattering coefficient is sufficient above
600 nm, where absorption is well known and independent of impurities [7, 8].
Physical scattering length $l_{\mathrm{sca}}$ has to be independent of
wavelength, and the wavelength-dependence of the scattering angle distribution
(asymmetry factor) can be modeled through Mie scattering theory, which can
serve as a waypoint to extrapolate scattering coefficients measured above 600
nm to short wavelengths.

Even though information for large scaling parameters is hard to retrieve with
the presented method, the small experimental effort necessary makes it an
attractive measurement regardless of its limitations. Only a diode laser
source, detector and collection optics are required.

###  IV.2 Time-of-flight measurements

Time-resolved measurements can potentially provide more information, but at
the cost of lower counting rates, since the signal will be distributed in many
discrete time bins. We first assume that the volume depth is either known or
can be assumed sufficiently large not to play a role. Using these assumptions,
we apply a least-square fit to the simulated arrival time histograms and
attempt to retrieve the simulation parameters for scattering and absorption
length.

Figure 6: Simulated arrival time histogram of detector counts per injected
photon for 100cm2 detector area, fitted with diffusion theory curves. The
dataset combines nine different combinations of scattering and absorption
length in the range of $l_{\mathrm{eff}}=0.5...1.5m$ and
$l_{\mathrm{abs}}=25...250m$. (a) Reconstructed vs simulated effective
scattering length. (b) Reconstructed vs simulated absorption length. (c)
Example on a simulated arrival time histogram at a detector position 15m away
from the source, scattering length of $l_{\mathrm{eff}}=1$m and absorption
length $l_{\mathrm{abs}}=100$m. The least-square fit retrieves parameters
$l_{\mathrm{eff,fit}}=1.14(\pm 0.03)$m and $l_{\mathrm{abs,fit}}=88(\pm 5)$m.

The results for several combinations of scattering and absorption length are
shown in Figure 6. The ratio between retrieved and simulated parameter is
shown as a function of detector position. It can be seen that for large
detector offset, both parameters asymptotically reproduce the simulated
parameters within the variance of the dataset. Achievable relative error is
below 20% for detector positions at least 15m away from the source. The
asymptotic behavior can be explained by the fact that the diffusion
approximation only holds for large detector separations. This appears to be
much more apparent for time-resolved measurements than it is for time-
integrated ones. Figure 6c shows an example of a fit. In general, these fits
are much more robust than the ones of surface distribution. This is due to the
separate influence of scattering and absorption length on different sections
of the curve, but also due to the smaller dynamic range of 1-2 orders of
magnitude. We found that fit solutions are largely independent of starting
values.

For the evaluation of the simulated arrival time distributions we chose a bin
width of 50 ns, which yields a reasonable compromise between sampling the
distribution with enough points and producing sufficient number of counts per
bin. Using this temporal resolution of 50 ns we estimate return loss of the
order of $10^{-10}$ for the bins in the tail of the distribution. This is well
above the temporal resolution of modern photon counting detection systems,
which is typically of the order of 100ps for combined response of the detector
and electronics. Since optical parameters can be retrieved using such bin
widths, the use of inexpensive nanosecond-pulsed lasers is possible. These
lasers typically provide pulse energies of the order of 1nJ ($10^{9}$ photons
per pulse), with pulse repetition rates of $10^{5}\mathrm{s}^{-1}$, limited by
the width of the arrival time distribution. Taking combined detector and
collection efficiency of $10^{-3}$ into account, we expect count rates for the
lowest counting bins of the order of $10s^{-1}$. This is reasonable for a
time-resolved measurement, where any background is constant. It is worth
noting that the technique gives a truly independent measurement of scattering
and absorption length, without assumptions or reference values used previously
[7, 8]

###  IV.3 Thickness

Figure 7: Simulated arrival time histogram of detector counts per injected
photon for 100cm2 detector area. (a) Reconstructed vs simulated volume depth.
Confidence intervals for unreliable fits are not shown. (b) Example on a
simulated arrival time histogram for a volume depth of $Z=25m$ with the
detector situated at $\rho=35$m. The recovered volume depth is $Z=24(\pm 9)m$.

Measuring ice thickness is the most challenging of the measurements
investigated here, since it is only far away from the source that the volume
depth causes a large modification to the diffusion model. Results for volume
depths up to 25m are shown in Figure 7. It can be seen that at small detector
separation the thickness is systematically overestimated. However, in those
cases the confidence intervals of the fits (not shown where the relative
confidence interval is larger than 100%) are not useful. Retrieved thicknesses
trend towards the simulated value at a detector position around 30 m from the
source. At larger separation, fits still converge, but do not converge to the
correct value. This can be explained in terms of the theoretical curves shown
earlier: At large separation, the signature of the truncated volume turns into
what is almost symmetrical around the peak, and is not strongly dependent on
the volume thickness anymore. To avoid mismeasurement it seems appropriate to
perform the measurement as a function of detector position, and choose the
closest separation where fits converge with reasonable confidence. If we
assume loss at the lower boundary, parameter retrieval is less reliable, but
still possible for moderate loss. We show results for reflectivity $R=0.9$ and
$R=0.75$ in Figure 8. For $R=0.9$, the fits retrieve thickness values that are
comparable to the ones shown previously in Figure 7 for the case of a
perfectly reflecting lower boundary. For $R=0.75$, it can be seen that
especially for small source-detector separation, the thickness retrieval
fails. As loss is introduced to the lower boundary, the portion of the
detectable signal stemming from the reflection is decreased, explaining the
loss of this signature. For shallow volumes, this effect is more pronounced
because multiple reflections are more likely. In addition, systematic
deviations from the simulated volume depth occur, with relative errors of the
order of 50%. We can therefore conclude that the idealized diffusion model can
tolerate small amounts of loss of the order of 10%, but further modeling is
needed to account for larger losses at the boundary.

Figure 8: Reconstructed vs simulated volume dept with a lossy lower boundary.
Confidence intervals for unreliable fits are not shown. (a) Lower boundary
reflectivity $R=0.9$. (b) Lower boundary reflectivity $R=0.75$.

While this approach leaves much to be desired and further investigation is
warranted, the results show that thickness retrieval is possible in principle.
Askebjer et al. [22] used a more sophisticated approach for parameter
retrieval, where they used curve fitting to only extract maximum and mean of
the distribution, and then used the optical parameters extracted from those
two figures as a starting point for a maximum likelihood estimation using
simulations. However, they did not apply their procedure to retrieval of
propagation distances, since those were known. Given that we assumed exact
knowledge of optical properties and unit reflectivity of the lower boundary,
we conclude that this measurement of ice thickness may be very limited under
realistic conditions, and that further research is needed to develop this
technique.

To judge their usefulness, these results need to be put into context. The
assumptions made about optical properties and access to bare ice definitely
restrict the described experiment to the ablation zone of temperate glaciers.
Here in Oregon, the thickness of glaciers in the Cascade range was last
measured in the 1980’s [27]. At that time, ice thickness rarely exceeded 50m,
while some of the glaciers under investigation have since disappeared
completely [28]. The accuracy of depth retrieval demonstrated here is of the
orders of 10%, on par with comparable techniques using radio waves.

###  IV.4 Experimental considerations

Measuring glacier thickness is the most challenging experiment proposed here.
The results presented here show that a signature of volume thickness can be
extracted in principle and under ideal conditions. While we have also shown
that optical properties, which need to be known as precisely as possible to
retrieve thickness, can be measured independently, further modeling is needed
in order to test the robustness of parameter retrieval. This is particularly
true in regard of the lower boundary and realistic structure of the ice.
Lastly, if a precise model is needed, the commonly used Henyey-Greenstein
formula as well as mixed phase functions used by others may be insufficient
since they do not reproduce the back-scattering peak of the Mie distribution
[25, 13].

Bedrock reflection can be as high as 70 % for sandstone [26], 50 % for granite
and as low as 10 % for andesite [29]. Additional issues such as air and water
between the ice and bedrock further complicate formulating an appropriate
lower boundary for the diffusion model, since they could either enhance or
reduce reflectivity, depending on geometry and optical properties. Therefore,
knowledge on local geology and glacier morphology is essential for successful
depth measurement. Again, establishing bounds on how robust the diffusion
model needs to be for the measurement to be successful will be essential.

All measurements outlined here will require to measure small quantities of
light on photon counting detectors, which will be subject to intrinsic
detector dark counts of the order of 10 $\mathrm{s}^{-1}$ to 1000
$\mathrm{s}^{-1}$. In addition, we have to assume a sizable amount of
background light, even at night, which can be somewhat suppressed using
spectral filtering around the employed laser wavelength. Both dark and
background counts are in principle constant. Using modulation, either by
modulating the laser or using a mechanical shutter, changes in those values on
scales longer than the modulation period can be accounted for, and detection
thresholds of 1$\mathrm{s}^{-1}$ are possible. In the case of time-resolved
measurements, a pulsed laser provides the periodic pulse train. We have shown
that detection bin widths of roughly 50 ns are sufficient to resolve the
expected arrival time distributions, which span several microseconds. If a
long pulse is used, the measured distribution will be the convolution of the
instantaneous source distribution with the laser pulse shape. As long as the
pulse duration is much shorter than the detection bin width, the broadening
cannot be resolved, and the instantaneous distribution can be used for
fitting.

Given the simple nature of these proposed experiments we believe that they can
be implemented at relatively low cost. With a photon counting detector, a ns-
pulsed diode laser, acquisition electronics and minimal optics for collection,
an adequate setup can be implemented with off-the-shelf components for under
$10,000. Stationary measurements are possible at much lower cost. In medical
imaging, diffuse optical tomography in the frequency-domain, where a
continuous wave laser is modulated at different frequencies, has recently been
gaining in popularity [30, 31, 32]. This has resulted in further reduction in
cost.

##  V Conclusion

We extended the diffusion model used in AMANDA [22] with boundary conditions
necessary to describe a semi-infinite volume of glacier ice. The modifications
from these boundary conditions make it possible to extract both optical
properties and the size of a shallow volume. The stationary fluence encodes
information on optical properties.

In the framework of the simplified toy model presented here, the feasibility
of several different measurement scenarios using laser sources and photon
counting detectors has been demonstrated. Measuring the intensity distribution
on the ice surface can offer some insight into scattering and absorption
coefficients, but they cannot be measured independently without reference.
However, the simplicity of this measurement makes it easy to deploy,
potentially enabling a vast increase in data collection that can help refine
models for radiative heating of ice and snow. A measurement apparatus can
potentially be lightweight and easily fit in a backpack. Especially in the
mountain ranges of the western United States, where many glaciers are located
in designated Wilderness where motorized access is prohibited, the portability
of the measurement apparatus could open new frontiers and enable repeat
measurements. While time-of-flight measurements come with the increased effort
of coincidence measurements, nanosecond temporal resolution, and high demands
on collection efficiency, the measurement apparatus would remain similarly
mobile, and can be realized with off-the-shelf components. Not requiring
boreholes and being non-invasive would make such an apparatus even easier to
deploy in the field. Using diffuse optics and time-of-flight measurements is,
to the best of the authors’ knowledge, the only viable option to obtain
independent measurements of scattering and absorption coefficients, without
reference values and assumptions about the nature of scattering processes, and
has been proven in the field for deep ice [12, 13]. Moreover, these
measurements require only small separation between source and detector of the
order of 10 meters, where such measurements appear to insensitive to the
volume size. The simple retrieval scheme used here showed systematic
overestimation of the scattering length by up to 20%, and further
investigation is needed to remedy this issue. Once optical parameters are
known, such measurements are able to reveal the depth of shallow bodies of ice
up to 25m in thickness. Since the diffusion description and Monte Carlo
simulations are both approximations, further research in the physical nature
of the scattering phase function may be needed to facilitate reliable ranging.
Future research should also focus both on range of realistic values for
bedrock reflectivity and additional modeling including appropriate boundary
conditions for absorptive interfaces. While our results imply that diffuse
optical ranging cannot compete with radar techniques in range [16], an optical
measurement device can be built using commercially available components at
small cost and size. For the shrinking temperate glaciers in the western
United States this technique could serve as an early warning system for dying
glaciers and improve hydrological forecasting. Measurement of optical
properties of ice in the ablation zone can help to improve radiative transfer
models and reveal valuable information on paleoclimate history and local ice
properties.

#### Acknowledgments.

We thank Anders Carlson and Aaron Hartz of the Oregon Glacier Institute,
Stephen G. Warren and Matthew G. Cooper for fruitful discussions that inspired
this work.

#### Disclosures.

The authors declare no conflicts of interest.

#### Data availability.

Data underlying the results presented in this paper are not publicly available
at this time but may be obtained from the authors upon reasonable request.

## References

  * Bay  _et al._ [2001] R. C. Bay, P. B. Price, G. D. Clow, and A. J. Gow, Climate logging with a new rapid optical technique at siple dome, Geophysical Research Letters 28, 4635 (2001), https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2001GL013763 .
  * Bramall  _et al._ [2013] N. E. Bramall, R. C. Bay, K. Woschnagg, R. A. Rohde, and P. B. Price, South pole glacial climate reconstruction from multi-borehole laser particulate stratigraphy, Journal of Glaciology 59, 1117–1128 (2013).
  * Weikusat  _et al._ [2017] I. Weikusat, D. Jansen, T. Binder, J. Eichler, S. H. Faria, F. Wilhelms, S. Kipfstuhl, S. Sheldon, H. Miller, D. Dahl-Jensen, and T. Kleiner, Physical analysis of an antarctic ice core-towards an integration of micro- and macrodynamics of polar ice, Philosophical transactions. Series A, Mathematical, physical, and engineering sciences 375, 20150347 (2017).
  * Rongen  _et al._ [2020] M. Rongen, R. C. Bay, and S. Blot, Observation of an optical anisotropy in the deep glacial ice at the geographic south pole using a laser dust logger, The Cryosphere 14, 2537 (2020).
  * Hellmann  _et al._ [2021] S. Hellmann, J. Kerch, I. Weikusat, A. Bauder, M. Grab, G. Jouvet, M. Schwikowski, and H. Maurer, Crystallographic analysis of temperate ice on rhonegletscher, swiss alps, The Cryosphere 15, 677 (2021).
  * Kokhanovsky and Zege [2004] A. A. Kokhanovsky and E. P. Zege, Scattering optics of snow, Appl. Opt. 43, 1589 (2004).
  * Warren  _et al._ [2006] S. G. Warren, R. E. Brandt, and T. C. Grenfell, Visible and near-ultraviolet absorption spectrum of ice from transmission of solar radiation into snow, Appl. Opt. 45, 5320 (2006).
  * Cooper  _et al._ [2021] M. G. Cooper, L. C. Smith, A. K. Rennermalm, M. Tedesco, R. Muthyala, S. Z. Leidman, S. E. Moustafa, and J. V. Fayne, Spectral attenuation coefficients from measurements of light transmission in bare ice on the greenland ice sheet, The Cryosphere 15, 1931 (2021).
  * Demuth [2011] M. N. Demuth, Lidar in glaciology, in  _Encyclopedia of Snow, Ice and Glaciers_, edited by V. P. Singh, P. Singh, and U. K. Haritashya (Springer Netherlands, 2011) pp. 713–722.
  * Price and Bergström [1997] P. B. Price and L. Bergström, Optical properties of deep ice at the south pole: scattering, Appl. Opt. 36, 4181 (1997).
  * Askebjer  _et al._ [1995] P. Askebjer, S. W. Barwick, L. Bergström, A. Bouchta, S. Carius, A. Coulthard, K. Engel, B. Erlandsson, A. Goobar, L. Gray, A. Hallgren, F. Halzen, P. O. Hulth, J. Jacobsen, S. Johansson, V. Kandhadai, I. Liubarsky, D. Lowder, T. Miller, P. C. Mock, R. Morse, R. Porrata, P. B. Price, A. Richards, H. Rubinstein, E. Schneider, Q. Sun, S. Tilav, C. Walck, and G. Yodh, Optical properties of the south pole ice at depths between 0.8 and 1 kilometer, Science 267, 1147 (1995), https://science.sciencemag.org/content/267/5201/1147.full.pdf .
  * Ackermann  _et al._ [2006] M. Ackermann, J. Ahrens, X. Bai, M. Bartelt, S. W. Barwick, R. C. Bay, T. Becka, J. K. Becker, K.-H. Becker, P. Berghaus, E. Bernardini, D. Bertrand, D. J. Boersma, S. Böser, O. Botner, A. Bouchta, O. Bouhali, C. Burgess, T. Burgess, T. Castermans, D. Chirkin, B. Collin, J. Conrad, J. Cooley, D. F. Cowen, A. Davour, C. De Clercq, C. P. de los Heros, P. Desiati, T. DeYoung, P. Ekström, T. Feser, T. K. Gaisser, R. Ganugapati, H. Geenen, L. Gerhardt, A. Goldschmidt, A. Groß, A. Hallgren, F. Halzen, K. Hanson, D. H. Hardtke, T. Harenberg, T. Hauschildt, K. Helbing, M. Hellwig, P. Herquet, G. C. Hill, J. Hodges, D. Hubert, B. Hughey, P. O. Hulth, K. Hultqvist, S. Hundertmark, J. Jacobsen, K. H. Kampert, A. Karle, M. Kestel, G. Kohnen, L. Köpke, M. Kowalski, K. Kuehn, R. Lang, H. Leich, M. Leuthold, I. Liubarsky, J. Lundberg, J. Madsen, P. Marciniewski, H. S. Matis, C. P. McParland, T. Messarius, Y. Minaeva, P. Miočinović, R. Morse, K. Münich, R. Nahnhauer, J. W. Nam, T. Neunhöffer, P. Niessen, D. R. Nygren, P. Olbrechts, A. C. Pohl, R. Porrata, P. B. Price, G. T. Przybylski, K. Rawlins, E. Resconi, W. Rhode, M. Ribordy, S. Richter, J. Rodríguez Martino, H.-G. Sander, S. Schlenstedt, D. Schneider, R. Schwarz, A. Silvestri, M. Solarz, G. M. Spiczak, C. Spiering, M. Stamatikos, D. Steele, P. Steffen, R. G. Stokstad, K.-H. Sulanke, I. Taboada, O. Tarasova, L. Thollander, S. Tilav, W. Wagner, C. Walck, M. Walter, Y.-R. Wang, C. H. Wiebusch, R. Wischnewski, H. Wissing, and K. Woschnagg, Optical properties of deep glacial ice at the south pole, Journal of Geophysical Research: Atmospheres 111 (2006), https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2005JD006687 .
  * Aartsen  _et al._ [2013] M. Aartsen, R. Abbasi, Y. Abdou, M. Ackermann, J. Adams, J. Aguilar, M. Ahlers, D. Altmann, J. Auffenberg, X. Bai, M. Baker, S. Barwick, V. Baum, R. Bay, J. Beatty, S. Bechet, J. Becker Tjus, K.-H. Becker, M. Bell, M. Benabderrahmane, S. BenZvi, J. Berdermann, P. Berghaus, D. Berley, E. Bernardini, A. Bernhard, D. Bertrand, D. Besson, G. Binder, D. Bindig, M. Bissok, E. Blaufuss, J. Blumenthal, D. Boersma, S. Bohaichuk, C. Bohm, D. Bose, S. Böser, O. Botner, L. Brayeur, A. Brown, R. Bruijn, J. Brunner, S. Buitink, M. Carson, J. Casey, M. Casier, D. Chirkin, B. Christy, K. Clark, F. Clevermann, S. Cohen, D. Cowen, A. Cruz Silva, M. Danninger, J. Daughhetee, J. Davis, C. De Clercq, S. De Ridder, P. Desiati, M. de With, T. DeYoung, J. Díaz-Vélez, M. Dunkman, R. Eagan, B. Eberhardt, J. Eisch, R. Ellsworth, S. Euler, P. Evenson, O. Fadiran, A. Fazely, A. Fedynitch, J. Feintzeig, T. Feusels, K. Filimonov, C. Finley, T. Fischer-Wasels, S. Flis, A. Franckowiak, R. Franke, K. Frantzen, T. Fuchs, T. Gaisser, J. Gallagher, L. Gerhardt, L. Gladstone, T. Glüsenkamp, A. Goldschmidt, G. Golup, J. Goodman, D. Góra, D. Grant, A. Groß, M. Gurtner, C. Ha, A. Haj Ismail, A. Hallgren, F. Halzen, K. Hanson, D. Heereman, P. Heimann, D. Heinen, K. Helbing, R. Hellauer, S. Hickford, G. Hill, K. Hoffman, R. Hoffmann, A. Homeier, K. Hoshina, W. Huelsnitz, P. Hulth, K. Hultqvist, S. Hussain, A. Ishihara, E. Jacobi, J. Jacobsen, G. Japaridze, K. Jero, O. Jlelati, B. Kaminsky, A. Kappes, T. Karg, A. Karle, J. Kelley, J. Kiryluk, F. Kislat, J. Kläs, S. Klein, J.-H. Köhne, G. Kohnen, H. Kolanoski, L. Köpke, C. Kopper, S. Kopper, D. Koskinen, M. Kowalski, M. Krasberg, G. Kroll, J. Kunnen, N. Kurahashi, T. Kuwabara, M. Labare, H. Landsman, M. Larson, M. Lesiak-Bzdak, J. Leute, J. Lünemann, J. Madsen, R. Maruyama, K. Mase, H. Matis, F. McNally, K. Meagher, M. Merck, P. Mészáros, T. Meures, S. Miarecki, E. Middell, N. Milke, J. Miller, L. Mohrmann, T. Montaruli, R. Morse, R. Nahnhauer, U. Naumann, H. Niederhausen, S. Nowicki, D. Nygren, A. Obertacke, S. Odrowski, A. Olivas, M. Olivo, A. O’Murchadha, L. Paul, J. Pepper, C. Pérez de los Heros, C. Pfendner, D. Pieloth, N. Pirk, J. Posselt, P. Price, G. Przybylski, L. Rädel, K. Rawlins, P. Redl, E. Resconi, W. Rhode, M. Ribordy, M. Richman, B. Riedel, J. Rodrigues, C. Rott, T. Ruhe, B. Ruzybayev, D. Ryckbosch, S. Saba, T. Salameh, H.-G. Sander, M. Santander, S. Sarkar, K. Schatto, M. Scheel, F. Scheriau, T. Schmidt, M. Schmitz, S. Schoenen, S. Schöneberg, L. Schönherr, A. Schönwald, A. Schukraft, L. Schulte, O. Schulz, D. Seckel, S. Seo, Y. Sestayo, S. Seunarine, C. Sheremata, M. Smith, M. Soiron, D. Soldin, G. Spiczak, C. Spiering, M. Stamatikos, T. Stanev, A. Stasik, T. Stezelberger, R. Stokstad, A. Stößl, E. Strahler, R. Ström, G. Sullivan, H. Taavola, I. Taboada, A. Tamburro, S. Ter-Antonyan, S. Tilav, P. Toale, S. Toscano, M. Usner, D. van der Drift, N. van Eijndhoven, A. Van Overloop, J. van Santen, M. Vehring, M. Voge, M. Vraeghe, C. Walck, T. Waldenmaier, M. Wallraff, R. Wasserman, C. Weaver, M. Wellons, C. Wendt, S. Westerhoff, N. Whitehorn, K. Wiebe, C. Wiebusch, D. Williams, H. Wissing, M. Wolf, T. Wood, C. Xu, D. Xu, X. Xu, J. Yanez, G. Yodh, S. Yoshida, P. Zarzhitsky, J. Ziemann, S. Zierke, A. Zilles, and M. Zoll, Measurement of south pole ice transparency with the icecube led calibration system, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 711, 73 (2013).
  * Durduran  _et al._ [2010] T. Durduran, R. Choe, W. B. Baker, and A. G. Yodh, Diffuse optics for tissue monitoring and tomography, Reports on Progress in Physics 73, 076701 (2010).
  * Haskell  _et al._ [1994] R. C. Haskell, L. O. Svaasand, T.-T. Tsay, T.-C. Feng, M. S. McAdams, and B. J. Tromberg, Boundary conditions for the diffusion equation in radiative transfer, J. Opt. Soc. Am. A 11, 2727 (1994).
  * Plewes and Hubbard [2001] L. A. Plewes and B. Hubbard, A review of the use of radio-echo sounding in glaciology, Progress in Physical Geography: Earth and Environment 25, 203 (2001), https://doi.org/10.1177/030913330102500203 .
  * Rutishauser  _et al._ [2016] A. Rutishauser, H. Maurer, and A. Bauder, Helicopter-borne ground-penetrating radar investigations on temperate alpine glaciers: A comparison of different systems and their abilities for bedrock mapping, Geophysics 81, WA119 (2016), https://pubs.geoscienceworld.org/geophysics/article-pdf/81/1/WA119/3246748/geo2015-0144.pdf .
  * Libois  _et al._ [2013] Q. Libois, G. Picard, J. L. France, L. Arnaud, M. Dumont, C. M. Carmagnola, and M. D. King, Influence of grain shape on light penetration in snow, The Cryosphere 7, 1803 (2013).
  * Warren [2019] S. G. Warren, Optical properties of ice and snow, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 377, 20180161 (2019), https://royalsocietypublishing.org/doi/pdf/10.1098/rsta.2018.0161 .
  * Beres  _et al._ [2020] N. D. Beres, M. Lapuerta, F. Cereceda-Balic, and H. Moosmüller, Snow surface albedo sensitivity to black carbon: Radiative transfer modelling, Atmosphere 11, 10.3390/atmos11101077 (2020).
  * Bryan [1890] G. H. Bryan, An application of the method of images to the conduction of heat, Proceedings of the London Mathematical Society s1-22, 424 (1890), https://londmathsoc.onlinelibrary.wiley.com/doi/pdf/10.1112/plms/s1-22.1.424 .
  * Askebjer  _et al._ [1997] P. Askebjer, S. W. Barwick, L. Bergström, A. Bouchta, S. Carius, E. Dalberg, K. Engel, B. Erlandsson, A. Goobar, L. Gray, A. Hallgren, F. Halzen, H. Heukenkamp, P. O. Hulth, S. Hundertmark, J. Jacobsen, A. Karle, V. Kandhadai, I. Liubarsky, D. Lowder, T. Miller, P. Mock, R. M. Morse, R. Porrata, P. B. Price, A. Richards, H. Rubinstein, E. Schneider, C. Spiering, O. Streicher, Q. Sun, T. Thon, S. Tilav, R. Wischnewski, C. Walck, and G. B. Yodh, Optical properties of deep ice at the south pole: absorption, Appl. Opt. 36, 4168 (1997).
  * Chandrasekhar [1943] S. Chandrasekhar, Stochastic problems in physics and astronomy, Rev. Mod. Phys. 15, 1 (1943).
  * Bohren [1987] C. F. Bohren, Multiple scattering of light and some of its observable consequences, American Journal of Physics 55, 524 (1987), https://doi.org/10.1119/1.15109 .
  * Toublanc [1996] D. Toublanc, Henyey–greenstein and mie phase functions in monte carlo radiative transfer computations, Appl. Opt. 35, 3270 (1996).
  * Watson [1971] R. D. Watson, Spectral reflectance and photometric properties of selected rocks, Remote Sensing of Environment 2, 95 (1971).
  * Driedger and Kennard [1986] C. L. Driedger and P. M. Kennard,  _Ice Volumes on Cascade Volcanoes: Mount Rainier, Mount Hood, Three Sisters, and Mount Shasta_, Tech. Rep. (Survey, U.S. Geological, 1986) report.
  * Hartz and Carlson [2020] A. Hartz and A. E. Carlson, Glacier disappearance in the high alpine of oregon, usa, in  _AGU Fall Meeting 2020_ (AGU, 2020).
  * Zhou and Wang [2017] K.-F. Zhou and S.-S. Wang, Spectral properties of weathered and fresh rock surfaces in the xiemisitai metallogenic belt, nw xinjiang, china, Open Geosciences 9, 322 (2017).
  * O’Sullivan  _et al._ [2012] T. D. O’Sullivan, A. E. Cerussi, B. J. Tromberg, and D. J. Cuccia, Diffuse optical imaging using spatially and temporally modulated light, Journal of Biomedical Optics 17, 1 (2012).
  * Applegate  _et al._ [2020] M. B. Applegate, R. E. Istfan, S. Spink, A. Tank, and D. Roblyer, Recent advances in high speed diffuse optical imaging in biomedicine, APL Photonics 5, 040802 (2020), https://doi.org/10.1063/1.5139647 .
  * Fantini and Sassaroli [2020] S. Fantini and A. Sassaroli, Frequency-domain techniques for cerebral and functional near-infrared spectroscopy, Frontiers in Neuroscience 14, 300 (2020).