## Abstract

Here, previous work using photon diffusion theory to describe radiative transfer through dense plane-parallel clouds at nonabsorbing wavelengths is extended. The focus is on the scaling of space- and time-domain moments for transmitted light with respect to cloud thickness *H* and optical depth *τ*; and the new results are as follows: accurate prefactors for asymptotic scaling, preasymptotic correction terms in closed form, 3D effects for internal variability in *τ,* and the rms transit time or pathlength. Mean pathlength is ∝*H* for dimensional reasons and, from random-walk theory, we already know that it is also ∝(1 – *g*)*τ* for large enough *τ* (*g* being the asymmetry factor). Here, it is shown that the prefactor is precisely 1/2 and that corrections are significant for (1 – *g*)*τ* < 10, which includes most actual boundary layer clouds. It is also shown that rms pathlength is not much larger than the mean for transmittance (its prefactor is ≈ 0.59); this proves that, in sharp contrast with reflection, pathlength distributions are quite narrow in transmission. If the light originates from a steady point source on a cloud boundary, a fuzzy spot is observed on the opposite boundary. This problem is formally mapped to the pulsed source problem, and it is shown that the rms radius of this spot slowly approaches *H* as *τ* increases; it is also shown that the transmitted spot shape has a flat top and an exponential tail. Because all preasymptotic corrections are computed here, the diffusion results are accurate when compared to Monte Carlo counterparts for *τ* ≥ 5, whereas the classic scaling relations apply only for *τ* ≥ 70, assuming *g* = 0.85. The temporal quantities shed light on observed absorption properties and optical lightning waveforms. The spatial quantity controls the three-dimensional radiative smoothing process in transmission, which was recently observed in spectral analyses of time series of zenith radiance at 725 nm. Opportunities in ground-based cloud remote sensing using the new developments are described and illustrated with simulations of 3D solar radiative transfer in realistic models of stratocumulus. Finally, since this analytical diffusion study applies only to weakly variable stratus layers, extensions to more complex cloud systems using anomalous diffusion theory are discussed.

## 1. Introduction

There is a growing interest in the pathlengths cumulated by photons while scattering around in clouds. The most popular observational approach is based on oxygen A-band spectroscopy, capitalizing on the equivalence theorem (Irvine 1964; Ivanov and Gutshabash 1974; van de Hulst 1980): variable absorption by a well-mixed gas yields the Laplace transform of the pathlength distribution. The empirical A-band studies of Pfeilsticker et al. (1998) and Min and Harrison (1999) use ground-based observations of transmitted light, as do those of Portman et al. (2001), although they use the *γ*-band. The theoretical results of Stephens and Heidinger (2000) and Heidinger and Stephens (2000, 2002) create tremendous expectation for forthcoming space-based observations from future National Aeronautics and Space Administration (NASA) platforms.

Cloud remote sensing in the solar spectrum is also benefiting from recent pathlength studies and associated spatial considerations such as Marshak et al.'s (1995) concept of radiative smoothing. The radiative smoothing scale is a diagnostic determining the pixel size at which the operational assumption of plane-parallel theory breaks down (Davis et al. 1997a), although not always irretrievably (Marshak et al. 1998). The conservative scattering results of these authors for optical depth retrievals have now been extended by Platnick (2001a,b) to wavelengths where liquid water is absorbing, hence to effective droplet radius retrievals. Davis and Marshak (2001) also address absorbing media in a more physical than remote sensing frame of mind, by varying the single scattering albedo in closed-form analytical results.

Another approach to pathlengths uses pulsed sources, both natural and artificial. Following the historical development of 3D solar radiative transfer, lightning studies have used Monte Carlo (Thomason and Krider 1982) and diffusion (Koshak et al. 1994) methods to characterize the “waveforms” of transmitted pulses. These waveforms are now being routinely measured by spaceborne instruments (Suszcynski et al. 2000) so there is a renewed pressure to refine the theory (Light et al. 2001). New developments in cloud lidar described by Davis et al. (1997b, 1999, 2001a) make direct use of laser pulses returned from dense clouds; the reflected pulse is observed to be vastly broadened by the multiple scattering, also the area from which it comes is commensurate with cloud thickness.

With lidar applications in mind, Davis et al. (1999) derive rigorous diffusion-based results for the mean and rms dwelling times or pathlengths of photons reflected from dense plane-parallel clouds as functions of their physical thickness and optical depth. In the present paper, we complement these results for transmission in view of emerging applications to ground-based A-band spectroscopy and space-based lightning observations. As did Davis et al. (1999), we validate the analytical diffusion theory with high-accuracy numerical results from Monte Carlo simulations.

In section 2, we use space–time moments of the transmitted Green's function to define characteristic radiative scales and survey the literature. In section 3, the Green's function theory is worked out in the diffusion approximation with Fourier and/or Laplace transforms. Closed-form expressions for the moments are derived and numerically validated in section 4. Remote sensing applications are discussed in section 5, while extensions to clouds with more or less internal variability is covered in section 6. In section 7, we summarize and outline future work.

## 2. Moment-based definitions of space- and timescales

### a. The Green's function for transmitted light

In time-dependent radiative transfer, we use the linear transport equation

subject to boundary/initial conditions to determine the time-dependent radiance field *I*(*t, ***r**, **Ω**) at instant *t* ≥ 0, and position **r** ∈ M = {(*x, **y, **z*)^{T} ∈ ℜ^{3} | 0 ≤ *z* ≤ *H*}, that is, in a plane-parallel slab medium and with propagation direction **Ω** ∈ Ξ = {(Ω_{x}, Ω_{y}, Ω_{z})^{T} ∈ ℜ^{3} | Ω^{2}_{x} + Ω^{2}_{y} + Ω^{2}_{z} = 1}; we use here a superscript “T” to denote transpose. Recall that

using standard notations. In the nonstationary radiative transfer equation in (1), we have introduced the following notations: *c* is the speed of light; *σ*(**r**) ≥ 0 is the (generally position-dependent) extinction coefficient; *σ*_{s}(**r**) ≥ 0 is the scattering coefficient; and *p*(**Ω**′ · **Ω**) is the phase function, assumed azimuthally symmetric, where **Ω**′ · **Ω** = cos*θ*_{s} = *μ*_{s} is the cosine of the scattering angle (cf. Fig. 1).

In this study, we use primarily boundary/initial conditions for an isotropic pulsed point source at cloud top:

and an absorbing condition at cloud base:

Note that the factor *π*^{–1} in (3a) is replaced by *δ*(*μ* + 1)/2*π* when modeling a collimated beam under normal incidence, and by *δ*(*μ* + *μ*_{0})*δ*(*ϕ*) for a slant beam with direction **Ω**(–*μ*_{0}, 0), where 0 < *μ*_{0} < 1.

Figure 1 describes the Green's function (GF) problem of interest here for a pulsed boundary source. In this study, we briefly discuss properties of the reflected flux field

prevailing at cloud top, and focus on those of the transmitted flux field

at cloud base.

### b. Background

The steady-state incarnation of these problems, known as “pencil beam” problems, were addressed analytically and numerically in the limit *H* → ∞ by many authors, going at least back to Chandrasekhar (1958). In this case of a semiinfinite medium, there is no transmission per se, so the internal and reflected light fields are of interest. Definitive answers for this problem with isotropic scattering were obtained only quite recently by Ganapol et al. (1994) for *μ*_{0} = 1, followed by Kornreich and Ganapol (1997) for *μ*_{0} < 1, who propose these problems as benchmarks for 3D numerical transport codes, primarily in nuclear engineering applications.

As far as we know, Romanova (1968b) was the first to relax the assumption of isotropic scattering, being motivated by problems in oceanic and atmospheric optics. Assuming *μ*_{0} = 1, Romanova (1968a) allows for media with finite thickness but treats the forward-peaked scattering in the small-angle approximation. Under identical assumptions, Weinman and Shipley (1972) consider the temporal ramifications in transmission for a pulsed source, followed by Weinman (1976) for reflection. Romanova (1971b) generalizes her (1968a) steady-state theory to capture scattering through arbitrary angles and, still investigating finite anisotropically scattering media, Romanova (1971a) reconsiders the possibility of *μ*_{0} < 1. In these last two papers by Romanova, the spatial component of the moment-based characterization of the light fields used below is first introduced.

In the frame of diffusion theory, which we use extensively further on, Richards (1956) studied the case of a steady isoptropic source observed in transmission through a slab (here, the spatial Green's function problem). Finally, Weinman and Masutani (1987) extended this early theoretical work and applied their analysis to the appearance of city lights through clouds in nighttime satellite images.

Temporal aspects of the pulsed/isotropic source problem have been investigated independently by at least three distinct communities. In all cases, space and time are eventually combined in improved diagnostics.

In observational astrophysics, Predehl et al. (2000, and references therein) obtain the distance to the binary Cyg X-3 source using the scattering “halo” due to interstellar dust, as imaged in the X-ray spectrum by the Chandra satellite.

In the medical optics literature, Patterson et al. (1989) solve the pulsed source problem in diffusion theory for a slab, thus starting a new trend in noninvasive diagnostics of soft tissue (Yodh and Chance 1995).

In the atmospheric optics literature, Davis et al. (1997b, 2001a), Davis (1999), as well as Winker (1997) and Miller and Stephens (1999), were motivated by multiple scattering effects in Lidar In-Space Technology Experiment (LITE, Winker et al. 1996) data for clouds and by new instrumental developments in off-beam lidar observation of clouds (Davis et al. 1999; Love et al. 2001).

The former study uses a single scattering estimate of the transmitted field, and an eclipse rather than a pulse. In the latter cases, multiple scattering is essential, and the emphasis shifts from transmission to reflection. As we will see in the following, these two fields on opposite sides of the medium are intimately related when multiple scattering dominates.

### c. Probabilistic interpretation and Monte Carlo results

In the following, we systematically exploit the probabilistic interpretation of the radiance field: given the photon source distribution in space–time, *I*(*t, ***r**, **Ω**) ≥ 0 is the probability of finding a photon in state (**r**, **Ω**) at time *t.* So *t, ***r**, and **Ω** are all random variables in this picture. In this framework, *G*_{T}(*t, **x, **y*) in (4b) describes the probability and space–time variables for photons in transmission (i.e., escape with Ω_{z} < 0). Specifically, we can obtain the following conditional probability measure

or probability density function (PDF)

For reasons of symmetry, we are interested in PDFs dependent only on horizontal photon displacement or “lateral transport” *ρ* = (*x*^{2} + *y*^{2})^{1/2} (cf. Fig. 1). Moreover, it is convenient to use “pathlength” units for time, *λ* = *ct* (cf. Fig. 1). Hence, the element of probability

where *c*^{–1}*ρdρdϕdλ*/*dxdydt* is the Jacobian *D*(*ρ, **ϕ, **λ*)/*D*(*x, **y, **t*) associated with this change of variables. Returning to our definition of flux at the lower boundary in (4b), we can express the PDF in (7) as the “normalized” transmitted Green's function

where

From there, we can define various moments of *λ* and *ρ*:

Figure 2 shows the results from a pathlength-tracking Monte Carlo simulation of radiative transfer in a plane-parallel cloud with uniform extinction; physical thickness is *H* = 0.3 km and (dimensionless) optical thickness is *τ* = *σH* = 16. The scattering is conservative with a Henyey–Greenstein (1941) phase function for *g* = 0.85, hence rescaled optical depth (1 – *g*)*τ* = 2.4. This problem is 3D only because of the nonuniform illumination at a single point; in this case, and because of the isotropic illumination pattern, there is no loss of spatial information during the azimuthal integration in (7). The source is transient, with all the light emitted at once, so the problem is also nonstationary. Figure 2a shows the reflected Green's function while Fig. 2b is for transmission. Space–time integrals of these Green's functions are, respectively, *R* = 0.557 and *T* = 1 – *R* = 0.443.

The close resemblance of these space–time fields at large distances and long times is not accidental. Indeed, since *T* ≈ *R* [as expected when (1 – *g*)*τ* ≈ 2] a significant fraction of the light will reach the middle of the cloud via multiple scattering; from there it has about an equal chance of leaving the cloud in reflection and in transmission at any given *ρ* and *λ.* For denser (more reflecting) clouds, the same is true but the fraction involved in this equipartition decreases, being asymptotically equal to the transmitted light. The only difference between *G*_{R}(*λ, **ρ*) and *G*_{T}(*λ, **ρ*) is at early times and small distances. Note the linear shape of the ballistic (*ρ*∝*λ*) envelope of *G*_{R}(*λ, **ρ*) due to a significant number of contributions from low orders of scattering while the envelope of *G*_{T}(*λ, **ρ*) is parabolic with *ρ* roughly proportional to (*λ* – *H*)^{1/2}.

The symbols in Fig. 3 show numerical results for the quantities in Eqs. (9a)–(9c) as functions of *H* and *τ.* We refer the reader to other numerical (Platnick 2001a,b) and analytical (Davis and Marshak 2001) studies for the effects of absorption. To model the scattering details, we again used the Henyey–Greenstein phase function with asymmetry factors *g* = 0, 0.85. Numerical (Monte Carlo) results are plotted versus rescaled optical depth *τ*_{t} = (1 – *g*)*τ* to show the good collapse one expects in diffusion-dominated transport. Analytical (diffusion theory) results are derived in the next two sections.

## 3. Green's function theory in the diffusion regime

### a. Definitions

Diffusion theory is based on a highly smoothed version of the local radiance field *I*( · , **Ω**). One key quantity in diffusion theory is *radiant energy density* or simply *photon density*

which can also be written as (4*π*/*c*) × mean radiance, or scalar or actinic flux/*c.* The other important quantity is radiant energy flux vector or photon current density

Here, *U*(*t, ***r**) and **F**(*t, ***r**) obey two independent constraints. First, they enter the (exact) law of radiant energy conservation (Case and Zweifel 1967):

where *σ*_{a} = *σ* – *σ*_{s} is the absorption coefficient; Eq. (12) follows from the radiative transfer equation (RTE) in (1) by angular integration using the above definitions. Second, they are related by Fick's law for photon diffusion, operating as a constitutive relation that “closes” the transport problem in Eq. (12),

where *D*(**r**) is (radiative) diffusivity. We know from kinetic theory (e.g., Reif 1965) that

where *l*_{t}(**r**) is the transport mean free path

where ϖ_{0} = *σ*_{s}(**r**)/*σ*(**r**) = 1 – *σ*_{a}(**r**)/*σ*(**r**) ≤ 1 is the single scattering albedo, assumed constant, and *g* = 2*π* ∫ *μ*_{s}*p*(*μ*_{s}) *dμ*_{s} is the asymmetry factor. In essence, *l*_{t} is the “effective” mean free path for isotropic scattering in the following sense. After a single step, the photon propagates a distance 1/*σ* on average, then undergoes (say) a Mie scattering. After a large number of such forward-peaked scatterings, the photon has propagated on average (1 – ϖ_{0}*g*)^{–1} = 1 + ϖ_{0}*g* + (ϖ_{0}*g*)^{2} + (ϖ_{0}*g*)^{3} + · · · times further in the original direction, and has in the process all but “forgotten” this original direction; see Davis and Marshak (1997) for the ϖ_{0} = 1 case, and Platnick (2001b) for the general case.

### b. Green's functions for homogeneous clouds

For the GF problem, we substitute *U* in Eqs. (10) and (12)–(13) with *cG.* For **r**-independent extinction coefficients, Eqs. (12) and (13) can be combined into a standard parabolic partial differential equation (PDE),

The boundary/initial conditions are, however, less standard. Indeed, what needs to be matched at the boundaries are hemispheric fluxes along the vertical axis, which cannot be computed exactly by knowing only *G* and its partial derivatives. Following Case and Zweifel (1967), we introduce the “extrapolation length” *χl*_{t} into the resulting (mixed) boundary conditions

A numerical constant *χ* is essentially a free parameter determined by matching diffusion-based results to detailed numerical computations.

We now seek

Following standard practice in mathematical physics, this PDE problem can be (horizontal) Fourier- and (time) Laplace-transformed into an ordinary differential equation (ODE) where *z* is the only independent variable. So we introduce

where the semicolon separates the variable and the parameters (“constants”) of the resulting 1D problem. The ODE is formally equivalent to that of the homogeneous/steady-state two-stream problem in 1D with an absorption term:

where the single coefficient contains the Fourier–Laplace conjugate variables *k*^{2} = ‖ **k** ‖ ^{2} = *k*^{2}_{x} + *k*^{2}_{y} and *s,* as well as *σ*_{a} and *D*:

In the classic uniform (*k* = 0) and steady (*s* = 0) source problem used in solar studies, we have *L* = (*D*/*cσ*_{a})^{1/2} = (*l*_{t}/3*σ*_{a})^{1/2} = *σ*^{–1}[3(1 – ϖ_{0})(1 – ϖ_{0}*g*)]^{–1/2}, which is known as the “diffusion length.” One can easily see in Eq. (21) why Chandrasekhar (1958) perceived the steady (*s* = 0) pencil-beam problem (the radiative transfer version of a point source in diffusion in Fig. 1) as an interesting combination of all possible absorption problems. Equation (21) also illustrates the &ldquo=uivalence theorem” (e.g., van de Hulst 1980) that shows how one can map time-dependent problems to steady-state absorption problems for any given illumination pattern.

for the ODE in Eq. (20). This defines a two-point boundary value problem for *G̃* that is easily solved (cf. Meador and Weaver 1980) for the generic absorption case 0 ≤ *σ*_{a} < ∞ with *k* = *s* = 0 in Eq. (21). Apart from *z, **G̃* depends on *L*(*s*/*D, **k*^{2}; *cσ*_{a}/*D*) that in turn depends on two local length scales *l*_{t} and 1/*σ*_{a}, and on two other boundary-related length scales, namely, *χl*_{t} and *H.* Explicitly, we have

Fourier–Laplace transforms of the surface flux field in (18a) and (18b) are the formal counterparts of transmittance and reflectance in the associated two-stream problem; using (22a) and (22b), we find

Substituting Eq. (23) here leads to

### c. Cloud modulation transfer and point spread functions

For simplicity, we consider only the spatial Green's function problem by setting *s* = 0 (along with *σ*_{a} = 0) in Eqs. (25a) and (25b) and rewriting them as follows:

In optical engineering, the distribution of intensity at the focal plane of an imaging system for a point source is called the point spread function (PSF). That is the analog of our spatial Green's function, especially of *G*_{T}(*ρ*) for transmission with an assembly of lenses in mind, but also for *G*_{R}(*ρ*) since imaging systems often contain one or more mirrors. Fourier optics is the analog of our Green's function theory so far. Indeed, it is generally easier to compute the 2D Fourier transform of the PSF, otherwise known as the modulation transfer function (MTF). Equations (26a) and (26b) thus represent the cloud's MTFs for reflection and transmission, respectively.

The cloud MTF *G̃*_{T}(0, *k*; · ) in (26b) is only a function of *Hk* (a dimensionless wavenumber) and of the ratio *χ*/*τ*_{t} = *χl*_{t}/*H,* hence of *τ*_{t} = (1 – *g*)*τ* for a given value of the extrapolation factor *χ. *Figure 4 shows normalized MTFs, *G̃*_{T}(0, *k*; · )/*G̃*_{T}(0, 0; · ), plotted versus *Hk* for different cloud opacities (parameterized by *τ*_{t}). Notice the collapse of the curves at *k* = 0 for large values of *τ*_{t}; this translates to the near constancy of 〈*ρ*^{2}〉_{T}/*H*^{2} already shown in Fig. 3 and derived in section 4b [cf. Eq. (35)]. Cloud MTFs and impulse responses, the temporal counterpart of the MTF, for reflection are discussed by Davis et al. (2001b).

To compute the Green's function in physical space, we need the inverse of the Fourier part of (19) expressed for the time-integrated boundary fields, namely,

for subscript *F* = *T, **R.* The standard change of variables for azimuthally symmetric (or, at least, separable) functions leads to the (inverse) Hankel transform, still assuming steady-state illumination (*s* = 0):

where *J*_{0}(*x*) is the 0th-order Bessel function of the first kind.

Figure 5 shows the result of a numerical computation of (28) for *F* = *T* using an efficient random quadrature procedure for the transmitted Green's function in (26b) after normalization with *T.* The cloud parameters were as used in Fig. 2 (*H* = 0.3 km, *τ* = 16, *g* = 0.85), as well as for isotropic scattering (*g* = 0). For comparison, we also plot the Monte Carlo results from the corresponding simulation in Fig. 2 and another one with isotropic scattering (yielding *T* = 0.0765); both cases are highlighted in Fig. 3. Agreement between the diffusion approximation and radiative transfer theory is remarkably good. The far field behavior of *G*_{T}(*ρ*), and also for *G*_{R}(*ρ*) (not illustrated), is exponential with a decay rate that clearly depends on *H* for reasons of dimensionality, but also on *g,* through (1 – *g*)*τ.* The theory for this decay calls for an asymptotic expansion of the integral in (28) with (26b) and is out of the scope of the present study.

## 4. Expressions for moments of transmitted light

### a. Moments and the characteristic function

We now return to the probabilistic interpretation of radiative Green's functions reviewed in section 2c and compute statistical moments of the random variables pathlength *λ* (in lieu of time *t*) and magnitude *ρ* of the lateral transport vector (*x, **y*)^{T}.

Depending on their support, Fourier and/or Laplace transforms of PDFs are known as characteristic functions. Consider the multivariate Taylor expansion of the characteristic function of the axisymmetric PDF in Eq. (8), but unconditioned by *T*:

Coefficients of each term in *s* and/or *k* are related to moments; in particular, for the moments defined in Eqs. (9a)–(9c), we have:

So, by substitution of (25b) into (29), and obtaining the Taylor expansion by successive partial derivation, we obtain these space or time moments as functions of the cloud optical parameters.

### b. Results for conservative diffusion theory

We can now compute the integrals in (9a)–(9c) using the derivatives prescribed in (30a)–(30c) functions of the local quantities *g* and *σ,* as well as the bulk property *H.* In lieu of *σ,* we will continue to use *τ* = *σH.* In addition, we have to carry the boundary condition parameter *χ.* Some (computer assisted) algebra starting with (25b) and using the recipes in (30a)–(30c) leads to the following results.

At 0th-order, we retrieve Schuster's (1905) simple transmission formula from (30a):

Although more accurate formulas exist (cf. Meador and Weaver 1980), this shows that (1 – *g*)*τ* can be retrieved from the observed value of *T* at some nonabsorbing wavelength (e.g., Min and Harrison 1996). This of course calls for a radiometrically calibrated instrument in order to compare the measured flux at ground level with the known incoming flux at the top of the atmosphere.

At first-order (*q* = 1) in (30b), we find

where

with *R* = 1 – *T* denoting reflectance.

The linear term in *k* in (29) vanishes for (26b); this expresses the fact that the lateral displacement vector has a vanishing average, by symmetry. To estimate the size of the spot of transmitted light, we therefore use the quadratic term leading to (30c). Now there is a natural separation of variables in (25b); it is function of *k*^{2} + *s*/*D,* through *L.* So we can use the Jacobian

in (32) to obtain

remembering to account for the 1/2 coefficient of the *k*^{2} term in (29).

Finally, from (30b) at second order (*q* = 2), we obtain

The diffusion results in (31)–(36) were previously plotted in Fig. 3, showing excellent agreement with the radiative transfer numerics for *χ* = 0.7014, the value recommended by Case and Zweifel (1967) and others based on the (half space) Milne problem. However, as noted by Davis et al. (2001a) working on reflection properties, *χ* = 0.57 works better for *T* itself, equivalently 1 – *R.* Note that *τ* and *g* appear always in the combination *τ*_{t} = (1 – *g*)*τ,* as usual in diffusion theory (without the separation of diffuse and direct contributions). In view of the excellent agreement of the PDFs themselves (cf. Fig. 5), agreement in low-order moments is not surprising.

The leading terms in (32) and (35) were obtained by Davis et al. (1997b) using simple scaling arguments, but not the prefactors, nor the important correction terms for the most commonly observed optical depths between ≈5 and 70+. In particular, 〈*λ*〉_{T}∝(1 – *g*)*τH* in (32) follows directly from the well-known fact that the (mean) number of scatterings in transmitted light 〈*n*〉_{T} goes as *τ*^{2} (implicitly for *g* = 0). To get (32) from there, use *τ*_{t} = (1 – *g*)*τ* and recall that *λ* = *n* × mean free path (MFP), with the transport MFP being *l*_{t} = *H*/*τ*_{t}. However, we note in Fig. 2 that this simple relation does not actually occur until (1 – *g*)*τ* ≥ 10, hence *τ* ≥ 70 for *g* = 0.85. So the new correction terms are important.

The result in (36) is new although the leading term is not a surprise since the narrowness of the distribution we find for *λ* in transmission is often assumed, correctly as it turns out, at least for uniform clouds. In particular, Davis et al. (1997a,b) use this assumption in their scaling arguments leading to characteristic times and scales in reflected light. Specifically, we find 〈*λ*^{2}〉^{1/2}_{T}/〈*λ*〉_{T} → (7/5)^{1/2} ≈ 1.18 as *τ* increases, but we see in Fig. 2 that this remains approximately true for (1 – *g*)*τ* ≥ 1. A similar computation as here was done by Davis et al. (1999) for reflected quantities, with lidar applications in mind. In sharp contrast with Eqs. (32) and (36), they find different scalings for 〈*λ*〉_{R}∝*H* (coming from the well-know relation 〈*n*〉_{R}∝*τ* for reflected light), and for 〈*λ*^{2}〉^{1/2}_{R}∝[(1 – *g*)*τ*]^{1/2}*H.* This is traceable to the very broad distribution for *λ* or *n* in reflection from finite but optically thick media (Davis 1999).

## 5. Application to ground-based cloud remote sensing

Although the theory presented here is for homogeneous clouds illuminated at a point, it is closely related to the 3D phenomenon of radiative smoothing of sunlight, that is, the observed deficit of variance at small scales with respect to the overall turbulence-like variability of cloud structure. For instance, cloud optical depth (or liquid water content or liquid water path) has a power-law “*k*^{–5/3}” wavenumber spectrum; so does reflected radiance, but only down to a scale *η*_{R} well above the scaling limit for *τ.* Radiative smoothing was investigated theoretically by Marshak et al. (1995), and then empirically by Davis et al. (1997a) using Landsat data, demonstrating the key role of multiple scattering in 3D. These studies also showed numerically that the leading terms of the counterparts of (32) and (35) for reflectance hold well for plane-parallel stratocumulus models with a fractal internal structure but with somewhat different prefactors, especially when (remotely observable) radiance is used in lieu of cloud-top flux. However, both studies were only in reflection where smoothing occurs on scales ≲〈*ρ*^{2}〉^{1/2}_{R}∝*H*/[(1 – *g*)*τ*]^{1/2}. Recently, Savigny et al. (1999) observed the same phenomenon in ground-based observations of transmitted light—more precisely, zenith radiance—and the characteristic smoothing scale is found empirically to be ≈〈*ρ*^{2}〉^{1/2}_{T}∝*H,* as predicted in Eq. (35). The only operational difference with the reflection studies is that advection by the mean wind is responsible for the “scanning” of cloud base, and its value is needed to convert time to space.

Figure 6 illustrates radiative smoothing in transmission with an ensemble of fractal stratus cloud models using one-dimensional bounded cascades (Cahalan et al. 1994) to generate the random internal structure with variance and correlations that match typical observations of *τ*(*x*). A forward Monte Carlo scheme was used to generate transmitted flux *T*(*x*) and zenith radiance *I*_{zen}(*x*) fields where *x* = 1, … , 1024 measured in 25-m pixels. A simple and convenient way to measure signal smoothness is the first-order structure function

where the overscore means ensemble averaging (i.e., over *x* and over realizations). In theory, the exponent *ζ*_{f}(1) is = 1 for smooth (differentiable) nonstationary signals, <1 for rough (e.g., topography, turbulence) signals, and =0 for stationary signals (such as white noise). From bounded cascade theory (Marshak et al. 1994), optical depth field has *ζ*_{τ}(1) = min {*h,* 1} at all scales (1 ≤ *r* < 1024) where the scaling parameter *h* is set to 1/3. At the smallest scales, we find *ζ*_{T}(1) ≈ *ζ*_{I}(1) ≈ 1 for *r* ≲ *η*_{I} ≤ *η*_{T} ≤ *H* = 0.3 km in this case because lateral (cross pixel) photon diffusion smears out the internal structure that would otherwise be visible. At larger scales, 1D radiative transfer using locally averaged optical depth—the “local” independent pixel approximation (IPA)—describes the 3D radiative transfer well enough, so we find *ζ*_{T}(1) ≈ *ζ*_{I}(1) ≈ *ζ*_{τ}(1) ≈ 1/3, because the IPA is just a one-to-one mapping of *τ* to *T* or *I*_{zen}. Note that the mappings are nonlinear, and *I*_{zen}(*τ*) is not even monotonic across all possible optical depths; but that does not affect the scaling in (37)—only the prefactor, because of the absolute value.

Cloud thickness *H* is hard to obtain reliably by remote sensing methods since milimeter-radars do not always agree with other instruments about cloud base and cloud top (Clothiaux et al. 1995). In the absence of a priori information on cloud thickness, *H* ≈ *η*_{I} ≈ might as well be considered a reasonable first-order estimate of cloud thickness *H* from remote observations since the *O*(1) ratio *η*_{I}/*H* can be determined numerically. Furthermore, the instrumentation used by Savigny et al. (1999) is extremely simple and does not call for calibration.

As mentioned earlier, *T* is used routinely to estimate cloud optical depth *τ,* assuming *g* = 0.85 for dense overcast liquid clouds, but a well-calibrated radiometer is required (Min and Harrison 1996), which is expensive to purchase and to maintain. With Eq. (32), including correction terms, another way of inferring *τ* is now available, knowing *H* and 〈*λ*〉_{T} (as well as *g,* but it varies little in warm clouds). The time domain quantity 〈*λ*〉_{T} can in fact be obtained from oxygen-band spectroscopy at high (Pfeilsticker et al. 1998), intermediate (Portman et al. 2001), or even low (Min and Harrison 1999) resolution; furthermore, there is no need for instrumental calibration since only differential absorption is used.

## 6. Extension to variable clouds and cloud systems

### a. Weak variability: Stratocumulus

In spite of it successes, the analytical diffusion theory presented here for transmission, and by Davis et al. (1999) for reflection, needs refinement. Using multiply scattered lidar data collected in space from a marine stratocumulus (Sc) deck, Davis et al. (2001a) show that it is important to account for stratification of liquid water content, hence of extinction. Furthermore, the diffusion approximation is always improved for realistic (Mie based) phase functions by separating the direct and diffuse components of the radiance field and applying the approximation only to the latter part, using the *δ*-Eddington rescaling (Joseph et al. 1976) of all the relevant optical properties. Finally, one could use the rescaling of local optical properties proposed by Cairns et al. (2000) to account for 3D optical depth variability.

Alternatively, internal variability can be modeled using the parameterization of horizontal optical depth variability with Gamma distributions proposed by Barker et al. (1996), originally for transmission in Eq. (32) and other steady-state quantities with and without absorption. Gamma distributions have just one variability parameter *a* = ( – 1)^{–1} > 0 beyond the mean *τ*:

where Γ(*a*) is Euler's gamma function. The homogeneous model is retrieved in the limit *a* → ∞ and Barker et al. (1996) typically find *a* to be in the range 3–5 for marine Sc. By thus randomizing *τ,* the distribution in (38) could be used to estimate the effects of variability on the new observables in Eqs. (35) and (36), in the same spirit of a “global” IPA (Cahalan et al. 1994). This enhanced analytical approach yields more complicated versions of formulas (32) and (33), and (35) and (36): best derived by computer-assisted algebra, with *τ* replaced by *τ* and the extra parameter *a.*

As an illustration, we use 〈*ρ*^{2}〉_{T} in (35). Its leading term, (2/3)*H*^{2}, is not affected by the variability in *τ* but its correction term is. The new correction term is obtained by averaging *C* over the PDF in (38), using *T*(*τ*) from (31) as a weighting term; in other words, we compute the average of *TC* and divide it by that of *T*. This yields

for *a* > 2 (which is consistent with the weak variability assumption), where Γ(*a*, *z*) = ∫^{∞}_{z }*t*^{a–1}*e*^{–t }*dt* is the incomplete gamma function and

Notice that, in the variable *τ* scenario, ɛ in (33) is a *random variable* while ɛ′ is a *fixed parameter*. As soon as *a* < ∞, the new correction term in (39) is larger than the old one (using ɛ′) in (32) and (35) because of Jensen's (1906) inequality since *C* is a concave function of *τ*.

The ratio *C*^{′} (ɛ, *a*)/*C* (ɛ′) is maximum for ɛ′ = 0, but the correction itself is vanishingly small in this region of parameter space (very large *τ*); the ratio is minimum for ɛ′ = ∞, but it becomes irrelevant as ɛ′ gets much greater than 1 (very small *τ*) anyway. Barker et al.'s (1996) recommendation for Sc being *a* ≈ 4.5, the corresponding ratio of correction terms using (39) and (32) decreases monotonically from ≈1.8 to ≈1.3 as ɛ′ increases from 0 to ∞. This refinement therefore captures the small (≈10%) systematic increase in 〈*ρ*^{2}〉_{T}^{1/2} due to optical-depth variability at intermediate values for the mean (say, *τ* ≈ 10) observed by Marshak et al. (1995, Fig. 6 for 〈*ρ*〉_{T}) and Davis et al. (1997a, Fig. 10 for 〈*ρ*^{2}〉_{T}) when comparing Monte Carlo results for uniform and fractal (e.g., Cahalan et al. 1994) cloud models. Indeed, taking *τ* = 10 (ɛ′ ≈ 1.0) yields *C* ≈ 1.7 and *C*^{′} (ɛ′, *a*) ≈ 2.4 (ratio ≈ 1.4) for *a* = 4.5, hence an increase in 〈*ρ*^{2}〉_{T}^{1/2} of – 1 ≈ 0.12.

The 3D effects we just estimated on 〈*ρ*^{2}〉_{T} are rather small because its asymptotic scaling in (35) does not depend on the variable quantity, *τ*. Turning to 〈*λ*〉* _{T}*, we do not anticipate much more because the asymptotic scaling in (32) is linear in

*τ*, so the first-order variability effects cancel, as long as

*τ*is used instead of

*τ*. Furthermore, the residual effect will be comparable in magnitude and sign with the refinement in (39)–(40) since

*C*

^{(1)}

_{λT}(ɛ) =

*C*(ɛ), as stated in (35).

In the previous section, we outlined a viable exercise in ground-based remote sensing using observables 〈*ρ*^{2}〉^{1/2}_{T} and 〈*λ*〉_{T} to determine cloud properties *H* and *τ* or, rather, the “effective” *τ* for the homogeneous cloud model, which is systematically smaller than *τ* due to IPA-like variability effects. Conceivably, one can use {*T, *〈*ρ*^{2}〉^{1/2}_{T}, 〈*λ*〉_{T}} data to infer cloud parameters {*H, **τ*, *a*}; if accuracy does not allow a sufficiently independent and reliable determination of *a,* at least its value can be prescribed (as for *g* and *χ*) in order to avoid the systematic bias in estimating *τ*.

Looking towards future oxygen A-band observations from space, Heidinger and Stephens (2002) come to similar conclusions about spatial variability. They suggest using joint estimates of cloud albedo *R* and 〈*λ*〉_{R} to determine at once a (less biased) mean optical depth and a measure of its (unresolved) variability.

### b. Strong variability: Broken and/or multiple layers

Both Pfeilsticker (1999) and Min and Harrison (1999) find that the leading term in (32) applies to their oxygen A-band spectroscopy data but only when a single dense cloud layer is present. This means that the homogeneous theory developed here applies reasonably well to dense unbroken clouds such as stratus and stratocumulus. In turn, this is a confirmation of what is already known about photon diffusion dominating the radiative transfer in such clouds from the radiance studies of Davis et al. (1997a) in reflectance, those of Min and Harrison (1996) and of Savigny et al. (1999) in transmittance, and going back to those of King et al. (1990) using in situ radiometry.

In more complex 3D situations with broken clouds and/or multiple layers, Pfeilsticker (1999) finds (in different notations) that an alternative scaling applies:

where 1 ≤ *α* ≤ 2 is the Lévy index and is related to the degree of internal and external variability of the clouds. In principle, one should be able to predict *α* from cloud structure but, so far, the model is used only as a diagnostic (*α* is inferred from radiative observations). Note that the leading term in (32) is retrieved in the limit *α* → 2. At any rate, *H* is now considered to be the total thickness of the cloud system, possibly encompassing the entire troposphere. The opposite limit *α* → 1 corresponds to almost clear-sky situations (the prefactor is then 1/*μ*_{0}) and *H* is a scale height; the left-hand side of (41) is then precisely the cumulated airmass. These findings are confirmed by Min et al.'s (2001) investigation of the joint statistics of 〈*λ*〉_{T} and *τ* at the U.S. Department of Energy's (DOE) climate observatory in Oklahoma, established as part of the Atmospheric Radiation Measurement (ARM) program.

The new scaling for 〈*λ*〉_{T} in (41) as was originally proposed by Davis and Marshak (1997) when applying Lévy flight theory (e.g., Mandelbrot 1982) to photon transport. Their stated goal was actually to explain the observed albedo reduction (transmission enhancement) by 3D effects, and they indeed showed that *T*∝[(1 – *g*)*τ*]^{–α/2}, which is minimal for *α* = 2 at given (1 – *g*)*τ. *Davis et al. (2000) recently uncovered another form of empirical validation of the “anomalous” scaling in (41) using lightning waveforms captured in space by the *Fast On-Orbit Transient Experiment* (*FORTÉ*) satellite, a DOE nonproliferation satellite technology mission, and analyzed by Suszcynski et al. (2000).

The basic result that 〈*ρ*^{2}〉^{1/2}_{T} = *H* probably carries over to anomalous diffusion because it translates the simple fact that a spherical “wave” of diffusing photons is being intercepted by an absorbing plane surface at distance *H* from the source. However, the fate of the relation between 〈*λ*^{2}〉^{1/2}_{F} and 〈*λ*〉_{F} for *F* = *R, **T* as well as that of all the correction terms in this generalization are open questions. An analytical formulation of anomalous diffusion problems in finite slabs was recently obtained by Buldyrev et al. (2001) for an entirely different application; it is expected to bring specific answers to these questions. For the time being, the Lévy flight theory outlined above operates as a diagnostic for shortwave radiative transfer in the whole atmospheric column that quantities the overall spatial variability in cloudiness.

From the standpoint of ground-based optical remote sensing, broken clouds pose a real challenge. To the best of our knowledge, only a couple of schemes are currently under consideration. Marshak et al. (2000) proposed the Normalized Difference Cloud Index (NDCI) method, recently enhanced by Barker and Marshak (2001), which uses radiative Green's function theory in conjunction with the strong difference in surface reflectance on either side of the “chlorophyll edge” at the visible/near-IR transition. Davis (2002) proposes another method based on photon diffusion theory when both sun-lit and shaded sides of a cloud are in view.

## 7. Summary and outlook

We have applied photon diffusion theory to the characterization of radiative transfer through dense clouds in both space and time. Low-order moments of radiative Green's functions are used to define the scales of interest in terms of cloud thickness and optical depth. The closed-form results include preasymptotic corrections and are accurate over the full range of observed optical depths using Monte Carlo results as a benchmark. The present treatment refines and gives credence to the heuristic arguments used previously to partially explain recent path length measurements from ground-based oxygen A-band spectroscopy (Pfeilsticker et al. 1998; Min and Harrison 1999; Min et al. 2001). The recent observations of radiative smoothing in transmission by Savigny et al. (1999) are also explained more rigorously.

We have emphasized that simultaneous observation of the spatial property in Eq. (35) and the temporal quantity in Eq. (32) can be used to infer the optical depth and physical thickness of the cloud layer from ground. This only calls for an oxygen A-band spectrometric instrument using zenith radiance. Indeed, such a device naturally provides the required pathlength, and a straightforward correlation analysis of the time series of any continuum value provides the smoothing scale (Savigny et al. 1999), provided it is sampled frequently enough during the spectroscopic integration time. Moreover, neither of the observables in Eqs. (32) and (35) call for calibration. If one adds a (necessarily calibrated) measurement of diffusely transmitted flux, then, with refined theory à la Barker (1996) containing one variability parameter in Eq. (39), we may also be able to assess the degree of variability inside the cloud layer; if not, the refined theory with a prescribed variability parameter will still be an improvement. Finally, if the variability is too strong for standard diffusion theory, even refined, then there is ample empirical evidence that anomalous diffusion theory is a viable alternative, although many of the details remain to be worked out.

Transmission Green's functions for dense clouds have other applications worth mentioning in characterization of opaque media for known sources as well as in characterization of unknown sources embedded in clouds, fogs, plumes, etc., where sources are either steady or time-dependent. In the former type of problem, we discussed only passive (solar source) techniques. Active approaches are also possible (bistatic lidar) and already used, for instance, in medical diagnostics (Patterson et al. 1989). A noteworthy example of the later type of problem is optical detection of lightning from space currently being pioneered by NASA's *Tropical Rain Monitoring Mission* (*TRMM*) and the DOE's *FORTÉ* satellite. Because of the need to search a vast parameter space for a best fit to data, both kinds of problem will eventually call for fast yet accurate solutions for forward 3D radiative transfer. A promising candidate here is adjoint perturbation theory (Box et al. 1988) where the base case can be taken as homogeneous plane-parallel and the critical 3D information is contained in the Green's function for the same base case. In this context, the Green's function source models the detector position and angular acceptance; boundary source Green's functions for transmission (introduced here) and for reflection (already quite extensively used elsewhere) are needed, preferably in the analytical form used here. Green's functions for internal sources are needed for higher-order perturbation terms, a more challenging problem we will be addressing in the near future.

## Acknowledgments

This work was supported in part by the Environmental Sciences Division of U.S. Department of Energy (under Grant DE-A105-90ER61069 to NASA's Goddard Space Flight Center for A.M.), as part of the Atmospheric Radiation Measurement (ARM) program. The authors thank H. Barker, R. Cahalan, P. Chylek, H. Frisch, Y. Knyazikhin, D. Kornreich, T. Light, S. Love, L. Oreopoulos, K. Pfeilsticker, I. Polonsky, C. v. Savigny, G. Stephens, D. Suszcynski, D. Winker, J. Weinman, and W. Wiscombe for many fruitful discussions.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

_{1}theory.

**,**

**,**

**,**

**,**

**,**

_{2}A-band.

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

## Footnotes

*Corresponding author address:* Anthony B. Davis, LANL/NIS-2, P.O. Box 1663 (MS C-323), Los Alamos, NM 87545. Email: adavis@lanl.gov