A Generalized Linear Transport Model for Spatially Correlated Stochastic Media

Abstract

We formulate a new model for transport in stochastic media with long-range spatial correlations where exponential attenuation (controlling the propagation part of the transport) becomes power law. Direct transmission over optical distance τ(s), for fixed physical distance s, thus becomes (1 + τ(s)/a)^{− a}, with standard exponential decay recovered when a → ∞. Atmospheric turbulence phenomenology for fluctuating optical properties rationalizes this switch. Foundational equations for this generalized transport model are stated in integral form for d = 1, 2, 3 spatial dimensions. A deterministic numerical solution is developed in d = 1 using Markov Chain formalism, verified with Monte Carlo, and used to investigate internal radiation fields. Standard two-stream theory, where diffusion is exact, is recovered when a = ∞. Differential diffusion equations are not presently known when a < ∞, nor is the integro-differential form of the generalized transport equation. Monte Carlo simulations are performed in d = 2, as a model for transport on random surfaces, to explore scaling behavior of transmittance T when transport optical thickness τ_t ≫ 1. Random walk theory correctly predicts T∝τ^{− min {1, a/2}}_t in the absence of absorption. Finally, single scattering theory in d = 3 highlights the model’s violation of angular reciprocity when a < ∞, a desirable property at least in atmospheric applications. This violation is traced back to a key trait of generalized transport theory, namely, that we must distinguish more carefully between two kinds of propagation: one that ends in a virtual or actual detection and the other in a transition from one position to another in the medium.

Keywords:

• linear transport theory
• radiative transfer
• stochastic optical media
• turbulence
• clouds
• multiple scattering
• Markov chain formalism
• Monte Carlo
• propagation kernel
• angular reciprocity
• non-exponential extinction laws

APPENDIX A: MARKOV CHAIN FORMALISM FOR GENERALIZED RADIATIVE TRANSFER

For “literal” 1D RT (“d = 1” in the main text), the plane-parallel medium becomes a finite line segment [0,τ*] and each layer is bounded by two adjacent “nodes” on such a line. Accordingly, the phase function describes light scattering in only two possible directions: forward (with probability p₊) and backward (with probability p_– = 1–p₊). These probabilities are conveniently parameterized in Table 1 as p_± = (1±g)/2 where g is the mean cosine of the scattering angle θ_s, which is either +1 (θ_s = 0) or –1 (θ_s = π). Thus, g = (+1)×p₊+(–1)×p_– is between –1 and +1, and g = 0 (isotropic scattering) leads to p₊ = p_– = 1/2.

To construct a computational Markov chain model, we discretize the medium into N “layers” (subsegments) so that each one has an optical thickness Δτ_n = τ_n − τ_n−1 (1 ≤ n ≤ N) with τ₀ = 0 and τ_N = τ*. Moreover, we set the positive τ-axis downward (or rightward, as in Figure 3); that way, the solar source is at the upper (or left, as in Figure 3) boundary at τ = τ₀ = 0.

Markov chain formalism, applied to computational RT in particular, describes a spatio-directional distribution of “particles” executing random walks with no memory beyond the present state. It has two key ingredients: (i) the initial distribution across all possible “states” (position on grid and direction of propagation), and (ii) the transition matrix, which describes the probability of a particle to jump from any given state to any other.

A.1. Initial Light Distribution

The initial amount of light in the nth layer moving in direction i = 0 (downward) or i = 1 (upward) is created by a single scattering of the solar light that directly propagates from the upper boundary to the layer without suffering a collision. Assuming a uniform source of diffuse light in the layer, the first scattering redistributes (according to Bernoulli trial probabilities p_±) the solar light stopped in the layer through scattering.

Adopting conventional notation from Markov chain theory, we thus have the following fluxes in each direction (i = 0,1) per unit of optical distance: (A.1) where F₀ is the incident flux at τ = 0 (which we can assume is unity without loss of generality), ω is the single scattering albedo, p(Δθ_i0) is the phase function, with Δθ_i0 denoting the scattering angle formed between the incidence angle θ₀ = 0 and the internal source angle θ_i: Δθ_i0 = |θ_i − θ₀|. For downward propagation (i = 0), θ_i = 0 and Δθ_i0 = 0 resulting in p(Δθ_i0) = p₊, while for upward propagation (i = 1), θ_i = π and Δθ_i0 = π resulting in p(Δθ_i0) = p_–. Lastly, T_a(τ) denotes here the generalized transmission law in Equation (34(34) ) of the main text, assuming unitary mean extinction ( = 1), while –(τ) > 0 is the absolute value of its derivative, denoted in the main text. Recall that T_a(τ) = –(τ) = exp(–τ) when a = ∞.

A.2. Transition Matrix

The transition matrix Q (continuing to follow conventional notation from Markov chain theory) describes a particle's change from state (n,i), meaning in layer n and going in direction i, into another state (n′,j). Imagine light traveling from the nth to n′th layer, extinguished in the n′th layer in its original direction i, and is then scattered into the new direction j. Each element Q_{(n′,j),(n,i)} of Q quantifies the probability of such a change: (A.2) where Δθ_ji = |θ_j − θ_i| is the scattering angle and the average escaping probability of a particle leaving nth layer is (A.3)

Note in (A.2(A.2) ) the finite probability in the discrete world of the particle remaining in the same layer. Equation (A.2(A.2) ) is evaluated analytically. For n > n′: (A.4) where (A.5) including –exp[–(τ–x)] when a = ∞. For n < n′: (A.6)

Finally, for n = n′: (A.7) where F(0, 0) = 0 for a = 1, and −a/(a−1) for a ≠ 1 (–1 when a = ∞).

A.3. Markov Chain Model

With the initial light distribution vector Π₀ and transition matrix Q, multiple light scattering processes in the spatially correlated stochastic medium with a general (exponential or not) transmission law or propagation kernel can be expressed in the form of a matrix series, namely, (A.8) where E is the identity matrix, and ₀, Q₀, QQ₀, … represent the contributions from first, second, third, and higher orders of scattering, respectively.

The total diffusely reflected (R_dif) and transmitted (T_dif) light are contributed by different layers, namely, (A.9) (A.10) at the top and bottom of medium, respectively. For consistency with the assumption of uniform source distribution in every layer, an average transmission is used for the particles leaving the nth layer for a given location x in the medium, namely, (A.11)

Invoking Equation (A.5(A.5) ), the expressions become (A.12)

Finally, to compute the diffuse intensity field at the top and bottom of the nth layer, (A.9(A.9) ) and (A.10(A.10) ) should be generalized to (A.13) (A.14)

Sums J and differences F of these quantities, as defined from (45(45) ) and (46(46) ), are plotted in Figure 3 for a variety of parameters {ω,g,a} when τ* = 10.

A.4. Energy Conservation

We have the following energy conservation when ω = 1 (no absorption), after including the directly transmitted light T_a(τ*): (A.15)

This identity was used as a first verification test for the Markov chain code. Other tests are described in §4.3 of the main text.

A.5. Relation to Quantities and Methods Used in Main Text

The way the formulation of the literal 1D RT problem is derived from the first principles of Markov chain theory is interesting because we arrive immediately at the Neumann series solution of a discrete version of the ancillary integral form of the RTE for d = 1, namely, Equation (53(53) ) in the main text.

Indeed, the large-but-finite state vector _tot is nothing more than a discrete version of the source function S_±(τ) as can be seen already by comparing (A.1(A.1) ) for its initial value ₀ and the source term Q_S±(τ) in Equation (54(54) ), paying attention to definitions rather than to notations. Consequently, the large system of coupled linear equations that is solved in (A.8(A.8) ), namely, _tot = Q_tot + ₀, is nothing more than a discrete version of the generalized ancillary integral RTE in Equation (53(53) ), solved formally in Equation (55(55) ). The transition matrix Q is therefore just a discretized version of the kernel K_S(τ’,±’;τ,±) for that integral equation, written out explicitly in Equation (53(53) ). Finally, boundary-leaving and internal radiances are obtained from the known _tot with (A.9(A.9) )–(A.10(A.10) ) and (A.13(A.13) )–(A.14(A.14) ) respectively. These are simply the discrete-space counterparts of the formal solution of the 1D RTE, as used in generalized RT to compute I_±(τ) from S_±(τ), once it is a known quantity; its expression for d = 1 can be seen in Equation (56(56) ).

Before closing, it is important to recall that the Monte Carlo method in linear transport/RT theory is put on a solid mathematical footing using the same Markov chain concepts as used here (Marchuk et al., 1980). Finally, we recall that the key assumption in our generalized RT theory is to modify all the spatial integral expressions without requiring that there be a differential formulation from which they are normally derived.

Notes

Following many others, we borrow here the terminology of time-series analysis since the proper language of statistical “homogeneity” might be confused with structural homogeneity, a usage we have already introduced.

²This limit is to be understood physically as to the scale where noise-like fluctuations occur, which is at least the inter-particle distance in a cloud but could be larger (Davis et al., 1999; Gerber et al., 2001).