# Moments of Sunlight Pathlength in Water and Aerosol Clouds from O2
Spectroscopy: Exploitable Parameter Sensitivities

Anthony B. Davis  Quentin Libois  Nicolas Ferlay  Alexander Marshak

###### Abstract

SYNOPSIS: According to the latest IPCC assessment, clouds and aerosols remain
a major challenge in future climate prediction using Global Climate Models
(GCMs). Consequently, NASA’s 2017 Decadal Survey has made them, along with
convection and precipitation, a high priority by calling them out as
Designated Observables for future coordinated missions under the “ACCP”
banner, now known as the Atmospheric Observing System. Atmospheric science is
now more than ever driving a renewal in remote sensing techniques that can
probe clouds and aerosols more accurately and with improved sampling. Herein,
we bring new theoretical results that support a special form of Differential
Optical Absorption Spectroscopy (DOAS) that uses oxygen absorption features in
the visible/near-IR (VNIR) spectrum where there is an abundance of sunlight.
Having known concentration and cross-section, the remaining unknown in
${{\mathrm{O}}_{2}}$ DOAS (DO2AS) is the path that the light has followed
through the absorbing gas. In the presence of scattering by cloud and aerosol
particles, that path is broken at each interaction. Cumulative pathlength
through the well-mixed ${{\mathrm{O}}_{2}}$ gas thus becomes a continuous
random variable, and its probability distribution function (PDF) will contain
desirable information about the clouds and aerosols such as top and bottom
altitudes of the layer they occupy as well as its optical thickness. In this
study, we build on previous work to compute statistical moments of the
pathlength PDF. Specifically, we show that mean and variance of the pathlength
inside the scattering particulate medium convey different pieces of
information, namely, size and opacity of the optical medium, hence geometric
and optical thicknesses in the case of plane-parallel layers. We also extend,
with dense aerosol plumes in mind (e.g., wildfire smoke, volcanic ash,
elevated dust), the close connection between geometric thickness and mean
pathlength to the case of absorption by the scattering particles themselves.
In summary, we view pathlength moments as intermediate DO2AS products that can
be obtained from ${{\mathrm{O}}_{2}}$ spectroscopic data, and we describe a
preliminary algorithm to do just that. In turn, these moments yield cloud or
aerosol profile parameters of high interest: layer geometric and optical
thicknesses. One normally thinks about active sensors, radars and lidars, to
do such atmospheric profiling. In this case, the profile is parameterized and
representative of some kind of horizontal average determined by the multiple
scattering. However, there is a clear advantage in using passive over active
instrumentation, starting with the possibility of imaging over a large swath.

###### keywords:

oxygen A-band spectroscopy , DOAS , clouds , aerosols , pathlength
distribution , pathlength moments

††journal: JQSRT

[inst1] organization=NASA Jet Propulsion Laboratory, California Institute of
Technology, city=Pasadena, state=CA, country=United States

[inst2] organization=Centre National de Recherche en Meteorologie CNRS,
Universite de Toulouse, Meteo-France, city=Toulouse, country=France

[inst3] organization=Laboratoire d’Optique Atmospherique CNRS, Universite
Lille-1, city=Villeneuve d’Ascq, country=France

[inst4] organization=NASA Goddard Space Flight Center, Climate and Radiation
Laboratory, city=Greenbelt, state=MD, country=United States

New derivation of the “4V/S” invariance property of mean pathlength in
arbitrarily-shaped optical media from first principles of radiative transfer.

Verification that, while mean pathlength conveys direct information about the
size of the optical medium, pathlength variance increases with its opacity
and, unlike reflectivity, without bound.

First prediction of how absorption by the scattering particles affects mean
pathlength in plane-parallel media.

Preliminary algorithm for extracting pathlength moments from spectroscopic
data.

## Outline

In §1, we set up a theoretical and observational framework for relating DO2AS
at any spectral resolution to cloud properties through moments of the in-cloud
pathlength distribution that is grounded in time-dependent 3D radiative
transfer (RT). In §2, we present a new derivation for the mean pathlength in
arbitrary cloud geometry. In §3, new results for pathlength mean and variance
in plane-parallel geometry are obtained. On the one hand, we confirm that the
mean is all about cloud size, variance informs us directly about cloud opacity
in a regime where operational determinations of cloud optical thickness (COT)
loose sensitivity. An alternative approach for connecting 1st- and 2nd-order
pathlength moments to cloud properties based on discrete-time random walk
theory is presented in A, thus reinforcing our findings. On the other hand,
with aerosol clouds in mind, we wonder how absorption by the scattering
particles affects pathlength moments and quantify its impact on the mean. In
short, pathlength moments yield cloud properties, but it is not obvious how to
go from DO2AS observations to pathlength moments, so we describe a tentative
algorithm in B.

##  1 Time-Dependent RT and Pathlength Moments

Let $I(t,\mathbf{x},\boldsymbol{\Omega})$ be the time-dependent specific
intensity for a light beam intercepted at point $\mathbf{x}$ propagating in
direction $\boldsymbol{\Omega}$ expressed in [W/m2/sr]. Rather than time $t>0$
per se, we can use pathlength $ct$ cumulated since _everywhere simultaneous_
entrance into the arbitrarily-shaped optical medium M $\subset\mathbb{R}^{3}$
at its boundary $\partial$M, with $c$ denoting the speed of light in the
optical medium. Note that the physical units of
$I(ct,\mathbf{x},\boldsymbol{\Omega})=I(t,\mathbf{x},\boldsymbol{\Omega})/c$
are [J/m3/sr]. $I(ct,\mathbf{x},\boldsymbol{\Omega})$ satisfies the integro-
differential 3+1D RT equation (RTE):

$\left[\frac{\partial\;}{\partial
ct}+\boldsymbol{\Omega}\cdot\boldsymbol{\nabla}+\sigma(\mathbf{x})\right]\,I=\sigma_{\rm
s}(\mathbf{x})\int_{4\pi}P(\boldsymbol{\Omega}\cdot\boldsymbol{\Omega}^{\prime})I(ct,\mathbf{x},\boldsymbol{\Omega}^{\prime}){\rm
d}\boldsymbol{\Omega}^{\prime},$ (1)

where $\sigma(\mathbf{x})$ is the extinction coefficient in [1/m], while
$\sigma_{\rm s}(\mathbf{x})$ is the scattering coefficient, with
$P(\boldsymbol{\Omega}\cdot\boldsymbol{\Omega}^{\prime})$ being the scattering
phase function (PF) in [1/sr]. The PF is normalized as
$2\pi\int_{-1}^{+1}P(\mu_{\rm s}){\rm d}\mu_{\rm s}=1$, where $\mu_{\rm
s}=\boldsymbol{\Omega}\cdot\boldsymbol{\Omega}^{\prime}$ is the cosine of the
scattering angle. Sinks of radiant energy at point
$(ct,\mathbf{x},\boldsymbol{\Omega})$ in transport space are on the left-hand
side of (1), namely advection and extinction, while sources are on the right-
hand side, but only in-scattering of ambient light is considered in the
present study. Finally, a statement of initial ($t=0$) and boundary
($\mathbf{x}\in\partial{\rm M}$) conditions for incoming radiation is required
to define $I(t,\mathbf{x},\boldsymbol{\Omega})$ uniquely everywhere and at all
times.

We set the initial and boundary conditions (BCs) for the time-dependent 3D RTE
(1) to

$I(ct,\mathbf{x},\boldsymbol{\Omega})=\delta(ct)E_{0}/\pi,\text{ for
}\mathbf{x}\in\partial{\rm M}\text{ and
}\boldsymbol{\Omega}\cdot\mathbf{n}(\mathbf{x})<0,$ (2)

where $\delta(\cdot)$ is Dirac’s delta function, $E_{0}$ is a uniform amount
of radiant energy released into M per unit of surface on $\partial$M in
[J/m2], and division by $\pi$ accounts for the isotropic distribution of the
incoming light. Equation (2) describes a sudden, spatially and angularly
uniform injection of light into the medium M.

Because of the multiplicity of possible paths, starting at the boundary,
changing direction at every scattering, ending at
$(\mathbf{x},\boldsymbol{\Omega})$, after covering a distance $ct$, this
pathlength can be viewed as a random variable. Its probability distribution
function (PDF) in [1/m] is defined as

$f(ct,\mathbf{x},\boldsymbol{\Omega})=I(ct,\mathbf{x},\boldsymbol{\Omega})\left/\int_{0}^{\infty}I(ct,\mathbf{x},\boldsymbol{\Omega}){\rm
d}ct\right.,$ (3)

which is dependent on position in transport space
$(\mathbf{x},\boldsymbol{\Omega})$. From there, we can compute the
$q^{\text{th}}$-order statistical moment of the random variable $ct$:

$\Big{\langle}{(ct)^{q}}\Big{\rangle}(\mathbf{x},\boldsymbol{\Omega})=\int_{0}^{\infty}(ct)^{q}f(ct,\mathbf{x},\boldsymbol{\Omega})\,{\rm
d}ct,$ (4)

where $q\in\mathbb{N}$ might need to have an upper bound to ensure
convergence. We are particularly interested in $q$ = 1,2.

Consider the Laplace transform of $f(ct,\mathbf{x},\boldsymbol{\Omega})$:

$\tilde{f}(s,\mathbf{x},\boldsymbol{\Omega})=\int_{0}^{\infty}\exp(-s\,ct)\,f(ct,\mathbf{x},\boldsymbol{\Omega})\,{\rm
d}ct=\langle{\exp(-s\,ct)}\rangle(\mathbf{x},\boldsymbol{\Omega}),$ (5)

which is non-dimensional. This is known as the characteristic (or “moment-
generating”) function of a PDF supported by $\mathbb{R}^{+}$, and it can be
used to compute statistical moments using its successive derivatives at $s=0$:

$\displaystyle\Big{\langle}{(ct)^{q}}\Big{\rangle}(\mathbf{x},\boldsymbol{\Omega})={}$
$\displaystyle\left(-\frac{\partial\;}{\partial
s}\right)^{q}\left.\tilde{f}(s,\mathbf{x},\boldsymbol{\Omega})\right|_{s=0}$
$\displaystyle={}$ $\displaystyle\left(-\frac{\partial\;}{\partial
s}\right)^{q}\left.\left.\tilde{I}(s,\mathbf{x},\boldsymbol{\Omega})\right|_{s=0}\right/\tilde{I}(0,\mathbf{x},\boldsymbol{\Omega}).$
(6)

These are _non_ -centered moments. There is also the “cumulant-generating”
function $\log\tilde{f}(s,\mathbf{x},\boldsymbol{\Omega})$ that is used to
compute cumulants, i.e., the centered moments that are additive for sums of
random variables. Of particular interest here are the 1${}^{\text{st}}$\- and
2${}^{\text{nd}}$-order cumulants, namely, the mean and variance of pathlength
$ct$, respectively from $q=1,2$:

$\displaystyle\langle{ct}\rangle(\mathbf{x},\boldsymbol{\Omega})={}$
$\displaystyle\left.\left(-\frac{\partial\;}{\partial
s}\right)\log\tilde{f}(s,\mathbf{x},\boldsymbol{\Omega})\right|_{s=0};$ (7)
$\displaystyle\text{Var}[ct](\mathbf{x},\boldsymbol{\Omega})=\langle{(ct)^{2}}\rangle(\mathbf{x},\boldsymbol{\Omega})-\langle{ct}\rangle^{2}(\mathbf{x},\boldsymbol{\Omega})={}$
$\displaystyle\left.\left(-\frac{\partial\;}{\partial
s}\right)^{2}\log\tilde{f}(s,\mathbf{x},\boldsymbol{\Omega})\right|_{s=0}.$
(8)

By Laplace transforming the 3+1D RTE in (1) term-by-term, we see that

$\tilde{I}(s,\mathbf{x},\boldsymbol{\Omega})\propto\int_{0}^{\infty}\exp(-s\,ct)\,I(ct,\mathbf{x},\boldsymbol{\Omega})\,{\rm
d}ct$ (9)

obeys a steady-state 3D RTE:

$\left[\boldsymbol{\Omega}\cdot\boldsymbol{\nabla}+(\sigma(\mathbf{x})+s)\right]\,\tilde{I}=\sigma_{\rm
s}(\mathbf{x})\int_{4\pi}P(\boldsymbol{\Omega}\cdot\boldsymbol{\Omega}^{\prime})\tilde{I}(s,\mathbf{x},\boldsymbol{\Omega}^{\prime}){\rm
d}\boldsymbol{\Omega}^{\prime},$ (10)

irrespective of the proportionality constant that may be used in (9). Note
that the right-hand side of (9) is in [J/m2/sr]. For instance, at $s=0$, it
expresses what portion of the original radiant energy release $E_{0}$ [J/m2]
in (2) ends up at $(\mathbf{x},\boldsymbol{\Omega})$, irrespective of the path
$ct$ covered to get there. At the same time, the steady-state radiance on the
left-hand side of (9) is conventionally expressed in [W/m2/sr] because of the
steady sources. This minor conundrum is resolved when we set the BCs for the
steady-state problem further on.

In the above time-dependent 3D RTE (1), we can now identify $ct$ with the
total pathlength $L$, and the Laplace conjugate variable $s$ in (10) with the
spectral absorption coefficient $k_{\lambda}$ in [1/m] of an absorbing gas
that uniformly fills the scattering medium M. In that form, (9) becomes a
predictive model for the spectral radiance
$I_{\lambda}(\mathbf{x},\boldsymbol{\Omega})\equiv\tilde{I}(k_{\lambda},\mathbf{x},\boldsymbol{\Omega})$
based on the radiance field $I(L,\mathbf{x},\boldsymbol{\Omega})$ viewed as
steady-state, but decomposed according to $L$, hence with units of
[(W/m2/sr)/m]. In that guise, Eq. (9) has been described as the “equivalence”
theorem [1, 2]. We finally note that in (10), now viewed as a spectral RT
problem, with $s\mapsto k_{\lambda}$, that quantity is no longer an
independent variable. Rather, $k_{\lambda}$ is a parameter—a given property of
the optical medium.

##  2 Universality of $\langle{ct}\rangle$ for Arbitrary Conservative Optical
Media

Returning to the time-dependent 3+1D RTE in (1), we integrate over
$\boldsymbol{\Omega}$ in all $4\pi$ directions to obtain the continuity
equation

$\frac{\partial J}{\partial
ct}+\boldsymbol{\nabla}\cdot\mathbf{F}=-\sigma_{\rm a}(\mathbf{x})\,J,$ (11)

where $\sigma_{\rm a}(\mathbf{x})=\sigma(\mathbf{x})-\sigma_{\rm
s}(\mathbf{x})$ is the (particulate) absorption coefficient, or

$\boldsymbol{\nabla}\cdot\tilde{\mathbf{F}}=-\left(\sigma_{\rm
a}(\mathbf{x})+s\right)\,\tilde{J}$ (12)

in Laplace space, from (10), where

$\begin{array}[]{ccc}J(ct,\mathbf{x})=\int_{4\pi}I(ct,\mathbf{x},\boldsymbol{\Omega})\,{\rm
d}\boldsymbol{\Omega},&\mathbf{F}(ct,\mathbf{x})=\int_{4\pi}\boldsymbol{\Omega}\,I(ct,\mathbf{x},\boldsymbol{\Omega})\,{\rm
d}\boldsymbol{\Omega},&\text{in }[\text{J/m}^{3}]\\\
\tilde{J}(s,\mathbf{x})=\int_{4\pi}\tilde{I}(s,\mathbf{x},\boldsymbol{\Omega})\,{\rm
d}\boldsymbol{\Omega},&\tilde{\mathbf{F}}(s,\mathbf{x})=\int_{4\pi}\boldsymbol{\Omega}\,\tilde{I}(s,\mathbf{x},\boldsymbol{\Omega})\,{\rm
d}\boldsymbol{\Omega},&\text{in }[\text{W/m}^{2}]\\\ \end{array}$ (13)

in (11) and (12), respectively. $J(ct,\mathbf{x})$ is the instantaneous
_scalar_ (or actinic) flux, closely related to radiant energy _density_
$u(t,\mathbf{x})=J(ct,\mathbf{x})/c$, and $\mathbf{F}(ct,\mathbf{x})$ is its
_vector_ flux counterpart, a.k.a. radiant energy _current_. Their Laplace
transforms, $\tilde{J}(s,\mathbf{x})$ and $\tilde{\mathbf{F}}(s,\mathbf{x})$,
behave like steady-state versions with variable amounts of gaseous absorption.
Equations (11) and (12) express the conservation of radiant energy,
respectively, in the time- and Laplace domains.

We now apply the divergence theorem to the differentiable vector field
$\tilde{\mathbf{F}}(s,\mathbf{x})$:

$\int_{\rm M}\boldsymbol{\nabla}\cdot\tilde{\mathbf{F}}\,{\rm
d}V(\mathbf{x})=\int_{\partial{\rm
M}}\tilde{\mathbf{F}}\cdot\mathbf{n}(\mathbf{x})\,{\rm d}S(\mathbf{x}),$ (14)

where $\mathbf{n}(\mathbf{x})$ is the unitary outgoing normal vector at point
$\mathbf{x}$ on $\partial$M. As is standard practice in radiative heating rate
estimation, the integrand on the left-hand side can be replaced by
$-\left(\sigma_{\rm a}(\mathbf{x})+s\right)\tilde{J}(\mathbf{x})$, per radiant
energy conservation expressed in (12).

The Laplace transform of a delta function at the origin is unity. Therefore,
(2) transforms into the (steady-state) BC

$\tilde{I}(s,\mathbf{x},\boldsymbol{\Omega})=F_{0}/\pi,\text{ for
}\mathbf{x}\in\partial{\rm M}\text{ and
}\boldsymbol{\Omega}\cdot\mathbf{n}(\mathbf{x})<0,$ (15)

where $F_{0}$ is a steady and uniform irradiance over $\partial$M expressed in
the usual units of [W/m2]. We thus ensure that
$\tilde{I}(s,\mathbf{x},\boldsymbol{\Omega})$ has the usual units of
[W/m2/sr]. In other words, the so-far unspecified proportionality constant in
(9) is set to $F_{0}/E_{0}$ in [1/s].

The right-hand side of (14) can be used to compute the mean in-cloud
pathlength, spatially and angularly-averaged over $\partial$M and the 2$\pi$
solid angle for radiance contributing to the outgoing flux. Formally,

$\langle{ct}\rangle=\left.\int_{\partial{\rm
M}}\int_{\boldsymbol{\Omega}\cdot\mathbf{n}(\mathbf{x})>0}\langle{ct}\rangle(\mathbf{x},\boldsymbol{\Omega})\times\boldsymbol{\Omega}\cdot\mathbf{n}(\mathbf{x})\,{\rm
d}\boldsymbol{\Omega}\,{\rm d}S(\mathbf{x})\right/(\pi S)$ (16)

where the normalization uses
$\int_{\boldsymbol{\Omega}\cdot\mathbf{n}(\mathbf{x})>0}\boldsymbol{\Omega}\cdot\mathbf{n}(\mathbf{x})\,{\rm
d}\boldsymbol{\Omega}=\pi$, for any $\mathbf{x}\in\partial$M, and
$S=\int_{\partial{\rm M}}{\rm d}S(\mathbf{x})$, which is the surface of the
optical medium M. Indeed, from (1) with $q=1$, but for an outgoing flux-based
PDF, hence normalization by the uniform boundary irradiation $F_{0}$ from
(15), we have

$\displaystyle\langle{ct}\rangle={}$
$\displaystyle\left.-\frac{\partial\;}{\partial s}\left.\int_{\partial{\rm
M}}\tilde{\mathbf{F}}(s,\mathbf{x})\cdot\mathbf{n}(\mathbf{x})\,{\rm
d}S(\mathbf{x})\right|_{s=0}\right/(F_{0}S)$ $\displaystyle={}$
$\displaystyle\left.\frac{\partial\;}{\partial s}\left.\int_{{\rm
M}}\left(\sigma_{\rm a}(\mathbf{x})+s\right)\,\tilde{J}(s,\mathbf{x})\,{\rm
d}V(\mathbf{x})\right|_{s=0}\right/(F_{0}S).$ (17)

The 3D RT problem defined in (10) with BCs (15) can in fact be solved easily
for conservatively scattering media where $\sigma_{\rm a}(\mathbf{x})\equiv 0$
(no absorption by particulates, hence $\sigma_{\rm
s}(\mathbf{x})\equiv\sigma(\mathbf{x})$) and in steady state, i.e., $s=0$
(equivalently, no gaseous absorption, $k_{\lambda}=0$). Indeed, a uniform and
isotropic radiance field $\tilde{I}(0,\mathbf{x},\boldsymbol{\Omega})\equiv
I_{0}=F_{0}/\pi$ from (15), leading to $\tilde{J}(0,\mathbf{x})\equiv
J_{0}=4\pi I_{0}=4F_{0}$ and
$\tilde{\mathbf{F}}(0,\mathbf{x})\equiv\mathbf{0}$ verifies (10), hence (12).
Moreover, this solution is valid not only for media with arbitrary outer
shapes, but also any internal structure as defined by the extinction field
$\sigma(\mathbf{x})$, including an optical void (where
$\sigma(\mathbf{x})\equiv 0$).

In this special case, the numerator in (2) simplifies to $4F_{0}V$, where
$V=\int_{\rm M}{\rm d}V(\mathbf{x})$ is the volume of M
$\subset\mathbb{R}^{3}$. From there, we recover the known universal result
that

$\langle{ct}\rangle=4V/S$ (18)

in 3D, as first proven most elegantly by [3]. In 2D RT, all physical units
loose a per-m and per-sr becomes per-rad; the isotropic BC in (15) becomes
$I_{0}=F_{0}/2$ and $J_{0}=2\pi I_{0}=\pi F_{0}$, hence
$\langle{ct}\rangle=\pi A/P$ where $A$ is the area of M
$\subset\mathbb{R}^{2}$ and $P$ is its perimeter. In _literal_ 1D RT, where
light can only travel up and down the $z$-axis (as in two-stream theory),
physical units all loose another per-m and per-rad goes away; the BC in (15)
becomes $I_{0}=F_{0}$ and $J_{0}=2I_{0}=2F_{0}$, hence $\langle{ct}\rangle=2H$
where $H$ is the length of M $\subset\mathbb{R}$. In short, the mean
pathlength $\langle{ct}\rangle$ is invariably a direct measure of the bulk
size of the medium.

##  3 Moments of $ct$ for Uniform Optically Thick Plane-Parallel Media

###  3.1 Diffusion-Theoretical Model Formulation

To explore the cloud information content of pathlength moments beyond
$\langle{ct}\rangle$, starting with $\langle{(ct)^{2}}\rangle$ or
$\text{Var}[ct]=\langle{(ct)^{2}}\rangle-\langle{ct}\rangle^{2}$, we consider
homogeneous plane-parallel scattering media that are sufficiently opaque for
the RT to be treated in the diffusion limit. In the diffusion regime, the
_exact_ (radiant energy) conservation law in (12) is complemented with Fick’s
_approximate_ constitutive law:

$\tilde{\mathbf{F}}=-\,\frac{1}{3\sigma_{\rm t}}\boldsymbol{\nabla}\tilde{J},$
(19)

where

$\sigma_{\rm t}=(1-g)\sigma_{\rm s}+\sigma_{\rm a}$ (20)

is the scaled or “transport” extinction coefficient, with $g$ being the
scattering PF’s asymmetry factor, $2\pi\int_{-1}^{+1}\mu_{\rm s}P(\mu_{\rm
s}){\rm d}\mu_{\rm s}$. Noting that $\tilde{F}_{x}=\tilde{F}_{y}=0$ in
uniformly illuminated plane-parallel, we only need to track $\tilde{F}_{z}$
(where we will drop the subscript for simplicity).

We have thus reduced the 1D (i.e., homogeneous plane-parallel) RT problem form
(12) and (19) to the solution of two coupled ordinary differential equations
(ODEs),

$\displaystyle{\rm d}\tilde{F}/{\rm d}z=$ $\displaystyle{}-(\sigma_{\rm
a}+s)\tilde{J}(s,z)\text{ and}$ (21) $\displaystyle{\rm d}\tilde{J}/{\rm d}z=$
$\displaystyle{}-3\sigma_{\rm t}\tilde{F}(s,z),$ (22)

subject to BCs

$\displaystyle\left.(\tilde{J}+3\chi\tilde{F})\right|_{z=0}$
$\displaystyle{}=4F_{0}\text{ and}$ (23)
$\displaystyle\left.(\tilde{J}-3\chi\tilde{F})\right|_{z=H}$
$\displaystyle{}=0,$ (24)

where $H$ is the physical thickness of the medium, and $\chi$ is the
“extrapolation” length factor (typically set to 2/3) that controls the mix of
density and current contributions in the above 3${}^{\text{rd}}$-type (a.k.a.
Robin) BCs. We wish to compute the outgoing boundary fluxes

$\displaystyle R(s)=$
$\displaystyle{}\left.\left.(\tilde{J}-3\chi\tilde{F})\right|_{z=0}\right/(4F_{0})=\tilde{J}(s,0)/(2F_{0})-1\text{
and}$ (25) $\displaystyle T(s)=$
$\displaystyle{}\left.\left.(\tilde{J}+3\chi\tilde{F})\right|_{z=H}\right/(4F_{0})=\tilde{J}(s,H)/(2F_{0}),$
(26)

respectively for diffuse reflection and total transmission. We refer to [4]
for derivations. In the remainder, we assume without loss of generality that
$F_{0}$ is unity.

Those authors also conducted detailed separate investigations of $R(s)$ and
$T(s)$ for the conservative scattering case ($\sigma_{\rm a}=0$) in the limit
$s\to 0$ as well as of $\langle{ct}\rangle_{R,T}$ and
$\langle{(ct)^{2}}\rangle_{R,T}$, from the behavior of $R(s)$ and $T(s)$ as
$s\to 0$. Specifically, considering these quantities as observables (obtained
from either active multiple-scattering cloud lidar or passive DO2AS), their
dependence on inherent cloud properties was examined from the standpoint of
remote sensing. Key cloud properties targeted are geometric and (scaled)
optical thicknesses: $H$ and $\tau_{\rm t}=\sigma_{\rm t}H=(1-g)\sigma H$. As
is well-known, $\tau_{\rm t}$ (hence COT $\tau=\sigma H$, for given $g$) can
be derived from $R(0)$ or $T(0)$, respectively from above and below, but
sensitivity deteriorates at large $\tau_{\rm t}$. Specifically,

$R(0)=1-T(0)=\frac{\tau_{\rm t}}{2\chi+\tau_{\rm t}},$ (27)

approaches unity, and $T(0)=1-R(0)$ vanishes, as $\tau_{\rm t}$ increases
without bound. The main new finding is that $\langle{ct}\rangle_{R}$ or
$\langle{ct}\rangle_{T}$ can be used to infer $H$, knowing $\tau_{\rm t}$ from
(27) or other determinations of COT $\tau$ and $g$. At large COT (and moderate
$\tau_{\rm t}$ if $g=0.85$, the canonical value for liquid clouds), $R(0)$
looses its sensitivity. [4] went on to show, however, that large COT is
precisely when $\langle{(ct)^{2}}\rangle_{R}$ becomes highly sensitive to it,
and COT can thus be inferred from this 2${}^{\text{nd}}$-order pathlength
moment, knowing $H$ (e.g., from $\langle{ct}\rangle_{R}$) and $g$. In a sharp
contrast, $\langle{(ct)^{2}}\rangle_{T}$ contains no new cloud information
beyond what can be learned from ground-based sensing of $T$ and
$\langle{ct}\rangle_{T}$.

###  3.2 Mean _and_ Variance of $ct$: Non-Absorbing Case

Here, we are interested in pathlength properties derived from
$R(s)+T(s)=1-A(s)$, where $A(k_{\lambda})$ is the bulk gaseous absorption by
the medium. That way, we are effectively assuming that both top and bottom of
the slab are irradiated with diffuse light.

The first step is to solve in closed-form the coupled ODE boundary-value
problem in (21)–(24). Then, knowing $\tilde{J}(s,z)$, one computes

$R(s)+T(s)=\frac{\tilde{J}(s,0)+\tilde{J}(s,H)}{2}-1,$ (28)

from (25)–(26), leading to:

$[R+T](s)=\frac{\sigma_{\rm t}-3\chi^{2}s+2\chi\sqrt{3s\sigma_{\rm
t}}/\sinh(H\sqrt{3s\sigma_{\rm t}})}{\sigma_{\rm
t}+3\chi^{2}s+2\chi\sqrt{3s\sigma_{\rm t}}/\tanh(H\sqrt{3s\sigma_{\rm t}})}.$
(29)

Figure 1 shows $\log[R+T](s,\sigma_{\rm t})$ for $\chi=2/3.$ We also set $H=1$
without loss of generality since we can then identify $\sigma_{\rm t}$ with
the non-dimensional $\tau_{\rm t}$. Because of its connection to pathlength
moments, we are especially interested in the behavior of
$\log[R+T](s,\sigma_{\rm t})$ for small values of $s$ and arbitrary—but not
too small—values of $\sigma_{\rm t}$ ($=\tau_{\rm t}$), to ensure the validity
of the diffusion approximation in (22).

Figure 1:  $\log[R+T]$ from (29) is plotted as a 2D function of
$(s;\sigma_{\rm t})$ after setting $\chi=2/3$ and $H=1$, which can be done
without loss of generality. We can, in that case, simply read $\sigma_{\rm t}$
as the scaled COT $\tau_{\rm t}$. We note that the slope in $s$ at $s=0$ is
$-$2, independently of the scaled extinction $\sigma_{\rm t}$, that is, the
opacity of the medium. In (30) and (31), we see that this follows directly
from the universal law in (18) for $\langle{ct}\rangle=4V/S$ (= 2$H$ in plane-
parallel media). Curvature of $\log[R+T]$ at $s=0$ does depend on $\sigma_{\rm
t}$, as we will see in (31)–(33).

The next step is to expand $\log[R+T](s)$ in Fig. 1 into a two-term Taylor
series at $s=0$. Based on (29), we find

$\log[R+T](s)=-3\chi\,H\,s+\frac{3}{4}\,\chi\,(\sigma_{\rm
t}H)\,H^{2}\,s^{2}+\mathcal{O}(s^{5/2}).$ (30)

We can identify the terms of that series with those of

$\log[R+T](s)=-\langle{ct}\rangle\,s+\frac{1}{2}\,\text{Var}[ct]\,s^{2}+\mathcal{O}(s^{5/2}),$
(31)

where we used (7) and (8) in the 1${}^{\text{st}}$\- and
2${}^{\text{nd}}$-order terms, respectively. we note that the
0${}^{\text{th}}$-order term vanishes since $[R+T](0)=1$.

From the 1${}^{\text{st}}$ terms in (30) and (31), we see that

$\langle{ct}\rangle=2H$ (32)

for $\chi=2/3$. This is as expected from the powerful prediction by [3] that
$\langle{ct}\rangle=4V/S$ for arbitrary _finite-sized_ media, as re-derived in
the above. Indeed, starting for a finite parallelepiped of dimensions $W\times
W\times H$, its volume is $V=W^{2}H$ and its surface is $S=2W^{2}+4WH$; so the
limit of $4V/S$ as $W\to\infty$ is indeed $2H$, which is also the above
prediction for $\langle{ct}\rangle$ in literal 1D RT theory using the
divergence theorem.

From the 2${}^{\text{nd}}$ terms in (30) and (31), we see that

$\text{Var}[ct]=\frac{3}{2}\,\chi\,(\sigma_{\rm t}H)\,H^{2}=\tau_{\rm
t}\,H^{2},$ (33)

again for $\chi=2/3$. It follows from (33) that
$\langle{(ct)^{2}}\rangle=\text{Var}[ct]+\langle{ct}\rangle^{2}=(\tau_{\rm
t}+4)\,H^{2}\sim\tau_{\rm t}$ when $\tau_{\rm t}\gg 1$. That is indeed the
predicted scaling by [5], $\langle{(ct)^{q}}\rangle\sim\tau_{\rm
t}^{q-1}H^{q}$ (in our notations) for arbitrary optically thick media, in the
$q=2$ case. In A, the scaling of $ct$ moments of arbitrary order for reflected
light is derived from discrete-time random walk theory, thus reinforcing the
robustness of the $ct$-moments/cloud-properties connection. Moreover, their
statement about the prefactor
$\partial\langle{(ct)^{2}}\rangle/\partial\tau_{\rm t}=H^{2}$ agrees with our
result for $\chi$ = 2/3.

We thus confirm that knowledge of $\langle{ct}\rangle$ yields $H$ while
knowledge of Var[$ct$] yields $\tau_{\rm t}=(1-g)\tau$, knowing $H$ (from
$\langle{ct}\rangle$). We also confirm that the sensitivity of Var[$ct$] to
COT $\tau$ is not affected by its value, no matter how large. The conclusion
is the same as obtained by [4] because, at large COTs where the diffusion
approximation works best, $R$ dominates in the sum $[R+T](s)$ at small $s$.
However, the asymptotic ($\tau\to\infty$) partition in the pathlength moments
between $R$ and $T$ contributions varies. Based on individual expressions by
[4], the $R$-to-$T$ is 2-to-1 for $\langle{ct}\rangle$ and 8-to-7 for
$\langle{(ct)^{2}}\rangle$. This rapid trend toward unity makes sense since
the higher the moment order, the more it is dominated by the longest paths
that typically end in transmission, which compensates for the dwindling weight
of the transmitted light, in $2\chi/\tau_{\rm t}$, cf. (27).

Numerical validation of the diffusion-theoretical predictions for
$\langle{ct}\rangle$ and Var[$ct$] in plane-parallel media is in order. Figure
2 shows Monte Carlo simulation results for two choices of scattering PF:
isotropic and [6] (with $g$ = 0.85, as observed in liquid clouds). COT $\tau$
was varied by factors of 2 from 0.5 to 256 for the isotropic ($g$ = 0) PF, and
from 4 to 512 for the Henyey–Greenstein (H–G) PF. As expected, we see that the
$\langle{ct}\rangle/H=2$ prediction is followed exactly to numerical
precision. Also as expected, the Var[$ct$]$/H^{2}=\tau_{\rm t}=(1-g)\tau$
prediction is asymptotically exact, and already very good when $\tau_{\rm t}$
exceeds $\sim$4\. This is consistent with the numerical simulations for plane-
parallel slabs by [5] using a linearized Monte Carlo code to compute
derivatives with respect to extinction, as displayed in their Table 1.

Equations (32) and (33), or their more complicated but more relevant
counterparts derived by [4] for reflected light only, can be used to infer
cloud properties $(H,\tau_{\rm t})$ from the mean and variance of $ct$.
However, it is not obvious how to get these pathlength moments from raw DO2AS
data. We therefore describe in B a tentative algorithm for consideration in
existing or future missions.

Figure 2:  Numerical validation of predictions for $\langle{ct}\rangle$ and
Var[$ct$] based on radiative diffusion theory in plane-parallel cloud models.
The original prediction that $\langle{ct}\rangle=2H$ by [3] for any COT is
verified, as a special case of their general result that
$\langle{ct}\rangle=4V/S$ for arbitrary optical media. Furthermore, our
prediction in (33) that Var[$ct$] = $\tau_{\rm t}\,H^{2}$ is verified in the
asymptotic regime ($\tau_{\rm t}\gg 1$).

###  3.3 Mean of $ct$: Absorbing Case

The physical principles of the retrieval of $(H,\tau_{\rm t})$, a
parameterized cloud profile, from reflection DO2AS are very general. In
essence, the light reflected from an optically thick medium is a combination
of short and long paths. The later dominate the mean pathlength
$\langle{ct}\rangle$ and the associated light has penetrated the whole cloud,
which therefore contains information about the geometrical cloud thickness
$H$. In contrast, pathlength variance Var[$ct$] is determined by the spread
between long and short paths. If the cloud is optically thin, i.e., the
_transport_ mean-free-path (tMFP) between (effectively isotropic) scatterings,
$\ell_{\rm t}=1/\sigma_{\rm t}$ is commensurate with $H$ and there is not
enough space in the medium M to build much diversity in pathlength before the
light escapes the medium, so Var[$ct$] is small. If the cloud is optically
thick, $1/\sigma_{\rm t}\ll H$ then the combination of short and long paths is
more balanced. Consequently, Var[$ct$] is large and will increase as
$\sigma_{\rm t}$ or $H$ increase since these two parameters respectively make
the short paths shorter and the long paths longer. It is therefore no surprise
to find that Var[$ct$] is $\propto\sigma_{\rm t}H=\tau_{\rm t}$.

Cloud particles are non-absorbing at ${{\mathrm{O}}_{2}}$ A- and B-bands, both
in- and out-of-band. Indeed, only the cloud boundaries are absorbing in the RT
problem investigated so far. However, since the physics behind the inference
of $(H,\tau_{\rm t})$ for an optically thick layer from DO2AS observations is
so general, we wonder how much absorption can be tolerated for the retrieval
method to still work. The same question arises for the inference of just $H$
at any optical thickness. We are motivated here by the potential for DO2AS to
be applied to primarily opaque aerosol “clouds” of desert dust, wildfire smoke
or volcano ash near their sources, where the AOT/COT can be very big. Can we
use DO2AS to go beyond aerosol (layer) top height retrieval, and also obtain
its geometrical thickness?

We can begin to address this question theoretically with the same diffusion
model used for water clouds as long as we keep absorption small compared to
scattering, i.e., $0<\sigma_{\rm a}\ll(1-g)\sigma_{\rm s}$ in the make up of
$\sigma_{\rm t}$ in (20). With that caveat in mind, we can return to the model
laid out in (21)–(26). The quantity of interest from that model, $[R+T](s)$,
is obtained in (29) under the assumption that $\sigma_{\rm a}=0$ (no
particulate absorption). By examining (21), we see that $[R+T](s)$ for the
case $\sigma_{\rm a}>0$ is the same expression, but with $s\mapsto
s+\sigma_{\rm a}$. Therefore, to compute the successive pathlength moments for
media with weakly absorbing particles, we only need to expand $\log[R+T](s)$
as it stands in (29) into a Taylor series at $s=\sigma_{\rm a}$.

Truncated at 1st-order, the Taylor series of $\log[R+T](s)$ at $\sigma_{\rm
a}$ yields (29) with $s=\sigma_{\rm a}$ at order 0 and, at 1st-order, an
expression for $\langle{ct}\rangle(\sigma_{\rm a})$ that is too complicated to
be illuminating. Figure 3 (left panel) shows $\langle{ct}\rangle/H$ as a
function of the single-scattering albedo (SSA) $\omega=\sigma_{\rm s}/\sigma$
and COT $\tau=\sigma H$ for $g=0.7$ and $\chi=2/3$. The SSA range is limited
to [0.95,1] and the aerosol COT varies from 0 to 30, although the diffusion
model becomes accurate only past $1/(1-g)$. A typical value of $g$ for
absorbing aerosols is $\approx$0.7, and can be less for fine smoke particles.
The SSA range is limited to [0.95,1] and the aerosol COT varies from 0 to 30,
although the diffusion model becomes accurate only somewhat past
$1/(1-g)\approx 3.3$ (i.e., $\tau_{\rm t}\approx 1$, about
1/10${}^{\text{th}}$ of the $\tau$-range in Fig. 3). The maximum stress for
the diffusion model is thus for $\sigma_{\rm a}/\sigma_{\rm
s}=1/\omega-1\approx 0.05$, which is indeed quite small compared to $1-g=0.3$.
We see in the left panel of Fig. 3 that $\langle{ct}\rangle=2H$ at $\tau=0$
for any $\omega$ as well as for $\omega=1$ for any $\tau$, as shown in the
above. In all other cases, $\langle{ct}\rangle<2H$ since both absorption and
opacity reduce pathlengths. In the left panel of Fig. 3, we have highlighted
the line in $(\tau,\omega)$-space where $\langle{ct}\rangle=H$, i.e., half of
the pathlength is gone on average. The lower the SSA, the quicker this happens
as COT increases.

Figure 3:  Left: Mean pathlength $\langle{ct}\rangle$ in units of layer
thickness $H$ for conservative ($\omega=1$) and absorbing $0.95\leq\omega<1$
aerosol clouds with varying COT, assuming $g=0.7$ and $\chi=2/3$. The
$\langle{ct}\rangle=H$ level (1/2 of the $\omega=1$ case) is plotted for
visual reference. Right: Reflectivity $R(\omega,\tau)$ and transmittivity
$T(\omega,\tau)$ over the same parameter domain.

An important question arises, from the standpoint of aerosol plume or cloud
remote sensing: Can we still recover layer thickness $H$ from a DO2AS-derived
estimate of in-cloud mean pathlength $\langle{ct}\rangle$? Based on the left
panel of Fig. 3, the answer is affirmative, but the algorithm will need to be
informed from other sources about the values of SSA and COT to get the right
$\langle{ct}\rangle/H$ ratio.

That said, we need to be aware that overall light levels decrease rapidly as
the co-SSA $1-\omega$ increases, and that is before adding the gaseous
absorption. The right-hand panel in Fig. 3 quantifies this issue. Therein, we
have have plotted $R(\omega,\tau)$ and $T(\omega,\tau)$ for $g=0.7$ (most
representative of aerosols) and $\chi=2/3$, from expressions in [4]. When
$\omega<1$, $T(\omega,\tau)$ decreases exponentially with $\tau$ at a rate

$L_{\text{dif}}(\omega,g)=\frac{1}{\sqrt{3(1-\omega)(1-\omega g)}},$

the so-called diffusion length scale in units where MFP = $1/\sigma$ is unity.
It rapidly goes to $\infty$ as $\omega$ approaches unity, and
$T(\omega,\tau)$’s downward trend becomes $\propto$$1/\tau$; for $g=0.7$ and
$\omega=0.95$, the diffusion scale is only $\approx$4.5 MFPs,
$\approx$1/6${}^{\text{th}}$ of the whole $\tau$-range in Fig. 3. Meanwhile,
$R(\omega,\tau)$ increases at the same exponential rate, asymptoting at

$R(\omega,\infty)=\frac{1-\chi S(\omega,g)}{1+\chi S(\omega,g)},\text{ where
}S(\omega,g)=\sqrt{\frac{3(1-\omega)}{1-\omega g}}$

is the so-called similarity factor [7]. For $g=0.7$ and $\chi=2/3$,
$R(0.95,\infty)$ is only $\approx$0.38, but the reflectivity of a moderately
opaque aerosol layer would be even less. In DO2AS observations, even more
light is removed by gaseous absorption but, as shown in Appendix B, it is
fortunately the low values of $k_{\lambda}$ that are the most informative, at
least for the extraction of $\langle{ct}\rangle$ and, from there, $H$.

##  Appendix A Scaling of Pathlength Moments from Discrete-Time Random Walk
Theory at Any Order

It is well-known from the theory of _continuous_ -time random walks (RWs),
a.k.a. Brownian motion, that the variance $\langle{\mathbf{x}^{2}}\rangle$ of
the RW position increases linearly with time $t$. Specifically,
$\langle{\mathbf{x}^{2}}\rangle(t)=6Dt$, resulting from the Gaussian solution
of $\partial_{t}U=D\boldsymbol{\nabla}^{2}U$, subject to
$U(0,\mathbf{x})=\delta(\mathbf{x})$, where $D=c\ell_{\rm t}/3$ is the
diffusivity constant in 3D space. Introducing pathlength $ct$, we have
$\langle{\mathbf{x}^{2}}\rangle=(2\ell_{\rm t})ct$. In _discrete_ -time RW
theory, we investigate the random position $\mathbf{x}_{n}$ of the RWing
particle after $n\geq 0$ steps. Since every discrete step corresponds, on
average, to a time-step $\ell_{\rm t}/c$, or a pathlength step $\ell_{\rm t}$,
we have $ct=\ell_{\rm t}n$, hence

$\langle{\mathbf{x}_{n}^{2}}\rangle=(2\ell_{\rm t}^{2})n.$ (A-1)

This results directly from the additivity of means and variances over $n$
independent identically-distributed (iid) random variables. Indeed, the RW is
composed of $n$ statistically-uncorrelated (i.e., randomly-oriented) steps,
$\mathbf{x}_{i}-\mathbf{x}_{i-1}$, each with variance $2\ell_{\rm t}^{2}$,
starting at $\mathbf{x}_{0}=\mathbf{0}$. Thus,
$\mathbf{x}_{n}=\sum_{i=1}^{n}(\mathbf{x}_{i}-\mathbf{x}_{i-1})$, where
$\langle{\mathbf{x}_{n}}\rangle=\langle{\mathbf{x}_{n-1}}\rangle=\dots=\mathbf{x}_{0}=\mathbf{0}$,
and (A-1) follows for
$\langle{\mathbf{x}_{n}^{2}}\rangle=\sum_{i=1}^{n}\langle{(\mathbf{x}_{i}-\mathbf{x}_{i-1})^{2}}\rangle$.

So far, $\mathbf{x}_{n}$ is the random variable and $n$ is the given quantity
and, so far, the RW unfolds in unbounded 3D space. Turning the tables, we ask:
On average, what is the number of steps at which a RW starting at
$\mathbf{x}_{0}=\mathbf{0}$ first reaches the $z=H$ plane? If the RW is then
terminated, it contributes a transmission event in a RW problem over a finite
domain

$\text{M}(H)=\\{\mathbf{x}\in\mathbb{R}^{3};0<z<H\\},$ (A-2)

so we denote the answer as $\langle{n}\rangle_{T}$. In the spirit of
dimensional analysis, we posit that, based on (A-1), this number surely scales
as

$\langle{n}\rangle_{T}\sim(H/\ell_{\rm t})^{2}=\tau_{\rm t}^{2},$ (A-3)

the scale COT from the main text, squared.

How do we now translate this scaling law into information about RWs that end
in reflection at $z=0$? We need to invoke the lesser-known _law of first
returns_ [8], which assumes the RW is on a half-space, i.e., M($\infty$) in
(A-2), and that the first step is in the positive-$z$ direction. Let $n$ be
the number of steps taken by the RWer before it crosses into the $z<0$ half-
space. Then Prob$\\{n\geq N\\}\sim 1/\sqrt{N}$ [9], equivalently,

$p_{n}\sim n^{-3/2},$ (A-4)

since Prob$\\{n\geq N\\}=\sum_{n=N}^{\infty}p_{n}$. Note that both mean and
variance of $n$, let alone higher moments, are divergent due to the frequent
occurrence of RWs that go deep into the positive half-space, even though these
RWs will eventually return to $z=0$ with probability 1. We also note that
(A-4) only required that the steps be isotropic [e.g., 10]; their variance can
be infinite, as in so-called Lévy flights [e.g., 11].

Now, returning to the RWs in the finite slab M($H$), in principle, with
$H\gg\ell_{t}$, we simply truncate the discrete PDF in (A-4) for
$n>\langle{n}\rangle_{T}$ from (A-3). In that case,

$\langle{n^{q}}\rangle_{R}=\left.\sum_{n=1}^{[\langle{n}\rangle_{T}]}n^{q}p_{n}\right/\sum_{n=1}^{[\langle{n}\rangle_{T}]}p_{n},$
(A-5)

where $q\geq 0$ and $[\cdot]$ designates the closest integer value. For
$\langle{n}\rangle_{T}$, that is, $\tau_{\rm t}^{2}$ from (A-3), sufficiently
large, the denominator in (A-5) is close to unity since we assume that the PDF
in (A-4) is normalized. Consequently,

$\langle{n^{q}}\rangle_{R}\approx\sum_{n=1}^{[\tau_{\rm
t}^{2}]}n^{q-3/2}\sim(\tau_{\rm t}^{2})^{q-1/2}=\tau_{\rm t}^{2q-1}.$ (A-6)

From there, we can derive the anticipated scaling for

$\langle{(ct)^{q}}\rangle_{R}\sim\tau_{\rm
t}^{2q-1}\,\ell_{r}t^{q}=H^{q}\times\tau_{\rm t}^{q-1},$ (A-7)

using $\ell_{r}t=H/\tau_{\rm t}$ Equation (A-7) is exactly as predicted by
[5], but in our notations. In particular, we retrieve
$\langle{ct}\rangle_{R}\sim H$ and $\langle{(ct)^{2}}\rangle_{R}\sim
H^{2}\,\tau_{\rm t}$.

##  Appendix B Algorithm for Deriving Low-Order Pathlength Moments from DO2AS
Data

In the presence of multiple scattering, spectroscopic observations with real-
world instruments involve not one but two averages at a miniumum. Even in the
line-by-line limit (i.e., ultra-fine spectral resolution), one must average
over all possible pathlengths $L=ct$, as in (9), from the main text, where
$s=k_{\lambda}$. The underlaying pathlength distribution in (9) is precisely
the target of the study since it contains information about the cloud.

Real spectrometers have finite spectral resolution $\Delta\lambda$, and inside
a spectral bin $[\lambda-\Delta\lambda/2,\lambda+\Delta\lambda/2)$ one can
potentially find many values of the molecular absorption coefficient
$k_{\lambda}$. The instrument is therefore performing a de facto integration
over all wavelengths $\lambda$ in that interval, which is equivalent to an
integration over a whole distribution of $k_{\lambda}$ values. In RT
simulations, this last spectral integration is commonly approximated by a
weighted sum,
$\int_{\lambda-\Delta\lambda/2}^{\lambda+\Delta\lambda/2}\Phi(k_{\lambda}){\rm
d}\lambda\approx\sum_{n}g_{n}(\lambda,\Delta\lambda)\Phi(k_{\lambda})$, where
the number of terms as well as the weights will depend on the diversity of
$k_{\lambda}$ values in $[\lambda-\Delta\lambda/2,\lambda+\Delta\lambda/2)$.

Suppose now that we are given DO2AS ratio data in the form of

$r_{\text{DO2AS}}(\lambda,\Delta\lambda;\cdots)=\frac{1}{\tilde{I}(0,\cdots)}\sum\limits_{n}g_{n}(\lambda,\Delta\lambda)\tilde{I}(k_{\lambda},\cdots),$
(B-1)

for $\lambda_{\text{min}}<\lambda<\lambda_{\text{max}}$, hopefully with
several if not many spectral samples. It will likely be based on some
averaging and/or sampling of the transport space coordinates
$(\mathbf{x},\boldsymbol{\Omega})$, e.g., $\mathbf{x}$ spans an imager pixel
and $\boldsymbol{\Omega}$ is some viewing direction; we have represented this
averaging/sampling here with “$\cdots$” as a placeholder. The spectral signal
in (B-1) fluctuates up and down, starting at the continuum where
$r_{\text{DO2AS}}(\lambda,\Delta\lambda;\cdots)$ is unity because
$k_{\lambda}$ vanishes there by definition. The first order-of-business is to
produce a scatter plot the $r_{\text{DO2AS}}(\lambda)$ data versus the known
values of $k_{\lambda}$. This scatter plot starts a unity (with $k_{\lambda}$
= 0) and trends downward to a minimum associated with the largest value of
$k_{\lambda}$ in the spectral window
$[\lambda_{\text{min}},\lambda_{\text{max}}]$.

As previously stated, the individual spectral channels will in general be wide
enough to contain a possibly large range of $k_{\lambda}$ values in the line-
by-line sense. It is therefore key to assign a meaningful _effective_ value
$k_{\lambda}^{\text{eff}}$ to each spectral channel. There is no unique way of
doing this, but a rational approach is to use the same weights as in (B-1):
$k_{\lambda}^{\text{eff}}(\Delta\lambda)=\sum_{n}g_{n}(\lambda,\Delta\lambda)k_{\lambda}$.
Based on (1) in the main text, the forward model for the DO2AS ratio data in
the form of a truncated Taylor series expansion at $k_{\lambda}=0$ is then

$r_{\text{DO2AS}}(\lambda,\Delta\lambda;\cdots)=1-\langle{L}\rangle\,k_{\lambda}^{\text{eff}}(\Delta\lambda)+\frac{\langle{L^{2}}\rangle}{2}\,\left(k_{\lambda}^{\text{eff}}(\Delta\lambda)\right)^{2}+\dots$
(B-2)

Based on (B-2), the algorithm to extract $\langle{L}\rangle$ and
$\langle{L^{2}}\rangle$ from DO2AS observations is as follows.

  * 1.

Generate scatter plot of
$\left(k_{\lambda}^{\text{eff}}(\Delta\lambda),r_{\text{DO2AS}}(\lambda,\Delta\lambda;\cdots)\right)$
pairs where at least one data point is out-of-band (i.e.,
$k_{\lambda}^{\text{eff}}(\Delta\lambda)$ = 0 and
$r_{\text{DO2AS}}(\lambda,\Delta\lambda;\cdots)$ = 1).

  * 2.

Apply to these data a constrained least-squares polynomial fit of order at
least order 2, the constraint being to go through the (0,1) point.

  * 3.

Identify the 1st-order coefficient with $\langle{L}\rangle$.

  * 4.

Identify the 2nd-order coefficient with $\langle{L^{2}}\rangle$/2.

Depending on the total number of available data points, hopefully more than 3,
one can fit a polynomial of order $>$2 if it helps to reduce the uncertainty
on the first two coefficients. To obtain cloud properties from these two
pathlength moments, we can use (32) and (33) from the main text. However,
their simplicity comes from the simplifying assumption of uniform (top _and_
bottom) illumination, hence a mix of reflected and transmitted light. [4]
offer more involved formulas for the two pathlength moments that are more
realistic since derived only for the observed (namely, reflected) light and,
moreover, account for important parameters such as the solar zenith angle.

## References

  * [1] W. Irvine, The formation of absorption bands and the distribution of photon optical paths in a scattering atmosphere, Bull. Astron. Inst. Neth. 17 (1964) 226–279.
  * [2] V. Ivanov, S. A. Sabashvili, Transfer of resonance radiation and photon random walks, Astrophysics and Space Science 17 (1972) 13–22.
  * [3] S. Blanco, R. Fournier, An invariance property of diffusive random walks, EPL (Europhysics Letters) 61 (2003) 168.
  * [4] A. Davis, I. Polonsky, A. Marshak, Space-time Green functions for diffusive radiation transport, in application to active and passive cloud probing, in: A. Kokhanovsky (Ed.), Light Scattering Reviews, Vol. 4, Springer-Praxis, Heidelberg (Germany), 2009, pp. 169–292.
  * [5] S. Blanco, R. Fournier, Short-path statistics and the diffusion approximation, Phys. Rev. Lett. 97 (2006) 230604.
  * [6] L. Henyey, J. Greenstein, Diffuse radiation in the galaxy, Astrophys. J. 93 (1941) 70–83.
  * [7] M. King, Determination of the scaled optical thickness of clouds from reflected solar radiation measurements, J. Atmos. Sci. 44 (1987) 1734–1751.
  * [8] S. Redner, A Guide to First-Passage Processes, Cambridge University Press, Cambridge, UK, 2001.
  * [9] E. Sparre Anderson, On the fluctuations of sums of random variables, Math. Scand. 1 (1953) 236–285.
  * [10] U. Frisch, H. Frisch, Universality in escape from half space of symmetrical random walks, in: M. Shlesinger, G. Zaslavsky, U. Frisch (Eds.), Lévy Flights and Related Topics in Physics, Springer-Verlag, New York (NY), 1995, pp. 262–268.
  * [11] A. Davis, A. Marshak, Lévy kinetics in slab geometry: Scaling of transmission probability, in: M. Novak, T. Dewey (Eds.), Fractal Frontiers, World Scientific, Singapore, 1997, pp. 63–72.

©2021. All rights reserved.