Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting layered media

Knut Stamnes; S-Chee Tsay; Warren Wiscombe; Kolf Jayaweera

doi:10.1364/AO.27.002502

I. Introduction

The discrete ordinate method for radiative transfer is commonly ascribed to Chandrasekhar.[1] Computer implementations of that method were, however, plagued by numerical difficulties[2] to such an extent that researchers made little use of it. The purpose of this paper is to alert the community to a new numerical implementation of the discrete ordinate method for vertically inhomogeneous layered media which is free of these difficulties and to give a summary of its equations and its various advanced features. The resulting computer code represents the culmination of years of effort[3]–[9] on the part of ourselves and our collaborators to make it the finest algorithm available. Our intent is that the code be so well documented, so versatile, and so error-free that other researchers can easily and safely use it both in data analysis and as a component of large models.

The problem to be solved is the transfer of monochromatic radiation in a scattering, absorbing, and emitting plane-parallel medium with a specified bidirectional reflectivity at the lower boundary. Section II summarizes the equations and boundary conditions. Section III discusses the numerical implementation of the theory, in particular, (a) how to compute eigenvalues and eigenvectors reliably and efficiently and (b) how to avoid fatal overflows and ill-conditioning in the matrix inversion needed to determine the constants of integration. Reference [9] provides a more detailed account as well as documentation, test problems, and a listing of the fortran code.

II. Theory

A. Basic Equations

The purpose of this section is to present the basic radiative transfer formulas with a minimum of definition and explanation to establish notation and conventions. (For more comprehensive discussions of the subject, the reader is referred to textbooks[1],[10] and review articles.[11]–[15])

The equation describing the transfer of monochromatic radiation at frequency ν through a plane-parallel medium is given by[1]

μ \frac{d u_{ν} (τ_{ν}, μ, ϕ)}{d τ} = u_{ν} (τ_{ν}, μ, ϕ) - S_{ν} (τ_{ν}, μ, ϕ),

where u_ν(τ_ν,μ,ϕ) is the specific intensity along direction μ,ϕ at optical depth τ_ν measured perpendicular to the surface of the medium (ϕ is the azimuthal angle, and μ is the cosine of the polar angle). S_ν is the source function given by

S_{ν} (τ_{ν}, μ, ϕ) = \frac{ω_{ν} (τ_{ν})}{4 π} \int_{0}^{2 π} d ϕ^{'} \int_{- 1}^{1} d μ^{'} P_{ν} (τ_{ν}, μ, ϕ; μ^{'}, ϕ^{'}) \times u_{ν} (τ_{ν}, μ^{'}, ϕ^{'}) + Q_{ν} (τ_{ν}, μ, ϕ),

where ω_ν(π_ν) is the single-scattering albedo, and P_ν(τ_ν,μ,ϕ;μ′,ϕ′) is the phase function. For thermal emission in local thermodynamic equilibrium (LTE), the source term Q_ν is

Q_{ν}^{(thermal)} (τ_{ν}) = [1 - ω_{ν} (τ_{ν})] B_{ν} [T (τ_{ν})],

where B_ν(T) is the Planck function at frequency ν and temperature T. If the usual diffuse–direct distinction is made (Ref. [1], p. 22), so that u_ν in Eqs. (1) and (2) describes the diffuse radiation only, then for a parallel beam incident in direction μ₀,ϕ₀ on a nonemitting medium,

Q_{ν}^{(beam)} (τ_{ν}, μ, ϕ) = \frac{ω_{ν} (τ_{ν}) I_{0}}{4 π} P_{ν} (τ_{ν}, μ, ϕ; - μ_{0}, ϕ_{0}) exp (- τ_{ν} / μ_{0}),

where μ₀I₀ is the incident flux. We take

Q_{ν} (τ_{ν}, μ, ϕ) = Q_{ν}^{(thermal)} (τ_{ν}) + Q_{ν}^{(beam)} (τ_{ν}, μ, ϕ) .

(Many radiative transfer models consider one or the other source term but not both.)

To simplify the notation we omit the ν subscript in the following. By expanding the phase function P(τ, cosθ) in a series of 2N Legendre polynomials and the intensity in a Fourier cosine series,[1],[9]

u (τ, μ, ϕ) = \sum_{m = 0}^{2 N - 1} u^{m} (τ, μ) cos m (ϕ_{0} - ϕ),

we find that Eqs. (1) and (2) are replaced by 2N independent equations (one for each Fourier component)

μ \frac{d u^{m} (τ, μ)}{d τ} = u^{m} (τ, μ) - \int_{- 1}^{1} D^{m} (τ, μ, μ^{'}) u^{m} (τ, μ^{'}) d μ^{'} - Q^{m} (τ, μ) (m = 0, 1, 2, \dots, 2 N - 1),

where

D^{m} (τ, μ, μ^{'}) = \frac{ω (τ)}{2} \sum_{l = m}^{2 N - 1} (2 l + 1) g_{l}^{m} (τ) P_{l}^{m} (μ) P_{l}^{m} (μ^{'})

Q^{m} (τ, μ) = X_{0}^{m} (τ, μ) exp (- τ / μ_{0}) + δ_{m 0} Q^{(thermal)} (τ),

X_{0}^{m} (τ, μ) = \frac{ω (τ) I_{0}}{4 π} (2 - δ_{m 0}) \sum_{l = 0}^{2 N - 1} {(- 1)}^{l + m} (2 l + 1) \times g_{l}^{m} (τ) P_{l}^{m} (μ) P_{l}^{m} (μ_{0}),

\begin{array}{l} δ_{m 0} = 1 if m = 0 (0 otherwise), g_{l}^{m} (τ) = g_{l} (τ) \frac{(l - m)!}{(l + m)!}, \\ g_{l} (τ) = \frac{1}{2} \int_{- 1}^{1} P_{l} (cos θ) P (τ, cos θ) d cos θ . \end{array}

P_l(cosθ) is the Legendre polynomial, $P_{l}^{m} (μ)$ is the associated Legendre polynomial, and θ is the angle between the direction vectors before and after scattering. Solutions to Eq. (6) give the azimuthal components, and then Eq. (5) gives the complete azimuthal dependence of the intensity. [When there is no beam source, the sum in Eq. (5) reduces to the m = 0 term, and μ₀, ϕ₀ are irrelevant.]

B. Discrete Ordinate Approximation—Matrix Formulation

The discrete ordinate approximation to (6) can be written[3],[4]

μ_{i} \frac{d u^{m} (τ, μ_{i})}{d τ} = u^{m} (τ, μ_{i}) - \sum_{\begin{array}{l} j = - N \\ j \neq 0 \end{array}}^{N} w_{j} D^{m} (τ, μ_{i}, μ_{j}) \times u^{m} (τ, μ_{j}) - Q^{m} (τ, μ_{i}) (i = \pm 1, \dots \pm N),

where μ_i and w_i are quadrature points and weights.

Since the single-scattering albedo ω(τ) and the phase function P(τ,μ,ϕ;μ′,ϕ′) are functions of τ in a vertically inhomogeneous medium, Eqs. (7) constitute (for each m) a system of 2N coupled differential equations with nonconstant coefficients for which analytic solutions do not exist. To obtain analytic solutions, the medium is assumed to consist of L adjacent homogeneous layers in which the single-scattering albedo and the phase function are taken to be constant within each layer (but allowed to vary from layer to layer, see Fig. 1), and the thermal source term is approximated by a polynomial in τ. For the time being, we shall consider a single homogeneous layer for which τ_p₋₁ ≤ τ ≤ τ_p. The τ arguments of D and X will be omitted since D and X are by assumption independent of τ in any one layer. We will also omit the p subscripts, which are implicit on all quantities, until later.

For any layer in Fig. 1, then, we write Eq. (7) in matrix form as

| \begin{matrix} \frac{d u^{+}}{d τ} \\ \frac{d u^{-}}{d τ} \end{matrix} | = | \begin{array}{r} - α & - β \\ β & α \end{array} | | \begin{matrix} u^{+} \\ u^{-} \end{matrix} | + | \begin{matrix} {\tilde{Q}}^{+} \\ {\tilde{Q}}^{-} \end{matrix} |,

where[3],[16]

\begin{array}{l} u^{\pm} = [u^{m} (τ, \pm μ_{i})] i = 1, \dots, N, \\ {\tilde{Q}}^{\pm} = M^{- 1} Q^{\pm}, \\ Q^{\pm} = Q^{m} (τ, \pm μ_{i}) i = 1, \dots, N, \\ M = {μ_{i} δ_{i j}} i, j = 1, \dots, N, \\ α = M^{- 1} (D^{+} W - I), \\ β = M^{- 1} D^{-} W, \\ W = (w_{i} δ_{i j}) i, j = 1, \dots, N, \\ D^{+} = [D^{m} (μ_{i}, μ_{j})] = [D^{m} (- μ_{i}, - μ_{j})] i, j = 1, \dots, N, \\ D^{-} = [D^{m} (- μ_{i}, μ_{j})] = [D^{m} (μ_{i}, - μ_{j})] i, j = 1, \dots, N . \end{array}

C. Quadrature Rule

The use of Gaussian quadrature is essential because it makes phase function renormalization unnecessary, i.e.,

\sum_{\begin{array}{l} j = - N \\ j \neq 0 \end{array}}^{N} w_{j} D^{0} (τ, μ_{i}, μ_{j}) = \sum_{\begin{array}{l} i = - N \\ i \neq 0 \end{array}}^{N} w_{i} D^{0} (τ, μ_{i}, μ_{j}) = ω (τ),

implying that energy is conserved in the computation.[17]

The quadrature points and weights of the double-Gauss scheme adopted here satisfy μ₋_j = −μ_j and w₋_j = w_j. Double-Gauss simply refers to a quadrature rule suggested by Sykes[18] in which the Gaussian formula is applied separately to the half-ranges −1 < μ < 0 and 0 < μ < 1. The main advantage of this double-Gauss scheme is that the quadrature points (in even orders) are distributed symmetrically around |μ| = 0.5 and clustered both toward |μ| = 1 and μ = 0, whereas in the Gaussian scheme for the complete range, −1 < μ < 1, they are clustered toward μ = ±1. The clustering toward μ = 0 will give superior results near the boundaries where the intensity varies rapidly around μ = 0. A half-range scheme is also preferred since the intensity is discontinuous at the boundaries. Another advantage is that upward and downward fluxes and average intensities are obtained immediately without further approximations.

D. Basic Solution

Equation (7b) is a system of 2N-coupled ordinary differential equations with constant coefficients. These coupled equations are linear, and our goal is to uncouple them by using well-known methods of linear algebra.

Seeking solutions to the homogeneous version ( $\tilde{Q}$ = O) of Eq. (7b) of the form

u^{\pm} = G^{\pm} exp (- k τ),

we find

| \begin{array}{r} α & β \\ - β & - α \end{array} | | \begin{matrix} G^{+} \\ G^{-} \end{matrix} | = k | \begin{matrix} G^{+} \\ G^{-} \end{matrix} |,

which is a standard algebraic eigenvalue problem of order 2N × 2N determining the eigenvalues k and the eigenvectors G^±.

Because of the special structure of the matrix in Eq. (8b), the eigenvalues occur in pairs (±k), and the order of the algebraic eigenvalue problem [Eq. (8b)] may be reduced by a factor of 2 as follows.[3] Breaking Eq. (8b) into

\begin{array}{l} α G^{+} + β G^{-} = k G^{+}, \\ β G^{+} + α G^{-} = - k G^{-}, \end{array}

and then adding and subtracing these two equations, we find

(α - β) (G^{+} - G^{-}) = k (G^{+} + G^{-}),

(α + β) (G^{+} + G^{-}) = k (G^{+} - G^{-}) .

Combining Eqs. (8c) and (8d), we obtain

(α - β) (α + β) (G^{+} + G^{-}) = k^{2} (G^{+} + G^{-}),

which completes the reduction of the order. Following Stamnes and Swanson[3] we solve Eq. (8e) to obtain eigenvalues and eigenvectors (G⁺ + G⁻). We then use Eq. (8d) to determine (G⁺ − G⁻).

For beam sources Q^(beam)(τ,μ) = X₀(μ) exp(−τ/μ₀) [Eq. (6d)], it is easily verified that a particular solution of Eq. (7) is (omitting m superscript)

u (τ, μ_{i}) = Z_{0} (μ_{i}) exp (- τ / μ_{0}) .

For thermal sources the emitted radiation is isotropic, so that Q^m = 0, m > 0, Q(τ) ≡ Q⁰(τ) = (1 − ω)B(τ). By approximating the Planck function for each layer by a polynomial in τ,[3],[19] i.e.,

B [T (τ)] = \sum_{l = 0}^{K} b_{l} τ^{l},

we find that the particular solution may also be expressed as a polynomial in τ:

u (τ, μ_{i}) = \sum_{l = 0}^{K} Y_{l} (μ_{i}) τ^{l} .

The coefficients Z₀(μ_i) and Y_l(μ_i) are readily determined by solving systems of linear algebraic equations (see Refs. [3] and [9] for details).

The general solution to Eq. (7a) consists of a linear combination, with coefficients C_j, of all the homogeneous solutions plus the particular solutions for beam and thermal emission sources (omitting the m superscript):

u (τ, μ_{i}) = \sum_{j = - N}^{N} C_{j} G_{j} (μ_{i}) exp (- k_{j} τ) + δ_{m 0} \sum_{l = 0}^{N} Y_{l} (μ_{i}) τ^{l} .

The k_j and the G_j(μ_i) are the eigenvalues and eigenvectors obtained as described in Sec. III.A, the μ_i the quadrature angles, and the C_j the constants of integration. We have for convenience included the beam solution in the sum by defining for j = 0

C_{0} G_{0} (μ_{i}) \equiv Z_{0} (μ_{i}),

k_{0} \equiv 1 / μ_{0} .

E. Boundary Conditions

Below we shall limit ourselves to a linear-in-optical-depth dependence of the Planck function within each layer, since this will give sufficient accuracy for most cases of practical interest. We shall discuss the inclusion of boundary conditions in some detail, since a proper formulation of this aspect of the problem (in terms of the bidirectional reflectivity) is rarely discussed in the open literature.

In general, Eq. (1) must be solved subject to the boundary conditions

u (τ = 0, - μ, ϕ) = u_{\infty} (μ, ϕ),

u (τ = τ_{L}, + μ, ϕ) = u_{g} (μ, ϕ),

where u_∞ and u_g are the intensities incident at the top and bottom boundaries, respectively, and τ_L is the total optical depth.

Quite general boundary conditions can be accommodated in the discrete ordinate method. To be specific, we shall assume that the medium is illuminated from above by known diffuse radiation (incident parallel beams are treated as pseudosources) and that the bottom boundary has a known bidirectional reflectivity. Equation (11b) may now be rewritten as

u (τ_{L}, + μ, ϕ) = ∊ (μ) B (T_{g}) + \frac{1}{π} \int_{0}^{2 π} d ϕ^{'} \int_{0}^{1} ρ_{d} (μ, ϕ; - μ^{'}, ϕ^{'}) u (τ_{L}, - μ^{'}, ϕ^{'}) μ^{'} d μ^{'} + \frac{μ_{0}}{π} I_{0} exp (- τ_{L} / μ_{0}) ρ_{d} (μ, ϕ; - μ_{0}, ϕ_{0}),

where ∊(μ) is the directional emissivity and T_g is the temperature of the bottom boundary; ρ_d(μ,ϕ;−μ′,ϕ′) is the bidirectional reflectivity, and I₀ is the incident beam intensity at the upper boundary. Kirchhoff’s law

∊ (μ) + \frac{1}{π} \int_{0}^{2 π} d ϕ^{'} \int_{0}^{1} ρ_{d} (μ, ϕ; - μ^{'}, ϕ^{'}) μ^{'} d μ^{'} = 1

is applied to calculate the emissivity from the reflectivity. Note that some authors include the (1/π) factor in ρ_d; we find our definition more useful, however, since when ρ_d is independent of angle it reduces to an albedo.

Below we will assume that the bidirectional reflectivity is a function only of the angle θ between the incident and the reflected radiation (i.e., there are no preferred directions at the lower boundary). Then it can be expanded in a series of 2N Legendre polynomials:

ρ_{d} (μ, ϕ; - μ^{'}, ϕ^{'}) = \sum_{l = 0}^{2 N - 1} (2 l + 1) h_{l} P_{l} (cos θ),

where

cos θ = - μ μ^{'} + {(1 - μ^{2})}^{1 / 2} {(1 - μ^{2})}^{1 / 2} cos (ϕ^{'} - ϕ) .

Using the addition theorem for spherical harmonics, we find

ρ_{d} (μ, ϕ; - μ^{'}, ϕ^{'}) = \sum_{m = 0}^{2 N - 1} ρ_{d}^{m} (μ, - μ^{'}) cos m (ϕ^{'} - ϕ),

ρ_{d}^{m} (μ, - μ^{'}) = (2 - δ_{0 m}) \sum_{l = m}^{2 N - 1} (2 l + 1) h_{l}^{m} P_{l}^{m} (- μ^{'}) P_{l}^{m} (μ),

where

\begin{array}{l} h_{l}^{m} = h_{l} \frac{(l - m)!}{(l + m)!}, \\ h_{l} = \frac{1}{2} \int_{- 1}^{1} P_{l} (cos θ) ρ_{d} (cos θ) d cos θ . \end{array}

Substituting Eq. (13c) into Eq. (12) and using Eq. (5), we find that each Fourier component must satisfy the bottom boundary condition

u^{m} (τ_{L}, + μ) = u_{g}^{m} (μ),

where

u_{g}^{m} (μ) = δ_{0 m} ∊ (μ) B (T_{g}) + (1 + δ_{0 m}) \int_{0}^{1} μ^{'} ρ_{d}^{m} (μ, - μ^{'}) u^{m} (τ_{L}, - μ^{'}) d μ^{'} + \frac{μ_{0}}{π} I_{0} exp (- τ_{L} / μ_{0}) ρ_{d}^{m} (μ, - μ_{0}) .

In a multilayered medium we must also require the intensity to be continuous across layer interfaces.[7] Thus Eq. (7) must satisfy boundary and continuity conditions as follows:

u_{1}^{m} (0, - μ_{i}) = u_{\infty}^{m} (- μ_{i}), i = 1, \dots, N,

\begin{matrix} u_{p}^{m} (τ_{p}, μ_{i}) = u_{p + 1}^{m} (τ_{p}, μ_{i}), i = \pm 1, \dots, \pm N; \\ p = 1, \dots, L - 1; \end{matrix}

u_{L}^{m} (τ_{L}, + μ_{i}) = u_{g}^{m} (μ_{i}), i = 1, \dots, N,

where

u_{g}^{m} (μ_{i}) = δ_{0 m} ∊ (μ_{i}) B (T_{g}) + (1 + δ_{0 m}) \sum_{j = 1}^{N} w_{j} μ_{j} ρ_{d}^{m} (μ_{i}, - μ_{j}) u^{m} (τ_{L}, - μ_{j}) + \frac{μ_{0}}{π} I_{0} exp (- τ_{L} / μ_{0}) ρ_{d}^{m} (μ_{i}, - μ_{0}) .

Evaluation of the integral in Eq. (16b) by numerical quadrature as in Eq. (17d) requires a quadrature rule which integrates separately over the downward and upward directions. As noted previously, the double-Gauss rule adopted here satisfies this requirement.

When discussing boundary conditions, it is convenient to write the discrete ordinate solution in the following form (k_jp > 0 and k₋_jp = −k_jp):

u_{p} (τ, μ_{i}) = \sum_{j = 1}^{N} [C_{j p} G_{j p} (μ_{i}) exp (- k_{j p} τ) + C_{- j p} G_{- j p} (μ_{i}) exp (+ k_{j p} τ)] + R_{p} (τ, μ_{i}),

where the sum contains the homogeneous solution involving the unknown coefficients (the C_jp) to be determined, and R_p(τ,μ_i) is the particular solution given by [Eq. (10a)]

R_{p} (τ, μ_{i}) = Z_{0} (μ_{i}) exp (- τ / μ_{0}) + δ_{m 0} [Y_{0} (μ_{i}) + Y_{1} τ] .

Insertion of Eq. (18) into Eqs. (17a)–(17c) yields (omitting the m superscript)

\sum_{j = 1}^{N} [C_{j 1} G_{j 1} (- μ_{i}) + C_{- j 1} G_{- j 1} (- μ_{i})] = u_{\infty} (- μ_{i}) - R_{p} (0, - μ_{i}), i = 1, \dots, N

\begin{array}{l} \sum_{j = 1}^{N} {C_{j p} G_{j p} (μ_{i}) exp (- k_{j p} τ_{p}) \\ + C_{- j p} G_{- j p} (μ_{i}) exp (k_{j p} τ_{p}) \\ - [C_{j, p + 1} G_{j, p + 1} (μ_{i}) exp (- k_{j, p + 1} τ_{p}) \\ + C_{- j, p + 1} G_{- j, p + 1} (μ_{i}) exp (k_{j, p + 1} τ_{p})]} \\ = [R_{p + 1} (τ_{p}, μ_{i}) - R_{p} (τ_{p}, μ_{i})], i = \pm 1, \dots \pm N; p = 1, \dots, L - 1 \end{array}

\sum_{j = 1}^{N} [C_{j L} r_{j L} (μ_{i}) exp (- k_{j L} τ_{L}) + C_{- j L} r_{- j L} (μ_{i}) exp (k_{j L} τ_{L})] = γ (τ_{L}, μ_{i}), i = 1, \dots, N,

where

γ (τ_{L}, μ_{i}) = δ_{0 m} ∊ (μ_{i}) B (T_{g}) - R_{L} (τ_{L}, + μ_{i}) + \frac{μ_{0}}{π} I_{0} exp (- τ_{L} / μ_{0}) ρ_{d} (μ_{i}, μ_{0}) + (1 + δ_{0 m}) \sum_{j = 1}^{N} ρ_{d} (μ_{i}, - μ_{j}) w_{j} μ_{j} R_{L} (τ_{L}, - μ_{j}),

r_{j L} (μ_{i}) = G_{j L} (μ_{i}) - (1 + δ_{0 m}) \times \sum_{n = 1}^{N} ρ_{d} (μ_{i}, - μ_{n}) w_{n} μ_{n} G_{j L} (- μ_{n}) .

Equations (20a)–(20c) constitute a (2N × L) × (2N × L) system of linear algebraic equations from which the 2N × L unknown coefficients, the C_jp (j = ±1,…,±N; p = 1,…,L) are determined. To solve this system of equations, we take advantage of the fact that the coefficient matrix is a (6N − 1) diagonal band matrix.[7]

III. Numerical Implementation

A. Computation of Eigenvalues and Eigenvectors

Stamnes and Swanson[3] obtained the eigenvalues and eigenvectors of Eq. (8e) using subroutines in eispack[20] for a general real asymmetric matrix. Since such matrices may have complex eigenvalues/eigenvectors, the eispack routines use complex arithmetic. Yet invariably the eigenvalues returned were purely real as they should be.

As discussed by Stamnes et al.,[8] the advantage of solving the eigenvalue problem involving the asymmetric matrix (α − β) (α + β) in Eq. (8e) is that only one matrix multiplication is necessary. Nakajima and Tanaka[16] showed that the eigenvalue problem can be reduced to finding eigenvalues and eigenvectors of the product of two symmetric matrices, one of which is positive definite. Stamnes et al.[8] used the Cholesky decomposition of this positive definite matrix to convert the eigenvalue problem into one involving only one symmetric matrix. But the transformation to symmetric matrices (whose eigensolutions can be found faster and by using real arithmetic) and subsequent solution procedures[8],[16] introduce matrix operations in which the effect of rounding error can seriously degrade the results. Stamnes et al.[8] compared the original Stamnes/Swanson method with the new procedure based on the Cholesky decomposition and with Nakajima and Tanaka’s symmetrizing procedure. They summarize the outcome of these comparisons as follows:

In single precision the Stamnes/Swanson procedure is more accurate than the Cholesky procedure, which in turn is slightly more accurate than the Nakajima/Tanaka procedure.
The Cholesky procedure is slightly faster than Nakajima/Tanaka’s, which is slightly faster than Stamnes/Swanson’s.

Based on the findings, they decided to speed up the Stamnes/Swanson procedure since it is the most accurate. By eliminating complex arithmetic in the eispack subroutines, the speed became comparable with that of the Cholesky approach. We have adopted this faster algorithm in our code.

B. Scaling Transformation

As discussed by Stamnes and Conklin,[7] to avoid numerical ill-conditioning, it is necessary to remove the positive exponentials in Eq. (20) (remember k_jp > 0 by convention). This is achieved by the scaling transformation[7]

C_{+ j p} = {\tilde{C}}_{+ j p} exp (k_{j p} τ_{p - 1}),

C_{- j p} = {\tilde{C}}_{- j p} exp (- k_{j p} τ_{p}) .

Inserting Eqs. (22) in Eqs. (20) and solving for the ${\tilde{C}}_{j p}$ instead of the C_jp, we find that all the exponential terms in the coefficient matrix have the form exp[−k_jp(τ_p − τ_p₋₁)];i.e., all the exponentials in the coefficient matrix have negative arguments (k_jp > 0, τ_p > τ_p₋₁). Consequently, numerical ill-conditioning is avoided implying that our numerical scheme for finding the ${\tilde{C}}_{j p}$ is unconditionally stable for arbitrary total optical depths and arbitrary individual layer thicknesses. [To solve Eqs. (20) for the ${\tilde{C}}_{j p}$ we use the linpack[21] routine sgbsl designed to solve a banded system of linear algebraic equations.]

C. Scaled Solutions—Angular Distributions

Since only the homogeneous solution is affected, we start by introducing the scaling transformation (22) into Eq. (18), i.e.,

u_{p} (τ, μ_{i}) = \sum_{j = 1}^{N} {{\tilde{C}}_{j p} G_{j p} (μ_{i}) exp [- k_{j p} (τ - τ_{p - 1})] + {\tilde{C}}_{- j p} G_{- j p} (μ_{i}) exp [- k_{j p} (τ_{p} - τ)]} .

Since k_jp > 0 and τ_p₋₁ ≤ τ ≤ τ_p, all exponentials in Eq. (23) have negative arguments as they should to avoid fatal overflows.

For a slab of thickness τ_L, we may solve Eq. (1) formally to obtain (μ > 0)

u (τ, + μ) = u (τ_{L}, + μ) exp [- (τ_{L} - τ) μ] + \int_{τ}^{τ_{L}} S (t, + μ) exp [- (t - τ) / μ] \frac{d t}{μ},

u (τ, - μ) = u (0, - μ) exp (- τ / μ) + \int_{0}^{τ} S (t, - μ) exp [- (τ - t) / μ] \frac{d t}{μ},

where we have again omitted the m superscript. These two equations show that if we know the source function S(t,±μ), we can find the intensity at arbitrary angles by integrating the source function. The discrete ordinate solutions are used to derive explicit expressions for the source function, which can be integrated analytically.[3]–[5],[9] This procedure is sometimes referred to as the iteration of the source-function technique, although as pointed out by Kourganoff[22] and discussed in some detail by Stamnes,[5] this procedure essentially amounts to an interpolation.

In a multilayered medium we may now integrate Eq. (24a) layer by layer[5] to obtain the following expression for the upward intensity at arbitrary τ and μ:

\begin{array}{l} u_{p} (τ, + μ) = u (τ_{L}, + μ) exp [- (τ_{L} - τ) / μ] \\ + \sum_{n = p}^{L} {\sum_{j = 1}^{N} [{\tilde{C}}_{- j n} \frac{G_{- j n} (+ μ)}{1 + k_{j n} μ} E_{- j n} (τ, + μ) \\ + \sum {\tilde{C}}_{+ j n} \frac{G_{+ j n} (+ μ)}{1 + k_{j n} μ} E_{+ j n} (τ, + μ)]] \\ + \frac{G_{0 n}}{1 + μ / μ_{0}} E_{0 n} (τ, + μ) \\ + δ_{m 0} [V_{0 n} (μ) F_{0 n} (τ, μ) + b_{1 n} F_{1 n} (τ, μ)]}, \end{array}

where the scaling given by Eqs. (22) has been incorporated and

E_{- j n} (τ, + μ) = exp [- (k_{j n} Δ τ_{n} - δ τ / μ)] - exp [- (τ_{n} - τ) / μ],

F_{0 n} (τ, + μ) = exp (- δ τ / μ) - exp [- (τ_{n} - τ) / μ],

F_{1 n} (τ, + μ) = (τ_{n - 1} + μ) exp (- δ τ / μ) - (τ_{n} + μ) exp [- (τ_{n} - τ) / μ],

with

\begin{array}{l} δ τ = τ_{n - 1} - τ for n > p and δ τ = 0 for n = p; \\ Δ τ_{n} = τ_{n} - τ_{n - 1} for n > p and Δ τ_{p} = τ_{p} - τ for n = p, \\ E_{+ j n} (τ, + μ) = exp [- (τ_{n - 1} - τ) / μ] - exp {- [k_{j n} (τ_{n} - τ_{n - 1}) - (τ - τ_{n}) / μ]}, \end{array}

for n > p and

E_{+ j p} (τ, + μ) = exp [- k_{j p} (τ - τ_{p - 1})] - exp {- [k_{j p} (τ_{p} - τ_{p - 1}) + (τ_{p} - τ) / μ]} .

Similarly, we find that the downward intensity becomes

\begin{array}{l} u_{p} (τ, - μ) = u (0, - μ) exp (- τ / μ) \\ + \sum_{n = 1}^{p} {\sum_{j = 1}^{N} [{\tilde{C}}_{- j n} \frac{G_{- j n} (- μ)}{1 + k_{j n} μ} E_{- j n} (τ, - μ) \\ + {\tilde{C}}_{+ j n} \frac{G_{+ j n} (- μ)}{1 - k_{j n} μ} E_{+ j n} (τ, - μ)] \\ + \frac{G_{0 n} (- μ s)}{1 - μ / μ_{0}} E_{0 n} (τ, - μ) \\ + δ_{m 0} [V_{0 n} (μ) F_{0 n} (τ, - μ) + b_{1 n} F_{1 n} (τ, - μ)]}, \end{array}

where

E_{+ j n} (τ, - μ) = exp [- (k_{j n} Δ τ_{n} + δ τ / μ)] - exp [- (τ - τ_{n - 1}) / μ]

F_{0 n} (τ, - μ) = exp (- δ τ / μ) - exp [- (τ - τ_{n - 1}) / μ]

F_{1 n} (τ, - μ) = (τ_{n} - μ) exp (- δ τ / μ) - (τ_{n - 1} - μ) exp [- (τ - τ_{n - 1}) / μ],

with

\begin{array}{l} δ τ = τ - τ_{n} for n < p and δ τ = 0 for n = p, \\ Δ τ_{n} = τ_{n} - τ_{n - 1} for n < p and Δ τ = τ_{p - 1} for n = p, \\ E_{- j n} (τ, - μ) = exp [- (τ - τ_{n}) / μ] - exp {- [k_{j n} (τ_{n} - τ_{n - 1}) + (τ - τ_{n - 1}) / μ]}, \end{array}

for n < p and

E_{- j p} (τ, - μ) = exp [- k_{j p} (τ_{p} - τ)] - exp {- [k_{j p} (τ_{p} - τ_{n - 1}) + (τ - τ_{p - 1}) / μ]} .

We see that all exponentials in Eqs. (25b)–(25f) and (26b)–(26f) have negative arguments as they should. This ensures that fatal overflow errors are avoided in the computations.

The G_jn terms and ${\bar{V}}_{0 n} (μ)$ in Eqs. (25a) and (26a) are given by the following explicit expressions (omitting the n subscript):

G_{j} (μ) = \sum_{\begin{array}{l} i = - N \\ i \neq 0 \end{array}}^{N} w_{i} D (μ, μ_{i}) G_{j} (μ_{i}) for j \neq 0;

G_{0} (μ) \equiv Z_{0} (μ) = \sum_{\begin{array}{l} i = - N \\ i \neq 0 \end{array}}^{N} w_{i} D (μ, μ_{i}) Z_{0} (μ_{i}) + X_{0} (μ);

V_{0} (μ) = \sum_{\begin{array}{l} i = - N \\ i \neq 0 \end{array}}^{N} w_{i} D^{0} (μ, μ_{i}) Y_{0} (μ_{i}) + (1 - ω) b_{0};

where b₀ is defined in Eq. (9b). Formulas (27a)–(27c) clearly reveal the interpolatory nature of Eqs. (25) and (26). Although we have generally omitted the m superscript above, we have explicitly written D⁰ in Eq. (27c) to indicate that thermal radiation contributes only to the azimuth-independent part of the intensity.

D. Computational Shortcut for Highly Absorbing Media

For highly absorbing media without thermal emission, very little radiation will reach levels for which the absorbing optical depth is greater than, say, 10. In such cases we ignore the medium below this depth and set the surface reflection to zero. This may result in enormous computation savings. For example, in an atmosphere with clouds it prevents unnecessary scattering computations at wavelengths where sunlight never reaches the clouds.

E. Simplified Albedo and Transmissivity Computations

By appealing to the reciprocity principle we may derive simple expressions for the plane albedo a(μ) and the transmissivity t(μ) of a vertically inhomogeneous plane-parallel medium lacking thermal sources. Thus

the plane albedo for a given angle of parallel incidence equals the azimuthally averaged reflected intensity u⁰(0,μ) for isotropic illumination of unit intensity incident at the top boundary;
the transmissivity equals the azimuthally averaged transmitted intensity u⁰*(0,μ) for isotropic illumination of unit intensity incident at the bottom boundary.

Consequently, for an inhomogeneous medium the following duality relations apply[6]:

a (μ) = u^{0} (0, μ); a^{*} (μ) = u^{0 *} (τ_{L}, - μ);

t (μ) = u^{0 *} (0, μ); t^{*} (μ) = u^{0} (τ_{L}, - μ);

where τ_L is to the optical thickness of the inhomogeneous medium, and the asterisk refers to illumination from below.

Equations (28) and (29) offer substantial computational advantages when only albedo and transmissivity are required. The conventional procedure requires computation of the angular distribution of the diffuse intensity, which must then be integrated over angle to obtain the desired result. This necessitates the computation of a particular solution, which can be quite costly in an inhomogeneous (multilayered) medium, since such a solution is required for each layer and every direction of incidence considered. In contrast, by applying an isotropic boundary condition and using Eqs. (28) and (29), one may solve for all desired angles of parallel beam incidence simultaneously and avoid entirely the computation of a particular solution. (Details about the numerical implementation of this procedure are provided in Refs. [6] and [9].)

IV. Summary

We have described an advanced discrete ordinate algorithm for radiative transfer calculations, including multiple scattering, in vertically inhomogeneous nonisothermal plane-parallel media. This algorithm is the product of over seven years of development by ourselves and our collaborators. It has been designed to be as versatile and general as possible and should find applications from the UV through to the radar region of the electromagnetic spectrum.

The algorithm has the following important features:

It is unconditionally stable for an arbitrarily large number of quadrature angles (streams) and arbitrarily large optical depths.
It can be forced by any combination of parallel beam or diffuse incidence and thermal sources in the medium and at the boundaries.
It allows for an arbitrary bidirectional reflectivity at the lower boundary; directional emissivity is computed from this reflectivity.
It offers rapid computation of bulk albedo and transmissivity when thermal sources are absent.
Unlike the popular doubling method, computing time for individual layers is independent of optical thickness.
Because the solution is analytic, and by using the iteration-of-the-source-function method, intensities can be returned at arbitrary angles and optical depths, unrelated to the computational meshes for these quantities.
It has been thoroughly tested against a wide variety of published solutions.
It is thoroughly documented both in Ref. [9] and in the code itself (with extensive references to published equation numbers).

Copies of the fortran program are available from the third author on microcomputer diskettes (IBM or MacIntosh).

This work was supported by the National Science Foundation through grant DPP 84-06093 to the University of Alaska. One of us (K.S.) would like to acknowledge partial support of this work from the Royal Norwegian Council on Science and Technology.

Figure

Fig. 1 Schematic illustration of a multilayered medium.

Download Full Size | PDF

References

1. S. Chandrasekhar, Radiative Transfer (Dover, New York, 1960).

2. K.-N. Liou, “A Numerical Experiment on Chandrasekhar’s Discrete-Ordinate Method for Radiative Transfer: Applications to Cloudy and Hazy Atmosphere,” J. Atmos. Sci. 30, 1303 (1973). [CrossRef]

3. K. Stamnes and R. A. Swanson, “A New Look at the Discrete Ordinate Method for Radiative Transfer Calculations in Anisotropically Scattering Atmospheres,” J. Atmos. Sci. 38, 387 (1981). [CrossRef]

4. K. Stamnes and H. Dale, “A New Look at the Discrete Ordinate Method for Radiative Transfer Calculation in Anisotropically Scattering Atmospheres, II: Intensity Computations,” J. Atmos. Sci. 38, 2696 (1981). [CrossRef]

5. K. Stamnes, “On the Computation of Angular Distributions of Radiation in Planetary Atmospheres,” J. Quant. Spectrosc. Radiat. Transfer 28, 47 (1982a). [CrossRef]

6. K. Stamnes, “Reflection and Transmission by a Vertically Inhomogeneous Planetary Atmosphere,” Planet. Space Sci. 30, 727 (1982b). [CrossRef]

7. K. Stamnes and P. Conklin, “A New Multi-Layer Discrete Ordinate Approach to Radiative Transfer in Vertically Inhomogeneous Atmospheres,” J. Quant. Spectrosc. Radiat. Transfer 31, 273 (1984). [CrossRef]

8. K. Stamnes, S.-C. Tsay, and T. Nakajima, “Computation of Eigenvaluesand Eigenvectors for Discrete Ordinate and Matrix Operator Method Radiative Transfer,” J. Quant. Spectrosc. Radiat. Transfer in press (1988). [CrossRef]

9. K. Stamnes, S.-C. Tsay, W. J. Wiscombe, and K. Jayaweera, “An Improved, Numerically Stable Computer Code for Discrete-Ordinate-Method Radiative Transfer in Scattering and Emitting Layered Meda,” NASA Report, in press (1988).

10. K. N. Liou, An Introduction to Atmospheric Radiation (Academic, Orlando, FL, 1980).

11. J. E. Hansen and L. D. Travis, “Light Scattering in Planetary Atmospheres,” Space Sci. Rev. 16, 527 (1974). [CrossRef]

12. W. M. Irvine, “Multiple Scattering in Planetary Atmospheres,” Icarus 25, 175 (1975). [CrossRef]

13. J. Lenoble, Ed., Radiative Transfer in Scattering and Absorbing Atmospheres: Standard Computational Procedures (A. Deepak, Hampton, VA, 1985).

14. W. J. Wiscombe, “Atmospheric Radiation: 1975–1983,” Rev. Geophys. 21, 957 (1983). [CrossRef]

15. K. Stamnes, “The Theory of Multiple Scattering of Radiation in Plane Parallel Atmospheres,” Rev. Geophys. 24, 299 (1986). [CrossRef]

16. T. Nakajima and M. Tanaka, “Matrix Formulations for the Transfer of Solar Radiation in a Plane-Parallel Atmosphere,” J. Quant. Spectrosc. Radiat. Transfer 35, 13 (1986). [CrossRef]

17. W. J. Wiscombe, “The Delta-M Method: Rapid Yet Accurate Radiative Flux Calculations for Strongly Asymmetric Phase Functions,” J. Atmos. Sci. 34, 1408 (1977). [CrossRef]

18. J. B. Sykes, “Approximate Integration of the Equation of Transfer,” Mon. Not. Roy. Astron. Soc. 11, 377 (1951).

19. W. J. Wiscombe, “Extension of the Doubling Method to Inhomogeneous Sources,” J. Quant. Spectrosc. Radiat. Transfer 16, 477 (1976). [CrossRef]

20. W. R. Cowell, Ed., Sources and Developments of Mathematical Software (Prentice Hall, Englewood Cliffs, NJ, 1980).

21. J. J. Dongarra, C. B. Moler, J. R. Bunch, and G. W. Stewart, linpack User’s Guide (SIAM, Philadelphia, 1979).

22. V. Kourganoff, Basic Methods in Transfer Problems (Dover, New York, 1963).

Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting layered media

Abstract

I. Introduction

II. Theory

A. Basic Equations

B. Discrete Ordinate Approximation—Matrix Formulation

C. Quadrature Rule

D. Basic Solution

E. Boundary Conditions

III. Numerical Implementation

A. Computation of Eigenvalues and Eigenvectors

B. Scaling Transformation

C. Scaled Solutions—Angular Distributions

D. Computational Shortcut for Highly Absorbing Media

E. Simplified Albedo and Transmissivity Computations

IV. Summary

Figure

References

Cited By

Figures (1)

Equations (71)

Applied Optics