# The 2 stream-exact single scattering (2S-ESS) radiative transfer model

V. Natraj, R. Spurr, A. Gao, T. Le, Z. C. Zeng, S. Fan, Y. L. Yung

► To cite this version:

V. Natraj, R. Spurr, A. Gao, T. Le, Z. C. Zeng, et al.. The 2 stream-exact single scattering (2S-ESS) radiative transfer model. Journal of Quantitative Spectroscopy and Radiative Transfer, 2023, 295, pp.108416. 10.1016/j.jqsrt.2022.108416. insu-03993951

HAL Id: insu-03993951

https://insu.hal.science/insu-03993951v1

Submitted on 15 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The 2 stream-exact single scattering (2S-ESS) radiative transfer model

V. Natraj<sup>a,*</sup>, R. Spurr<sup>b</sup>, A. Gao<sup>c</sup>, T. Le<sup>d</sup>, Z.C. Zeng<sup>c,d</sup>, S. Fan<sup>e</sup>, Y.L. Yung<sup>d</sup>

<figure><img src="image_2.png" /></figure>

## ARTICLE INFO

**Article history:**

Received 30 June 2022

Revised 31 October 2022

Accepted 7 November 2022

Available online 8 November 2022

**Keywords:**

Two stream

Exact single scattering

Radiative transfer

Spherical geometry

## ABSTRACT

The plane-parallel two-stream approximation is a popular radiative transfer approach for the calculation of fluxes and heating rates. However, it is typically not accurate enough for remote sensing applications involving the analysis of hyperspectral radiances. We present the 2 stream-exact single scattering (2S-ESS) radiative transfer model, which performs an exact calculation of single scattering in a spherically curved medium using an accurate treatment of the phase function and curved ray-tracing of the solar and line-of-sight paths, while approximating multiple scattering with the plane-parallel two-stream approach.

The 2S-ESS model has three important features. First, it can be deployed for calculations in vertically inhomogeneous atmospheres. Second, the sphericity capability makes it applicable to large solar and/or viewing zenith angle scenarios such as those encountered close to sunrise or sunset. Third, it is fully linearized: in addition to generating radiances, the model can also compute Jacobians analytically with respect to any atmospheric or surface property (e.g., trace gases, aerosols and surface reflectance). These features of the model are especially useful for remote sensing retrieval applications.

We examine the accuracy of this model for homogeneous slab and inhomogeneous multi-layer scenarios. The results show that this methodology introduces less than a few percent error in most situations, with the exception of heavy aerosol or cloud loading events, while providing three orders of magnitude improvement in computational efficiency. The code is publicly available along with documentation and test cases to assist the user.

© 2022 The Authors. Published by Elsevier Ltd.

This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/)

## 1. Introduction

Radiative transfer (RT) calculations are required for many remote sensing applications. However, they are often computationally expensive. Two-stream (2S) approximations — which represent the angular dependence of the intensity in terms of two functions of optical depth — are a popular alternative to numerically more precise models, especially for the calculation of fluxes and heating rates relevant to climate modeling. 2S methods avoid the complexity of numerical solutions to the RT equation (RTE), instead yielding closed-form analytical results.

A number of 2S approximations have appeared in the literature; for a review of 2S approaches, see Ref. [1]. Ref. [2] showed that the Eddington approximation can be used to compute fluxes in inhomogeneous atmospheres within a few percent. Ref. [3] derived analytic closed-form expressions for fluxes in the 2S discrete-ordinate

approximation. Ref. [4] introduced the delta-Eddington approximation, where the forward peak of the phase function is represented by a Dirac delta function and the rest of the scattering hemisphere is treated by a two-term expansion. This methodology was shown to provide more accurate flux estimates for highly asymmetric phase functions. Ref. [5] presented a framework to unify various 2S approximations. Ref. [6] investigated the accuracy of several 2S approaches. The 2S approximation was extended to inhomogeneous atmospheres by Ref. [7]. Ref. [8] introduced corrections for spherical geometry and added treatment of internal sources. The accuracy of the 2S approximation in vertically inhomogeneous media has also been investigated [8–10].

However, 2S models are not accurate enough when radiances, rather than fluxes or azimuthally averaged intensities, are the quantities of interest. Ref. [11] outlined an approach to improve computational efficiency by separating the radiation field into single and multiple scattering components, with single scattering (SS) computed accurately and multiple scattering (MS) calculated approximately. The logic behind this approach is that the SS calcu-

<figure><img src="image_3.png" /></figure>

<footnote><sup>a</sup>Jet Propulsion laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA</footnote>

<footnote><sup>b</sup>RT Solutions, Inc., 9 Channing Street, Cambridge, MA 02138, USA</footnote>

<footnote><sup>c</sup>Peking University, 5 Yiheyuan Rd, Beijing 100871, China</footnote>

<footnote><sup>d</sup>California Institute of Technology, 1200 E California Blvd, Pasadena, CA 91125, USA</footnote>

<footnote><sup>e</sup>LMD/IPSL, CNRS, PSL Research University, école Normale Supérieure, école Polytechnique, Sorbonne Université, Paris 75005, France</footnote>

<footnote>* Corresponding author.
E-mail address: vijay.natraj@jpl.nasa.gov (V. Natraj).</footnote>
lation, which often captures a large fraction of the angular dependence of the radiance, is fast; most of the computational effort is spent on calculating the MS contribution, which has smoother angular variation. It follows that, if MS is computed approximately using a fast approach, the overall RT calculation would be much faster. For this, the simplest method would be to use a 2S approximation for MS. Indeed, this is the technique employed by Ref. [12]. The combination of 2S MS and exact SS has been used in several studies [13–16].

This work has four major objectives. First, we provide a comprehensive theoretical basis for the 2 stream-exact single scattering (2S-ESS) framework for both solar beam and thermal emission forcing scenarios. Second, we perform extensive comparisons of the accuracy of the 2S-ESS model, both for the slab and multi-layer problems. Third, we describe the computational shortcuts used in the formulation and assess speed-ups relevant to remote sensing applications in the shortwave and thermal infrared spectral regions. Fourth, we provide adequate documentation so that the code, which will be publicly available, can be useful to a wide group of researchers.

The manuscript is organized as follows. Sections 2 and 3 provide a summary of the 2S and ESS theoretical basis for solar beam and emission sources, and for upwelling and downwelling radiance fields. We describe two approaches for solving the particular integral: the exponential substitution method and the Green's function method. In Section 4, we perform simulations over a large phase space of geometrical and optical parameters to assess the accuracy of the 2S-ESS model with respect to the numerically precise LIDORT [17] model. We also provide comparisons against a "Full-2S" model, which calculates not only the MS field but also a severely truncated form of the SS field. The multi-layer problem is analyzed in Section 5 for hyperspectral calculations over two spectral regions in the shortwave and thermal infrared. Comparisons against the test problems included with the popular DISORT code package [18] are considered in Section 6. Section 7 discusses computational considerations that make the 2S-ESS model appropriate for a variety of applications. The 2S approximation and the single scattering solution are well established and popularly used; differences between the established schemes and the methods employed in the 2S-ESS model are reviewed in Section 8.

The notation used in this manuscript is as follows. In addition to the use of Roman and Greek symbols in the mathematical exegeses, we use the “double struck” typeface to denote 2S RTE solutions in quadrature directions; thus, for example, $X^+$ and $X^-$ denote the downwelling (+) and upwelling (-) solutions to the homogeneous 2S RTE. For the Green’s function solutions, we make occasional use of the “script” typeface to preserve clarity.

## 2. 2S theoretical basis

In the 2S approximation, the cosine-azimuth expansion of the radiance field $I$ is restricted to the first two terms in the series:

$$I(\tau, \mu, \varphi) = I^0(\tau, \mu) + I^1(\tau, \mu)\cos\varphi, \quad (1)$$

where, $\tau$ is the optical depth coordinate, $\mu$ is the cosine of the polar angle, and $\varphi$ is the relative azimuthal angle.

### 2.1. Solar beam source

If there is no thermal emission, the RTE for Fourier component $m$ can be expressed for a single homogeneous layer as:

$$\mu \frac{dI^m(\tau, \mu)}{d\tau} = I^m(\tau, \mu) - \int_{-1}^{1} \Pi^m(\mu, \mu') I^m(\tau, \mu') d\mu' - F_0 J^m(\mu) \operatorname{Exp}(-\eta\tau) (m = 0, 1), \quad (2)$$

where,

$$\Pi^m(\mu, \mu') = \frac{\omega}{2} \sum_{l=1}^{m} \beta_l Y_l^m(\mu) Y_l^m(\mu'); \quad (3a)$$

$$J^m(\mu) = \frac{\omega}{4\pi} (2 - \delta_{m0}) \sum_{l=1}^{m} (-1)^{l+m} \beta_l Y_l^m(\mu) Y_l^m(\mu_0). \quad (3b)$$

We assume, for simplicity of exposition, that there is no surface-leaving radiance (e.g., from solar-induced fluorescence or ocean water-leaving terms) and no isotropic illumination at the top or bottom boundaries. However, the code does allow for these options. $F_0$ is the incoming solar flux, $\mu_0$ is the cosine of the solar zenith angle at the top of the atmosphere (TOA) and $\delta$ denotes the Kronecker delta function. Here, $\exp(-\eta\tau)$ is the in-layer transmittance expressed in terms of an “average-secant” $\eta$ [17]. Section 3 provides a brief discussion of the average-secant parameterization. For a plane-parallel atmosphere, $\eta = 1/\mu_0$. $T$ is the solar beam transmittance to the top boundary of the layer, and $\omega$ refers to the single scattering albedo of the layer. The quantities $Y_l^m$ are normalized associated Legendre polynomials and $\beta_l$ are the phase function expansion coefficients.

$Y_l^m$ are defined in terms of the regular associated Legendre polynomials $P_l^m$ as follows:

$$Y_l^m(\mu) = \sqrt{\frac{(l-m)!}{(l+m)!}} P_l^m(\mu). \quad (4)$$

The integral in Eq. (2) is approximated by a summation using double Gauss-Legendre quadrature, with one quadrature point (stream), $\bar{\mu}$, in each half-space. In the 2S case, superscript $m = 0$ or $1$. There are several choices for the average stream value $\bar{\mu}$; for a discussion on this issue, see Ref. [1]. The two most common choices are $\bar{\mu} = 0.5$, which is the result for a hemispherical one-point quadrature, and $\bar{\mu} = 1/\sqrt{3}$, which is the value for a two-point quadrature over the complete range of $\mu$ from -1 to +1. LIDORT and DISORT are both based on multi-stream Gaussian quadrature polar hemispheres, so it is necessary to set $\bar{\mu} = 0.5$ when validating 2S against these models.

Henceforth, we omit the superscript $m$, and in Eq. (2), we define:

$$\mathbb{I}^{\pm} = I(\tau, \pm\bar{\mu}); \quad \mathbb{J}^{\pm} = J(\pm\bar{\mu}). \quad (5)$$

Here, $\mathbb{I}^+$ and $\mathbb{J}^-$ refer to the downwelling and upwelling radiance fields, respectively.

Eq. (2) can then be written as:

$$\frac{d(\mathbb{I}^+)}{d\tau} = \begin{pmatrix} -\zeta & -\xi \\ \xi & \xi \end{pmatrix} (\mathbb{I}^+ - \mathbb{I}^-) + \frac{F_0 T}{\bar{\mu}} \exp(-\eta\tau) \begin{pmatrix} \mathbb{J}^+ \\ \mathbb{J}^- \end{pmatrix}. \quad (6)$$

where,

$$\zeta = \frac{1}{\bar{\mu}} (\mathbb{P}^+ - 1); \quad (7a)$$

$$\xi = \frac{1}{\bar{\mu}} \mathbb{P}^{-}; \quad (7b)$$

$$\mathbb{P}^{\pm} = \Pi(\bar{\mu}, \pm\bar{\mu}). \quad (7c)$$

Note that $\Pi(-\bar{\mu}, \bar{\mu}) = \Pi(\bar{\mu}, -\bar{\mu})$ and $\Pi(-\bar{\mu}, -\bar{\mu}) = \Pi(\bar{\mu}, \bar{\mu})$. Table 1 lists the values of $\mathbb{P}^{\pm}$ for $m = 0$ and $1$.

Making the transformation $\mathbb{M}^{\pm} = \mathbb{I}^{\pm} + \mathbb{I}^{-}$, we obtain:

$$\frac{d\begin{pmatrix} \mathbb{M}^+ \\ \mathbb{M}^- \end{pmatrix}}{d\tau} = \begin{pmatrix} 0 & -d \\ -s & 0 \end{pmatrix} \begin{pmatrix} \mathbb{M}^+ \\ \mathbb{M}^- \end{pmatrix} - \begin{pmatrix} \mathbb{Q}^- \\ \mathbb{Q}^+ \end{pmatrix} \operatorname{Exp}(-\eta\tau), \quad (8)$$
<table><caption>**Table 1**
Values of $P^{\pm}$ for different Fourier components in the 2S model. g refers to the asymmetry factor.</caption><thead><tr><th></th><th>m = 0</th><th>m = 1</th></tr></thead><tbody><tr><td>P<sup>+</sup></td><td>ω<sub>2</sub> (1 + 3gμ̃<sup>2</sup>)</td><td>3ω<sub>2</sub>g (1 - μ̃<sup>2</sup>)</td></tr><tr><td>P<sup>-</sup></td><td>ω<sub>2</sub> (1 - 3gμ̃<sup>2</sup>)</td><td>3ω<sub>2</sub>g (1 - μ̃<sup>2</sup>)</td></tr></tbody></table>

where,

$$s = \zeta + \xi; \qquad (9a)$$

$$d = \zeta - \xi; \qquad (9b)$$

$$Q^{\pm} = \frac{F_0}{\bar{\mu}} (J^{+} \pm J^{-}). \qquad (9c)$$

The eigenvalues of the homogeneous part of the system of Eq. (8) are $k^{\pm} = \pm\sqrt{sd}$, and the corresponding eigenvectors are $(\frac{1}{-s/k})$ and $(\frac{s/k}{1})$, respectively, where k refers to the positive eigenvalue.

The particular solution can be obtained by substituting the form $M^{\pm} = Z^{\pm}exp(-\eta\tau)$ into Eq. (8). We get:

$$Z^{+} = \frac{-\eta Q^{-} - dQ^{+}}{k^{2} - \eta^{2}}; \qquad (10a)$$

$$Z^{-} = \frac{-sQ^{-} - \eta Q^{+}}{k^{2} - \eta^{2}}. \qquad (10b)$$

Therefore, the general solution to Eq. (6) for a single optically homogeneous layer can be written as:

$$I^{\pm} = LX^{\pm}exp(-k\tau) + MX^{\mp}exp(-k(\Delta - \tau)) + W^{\pm}exp(-\eta\tau), \quad (11)$$

where,

$$X^{\pm} = \frac{1}{2}[1 \mp s/k]; \qquad (12a)$$

$$W^{\pm} = \frac{1}{2}T[Z^{+} \pm Z^{-}]. \qquad (12b)$$

L and M are constants of integration to be determined by solving the boundary value problem, and Δ is the layer optical depth. The exponential argument in the second term in Eq. (11) has been written as an optical depth difference to safeguard the numerical stability of the solution [19,20].

We remark here that the solutions $Z^{\pm}$ in Eq. (10) are invalid when $\eta = k$. Although an exact match between $\eta$ and $k$ is never realized in practice, the solutions for $Z^{\pm}$ could become numerically unstable in the limit $|\eta - k| \to 0$. Alternative expressions for $Z^{\pm}$ may be found using the infinite-medium Green's function solution to the 2S particular integral. The Green's function approach was introduced for discrete-ordinate RT in Ref. [21], and has been a standard feature of the LIDORT model for more than 20 years. Here, we write down the Green's function particular integral as a linear combination of the homogeneous solution vectors $X^{\pm}$, adapting the LIDORT formalism to the 2S case.

$$G^{+}(\tau) = \frac{1}{\mathcal{R}} [AC(\tau)X^{+} + BD(\tau)X^{-}]; \qquad (13a)$$

$$G^{-}(\tau) = \frac{1}{\mathcal{R}} [AC(\tau)X^{+} + BD(\tau)X^{-}]. \qquad (13b)$$

$$A = Q^{+}X^{+} + Q^{-}X^{-}; \qquad (14a)$$

$$B = Q^{-}X^{+} + Q^{+}X^{-}; \qquad (14b)$$

$$\mathcal{R} = \bar{\mu}[X^{+}X^{+} - X^{-}X^{-}]. \qquad (14c)$$

$$C(\tau) = T \frac{\exp(-k\tau) - \exp(-\eta\tau)}{\eta - k}; \qquad (15a)$$

$$\mathcal{D}(\tau) = T \frac{\exp(-\eta\tau) - \exp(-\eta\Delta)\exp(-[\Delta - \tau]k)}{\eta + k}. \qquad (15b)$$

In this solution, dependence on the optical depth (as measured from the top of the layer) comes through the multipliers $C(\tau)$ and $\mathcal{D}(\tau)$. We note in particular for Eq. (15a), that $C(\tau) = T\tau\exp(-k\tau) + O(\epsilon)$ in the limit as $\epsilon \equiv |\eta - k| \to 0$; the solution is always bounded.

### 2.1.1. Boundary value problem

We now look at the boundary value problem. The atmosphere is stratified into N layers, with $\tau_{n-1}$ and $\tau_n$ the optical depths at the top and bottom of layer n, respectively ($\tau_0 = 0$). In this case, $\tau$ and $\Delta - \tau$ in the homogeneous solution are replaced by $\tau - \tau_{n-1}$ and $\tau_n - \tau$, respectively, for layer n. The integration constants follow from the imposition of a set of three boundary conditions: (1) zero diffuse downwelling radiation at the upper boundary; (2) radiance continuity at intermediate layer boundaries; and (3) a surface reflection condition at the lower boundary. Using the exponential-substitution particular integral, the boundary value problem has the form:

$$X_{1}^{+}L_{1} + X_{1}^{-}Y_{1}M_{1} = -W_{1}^{+}; \qquad (16a)$$

$$X_{n}^{+}Y_{n}L_{n} + X_{n}^{-}M_{n} - X_{n+1}^{+}L_{n+1} - X_{n+1}^{-}Y_{n+1}M_{n+1} \\ = W_{n+1}^{+} - W_{n}^{+}\Lambda_{n}(n = 1, \dots, N-1); \qquad (16b)$$

$$X_{n}^{-}Y_{n}L_{n} + X_{n}^{+}M_{n} - X_{n+1}^{-}L_{n+1} - X_{n+1}^{+}Y_{n+1}M_{n+1} \\ = W_{n+1}^{-} - W_{n}^{-}\Lambda_{n}(n = 1, \dots, N-1); \qquad (16c)$$

$$V^{-}Y_{N}L_{N} + V^{+}M_{N} = -K\Lambda_{N} + D. \qquad (16d)$$

Here, we define the following quantities:

$$Y_{n} = \exp(-k_{n}\Delta_{n}); \qquad (17a)$$

$$\Lambda_{n} = \exp(-\eta_{n}\Delta_{n}); \qquad (17b)$$

$$V^{\pm} = X_{N}^{\pm} - (1 + \delta_{m0})\bar{\mu}\rho^{m}(\bar{\mu}, -\bar{\mu})X_{N}^{\mp}; \qquad (17c)$$

$$K = W_{N}^{-} - (1 + \delta_{m0})\bar{\mu}\rho^{m}(\bar{\mu}, -\bar{\mu})W_{N}^{+}; \qquad (17d)$$

$$D = \frac{\mu_{0}F_{0}}{(2 - \delta_{m0})\pi} \rho^{m}(\mu_{0}, -\bar{\mu})T_{N}\Lambda_{N}. \qquad (17e)$$

$\rho^m(\bar{\mu}, -\bar{\mu})$ refers to the mth Fourier component of the bidirectional reflectance function (BRDF) of the surface from the downwelling to upwelling quadrature direction; $\rho^m(\mu_0, -\bar{\mu})$ is the BRDF for reflection from the incoming solar direction at the surface to the upwelling quadrature direction. $T_N \Lambda_N$ is the solar beam attenuation to the surface. For a plane-parallel atmosphere and Lambertian surface with albedo $a$, we have the following expressions:

$$\rho^m(\bar{\mu}, -\bar{\mu}) = \rho^m(\mu_0, -\bar{\mu}) = \delta_{m0}a; \qquad (18a)$$

$$T_N \Lambda_N = \exp(-\tau_N / \mu_0). \qquad (18b)$$
<figure><img src="image_4.png" /><figcaption>Fig. 1. Radiance comparisons between LIDORT and (left) 2S-ESS, (right) Full-2S results for (top) upwelling at the TOA, and (bottom) downwelling at the surface. The incoming solar flux is assumed to be unity; the radiances are, therefore, unitless.</figcaption></figure>

In the 2S model, the boundary value problem is the only part requiring numerical solution techniques. However, the system of Eqs. (16a)-(16d) has a penta-diagonal form; therefore, specialized sparse matrix solvers can be used instead of traditional inversion routines. For this purpose, we have adapted the penta-diagonal solver taken from the website http://cococubed.asu.edu/codes/penta/matrix_pent.f. Note that it is possible to recast the linear algebra system in Eqs. (16a)-(16d) to read from the surface to the TOA; in effect, we would then solve the penta-diagonal system in the reverse direction. In practice, we have found that the reversed system is more stable numerically. Once found, the inverse matrix is stored and then applied repeatedly in the back substitutions required for solution of the linearized boundary value problem (not discussed here, but relevant to retrievals). The boundary value problem could also be cast in a tri-diagonal form (see, e.g., Ref. [7]). Both forms have similar computational cost.

For the Green's function solution, there will be a different right-hand-side in Eqs. (16a)-(16d); this takes the form:

$$-\mathbb{G}_{1}^{+}(0); \qquad (19a)$$

$$\mathbb{G}_{n+1}^{+}(0) - \mathbb{G}_{n}^{+}(\Delta_n)(n = 1, \dots, N-1); \qquad (19b)$$

$$\mathbb{G}_{n+1}^{-}(0) - \mathbb{G}_{n}^{-}(\Delta_n)(n = 1, \dots, N-1); \qquad (19c)$$

$$-\mathbb{H} + \mathbb{D}. \qquad (19d)$$

Here, the Green's function solutions are evaluated at the upper and lower layer-boundary optical depths, the direct reflectance term $\mathbb{D}$ has already been defined, and:

$$\mathbb{H} = \mathbb{G}_{N}^{-}(\Delta_{N}) - (1 + \delta_{m0})\bar{\mu}\rho^{m}(\bar{\mu}, -\bar{\mu})\mathbb{G}_{N}^{+}(\Delta_{N}). \qquad (20)$$

### 2.1.2. Post-processing

The discrete-ordinate solution gives the radiance at the Gaussian quadrature angle $\bar{\mu}$. Post-processing is required to obtain radiances at other viewing angles, and this is done by substituting the discrete-ordinate solution for the quadrature angle into the RTE, and integrating over the optical depth coordinate. In an inhomogeneous atmosphere, the source terms must be integrated on a layer-by-layer basis. The following is obtained for the upwelling radiance at the boundaries of layer n:

$$I(\tau_{n-1}, -\mu) = I(\tau_n, -\mu)e^{-\Delta_n/\mu} + \int_0^{\Delta_n} [\Pi(-\mu, \bar{\mu})\mathbb{I}^+ + \Pi(-\mu, -\bar{\mu})\mathbb{I}^-] \frac{e^{-\tau/\mu}}{\mu} d\tau. \qquad (21)$$

Note that we employ the 2S model in MS mode. Therefore, the post-processing result in Eq. (21) does not include the SS contribution. Eq. (21) is valid for both kinds of particular integral; the optical depth integrations can be performed analytically to obtain:

$$I(\tau_{n-1}, -\mu) = I(\tau_n, -\mu)e^{-\Delta_n/\mu} + L_n H_n^+(\mu)U_n^-(\mu) + M_n H_n^+(\mu)U_n^+(\mu) + \Psi_n^-(\mu), \qquad (22)$$
<table><caption>**Table 2**
2S solutions for the RTE with solar sources (quadrature and post-processed). Layer index has been omitted. Here, c = 3ωg.</caption><thead><tr><th>Solution type</th><th>m = 0</th><th>m = 1</th></tr></thead><tbody><tr><td>Homogeneous</td><td>s = <sup>1</sup>&frasl;<sub>&mu;</sub>(&omega; &minus; 1); d = <sup>1</sup>&frasl;<sub>&mu;</sub>(c&mu;&macr;<sup>2</sup> &minus; 1)</td><td>s = <sup>c(1 &minus; &mu;<sup>2</sup>) &minus; 2</sup>&frasl;<sub>2&mu;</sub>; d = &minus;<sup>1</sup>&frasl;<sub>&mu;</sub></td></tr><tr><td>Particular</td><td>Q<sup>+</sup> = <sup>&omega;&phi;<sub>n</sub></sup>&frasl;<sub>2&pi;&mu;</sub>; Q<sup>&minus;</sup> = &minus;<sup>c&phi;<sub>n</sub>&mu;<sub>0</sub></sup>&frasl;<sub>2&pi;</sub></td><td>Q<sup>+</sup> = <sup>c&phi;<sub>0</sub></sup>&frasl;<sub>2&pi;&mu;</sub>(1 &minus; &mu;&macr;<sup>2</sup>)<sup>1/2</sup>(1 &minus; &mu;<sub>0</sub><sup>2</sup>)<sup>1/2</sup>; Q<sup>&minus;</sup> = 0</td></tr><tr><td>Post-processed homogeneous</td><td>U<sup>&pm;</sup> = <sup>1</sup>&frasl;<sub>2</sub>(&omega; + c&mu;&macr;&mu;)<sup>X&pm;</sup> + <sup>1</sup>&frasl;<sub>2</sub>(&omega; &minus; c&mu;&macr;&mu;)<sup>X&mp;</sup></td><td>U<sup>&pm;</sup> = <sup>c</sup>&frasl;<sub>4</sub>(1 &minus; &mu;<sup>2</sup>)<sup>1/2</sup>(1 &minus; &mu;&macr;<sup>2</sup>)<sup>1/2</sup></td></tr><tr><td>Post-processed particular (substitution method)</td><td>G<sup>&pm;</sup> = <sup>1</sup>&frasl;<sub>2</sub>(&omega; + c&mu;&macr;&mu;)<sup>W&mp;</sup> + <sup>1</sup>&frasl;<sub>2</sub>(&omega; &minus; c&mu;&macr;&mu;)<sup>W&mp;</sup></td><td>G<sup>&pm;</sup> = <sup>c</sup>&frasl;<sub>4</sub>(1 &minus; &mu;<sup>2</sup>)<sup>1/2</sup>(1 &minus; &mu;&macr;<sup>2</sup>)<sup>1/2</sup>(W<sup>+</sup> + W<sup>&mp;</sup>)</td></tr></tbody></table>

<table><caption>**Table 3**
Parameters used for slab problem comparisons.</caption><thead><tr><th>Parameter</th><th>Values</th></tr></thead><tbody><tr><td>Solar zenith angle (SZA)</td><td>10, 20, ..., 80, 85, 89 (10 values)</td></tr><tr><td>Viewing zenith angle (VZA)</td><td>0, 10, ..., 80, 85, 89 (11 values)</td></tr><tr><td>Δ</td><td>0.001, 0.002, ..., 0.01, 0.02, ..., 0.1, 0.2, ..., 1, 2, ..., 10, 20, ..., 100, 200, ..., 1000 (55 values)</td></tr><tr><td>ω</td><td>0.001, 0.01, 0.05, 0.1, 0.2, ..., 0.9, 0.95, 0.99, 0.999 (15 values)</td></tr><tr><td>g</td><td>0.1, 0.5, 0.9, 0.95 (4 values)</td></tr><tr><td>a</td><td>0.1, 0.3, 0.5, 0.9 (4 values)</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Probability distribution functions for (left) surface downwelling and (right) TOA upwelling radiance. The incoming solar flux is assumed to be unity; the radiances are, therefore, unitless.</figcaption></figure>

where,

$$H_{n}^{-}(\mu) = \frac{1 - e^{-k_{n}\Delta_{n}}e^{-\Delta_{n}/\mu}}{1 + \mu k_{n}}; \qquad (23a)$$

$$U_{n}^{-}(\mu) = \frac{1}{2}[\Pi(-\mu, \bar{\mu})(1 - s_{n}/k_{n}) + \Pi(\mu, \bar{\mu})(1 + s_{n}/k_{n})]; \quad (23b)$$

$$H_{n}^{+}(\mu) = \frac{e^{-k_{n}\Delta_{n}} - e^{-\Delta_{n}/\mu}}{1 - \mu k_{n}}; \qquad (23c)$$

$$U_{n}^{+}(\mu) = \frac{1}{2}[\Pi(-\mu, \bar{\mu})(1 + s_{n}/k_{n}) + \Pi(\mu, \bar{\mu})(1 - s_{n}/k_{n})]. \quad (23d)$$

For the exponential-substitution particular integral, contribution $\Psi_n^-$ ($\mu$) takes the form:

$$\Psi_{n}^{-}(\mu) = E_{n}^{-}(\mu)G_{n}^{-}(\mu); \qquad (24a)$$

$$E_{n}^{-}(\mu) = \frac{1 - e^{-\eta_{n}\Delta_{n}}e^{-\Delta_{n}/\mu}}{1 + \mu\eta_{n}}; \qquad (24b)$$

$$G_{n}^{-}(\mu) = \Pi(-\mu, \bar{\mu})W_{n}^{+} + \Pi(\mu, \bar{\mu})W_{n}^{-}. \qquad (24c)$$

For the Green's function particular integral, contribution $\Psi_n^-$ ($\mu$) takes the form:

$$\Psi_{n}^{-}(\mu) = U_{n}^{+}(\mu)A_{n}(\mu) \frac{H_{n}^{-}(\mu) - E_{n}^{-}(\mu)}{\eta_{n} - k_{n}} + U_{n}^{-}(\mu)B_{n}(\mu) \frac{E_{n}^{-}(\mu) - H_{n}^{+}(\mu)\Delta_{n}}{\eta_{n} + k_{n}}. \qquad (25)$$

Given the expressions for $H_n^-$ ($\mu$) and $E_n^-$ ($\mu$), one can use Taylor series expansions to show that the first term on the right-hand side of Eq. (25) is finite in the limit $\eta_n \to k_n$; full expressions for these Taylor-series limits may be found in the literature (see, e.g., Ref. [22]).

The starting value for the recursion in Eq. (21) is the lower boundary source term $I(\tau_N, -\mu)$. This is given by:

$$I(\tau_N, -\mu) = (1 + \delta_{m0})\bar{\mu}\rho^m(\bar{\mu}, -\mu) \left[ L_N\left(\frac{1 - s_N/k_N}{2}\right) e^{-k_N \Delta_N} + M_N\left(\frac{1 + s_N/k_N}{2}\right) + \mathbb{E}_N \right] + \frac{\mu_0 F_0}{\pi} \rho^m(\eta_N^{-1}, -\mu) T_N \Lambda_N, \qquad (26)$$

where, $\mathbb{E}_N = W_N^+ \Lambda_N$ for the exponential-substitution particular integral, and $\mathbb{E}_N = G_N^+(\Delta_N)$ for the Green's function particular integral.

Similarly, for the downwelling radiance at the boundaries of layer n, we have:

$$I(\tau_n, \mu) = I(\tau_{n-1}, \mu)e^{-\Delta_n/\mu} + \int_0^{\Delta_n} [\Pi(\mu, \bar{\mu})\mathbb{I}^+ + \Pi(\mu, -\bar{\mu})\mathbb{I}^-] \frac{e^{-(\Delta_n - \tau)/\mu}}{\mu} d\tau. \quad (27)$$

Again, Eq. (27) is valid for both forms of the particular integral, and performing the optical depth integrations analytically, we obtain:
<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Radiance comparisons between LIDORT and (top) 2S-ESS, (bottom) Full-2S upwelling results for Case 1 from Table 5. The incoming solar flux is assumed to be unity; the radiances are, therefore, unitless. Positive (negative) relative errors imply over- (under-)estimation by the approximate models. Note the different color scales for the 2S-ESS and Full-2S errors.</figcaption></figure>

<table><caption>**Table 4**
Earth Mover's Distance comparison for 2S-ESS and Full-2S with respect to LIDORT.</caption><thead><tr><th></th><th>Model</th><th></th></tr><tr><th></th><th>2S-ESS</th><th>Full-2S</th></tr></thead><tbody><tr><td>Downwelling</td><td>0.03266</td><td>0.39348</td></tr><tr><td>Upwelling</td><td>0.00993</td><td>0.07218</td></tr></tbody></table>

$$I(\tau_n, \mu) = I(\tau_{n-1}, \mu)e^{-\Delta_n/\mu} + L_n H_n^+(\mu) U_n^+(\mu) \\ + M_n H_n^-(\mu) U_n^-(\mu) + \Psi_n^+(\mu). \qquad (28)$$

The starting value for the recursion is $I(0, \mu) = 0$.

For the exponential-substitution particular integral contribution $\Psi_n^+(\mu)$ is given by:

$$\Psi_n^+(\mu) = E_n^+(\mu)G_n^+(\mu); \qquad (29a)$$

$$E_n^+(\mu) = \frac{e^{-\eta_n \Delta_n} - e^{-\Delta_n/\mu}}{1 - \mu \eta_n}; \qquad (29b)$$

$$G_n^+(\mu) = \Pi(\mu, \tilde{\mu})W_n^+ + \Pi(-\mu, \tilde{\mu})W_n^-. \qquad (29c)$$

For the Green's function particular integral, we have:

$$\Psi_n^+(\mu) = U_n^-(\mu)A_n(\mu) \frac{H_n^+(\mu) - E_n^+(\mu)}{\eta_n - k_n} \\ + U_n^+(\mu)B_n(\mu) \frac{E_n^+(\mu) - H_n^-(\mu)\Delta_n}{\eta_n + k_n}. \qquad (30)$$

The first term in Eq. (30) will have a Taylor-series expansion in the limit $\eta_n \to k_n$. The starting value for the downwelling recursion is $I(0, \mu) = 0$. We remark that post-processed output can be found at any intermediate optical depth in any layer, as partial-layer integration over optical depth is also a straightforward procedure, given the exponential dependencies present in the homogeneous and particular integral solutions.

Expressions for the various solutions are written down in Table 2 for $m = 0$ and 1. We note that there is a singularity in the homogeneous solution for $m = 0$ for the case of conservative scattering. Indeed, the homogeneous solutions have an entirely different form for that scenario. This will be the topic of a subsequent paper; here, we handle that situation by toggling the single scattering albedo to make it slightly different from unity (see Ref. [23] for a more detailed discussion on this point).
<figure><img src="image_7.png" /><figcaption>Fig. 4. Same as Fig. 3, but for downwelling results.</figcaption></figure>

<table><caption>**Table 5**

Scenarios for radiance comparison plots. For all cases, *a* = 0.3.</caption><thead><tr><th>Case</th><th>Parameter Values</th><th>Gas Absorption</th><th>Scattering</th><th>Example Scatterer Type</th></tr></thead><tbody><tr><td>1</td><td>Δ = 0.001, ω = 0.999, g = 0.1</td><td>Very weak</td><td>Weak</td><td>Soot</td></tr><tr><td>2</td><td>Δ = 1, ω = 0.001, g = 0.95</td><td>Moderate</td><td>Weak</td><td>Sea salt</td></tr><tr><td>3</td><td>Δ = 1, ω = 0.999, g = 0.95</td><td>Weak</td><td>Moderate</td><td>Dust</td></tr><tr><td>4</td><td>Δ = 100, ω = 0.001, g = 0.95</td><td>Heavy</td><td>Moderate</td><td>Cirrus</td></tr><tr><td>5</td><td>Δ = 100, ω = 0.999, g = 0.95</td><td>Moderate</td><td>Heavy</td><td>Water cloud</td></tr></tbody></table>

<table><caption>**Table 6**

Error statistics (downwelling). The 2S-ESS and Full-2S values are listed outside and within parentheses, respectively.</caption><thead><tr><th rowspan="2">Scenario</th><th rowspan="2">Number of cases</th><th colspan="4">Percentage of cases with error less than</th></tr><tr><th>0.1%</th><th>0.5%</th><th>1%</th><th>5%</th></tr></thead><tbody><tr><td>Overall</td><td>1,452,000</td><td>22.72 (20.07)</td><td>35.56 (21.04)</td><td>43.83 (22.24)</td><td>62.77 (31.85)</td></tr><tr><td>SZA, VZA &lt; 85°</td><td>1,161,600</td><td>22.04 (19.80)</td><td>35.35 (20.75)</td><td>44.09 (21.92)</td><td>63.65 (32.05)</td></tr><tr><td>SZA, VZA &lt; 80°</td><td>1,016,400</td><td>21.60 (19.59)</td><td>34.76 (20.53)</td><td>43.57 (21.67)</td><td>63.46 (32.12)</td></tr><tr><td>Δ &lt; 0.1</td><td>475,200</td><td>12.78 (0.35)</td><td>41.94 (1.76)</td><td>58.87 (3.46)</td><td>85.50 (18.69)</td></tr><tr><td>0.1 ≤ Δ ≤ 10</td><td>501,600</td><td>4.67 (0.35)</td><td>13.80 (1.79)</td><td>21.52 (3.57)</td><td>49.59 (16.36)</td></tr><tr><td>Δ &gt; 10</td><td>475,200</td><td>51.71 (60.59)</td><td>52.15 (60.65)</td><td>52.34 (60.73)</td><td>53.96 (61.36)</td></tr><tr><td>ωΔ &lt; 0.1</td><td>617,760</td><td>12.98 (0.36)</td><td>40.45 (1.85)</td><td>56.57 (3.65)</td><td>83.77 (18.52)</td></tr><tr><td>0.1 ≤ ωΔ ≤ 10</td><td>469,920</td><td>10.59 (10.40)</td><td>14.12 (11.43)</td><td>18.43 (12.70)</td><td>40.77 (22.38)</td></tr><tr><td>ωΔ &gt; 10</td><td>364,320</td><td>51.44 (62.46)</td><td>51.49 (62.51)</td><td>51.56 (62.57)</td><td>52.12 (63.17)</td></tr><tr><td>g &lt; 0.95</td><td>1,089,000</td><td>23.01 (20.59)</td><td>36.21 (21.84)</td><td>44.98 (23.37)</td><td>65.08 (35.54)</td></tr><tr><td>g &lt; 0.9</td><td>726,000</td><td>11.28 (6.57)</td><td>26.66 (9.57)</td><td>37.11 (12.31)</td><td>61.21 (30.92)</td></tr></tbody></table>
<figure><img src="image_8.png" /><figcaption>Fig. 5. Same as Fig. 3, but for Case 2 from Table 5.</figcaption></figure>

## 2.2. Thermal emission source

If thermal emission is present, Eq. (6) can be written as (ignoring beam sources for the moment):

$$ \frac{d\left(\frac{\mathbb{I}^+}{\mathbb{I}^-}\right)}{d\tau} = \begin{pmatrix} -\zeta & -\xi \\ \xi & \zeta \end{pmatrix} \left(\frac{\mathbb{I}^+}{\mathbb{I}^-}\right) \mp \frac{1-\omega}{\bar{\mu}} \begin{pmatrix} B(\tau) \\ B(\tau) \end{pmatrix}, \quad (31) $$

where, $B(\tau)$ refers to the Planck function. The emitted radiation is isotropic, so that the thermal sources are present only for Fourier component $m = 0$. $B(\tau)$ is typically assumed to vary linearly with optical depth:

$$ B(\tau) = b_0 + b_1 \tau, \quad (32) $$

with $b_0$ and $b_1$ obtained from the known temperatures (and hence known Planck functions) at the layer boundaries.

The particular solution can be obtained by substituting the form $S^{\pm} = \mathbb{R}_{0}^{\pm} + \mathbb{R}_{1}^{\pm} \tau$ into Eq. (31). We get:

$$ \mathbb{R}_{0}^{+} = -\frac{1-\omega}{\bar{\mu}s}(b_0 + b_1/d); \quad (33a) $$

$$ \mathbb{R}_{0}^{+} = -\frac{1-\omega}{\bar{\mu}s}(b_0 - b_1/d); \quad (33b) $$

$$ \mathbb{R}_{1}^{+} = -\frac{(1-\omega)b_{1}}{\bar{\mu}s}; \quad (33c) $$

$$ \mathbb{R}_{1}^{-} = -\frac{(1-\omega)b_{1}}{\bar{\mu}s}. \quad (33d) $$

The general solution to Eq. (31) for a single optically homogeneous layer can be written as:

$$ \mathbb{I}^{\pm} = L\mathbb{X}^{\pm}\exp(-k\tau) + M\mathbb{X}^{\mp}\exp(-k(\Delta - \tau)) + \mathbb{R}_{0}^{\pm} + \mathbb{R}_{1}^{\pm}\tau. \quad (34) $$

### 2.2.1. Boundary value problem

The boundary value problem for the N-layer scenario takes the form:

$$ \mathbb{X}_{1}^{+}L_{1} + \mathbb{X}_{1}^{-}Y_{1}M_{1} = -\mathbb{R}_{0,1}^{+}; \quad (35a) $$

$$ \begin{aligned} & \mathbb{X}_{n}^{+}Y_{n}L_{n} + \mathbb{X}_{n}^{-}M_{n} - \mathbb{X}_{n+1}^{+}L_{n+1} - \mathbb{X}_{n-1}^{-}Y_{n+1}M_{n+1} \\ & = \mathbb{R}_{0,n+1}^{+} - \mathbb{R}_{0,n}^{-} - \mathbb{R}_{1,n}^{+}\Delta_{n}(n = 1, \dots, N-1); \end{aligned} \quad (35b) $$

$$ \begin{aligned} & \mathbb{X}_{n}^{-}Y_{n}L_{n} + \mathbb{X}_{n}^{+}M_{n} - \mathbb{X}_{n+1}^{-}L_{n+1} - \mathbb{X}_{n+1}^{+}Y_{n+1}M_{n+1} \\ & = \mathbb{R}_{0,n+1}^{-} - \mathbb{R}_{0,n}^{-} - \mathbb{R}_{1,n}^{-}\Delta_{n}(n = 1, \dots, N-1); \end{aligned} \quad (35c) $$

$$ \mathbb{V}^{-}Y_{N}L_{N} + \mathbb{V}^{+}M_{N} = -\mathbb{P} + \varepsilon B_{\text{surf}}. \quad (35d) $$

Here, $\varepsilon$ and $B_{surf}$ refer to the surface emissivity and Planck function, respectively, and $\mathbb{P}$ is defined as:

$$ \mathbb{P} = \mathbb{R}_{0,N}^{-} - \mathbb{R}_{1,N}^{-}\Delta_{N} - 2\bar{\mu}\rho^{0}(\bar{\mu}, -\bar{\mu})(\mathbb{R}_{0,N}^{+} - \mathbb{R}_{1,N}^{+}\Delta_{N}). \quad (36) $$
<figure><img src="image_9.png" /><figcaption>Fig. 6. Same as Fig. 5, but for downwelling results.</figcaption></figure>

The surface emissivity can be determined using Kirchhoff's law:

$$ \varepsilon = 1 - 2\bar{\mu}\rho^0(\bar{\mu}, -\bar{\mu}). \quad (37) $$

### 2.2.2. Post-processing

For upwelling radiance at user angles, we have the recursion:

$$ I(\tau_{n-1}, -\mu) = I(\tau_n, -\mu)e^{-\Delta_n/\mu} + L_n H_n^{-}(\mu)U_n^{-}(\mu) \\ + M_n H_n^{+}(\mu)U_n^{+}(\mu) + O_{n,0}^{-}(\mu)S_{n,0}^{-}(\mu) \\ + O_{n,1}^{-}(\mu)S_{n,1}^{-}(\mu), \quad (38) $$

where,

$$ O_{n,0}^{-}(\mu) = 1 - e^{-\Delta_n/\mu}; \quad (39a) $$

$$ S_{n,0}^{-}(\mu) = \delta_{m0}[\Pi(-\mu, \bar{\mu})\mathbb{R}_{0,n}^{+} + \Pi(\mu, \bar{\mu})\mathbb{R}_{0,n}^{-}]; \quad (39b) $$

$$ O_{n,1}^{-}(\mu) = \mu(1 - e^{-\Delta_n/\mu}) - \Delta_n e^{-\Delta_n/\mu}; \quad (39c) $$

$$ S_{n,1}^{-}(\mu) = \delta_{m0}[\Pi(-\mu, \bar{\mu})\mathbb{R}_{1,n}^{+} + \Pi(\mu, \bar{\mu})\mathbb{R}_{1,n}^{-}]. \quad (39d) $$

The starting value for the recursion is given by:

$$ I(\tau_N, -\mu) = (1 + \delta_{m0})\bar{\mu}\rho^m(\bar{\mu}, -\mu) \left[ L_N \left( \frac{1 - s_N/k_N}{2} \right) e^{-k_N \Delta_N} \right. \\ \left. + M_N \left( \frac{1 + s_N/k_N}{2} \right) + R_0^+ + R_1^+ \Delta_N \right]. \quad (40) $$

The downwelling radiance obeys the recursion:

$$ I(\tau_n, \mu) = I(\tau_{n-1}, \mu)e^{-\Delta_n/\mu} + L_n H_n^{+}(\mu)U_n^{+}(\mu) + M_n H_n^{-}(\mu)U_n^{-}(\mu) \\ + O_{n,0}^{+}(\mu)S_{n,0}^{+}(\mu) + O_{n,1}^{+}(\mu)S_{n,1}^{+}(\mu), \quad (41) $$

where,

$$ O_{n,0}^{+}(\mu) = 1 - e^{-\Delta_n/\mu}; \quad (42a) $$

$$ S_{n,0}^{+}(\mu) = \delta_{m0}[\Pi(\mu, \bar{\mu})\mathbb{R}_{0,n}^{+} + \Pi(-\mu, \bar{\mu})\mathbb{R}_{0,n}^{-}]; \quad (42b) $$

$$ O_{n,1}^{+}(\mu) = \Delta_n - \mu(1 - e^{-\Delta_n/\mu}); \quad (42c) $$

$$ S_{n,1}^{+}(\mu) = \delta_{m0}[\Pi(\mu, \bar{\mu})\mathbb{R}_{1,n}^{+} + \Pi(-\mu, \bar{\mu})\mathbb{R}_{1,n}^{-}]. \quad (42d) $$

As for the solar beam scenario, the starting value for the recursion is $I(0, \mu) = 0$.

The RTE with both beam and thermal sources can be handled in the same manner as described above for the solar- or emission-only scenarios. The Green's function particular integral is also available for the thermal regime [17]:

$$ T_n^{\pm}(\tau) = \frac{(1 - \omega_n)}{\mathcal{R}} [c_n^{(BB)}(\tau)X_n^{\mp} + D_n^{(BB)}(\tau)X_n^{\pm}]. \quad (43) $$

Here, $X_n^{\pm}$ and the normalization factor $\mathcal{R}$ have already been defined.
<figure><img src="image_10.png" /><figcaption>Fig. 7. Same as Fig. 3, but for Case 3 from Table 5.</figcaption></figure>

The Green's function thermal source multipliers $C_n^{(BB)}(\tau)$ and $D_n^{(BB)}(\tau)$ are given by:

$$C_n^{(BB)}(\tau) = e^{-\tau k_n} \int_0^\tau e^{+y k_n} B_n(y) dy; \quad (44a)$$

$$D_n^{(BB)}(\tau) = e^{+\tau k_n} \int_{\tau}^{\Delta_n} e^{-y k_n} B_n(y) dy. \quad (44b)$$

Once again, $\tau$ is the vertical optical depth, and we assume the usual piecewise-linear dependence of the Planck function on $\tau$.

We evaluate the integrals in Eq. (44) at $\tau = \Delta_n$ and $\tau = 0$, respectively, and express the results as linear combinations of the level boundary input Planck functions $B_{n-1}$ and $B_n$:

$$C_n^{(BB)}(\Delta_n) = P_n^+ B_{n-1} + Q_n^+ B_n; \quad (45a)$$

$$D_n^{(BB)}(0) = P_n^- B_{n-1} + Q_n^- B_n. \quad (45b)$$

Here we have:

$$P_n^+ = \Phi_n - \frac{1}{\Delta_n k_n}(\Delta_n - \Phi_n); \quad Q_n^+ = \frac{1}{\Delta_n k_n}(\Delta_n - \Phi_n); \quad (46a)$$

$$P_n^- = \Phi_n + \frac{1}{\Delta_n k_n}(\Delta_n Y_n - \Phi_n); \quad Q_n^- = -\frac{1}{\Delta_n k_n}(\Delta_n Y_n - \Phi_n); \quad (46b)$$

$$Y_n = e^{-\Delta_n k_n}; \quad \Phi_n = \frac{1}{k_n}(1 - Y_n). \quad (46c)$$

We note also that:

$$C_n^{(BB)}(0) = 0; \quad D_n^{(BB)}(\Delta_n) = 0. \quad (47)$$

We can now compute the particular integrals at layer boundaries, by substituting these results into Eq. (43); again, we write the answers as linear combinations of Planck functions:

$$T_n^{\pm}(0) = \frac{(1-\omega_n)S_n}{\mathcal{R}} [B_{n-1} P_n^{-} X_n^{\pm} + B_n Q_n^{-} X_n^{\pm}]; \quad (48a)$$

$$T_n^{\pm}(\Delta_n) = \frac{(1-\omega_n)S_n}{\mathcal{R}} [B_{n-1} P_n^{+} X_n^{\pm} + B_n Q_n^{+} X_n^{\pm}]. \quad (48b)$$

Although we have evaluated the multipliers in Eq. (44) at the layer boundaries, these integrals can be done readily enough for arbitrary optical depth $\tau$, and quantities $T_n^{\pm}(x)$ can then be derived analytically as linear combinations of Planck functions as in Eq. (48).

The boundary value problem including this Green's function thermal solution takes the same form as that in Eqs. (35a)-(35d), except that the right-hand-side is now:

$$-T_1^+(0); \quad (49a)$$

$$T_{n+1}^+(0) - T_n^+(\Delta_n)(n=1, \dots, N-1); \quad (49b)$$

$$T_{n+1}^-(0) - T_n^-(\Delta_n)(n=1, \dots, N-1); \quad (49c)$$

$$-P + \varepsilon B_{\text{surf}}. \quad (49d)$$
<figure><img src="image_11.png" /><figcaption>Fig. 8. Same as Fig. 7, but for downwelling results.</figcaption></figure>

<caption>**Table 7**</caption>

Same as Table 6, but for upwelling.

<table><thead><tr><th rowspan="2">Scenario</th><th rowspan="2">Number of cases</th><th colspan="4">Percentage of cases with error less than</th></tr><tr><th>0.1%</th><th>0.5%</th><th>1%</th><th>5%</th></tr></thead><tbody><tr><td>Overall</td><td>1,452,000</td><td>28.75 (16.35)</td><td>43.45 (25.11)</td><td>51.89 (30.45)</td><td>79.16 (48.68)</td></tr><tr><td>SZA, VZA &lt; 85°</td><td>1,161,600</td><td>32.31 (19.23)</td><td>47.67 (29.02)</td><td>56.23 (34.76)</td><td>83.24 (54.13)</td></tr><tr><td>SZA, VZA &lt; 80°</td><td>1,016,400</td><td>33.36 (20.66)</td><td>48.59 (30.72)</td><td>57.06 (36.67)</td><td>83.88 (56.54)</td></tr><tr><td>Δ &lt; 0.1</td><td>475,200</td><td>58.49 (44.65)</td><td>75.73 (62.43)</td><td>81.59 (69.00)</td><td>92.71 (81.96)</td></tr><tr><td>0.1 ≤ Δ ≤ 10</td><td>501,600</td><td>15.07 (4.35)</td><td>30.58 (10.72)</td><td>41.05 (16.44)</td><td>74.70 (38.78)</td></tr><tr><td>Δ &gt; 10</td><td>475,200</td><td>13.26 (0.56)</td><td>24.67 (2.94)</td><td>33.59 (6.67)</td><td>70.33 (25.84)</td></tr><tr><td>ωΔ &lt; 0.1</td><td>617,760</td><td>53.07 (34.91)</td><td>71.40 (51.60)</td><td>77.82 (58.64)</td><td>89.23 (72.83)</td></tr><tr><td>0.1 ≤ ωΔ ≤ 10</td><td>469,920</td><td>12.25 (0.88)</td><td>27.92 (4.53)</td><td>38.81 (9.32)</td><td>73.95 (30.96)</td></tr><tr><td>ωΔ &gt; 10</td><td>364,320</td><td>1.60 (0.50)</td><td>9.12 (2.58)</td><td>17.85 (5.74)</td><td>62.03 (26.54)</td></tr><tr><td>g &lt; 0.95</td><td>1,089,000</td><td>28.29 (16.25)</td><td>43.42 (25.40)</td><td>52.25 (31.42)</td><td>82.22 (52.66)</td></tr><tr><td>g &lt; 0.9</td><td>726,000</td><td>27.61 (16.20)</td><td>44.03 (26.31)</td><td>53.82 (33.54)</td><td>89.04 (60.90)</td></tr></tbody></table>

where,

$$P = T_{N}^{-}(\Delta_{N}) - (1 + \delta_{m0})\bar{\mu}\rho^{0}(\bar{\mu}, -\bar{\mu})T_{N}^{+}(\Delta_{N}). \quad (50)$$

For the post-processing, we proceed as before, substituting the Green's function result into the original RTE and recursively evaluating the optical depth integrals on a layer-by-layer basis to obtain the intensity at any boundary level. This is a straightforward if lengthy algebraic exercise, with a number of exponential-linear integrals of the form $f(a+bx) \exp[-kx]dx$ to determine. As with the solar source terms, it is possible to carry out the source function integration to any level in any layer, as optical depth integrals can be done readily enough for arbitrary optical depth $\tau$. These derivations are left as an exercise for the reader.

### 3. ESS theoretical basis

The formalism here is taken mostly from the exposition found in the (V)LIDORT review paper [22]. Consequently, the convention in the present work is to use a “-” sign for the downwelling polar angle direction, whereas a “+” sign is used for the same quantity in Section 2. Further, we note that the “down” and “up” arrows in the ESS formulation for the radiation field hemispheres correspond, respectively, to the signage “+” for downwelling and “-” for upwelling used in Section 2. The ESS formalism also includes the direct-beam surface reflection contribution.

Dealing first with SS in the upwelling direction $\Omega = (\mu, \varphi)$, for a given layer $n$ in an optically-stratified atmospheric model, the
<caption>**Table 8**</caption>

DISORT Test Problem 1 (Isotropic Scattering). For problems 1a, 1b, 1d, and 1e, the exact values (shown in parentheses) are taken from Table 12 of Ref. [37], except for the flux divergence where they are based on LIDORT calculations using 16 half-space streams (these agree with values calculated using Eqs. (70)–(72) and the “AVERAGE” intensity taken from Ref. [37]). For conservative scattering, the flux divergence is identically zero, as calculated using Eq. (71). The exact value of $F^+$ for these problems is obtained by subtracting $F_{dir}^+$, as calculated using Eq. (70), from the total downwelling flux tabulated in Ref. [37]. For problems 1c and 1f, the exact values are based on LIDORT calculations using 16 half-space streams. Intensities are multiplied by $\pi$ to be consistent with the definition in Ref. [37].

<table><thead><tr><th>τ</th><th>F<sup>+</sup><sub>dir</sub></th><th>F<sup>+</sup></th><th>F<sup>-</sup></th><th>dF/dτ</th><th>μ</th><th>I<sup>+</sup></th><th>I<sup>-</sup></th></tr></thead><tbody><tr><td rowspan="3">0</td><td rowspan="3">1 (1)</td><td rowspan="3">0 (0)</td><td>0.02622</td><td>Test Problem 1a</td><td>0.1</td><td>0 (0)</td><td>0.11689 (0.11779)</td></tr><tr><td>(0.02542)</td><td>8.04195</td><td>0.5</td><td>0 (0)</td><td>0.02622 (0.02642)</td></tr><tr><td></td><td>(8.08663)</td><td>1.0</td><td>0 (0)</td><td>0.01330 (0.01341)</td></tr><tr><td rowspan="3">0.03125</td><td rowspan="3">0.73162<br>(0.73162)</td><td rowspan="3">0.02613<br>(0.02526)</td><td rowspan="3">0 (0)</td><td>5.89474</td><td>0.1</td><td>0.11502 (0.11592)</td><td>0 (0)</td></tr><tr><td>(5.93618)</td><td>0.5</td><td>0.02614 (0.02634)</td><td>0 (0)</td></tr><tr><td></td><td>1.0</td><td>0.01328 (0.01339)</td><td>0 (0)</td></tr><tr><td rowspan="3">0</td><td rowspan="3">1 (1)</td><td rowspan="3">3.72944E-9<sup>-</sup> (0)</td><td>0.13440</td><td>Test Problem 1b</td><td>0.1</td><td>0 (0)</td><td>0.59912 (0.62350)</td></tr><tr><td>(0.13460)</td><td>1.0</td><td>0 (0)</td><td>0.13440 (0.13990)</td></tr><tr><td></td><td>(0)</td><td>1.0</td><td>0 (0)</td><td>0.06820 (0.07099)</td></tr><tr><td rowspan="3">0.03125</td><td rowspan="3">0.73162<br>(0.73162)</td><td rowspan="3">0.13398<br>(0.13378)</td><td rowspan="3">3.72980E-9<sup>-</sup> (0)</td><td>7.58412E-9<sup>-</sup></td><td>0.1</td><td>0.58977 (0.61404)</td><td>0 (0)</td></tr><tr><td>(0)</td><td>0.5</td><td>0.13398 (0.13947)</td><td>0 (0)</td></tr><tr><td></td><td>1.0</td><td>0.06809 (0.07088)</td><td>0 (0)</td></tr><tr><td rowspan="3">0</td><td rowspan="3">0 (0)</td><td rowspan="3">3.14159<br>(3.14159)</td><td>0.09422</td><td>Test Problem 1c</td><td>0.1</td><td>3.14159 (3.14159)</td><td>0.41789 (0.41834)</td></tr><tr><td>(0.09055)</td><td>0.5</td><td>3.14159 (3.14159)</td><td>0.09422 (0.09421)</td></tr><tr><td></td><td>(0.06665)</td><td>1.0</td><td>3.14159 (3.14159)</td><td>0.04784 (0.04782)</td></tr><tr><td rowspan="3">0.03125</td><td rowspan="3">0 (0)</td><td rowspan="3">3.04541<br>(3.04908)</td><td rowspan="3">0 (0)</td><td>0.06091</td><td>0.1</td><td>2.71501 (2.71421)</td><td>0 (0)</td></tr><tr><td>(0.05893)</td><td>0.5</td><td>3.04541 (3.04534)</td><td>0 (0)</td></tr><tr><td></td><td>1.0</td><td>3.09276 (3.09273)</td><td>0 (0)</td></tr><tr><td rowspan="3">0</td><td rowspan="3">1 (1)</td><td rowspan="3">0 (0)</td><td>0.08955</td><td>Test Problem 1d</td><td>0.1</td><td>0 (0)</td><td>0.25904 (0.26297)</td></tr><tr><td>(0.08266)</td><td>0.5</td><td>0 (0)</td><td>0.08955 (0.09070)</td></tr><tr><td></td><td>(0)</td><td>1.0</td><td>0 (0)</td><td>0.04977 (0.05029)</td></tr><tr><td rowspan="3">32</td><td rowspan="3">0 (0)</td><td rowspan="3">1.70407E-26 (0)</td><td rowspan="3">-1.40874E-43 (0)</td><td>2.72652E-26</td><td>0.1</td><td>2.11249E-27 (0)</td><td>0 (0)</td></tr><tr><td>(0)</td><td>0.5</td><td>1.70407E-26 (0)</td><td>0 (0)</td></tr><tr><td></td><td>1.0</td><td>9.07856E-16 (0)</td><td>0 (0)</td></tr><tr><td rowspan="3">0</td><td rowspan="3">1 (1)</td><td rowspan="3">-9.57183E-11<sup>-</sup> (0)</td><td>0.98182</td><td>Test Problem 1e</td><td>0.1</td><td>0 (0)</td><td>1.78909 (1.93322)</td></tr><tr><td>(0.97845)</td><td>0.5</td><td>0 (0)</td><td>0.98182 (1.02732)</td></tr><tr><td></td><td>(0)</td><td>1.0</td><td>0 (0)</td><td>0.79091 (0.79720)</td></tr><tr><td rowspan="3">32</td><td rowspan="3">0 (0)</td><td rowspan="3">0.01818<br>(0.02155)</td><td rowspan="3">-9.57183E-11<sup>-</sup> (0)</td><td>3.63636E-11<sup>-</sup></td><td>0.1</td><td>0.01091 (0.01164)</td><td>0 (0)</td></tr><tr><td>(0)</td><td>0.5</td><td>0.01818 (0.01878)</td><td>0 (0)</td></tr><tr><td></td><td>1.0</td><td>0.02727 (0.02713)</td><td>0 (0)</td></tr><tr><td rowspan="3">0</td><td rowspan="3">0 (0)</td><td rowspan="3">3.14159<br>(3.14159)</td><td>2.57039</td><td>Test Problem 1f</td><td>0.1</td><td>3.14159 (3.14159)</td><td>2.77199 (2.75678)</td></tr><tr><td>(2.49618)</td><td>0.5</td><td>3.14159 (3.14159)</td><td>2.57039 (2.56083)</td></tr><tr><td></td><td>(0)</td><td>1.0</td><td>3.14159 (3.14159)</td><td>2.35619 (2.36472)</td></tr><tr><td rowspan="3">32</td><td rowspan="3">0 (0)</td><td rowspan="3">0.00173<br>(0.00460)</td><td rowspan="3">0 (0)</td><td>3.45120E-5</td><td>0.1</td><td>0.00103 (0.00245)</td><td>0 (0)</td></tr><tr><td>(7.93633E-5)</td><td>0.5</td><td>0.00173 (0.00397)</td><td>0 (0)</td></tr><tr><td></td><td>1.0</td><td>0.00267 (0.00587)</td><td>0 (0)</td></tr></tbody></table>

single-scattering scalar RTE is (assuming non-refractive geometry):

$$ \frac{dI^{\uparrow}(x, \Omega)}{dx} = -I^{\uparrow}(x, \Omega) + \frac{\omega_n}{4\pi} A_n(x) \Phi_n^{\uparrow}(\Omega, \Omega_0) F_0. \quad (51) $$

Here, $x$ is the optical depth coordinate along the line of sight, $\Omega_0 = (-\mu_0, \varphi_0)$ is the solar beam direction, $\Phi_n^{\uparrow}(\Omega, \Omega_0)$ the scattering phase function, and $A_n(x)$ the solar beam attenuation to the point of scatter. $F_0$ and $\omega_n$ are the solar flux and layer single scattering albedo, respectively, as defined in Section 2. Note that, although the directional dyads $\{\Omega, \Omega_0\}$ will change along the line of sight, the angle of scatter remains constant in non-refractive geometry.

Layers are optically uniform, so $\omega_n$ and $\Phi_n^{\uparrow}(\Omega, \Omega_0)$ are constant for a given layer $n$, and $x = \rho_n s$ in terms of the (constant) layer extinction coefficient $\rho_n$, and path distance $s$. In terms of the layer vertical height difference $h_n$ and layer vertical optical depth $\Delta_n$, we have $\rho_n = \Delta_n/h_n$. Using path distance as the variable (more convenient, as this is independent of wavelength), we integrate along the line of sight from $s=0$ to $s=D_n$ (at the top boundary of the

layer):

$$ \begin{aligned} I^{\uparrow}(D_n, \Omega) &= I^{\uparrow}(0, \Omega)e^{-\rho_n D_n} + S_n^{\uparrow}; \quad S_n^{\uparrow} \\ &= \frac{\sigma_n}{4\pi} \Phi_n^{\uparrow}(\Omega, \Omega_0) F_0 \int_0^{D_n} A_n(s) e^{-\rho_n (D_n - s)} ds. \end{aligned} \quad (52) $$

This defines the upwelling layer source term $S_n^{\uparrow}$. In addition, we have $\sigma_n = \omega_n \rho_n$ (the scattering coefficient), and the solar-beam attenuation $A_n(s)$ is found by ray-tracing:

$$ A_n(s) = \exp \left[ - \sum_{k=1}^{n} d_{nk}(s) \rho_k \right], \quad (53) $$

where, $d_{nk}(s)$ are solar path distances through whole layers $k=1, \dots, n-1$ and $d_{nn}(s)$ is the partial distance in layer $k=n$ to the point of scatter.

To obtain the upwelling SS intensity at all levels in the atmosphere, Eq. (52) is applied recursively, starting at the surface; if we write $I_n^{\uparrow}$ to indicate the upwelling intensity at level boundary $n$ (lower boundary of the layer), then for instance at the TOA, we

<footnote>* Deviation from the exact value of 0 is a consequence of dithering the single scattering albedo for conservative scattering, and not due to the 2S approximation for MS.</footnote>
<figure><img src="image_12.png" /><figcaption>Fig. 9. Same as Fig. 3, but for Case 4 from Table 5.</figcaption></figure>

find:

$$I_0^\uparrow(\Omega) = C_N^\uparrow I_{\text{Surf}}^\uparrow(\Omega) + \sum_{k=1}^N C_{k-1}^\uparrow S_k^\uparrow(\Omega); \quad (54)$$

$$C_n^\uparrow = \prod_{k=1}^{n} e^{-\rho_k D_k}, \quad n = 1, \dots, N; \quad C_0^\uparrow = 1. \quad (55)$$

$\{C_n^\uparrow\}$, $n = 0, 1, \dots, N$ are the cumulative (upwelling) line of sight transmittances. Similar expressions pertain for output at other levels. The ESS model for the upwelling field does include the direct-bounce surface reflectance of the solar beam:

$$I_{\text{Surf}}^{\uparrow}(\Omega) = \frac{1}{\pi} A_{\text{Surf}} \mu_0 R(\Omega, \Omega_0) F_0, \quad (56)$$

with $A_{\text{Surf}}$ the solar beam attenuation to the surface, and $R(\Omega, \Omega_0)$ the bidirectional reflectance function of the surface.

For completeness, we also present the treatment for downwelling radiation. By analogy with Eqs. (51) and (52), we have:

$$-\frac{dI^{\downarrow}(x, \Omega)}{dx} = -I^{\downarrow}(x, \Omega) + \frac{\omega_n}{4\pi} A_n(x) \Phi_n^{\downarrow}(\Omega, \Omega_0) F_0. \quad (57)$$

$$I^{\downarrow}(D_n, \Omega) = I^{\downarrow}(0, \Omega)e^{-\rho_n D_n} + S_n^{\downarrow}; \quad S_n^{\downarrow} = \frac{\sigma_n}{4\pi} \Phi_n^{\downarrow}(\Omega, \Omega_0) F_0 \int_0^{D_n} A_n(s) e^{-\rho_n s} ds. \quad (58)$$

Here, the optical depth coordinate $x$ and the distance coordinate $s$ are measured from the top of the layer downwards. There

is also a cumulative recursion: if $I_n^{\downarrow}$ indicates the downwelling intensity at level boundary $n$, then the intensity at the bottom of the atmosphere is:

$$I_N^{\downarrow}(\Omega) = C_0^{\downarrow} I_0^{\downarrow}(\Omega) + \sum_{k=1}^{N} C_k^{\downarrow} S_k^{\downarrow}(\Omega). \quad (59)$$

$$C_n^{\downarrow} = \prod_{k=n+1}^{N} e^{-\rho_k D_k}, \quad n = 0, \dots, N-1; \quad C_N^{\downarrow} = 1. \quad (60)$$

$\{C_n^{\downarrow}\}$, $n = 0, 1, \dots, N$ are the cumulative (downwelling) line of sight transmittances.

In ESS, there are three options for computing path attenuations. The first option assumes the atmosphere is plane-parallel. In this case we write $s = z/\mu$ in terms of the vertical height difference $z$ from layer-top and line of sight cosine $\mu$, and then $A_n(s) = T_n \exp[-\rho_n z/\mu_0]$, where $T_n$ is the solar beam transmittance to layer-top. Then we obtain the familiar results for the upwelling and downwelling integrations in Eqs. (52) and (58):

$$\int_0^{D_n} A_n(s) e^{-\rho_n (D_n - s)} ds = \frac{T_n}{\mu \rho_n} \left[ \frac{1 - \exp(-\Delta_n \kappa_n)}{\kappa_n} \right]; \quad \kappa_n = \frac{1}{\mu_0} + \frac{1}{\mu}; \quad (61a)$$

$$\int_0^{D_n} A_n(s) e^{-\rho_n s} ds = \frac{T_n}{\mu \rho_n} \left[ \frac{1 - \exp(-\Delta_n \kappa_n)}{\kappa_n} \right]; \quad \kappa_n = \frac{1}{\mu_0} - \frac{1}{\mu}. \quad (61b)$$
<figure><img src="image_13.png" /><figcaption>Fig. 10. Same as Fig. 3, but for Case 5 from Table 5.</figcaption></figure>

The second option allows for curved ray-tracing of the solar beam, but scattering is still plane-parallel (this is the pseudospherical approximation). Thus, we still have $s = z/\mu$, so we could work with vertical coordinates. In the average-secant parameterization used for solar beam attenuation, $A_n(s) = T_n \exp[-\eta_n \rho_n z]$, where the average secant is given by $\eta_n = \ln[T_n/T_{n+1}]/\rho_n h_n$. In this case, the results in Eq. (61) still apply, only this time with $\kappa_n = \eta_n \pm \frac{1}{\mu}$.

The third choice for ESS is the fully-spherical situation in which both incoming and outgoing paths are treated for spherical geometry. In this case, ESS will evaluate attenuation integrals numerically, using a quadrature scheme:

$$ \int_{0}^{D_n} A_n(s) e^{-\rho_n(D_n-s)} ds \approx \sum_{j=1}^{N_q} c_j A_n(s_j) e^{-\rho_n(D_n-s_j)}; \quad (62a) $$

$$ \int_{0}^{D_n} A_n(s) e^{-\rho_n s} ds \approx \sum_{j=1}^{N_q} c_j A_n(s_j) e^{-\rho_n s_j}. \quad (62b) $$

Here, the quadrature is $\{s_j, c_j\}_{j=1, \dots, N_q}$ taken over the interval $[0, D_n]$. We have found that Gauss-Legendre integration yields accuracy at the $10^{-6}$ level in most situations, with $N_q$ varying from 5 to 12. However, there are occasions where care must be taken with this procedure. For an optically thick layer below clear skies, solar beam attenuation through such a layer may become vanishingly small. There is then no light reaching the lowest part of the LOS path, and quadrature over interval $[0, D_n]$ will

not work. The remedy for this is to introduce a cutoff distance $D_{nc}$ such that quadrature in Eq. (62) is over the interval $[D_{nc}, D_n]$. The cutoff is determined as follows: assuming that the layer attenuation $\exp[-d_{nn}(s)\rho_n]$ never falls below a certain threshold value $A_{crit}$, we use ray-tracing and numerical solution methods to solve for $D_{nc}$ through the equation $A_{crit} = \exp[-d_{nn}(D_{nc})\rho_n]$. We have found that $A_{crit} = 10^{-8}$ is a suitable threshold to maintain accuracy.

In ESS, the situation with direct thermal emission is simpler. The scalar RTE is then:

$$ \frac{dI^{\uparrow}(x, \Omega)}{dx} = -I^{\uparrow}(x, \Omega) + (1 - \omega_n)B_n(x). \quad (63) $$

Here, $x$ is now the vertical optical depth coordinate measured from the bottom of the layer, and $B_n(x)$ the Planck Function for temperatures in layer $n$. For the upwelling field, the solution corresponding to Eq. (63) is:

$$ I^{\uparrow}(D_n, \Omega) = I^{\uparrow}(0, \Omega)e^{-\rho_n D_n} + (1 - \omega_n) \int_{0}^{D_n} B_n(x)e^{-\rho_n(D_n-s)} ds. \quad (64) $$

There are two choices. The first is the plane-parallel option, with no sphericity along the line of sight; as before we have $x = -\rho_n z$; $s = z/\mu$. Using the linear parameterization $B_n(x) = a_n + b_n s$ in terms of the distance coordinate $s$, then Eq. (64) is then easy to integrate explicitly:

$$ \int_{0}^{D_n} B_n(s) e^{-\rho_n(D_n-s)} ds = \frac{(\rho_n a_n - b_n)(1 - e^{-\rho_n D_n})}{\rho_n^2} + \frac{b_n D_n}{\rho_n}. \quad (65) $$
<figure><img src="image_14.png" /><figcaption>Fig. 11. Same as Fig. 10, but for downwelling results.</figcaption></figure>

In the second choice, the line of sight passes through a curved atmosphere, and we must again use the distance quadrature as in Eq. (62) to perform the integral in Eq. (64).

Building the field recursively using Eq. (54) also applies in the thermal case, except that the surface term is now given by $I_{surf}^{\uparrow}(\Omega) = B_{surf}E(\mu)$, with emissivity $E(\mu)$ and surface Planck function $B_{surf}$.

Downwelling treatment for the thermal radiation unfolds in an analogous way. We have the following:

$$-\frac{dI^{\downarrow}(x, \Omega)}{dx} = -I^{\downarrow}(x, \Omega) + (1 - \omega_n)B_n(x). \qquad (66)$$

$$I^{\downarrow}(D_n, \Omega) = I^{\downarrow}(0, \Omega)e^{-\rho_n D_n} + (1 - \omega_n) \int_{0}^{D_n} B_n(x)e^{-\rho_n s} ds. \qquad (67)$$

Here, we take the optical depth coordinate $x$ and the distance coordinate $s$ to be measured from the top of the layer downwards; thus we now have $x = \mu \rho_n s$, and the Planck function linear param- eterization is now $B_n(x) = a'_n + b'_n s$ with coefficients $\{a'_n, b'_n\}$ simi- lar to $\{a_n, b_n\}$ in Eq. (65).

## 4. Accuracy of 2S-ESS: slab problems

We first evaluate the accuracy of the 2S-ESS model for a series of slab problems. Calculations were performed for a wide range of solar and viewing zenith angles, optical depths, single scattering albedos, phase function asymmetry parameters, and Lambertian

surface albedos. Table 3 lists the values assumed for these parame- ters. There are a total of 1,452,000 parameter combinations. Refer- ence calculations were performed using LIDORT with 64 discrete- ordinate streams in the half-space. Comparisons were also made with a “Full-2S” model, in which the complete radiance field (both SS and MS) is calculated using the 2S approximation. The atmo- sphere was assumed to be plane-parallel. The incoming solar flux $F_0$ was assumed to be unity.

The results are shown in Fig. 1. For both upwelling at the TOA and downwelling at the surface, 2S-ESS provides more accurate re- sults than Full-2S. In the Full-2S model, the SS contributions are based on severely truncated phase function information; the 2S- ESS model replaces these terms by accurate SS RT calculations us- ing the full phase function. Further, upwelling Full-2S results are more accurate than those for downwelling. This is because the ra- diation field is more isotropic for upwelling as a consequence of the presence of a Lambertian surface. It is appropriate here to com- pare the Pearson correlation coefficient ($r$) and the coefficient of determination ($R^2$). The former measures how close the predicted data is to the line of best fit with the observations. On the other hand, the latter measures how well the predicted values match the observed values. Thus, the coefficient of determination is a better indication of the accuracy of the model. Fig. 2 illustrates the accu- racy of the 2S-ESS model even more clearly. Table 4 lists the Earth Mover's Distance (EMD) values for the 2S-ESS and Full-2S mod- els in comparison with LIDORT. EMD is a metric to compare prob- ability distribution functions (PDFs). The EMD between two PDFs
<caption>**Table 9**</caption>

DISORT Test Problem 2 (Rayleigh Scattering, Beam Source). The exact values (shown in parentheses) are based on LIDORT calculations using 16 half-space streams. For conservative scattering, the flux divergence is identically zero, as calculated using Eq. (71).

<table><thead><tr><th>τ</th><th>F<sup>+</sup><sub>dir</sub></th><th>F<sup>+</sup></th><th>F<sup>-</sup></th><th>dF<br>dr</th><th>μ</th><th>I<sup>+</sup></th><th>I<sup>-</sup></th></tr></thead><tbody><tr><td rowspan="3">0</td><td rowspan="3">0.25272<br>(0.25272)</td><td rowspan="3">0 (0)</td><td rowspan="3">0.05286<br>(0.05351)</td><td colspan="4" style="text-align: center;">Test Problem 2a</td></tr><tr><td>1.62365<br>(1.66572)</td><td>0.018014</td><td>0 (0)</td><td>0.16467 (0.16184)</td></tr><tr><td></td><td>0.538263</td><td>0 (0)</td><td>0.02082 (0.02115)</td></tr><tr><td rowspan="3">0.2</td><td rowspan="3">0.02103<br>(0.02103)</td><td rowspan="3">0.04801<br>(0.04417)</td><td rowspan="3">0 (0)</td><td rowspan="3">0.19326<br>(0.18990)</td><td>0.981986</td><td>0 (0)</td><td>0.00770 (0.00787)</td></tr><tr><td>0.018014</td><td>0.02757 (0.02578)</td><td>0 (0)</td></tr><tr><td></td><td>0.538263</td><td>0.01993 (0.02008)</td><td>0 (0)</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.981986</td><td>0.00754 (0.00772)</td><td>0 (0)</td></tr><tr><td rowspan="3">0</td><td rowspan="3">0.25272<br>(0.25272)</td><td rowspan="3">-3.56301E-9<sup>*</sup> (0)</td><td rowspan="3">0.12050<br>(0.13460)</td><td colspan="4" style="text-align: center;">Test Problem 2b</td></tr><tr><td>3.38259E-9<sup>*</sup> (0)</td><td>0.018014</td><td>0 (0)</td><td>0.35895 (0.34786)</td></tr><tr><td></td><td>0.538263</td><td>0 (0)</td><td>0.04657 (0.04872)</td></tr><tr><td rowspan="3">0.2</td><td rowspan="3">0.02103<br>(0.02103)</td><td rowspan="3">0.11119<br>(0.10612)</td><td rowspan="3">-3.56301E-9<sup>*</sup> (0)</td><td rowspan="3">5.45743E-10<sup>*</sup> (0)</td><td>0.981986</td><td>0 (0)</td><td>0.01787 (0.01894)</td></tr><tr><td>0.018014</td><td>0.07508 (0.06784)</td><td>0 (0)</td></tr><tr><td></td><td>0.538263</td><td>0.04506 (0.04641)</td><td>0 (0)</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.981986</td><td>0.01749 (0.01861)</td><td>0 (0)</td></tr><tr><td rowspan="3">0</td><td rowspan="3">0.25272<br>(0.25272)</td><td rowspan="3">0 (0)</td><td rowspan="3">0.06310<br>(0.06247)</td><td colspan="4" style="text-align: center;">Test Problem 2c</td></tr><tr><td>1.63390<br>(1.67463)</td><td>0.018014</td><td>0 (0)</td><td>0.16548 (0.16261)</td></tr><tr><td></td><td>0.538263</td><td>0 (0)</td><td>0.02442 (0.02458)</td></tr><tr><td rowspan="3">5</td><td rowspan="3">2.56077E-28<br>(2.56077E-28)</td><td rowspan="3">6.03691E-5<br>(2.51685E-4)</td><td rowspan="3">0 (0)</td><td rowspan="3">6.03691E-5<br>(1.75465E-4)</td><td>0.981986</td><td>0 (0)</td><td>0.01016 (0.01015)</td></tr><tr><td>0.018014</td><td>4.72421E-6 (1.32511E-5)</td><td>0 (0)</td></tr><tr><td></td><td>0.538263</td><td>2.51169E-5 (3.97171E-5)</td><td>0 (0)</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.981986</td><td>1.59508E-4 (1.70006E-4)</td><td>0 (0)</td></tr><tr><td rowspan="3">0</td><td rowspan="3">0.25272<br>(0.25272)</td><td rowspan="3">-8.03962E-10<sup>*</sup> (0)</td><td rowspan="3">0.22775<br>(0.22591)</td><td colspan="4" style="text-align: center;">Test Problem 2d</td></tr><tr><td>3.59709E-9<sup>*</sup> (0)</td><td>0.018014</td><td>0 (0)</td><td>0.37668 (0.36419)</td></tr><tr><td></td><td>0.538263</td><td>0 (0)</td><td>0.08176 (0.08270)</td></tr><tr><td rowspan="3">5</td><td rowspan="3">2.56077E-28<br>(2.56077E-28)</td><td rowspan="3">0.024968<br>(0.026801)</td><td rowspan="3">-8.04319E-10<sup>*</sup> (0)</td><td rowspan="3">4.99360E-11<sup>*</sup> (0)</td><td>0.981986</td><td>0 (0)</td><td>0.04998 (0.04924)</td></tr><tr><td>0.018014</td><td>0.00412 (0.00380)</td><td>0 (0)</td></tr><tr><td></td><td>0.538263</td><td>0.00826 (0.00769)</td><td>0 (0)</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>0.981986</td><td>0.01147 (0.01060)</td><td>0 (0)</td></tr></tbody></table>

<caption>**Table 10**</caption>

DISORT Test Problem 3 (Henyey-Greenstein Scattering). The exact values (shown in parentheses) are taken from Table 37 of Ref. [38], except for the flux divergence where they are identically zero, as calculated using Eq. (71). The exact value of $F^+$ for these problems is obtained by subtracting $F^+_{dir}$, as calculated using Eq. (70), from the total downwelling flux tabulated in Ref. [38]. Intensities are multiplied by $\pi$ to be consistent with the definition in Ref. [38].

<table><thead><tr><th>τ</th><th>F<sup>+</sup><sub>dir</sub></th><th>F<sup>+</sup></th><th>F<sup>-</sup></th><th>dF<br>dr</th><th>μ</th><th>I<sup>+</sup></th><th>I<sup>-</sup></th></tr></thead><tbody><tr><td rowspan="3">0</td><td rowspan="3">1 (1)</td><td rowspan="3">-3.25192E-10<sup>*</sup> (0)</td><td rowspan="3">0.09230<br>(0.07874)</td><td colspan="4" style="text-align: center;">Test Problem 3a</td></tr><tr><td>1.18460E-9<sup>*</sup> (0)</td><td>0.1</td><td>0 (0)</td><td>0.17462 (0.15137)</td></tr><tr><td></td><td>0.5</td><td>0 (0)</td><td>0.08343 (0.10120)</td></tr><tr><td rowspan="3">1</td><td rowspan="3">0.36788<br>(0.36788)</td><td rowspan="3">0.53982<br>(0.55338)</td><td rowspan="3">-3.24835E-10<sup>*</sup> (0)</td><td rowspan="3">1.16975E-9<sup>*</sup> (0)</td><td>1.0</td><td>0 (0)</td><td>0.03880 (0.03909)</td></tr><tr><td>0.1</td><td>0.23188 (0.21380)</td><td>0 (0)</td></tr><tr><td></td><td>0.5</td><td>0.23032 (0.26663)</td><td>0 (0)</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>1.0</td><td>4.56915 (3.06524)</td><td>0 (0)</td></tr><tr><td rowspan="3">0</td><td rowspan="3">1 (1)</td><td rowspan="3">-1.53578E-10<sup>*</sup> (0)</td><td rowspan="3">0.56003<br>(0.50641)</td><td colspan="4" style="text-align: center;">Test Problem 3b</td></tr><tr><td>2.12006E-9<sup>*</sup> (0)</td><td>0.1</td><td>0 (0)</td><td>0.45756 (0.37997)</td></tr><tr><td></td><td>0.5</td><td>0 (0)</td><td>0.54790 (0.51971)</td></tr><tr><td rowspan="3">8</td><td rowspan="3">3.35463E-4<br>(3.35463E-4)</td><td rowspan="3">0.43964<br>(0.49325)</td><td rowspan="3">-1.52864E-10<sup>*</sup> (0)</td><td rowspan="3">8.49744E-10<sup>*</sup> (0)</td><td>1.0</td><td>0 (0)</td><td>0.53056 (0.49270)</td></tr><tr><td>0.1</td><td>0.24945 (0.23639)</td><td>0 (0)</td></tr><tr><td></td><td>0.5</td><td>0.40571 (0.42235)</td><td>0 (0)</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>1.0</td><td>2.17451 (0.67002)</td><td>0 (0)</td></tr></tbody></table>

can be thought of as the minimal work needed to transform one into the other. Two perfectly matching PDFs will have zero EMD; EMD increases with increasing dissimilarity between the two PDFs. Therefore, a 7–12 times smaller EMD, as indicated by Table 4, implies that, over the range of scenarios considered, the 2S-ESS results are 7–12 times closer to the LIDORT results than those from the Full-2S model.

Figs. 3–11 show comparative radiance results (for upwelling at the TOA and downwelling at the surface) for the 2S-ESS and Full-2S models with respect to LIDORT for the five scenarios described in Table 5. These scenarios encompass a range of gaseous absorption and a variety of scattering particulate types. Case 1 is a scenario with optically thin scattering, such that SS provides the dominant contribution to the radiance field. At the same time, the scat-

tering is fairly isotropic (low value of the asymmetry parameter g), so the phase function does not require many expansion coefficients to simulate accurately. Consequently, both the 2S-ESS and Full-2S models work well for this case, although the former is the more accurate of the two. Case 2 is similar, except that the scattering is much more anisotropic. 2S-ESS still provides accurate results, as it performs an exact SS calculation. On the other hand, Full-2S is unable to reproduce faithfully the scattering phase function, resulting in larger errors. Errors are in general larger for the downwelling field than those for the upwelling radiation, because the presence of a Lambertian surface reduces the anisotropy of the upwelling radiance field. Case 3 involves larger scattering optical depths and, therefore, more significant MS contributions. Both 2S-ESS and Full-2S produce large errors. However, 2S-ESS at least reproduces the

<footnote>* Deviation from the exact value of 0 is a consequence of dithering the single scattering albedo for conservative scattering, and not due to the 2S approximation for MS.</footnote>

<footnote>* Deviation from the exact value of 0 is a consequence of dithering the single scattering albedo for conservative scattering, and not due to the 2S approximation for MS.</footnote>
<caption>**Table 11**</caption>

DISORT Test Problem 4 (Haze-L Scattering, Beam Source). The exact values (shown in parentheses) are taken from Tables 12–16 and Table 22 of Ref. [39], except for the flux divergence where they are based on LIDORT calculations using 16 half-space streams. For conservative scattering, the flux divergence is identically zero, as calculated using Eq. (71). The exact value of F⁺ for these problems is obtained by subtracting F⁺ₗ, as calculated using Eq. (70), from the total downwelling flux tabulated in Ref. [39].

<table><thead><tr><th>τ</th><th>F<sup>+</sup><sub>dir</sub></th><th>F<sup>+</sup></th><th>F<sup>-</sup></th><th>dF/dτ</th><th>μ</th><th>I<sup>+</sup></th><th>I<sup>-</sup></th></tr></thead><tbody><tr><td colspan="8">Test Problem 4a</td></tr><tr><td rowspan="3">0</td><td rowspan="3">3.14159<br>(3.14159)</td><td rowspan="3">0 (0)</td><td rowspan="3">2.19464E-1<br>(1.73223E-1)</td><td rowspan="3">3.58052E-8<sup>*</sup> (0)</td><td>0.1</td><td>0 (0)</td><td>1.21702E-1 (3.61452E-2)</td></tr><tr><td>0.5</td><td>0 (0)</td><td>5.66859E-2 (6.60942E-2)</td></tr><tr><td>1.0</td><td>0 (0)</td><td>3.87987E-2 (9.29328E-2)</td></tr><tr><td rowspan="3">0.5</td><td rowspan="3">1.90547<br>(1.90547)</td><td rowspan="3">1.14705<br>(1.17401)</td><td rowspan="3">1.30390E-1<br>(1.11113E-1)</td><td rowspan="3">3.73373E-8<sup>*</sup> (0)</td><td>0.1</td><td>1.32194E-1 (1.35096E-1)</td><td>1.23318E-1 (1.76286E-2)</td></tr><tr><td>0.5</td><td>1.02873E-1 (1.19079E-1)</td><td>3.45556E-2 (4.02740E-2)</td></tr><tr><td>1.0</td><td>3.19717 (2.52255)</td><td>1.90276E-2 (1.24037E-1)</td></tr><tr><td rowspan="3">1</td><td rowspan="3">1.15573<br>(1.15573)</td><td rowspan="3">1.76640<br>(1.81264)</td><td rowspan="3">8.92894E-14<sup>*</sup> (0)</td><td rowspan="3">3.63904E-8<sup>*</sup> (0)</td><td>0.1</td><td>1.76497E-1 (1.56976E-1)</td><td>0 (0)</td></tr><tr><td>0.5</td><td>1.80539E-1 (2.19633E-1)</td><td>0 (0)</td></tr><tr><td>1.0</td><td>5.37349 (3.38091)</td><td>0 (0)</td></tr><tr><td colspan="8">Test Problem 4b</td></tr><tr><td rowspan="3">0</td><td rowspan="3">3.14159<br>(3.14159)</td><td rowspan="3">0 (0)</td><td rowspan="3">1.65964E-1<br>(1.23665E-1)</td><td rowspan="3">3.47352E-1<br>(3.43724E-1)</td><td>0.1</td><td>0 (0)</td><td>8.31380E-2 (6.55296E-2)</td></tr><tr><td>0.5</td><td>0 (0)</td><td>4.18697E-2 (4.56379E-2)</td></tr><tr><td>1.0</td><td>0 (0)</td><td>3.03338E-2 (2.79717E-2)</td></tr><tr><td rowspan="3">0.5</td><td rowspan="3">1.90547<br>(1.90547)</td><td rowspan="3">9.95297E-1<br>(1.01518)</td><td rowspan="3">9.93729E-2<br>(7.88688E-2)</td><td rowspan="3">3.45187E-1<br>(3.52390E-1)</td><td>0.1</td><td>9.69072E-2 (9.62543E-2)</td><td>8.56143E-2 (8.44628E-2)</td></tr><tr><td>0.5</td><td>8.56026E-2 (9.64183E-2)</td><td>2.58357E-2 (2.80086E-2)</td></tr><tr><td>1.0</td><td>2.78460 (2.24839)</td><td>1.49711E-2 (1.37519E-2)</td></tr><tr><td rowspan="3">1</td><td rowspan="3">1.15573<br>(1.15573)</td><td rowspan="3">1.47852<br>(1.51554)</td><td rowspan="3">-2.17992E-17 (0)</td><td rowspan="3">3.20125E-1<br>(3.19450E-1)</td><td>0.1</td><td>1.18598E-1 (1.08189E-1)</td><td>0 (0)</td></tr><tr><td>0.5</td><td>1.41469E-1 (1.67516E-1)</td><td>0 (0)</td></tr><tr><td>1.0</td><td>4.52924 (2.97766)</td><td>0 (0)</td></tr><tr><td colspan="8">Test Problem 4c</td></tr><tr><td rowspan="10">0</td><td rowspan="10">1.57080<br>(1.57080)</td><td rowspan="10">0 (0)</td><td rowspan="10">2.28960E-1<br>(2.25487E-1)</td><td rowspan="10">3.59951E-1<br>(3.85003E-1)</td><td>0.1<sup>1</sup></td><td>0 (0)</td><td>1.05786 (8.70325E-1)</td></tr><tr><td>0.5<sup>1</sup></td><td>0 (0)</td><td>1.25757E-1 (2.24768E-1)</td></tr><tr><td>1.0<sup>1</sup></td><td>0 (0)</td><td>2.09348E-2 (2.28190E-2)</td></tr><tr><td>0.1<sup>2</sup></td><td>0 (0)</td><td>8.32625E-2 (8.86327E-2)</td></tr><tr><td>0.5<sup>2</sup></td><td>0 (0)</td><td>3.94426E-2 (5.76960E-2)</td></tr><tr><td>1.0<sup>2</sup></td><td>0 (0)</td><td>2.09348E-2 (2.28190E-2)</td></tr><tr><td>0.1<sup>3</sup></td><td>0 (0)</td><td>6.96739E-2 (6.95600E-2)</td></tr><tr><td>0.5<sup>3</sup></td><td>0 (0)</td><td>3.86045E-2 (4.99464E-2)</td></tr><tr><td>1.0<sup>3</sup></td><td>0 (0)</td><td>2.09348E-2 (2.28190E-2)</td></tr><tr><td rowspan="10">0.5</td><td rowspan="10">5.77864E-1<br>(5.77864E-1)</td><td rowspan="10">7.16752E-1<br>(7.02764E-1)</td><td rowspan="10">1.12718E-1<br>(1.23848E-1)</td><td rowspan="10">2.81467E-1<br>(3.37317E-1)</td><td>0.1<sup>1</sup></td><td>2.13795 (1.41181)</td><td>6.66601E-1 (6.97514E-1)</td></tr><tr><td>0.5<sup>1</sup></td><td>8.19680 (3.03333)</td><td>6.14757E-2 (1.09068E-1)</td></tr><tr><td>1.0<sup>1</sup></td><td>4.23914E-2 (4.76083E-2)</td><td>8.75750E-3 (9.34719E-3)</td></tr><tr><td>0.1<sup>2</sup></td><td>8.11180E-2 (1.04622E-1)</td><td>6.76277E-2 (9.15413E-2)</td></tr><tr><td>0.5<sup>2</sup></td><td>6.56658E-2 (5.82338E-2)</td><td>2.05250E-2 (2.95696E-2)</td></tr><tr><td>1.0<sup>2</sup></td><td>4.23914E-2 (4.76083E-2)</td><td>8.75750E-3 (9.34719E-3)</td></tr><tr><td>0.1<sup>3</sup></td><td>5.26224E-2 (6.25009E-2)</td><td>5.35418E-2 (5.90302E-2)</td></tr><tr><td>0.5<sup>3</sup></td><td>2.65657E-2 (2.59261E-2)</td><td>1.88246E-2 (2.46595E-2)</td></tr><tr><td>1.0<sup>3</sup></td><td>4.23914E-2 (4.76083E-2)</td><td>8.75750E-3 (9.34719E-3)</td></tr><tr><td rowspan="10">1</td><td rowspan="10">2.12584E-1<br>(2.12584E-1)</td><td rowspan="10">8.46320E-1<br>(8.03294E-1)</td><td rowspan="10">0 (0)</td><td rowspan="10">2.11781E-1<br>(2.16403E-1)</td><td>0.1<sup>1</sup></td><td>1.69433 (8.76470E-1)</td><td>0 (0)</td></tr><tr><td>0.5<sup>1</sup></td><td>9.50815 (2.71055)</td><td>0 (0)</td></tr><tr><td>1.0<sup>1</sup></td><td>7.07344E-2 (8.37579E-2)</td><td>0 (0)</td></tr><tr><td>0.1<sup>2</sup></td><td>8.27138E-2 (8.95907E-2)</td><td>0 (0)</td></tr><tr><td>0.5<sup>2</sup></td><td>8.93454E-2 (9.43192E-2)</td><td>0 (0)</td></tr><tr><td>1.0<sup>2</sup></td><td>7.07344E-2 (8.37579E-2)</td><td>0 (0)</td></tr><tr><td>0.1<sup>3</sup></td><td>5.60230E-2 (4.66801E-2)</td><td>0 (0)</td></tr><tr><td>0.5<sup>3</sup></td><td>4.12978E-2 (3.99910E-2)</td><td>0 (0)</td></tr><tr><td>1.0<sup>3</sup></td><td>7.07344E-2 (8.37579E-2)</td><td>0 (0)</td></tr></tbody></table>

general angular variation, including the hot spots, while Full-2S does not. Case 4 provides an interesting illustration of the utility of 2S-ESS; it has optically thick gaseous absorption and moderately thick but highly anisotropic scattering, and the radiance field is dominated by SS. The Full-2S model gets this contribution completely wrong but the 2S-ESS model is much more accurate. Note that the downwelling radiance field at the surface is essentially zero due to the strong absorption, and is hence not plotted. The scattering is optically very thick in Case 5, and MS dominates the radiance field. Both 2S-ESS and Full-2S results are highly inaccurate, with errors larger than 100% in extreme cases.

The banded radiance distributions in Figs. 6 and 8 are interesting. Note that they occur when the SZA is equal to the VZA, and only for the downwelling radiance; in other words, when looking directly at the solar aureole. This is a consequence of the radiance

being significantly larger in the forward scattering direction compared with other scattering geometries, due to the extremely large asymmetry factor. Note also the maximum for large ω that is absent for ω = 0.001. This occurs because, in the former case, the SS contribution decreases with increasing SZA (due to increased attenuation) while the MS contribution increases up to SZA ~ 70° (due to the increased radiation path length) before decreasing for very large SZA as the attenuation rapidly increases. For small ω, on the other hand, both SS and MS decrease monotonically with increasing SZA. Finally, this banded feature is not seen for Case 1 because the phase function in this case is nearly isotropic, or for Cases 4 and 5 because here, the large optical depth makes the radiance distribution more isotropic.

In general, errors increase as scattering optical depth and/or anisotropy increases. For most scenarios, 2S-ESS provides results

<footnote>* Deviation from the exact value of 0 is a consequence of dithering the single scattering albedo for conservative scattering, and not due to the 2S approximation for MS.</footnote>

<footnote n="1">¹ φ = 0°</footnote>

<footnote n="2">² φ = 90°</footnote>

<footnote n="3">³ φ = 180°</footnote>
<caption>**Table 12**</caption>

DISORT Test Problem 5 (Cloud C₁ Scattering, Beam Source). The exact values (shown in parentheses) are taken from Tables 19, 20 and 22 of Ref. [39], except for the flux divergence where they are based on LIDORT calculations using 24 half-space streams. For conservative scattering, the flux divergence is identically zero, as calculated using Eq. (71). The exact value of F⁺ for these problems is obtained by subtracting Fₐᵦ, as calculated using Eq. (70), from the total downwelling flux tabulated in Ref. [39].

<table><thead><tr><th>τ</th><th>F<sup>+</sup><sub>dir</sub></th><th>F<sup>+</sup></th><th>F<sup>-</sup></th><th>dF/dτ</th><th>μ</th><th>I<sup>+</sup></th><th>I<sup>-</sup></th></tr></thead><tbody><tr><td colspan="8" style="text-align: center; font-style: italic;">Test Problem 5a</td></tr><tr><td rowspan="3">0</td><td rowspan="3">3.14159<br>(3.14159)</td><td rowspan="3">0 (0)</td><td rowspan="3">2.74460<br>(2.66174)</td><td rowspan="3">8.63080E-9<sup>*</sup> (0)</td><td>0.1</td><td>0 (0)</td><td>5.70503E-1 (4.58734E-1)</td></tr><tr><td>0.5</td><td>0 (0)</td><td>8.20285E-1 (7.72606E-1)</td></tr><tr><td>1.0</td><td>0 (0)</td><td>1.29400 (1.06370)</td></tr><tr><td rowspan="3">32</td><td rowspan="3">3.97856E-14<br>(3.97856E-14)</td><td rowspan="3">2.55344<br>(2.24768)</td><td rowspan="3">2.15646<br>(1.76783)</td><td rowspan="3">9.41734E-8<sup>*</sup> (0)</td><td>0.1</td><td>7.61795E-1 (6.50543E-1)</td><td>7.36690E-1 (6.27632E-1)</td></tr><tr><td>0.5</td><td>8.11937E-1 (6.96363E-1)</td><td>6.86380E-1 (5.81811E-1)</td></tr><tr><td>1.0</td><td>1.20995E+1 (7.53633E-1)</td><td>6.23669E-1 (5.24533E-1)</td></tr><tr><td rowspan="3">64</td><td rowspan="3">5.03852E-28<br>(5.03852E-28)</td><td rowspan="3">3.96984E-1<br>(4.79852E-1)</td><td rowspan="3">2.30199E-14<sup>*</sup> (0)</td><td rowspan="3">7.93967E-9<sup>*</sup> (0)</td><td>0.1</td><td>7.58181E-2 (7.20391E-2)</td><td>0 (0)</td></tr><tr><td>0.5</td><td>1.26363E-1 (1.31990E-1)</td><td>0 (0)</td></tr><tr><td>1.0</td><td>2.06822E-1 (1.95231E-1)</td><td>0 (0)</td></tr><tr><td colspan="8" style="text-align: center; font-style: italic;">Test Problem 5b</td></tr><tr><td rowspan="3">3.2</td><td rowspan="3">1.28058E-1<br>(1.28058E-1)</td><td rowspan="3">1.70558<br>(1.74767)</td><td rowspan="3">3.11605E-1<br>(2.70485E-1)</td><td rowspan="3">3.09166E-1<br>(3.10129E-1)</td><td>0.1</td><td>1.20840E-1 (1.36598E-1)</td><td>1.05195E-1 (1.14041E-1)</td></tr><tr><td>0.5</td><td>1.89284E-1 (2.20190E-1)</td><td>8.55910E-2 (8.72627E-2)</td></tr><tr><td>1.0</td><td>3.20094E2 (6.94794E1)</td><td>1.48962E-1 (8.66126E-2)</td></tr><tr><td rowspan="3">12.8</td><td rowspan="3">8.67322E-6<br>(8.67322E-6)</td><td rowspan="3">1.96903E-1<br>(2.33975E-1)</td><td rowspan="3">3.89658E-2<br>(3.74253E-2)</td><td rowspan="3">4.05141E-2<br>(4.52671E-2)</td><td>0.1</td><td>2.09990E-2 (2.65448E-2)</td><td>1.67235E-2 (2.02153E-2)</td></tr><tr><td>0.5</td><td>4.03506E-2 (4.92710E-2)</td><td>1.16477E-2 (1.29658E-2)</td></tr><tr><td>1.0</td><td>9.62796E1 (2.23085E-1)</td><td>1.24898E-2 (9.51105E-3)</td></tr><tr><td rowspan="3">48</td><td rowspan="3">4.47729E-21<br>(4.47729E-21)</td><td rowspan="3">1.09530E-5<br>(6.38345E-5)</td><td rowspan="3">2.34810E-6<br>(1.02904E-5)</td><td rowspan="3">2.49380E-6<br>(1.25021E-5)</td><td>0.1</td><td>1.43159E-6 (7.47040E-6)</td><td>1.10851E-6 (5.65602E-6)</td></tr><tr><td>0.5</td><td>2.92801E-6 (1.39916E-5)</td><td>7.28550E-7 (3.58246E-6)</td></tr><tr><td>1.0</td><td>9.60374E-3 (3.41286E-5)</td><td>5.79338E-7 (2.57858E-6)</td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Example TOA radiance spectra in the shortwave infrared (top) and thermal infrared (bottom). The radiance is in units of W/m²/sr/cm⁻¹.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 13.** Histogram of root mean square relative residual errors (RMSE) for the 2S-ESS model (compared to LIDORT with 16 half-space streams) in the (left) shortwave infrared and (right) thermal infrared.</figcaption></figure>

<footnote>* Deviation from the exact value of 0 is a consequence of dithering the single scattering albedo for conservative scattering, and not due to the 2S approximation for MS.</footnote>
<caption>**Table 13**</caption>

DISORT Test Problem 6 (No Scattering, Increasingly Complex Sources). The exact values (shown in parentheses) are based on LIDORT calculations using 16 (problems 6a–6c) and 80 (problems 6d–6h) half-space streams.

<table><thead><tr><th>τ</th><th>F<sup>+</sup><sub>dir</sub></th><th>F<sup>+</sup></th><th>F<sup>-</sup></th><th>dE/dτ</th><th>μ</th><th>I<sup>+</sup></th><th>I<sup>-</sup></th></tr></thead><tbody><tr><td rowspan="2">0</td><td rowspan="2">1 (1)</td><td rowspan="2">0 (0)</td><td rowspan="2">0 (0)</td><td>Test Problem 6a<br>200 (200)</td><td>0.1</td><td>0 (0)</td><td>0 (0)</td></tr><tr><td></td><td>1.0</td><td>0 (0)</td><td>0 (0)</td></tr><tr><td rowspan="2">0.5</td><td rowspan="2">100 (100)</td><td rowspan="2">0 (0)</td><td rowspan="2">0 (0)</td><td>Test Problem 6b<br>200 (200)</td><td>0.1</td><td>0 (0)</td><td>0 (0)</td></tr><tr><td></td><td>1.0</td><td>0 (0)</td><td>0 (0)</td></tr><tr><td rowspan="2">1</td><td>3.67879E1<br>(3.67879E1)</td><td>0 (0)</td><td>0 (0)</td><td>7.35759E1<br>(7.35759E1)</td><td>0.1</td><td>0 (0)</td><td>0 (0)</td></tr><tr><td>1.35335E1<br>(1.35335E1)</td><td>0 (0)</td><td>0 (0)</td><td>2.70671E1<br>(2.70671E1)</td><td>0.1</td><td>0 (0)</td><td>0 (0)</td></tr><tr><td rowspan="2">0</td><td rowspan="2">100 (100)</td><td rowspan="2">0 (0)</td><td rowspan="2">9.15782E-1<br>(1.48452)</td><td>Test Problem 6c<br>2.01832E2<br>(2.02010E2)</td><td>0.1</td><td>0 (0)</td><td>9.77882E-5 (9.77882E-5)</td></tr><tr><td></td><td>1.0</td><td>0 (0)</td><td>7.92386E-1 (7.92386E-1)</td></tr><tr><td rowspan="2">0.5</td><td>3.67879E1<br>(3.67879E1)</td><td>0 (0)</td><td>2.48935<br>(2.99909)</td><td>7.85546E1<br>(7.79965E1)</td><td>0.1</td><td>0 (0)</td><td>1.45131E-2 (1.45131E-2)</td></tr><tr><td>1.35335E1<br>(1.35335E1)</td><td>0 (0)</td><td>6.76676<br>(6.76676)</td><td>4.06006E1<br>(4.06006E1)</td><td>0.1</td><td>0 (0)</td><td>1.30642 (1.30642)</td></tr><tr><td rowspan="2">1</td><td rowspan="2">100 (100)</td><td rowspan="2">0 (0)</td><td rowspan="2">4.95269E-1<br>(6.70946E-1)</td><td>Test Problem 6d<br>2.00991E2<br>(2.00936E2)</td><td>0.1</td><td>0 (0)</td><td>6.80075E-5 (6.80075E-5)</td></tr><tr><td></td><td>1.0</td><td>0 (0)</td><td>3.15441E-1 (3.15441E-1)</td></tr><tr><td rowspan="2">0.5</td><td>3.67879E1<br>(3.67879E1)</td><td>0 (0)</td><td>1.34628<br>(1.39125)</td><td>7.62685E1<br>(7.57199E1)</td><td>0.1</td><td>0 (0)</td><td>1.00932E-2 (1.00932E-2)</td></tr><tr><td>1.35335E1<br>(1.35335E1)</td><td>0 (0)</td><td>3.65957<br>(3.31775)</td><td>3.43862E1<br>(3.45341E1)</td><td>0.1</td><td>0 (0)</td><td>5.20074E-1 (5.20074E-1)</td></tr><tr><td rowspan="2">0</td><td rowspan="2">100 (100)</td><td rowspan="2">0 (0)</td><td rowspan="2">4.74155E1<br>(7.95461E1)</td><td>Test Problem 6e<br>2.94831E2<br>(3.07078E2)</td><td>0.1</td><td>0 (0)</td><td>4.53789E-3 (4.53789E-3)</td></tr><tr><td></td><td>1.0</td><td>0 (0)</td><td>4.33773E1 (4.33773E1)</td></tr><tr><td rowspan="2">0.5</td><td>3.67879E1<br>(3.67879E1)</td><td>0 (0)</td><td>1.28889E2<br>(1.59898E2)</td><td>3.31353E2<br>(3.07078E2)</td><td>0.1</td><td>0 (0)</td><td>6.73483E-1 (6.73483E-1)</td></tr><tr><td>1.35335E1<br>(1.35335E1)</td><td>2.62901E-13 (0)</td><td>3.50356E2<br>(3.56409E2)</td><td>7.27779E2<br>(7.17450E2)</td><td>0.1</td><td>0 (0)</td><td>7.15170E1 (7.15170E1)</td></tr><tr><td rowspan="2">1</td><td rowspan="2">100 (100)</td><td rowspan="2">3.21497E2<br>(3.21497E2)</td><td rowspan="2">4.90078E1<br>(8.27883E1)</td><td>Test Problem 6f<br>9.41009E2<br>(9.54517E2)</td><td>0.1</td><td>1.02336E2 (1.02336E2)</td><td>4.76262E-2 (4.80532E-3)</td></tr><tr><td></td><td>1.0</td><td>1.02336E2 (1.02336E2)</td><td>4.43914E1 (4.50201E1)</td></tr><tr><td rowspan="2">0.5</td><td>3.67879E1<br>(3.67879E1)</td><td>1.18272E2<br>(1.42490E2)</td><td>1.33217E2<br>(1.66522E2)</td><td>5.76554E2<br>(5.27110E2)</td><td>0.1</td><td>6.89532E-1 (6.89532E-1)</td><td>7.06836E-1 (7.13172E-1)</td></tr><tr><td>1.35335E1<br>(1.35335E1)</td><td>4.35099E1<br>(7.05312E1)</td><td>3.62121E2<br>(3.71729E2)</td><td>8.38329E2<br>(8.45308E2)</td><td>0.1</td><td>6.20697E1 (6.20697E1)</td><td>7.31891E1 (7.42256E1)</td></tr><tr><td rowspan="2">1</td><td rowspan="2">100 (100)</td><td rowspan="2">3.21497E2<br>(3.21497E2)</td><td rowspan="2">3.23103E2<br>(3.35277E2)</td><td>Test Problem 6g<br>6.03212E2<br>(5.80375E2)</td><td>0.1</td><td>1.02336E2 (1.02336E2)</td><td>7.80735E1 (7.80733E1)</td></tr><tr><td></td><td>1.0</td><td>1.02336E2 (1.02336E2)</td><td>1.16571E2 (1.16440E2)</td></tr><tr><td rowspan="2">0.5</td><td>3.67879E1<br>(3.67879E1)</td><td>3.02026E2<br>(3.04774E2)</td><td>4.12287E2<br>(4.12513E2)</td><td>1.40616E2<br>(1.27083E2)</td><td>0.1</td><td>1.01048E2 (1.01048E2)</td><td>1.15842E2 (1.15819E2)</td></tr><tr><td>1.35335E1<br>(1.35335E1)</td><td>3.70021E2<br>(3.63633E2)</td><td>4.50413E2<br>(4.41067E2)</td><td>-1.69249E2<br>(-1.68107E2)</td><td>0.1</td><td>9.78748E1 (9.78748E1)</td><td>1.35197E2 (1.34982E2)</td></tr><tr><td rowspan="2">1</td><td rowspan="2">100 (100)</td><td rowspan="2">3.21497E2<br>(3.21497E2)</td><td rowspan="2">2.33387E2<br>(2.37350E2)</td><td>Test Problem 6h<br>4.23780E2<br>(4.23780E2)</td><td>0.1</td><td>1.02336E2 (1.02336E2)</td><td>7.12616E1 (7.12616E1)</td></tr><tr><td></td><td>1.0</td><td>1.02336E2 (1.02336E2)</td><td>7.80736E1 (7.80736E1)</td></tr><tr><td rowspan="2">1</td><td>1.35335E1<br>(1.35335E1)</td><td>2.48529E2<br>(2.55455E2)</td><td>2.57167E2<br>(2.61130E2)</td><td>5.73524E1<br>(6.19825E1)</td><td>0.1</td><td>7.73186E1 (7.73186E1)</td><td>7.88310E1 (7.88310E1)</td></tr><tr><td>2.06115E-7<br>(2.06115E-7)</td><td>4.47406E2<br>(4.43444E2)</td><td>4.67678E2<br>(4.55789E2)</td><td>-7.01499<br>(-3.12908E1)</td><td>0.1</td><td>8.49992E1 (8.49992E1)</td><td>8.56424E1 (8.56423E1)</td></tr><tr><td rowspan="2">10</td><td rowspan="2">100 (100)</td><td rowspan="2">3.21497E2<br>(3.21497E2)</td><td rowspan="2">2.33387E2<br>(2.37350E2)</td><td>Test Problem 6i<br>4.23780E2<br>(4.23780E2)</td><td>0.1</td><td>1.02336E2 (1.02336E2)</td><td>7.12616E1 (7.12616E1)</td></tr><tr><td></td><td>1.0</td><td>1.02336E2 (1.02336E2)</td><td>7.80736E1 (7.80736E1)</td></tr></tbody></table>

that are accurate to within a few percent. Overall, the 2S-ESS model is only unsuitable for heavy aerosol loading events (e.g., dust storms, wildfires or biomass burning events) and for scenarios with optically thick water clouds. For the former, 2S-ESS could be used as a starting point for a more accurate RT calculation, such as that obtained using a PCA-based approach [24–32]; for the latter, one can resort to asymptotic theory for analytic solutions (see, e.g., Ref. [33]).

We can further analyze the performance of 2S-ESS and Full-2S for various parameter thresholds. Tables 6 and 7 show error statistics for both models for these scenarios. Clearly, the 2S-ESS model vastly outperforms the Full-2S model.

## 5. Accuracy of 2S-ESS: multi-layer problem

2S-ESS model accuracy can also be evaluated for more realistic multi-layer problems. We use the Weather Research and Forecasting model coupled to Chemistry (WRF-Chem) [34–36] at 4 km spatial resolution to extract the trace gas column mixing ratio and aerosol concentration profiles at 115 atmospheric levels for 17 locations in the continental US at six local times for two days in July 2006. There are five aerosol types (black carbon, organic carbon, accumulation mode sea salt, coarse mode sea salt and sulfate) in varying amounts in each profile. The aerosol treatment is described in detail in Ref. [30]. The scenarios cover a range of SZA (16°–89°)
<caption>**Table 14**</caption>

DISORT Test Problem 7 (Absorbing/Scattering Medium, All Possible Sources). The exact values (shown in parentheses) are based on LIDORT calculations using 16 (problems 7a, 7b), 32 (problems 7c, 7d) and 64 (problem 7e) half-space streams, respectively.

<table><thead><tr><th>τ</th><th>F<sup>+</sup><sub>dir</sub></th><th>F<sup>+</sup></th><th>F<sup>-</sup></th><th>dF/dτ</th><th>μ</th><th>I<sup>+</sup></th><th>I<sup>-</sup></th></tr></thead><tbody><tr><td>0</td><td>0 (0)</td><td>0 (0)</td><td>9.41345E1<br>(8.62931E1)</td><td>-3.50918E1<br>(-5.13727E1)</td><td>-</td><td>-</td><td>-</td></tr><tr><td>1</td><td>0 (0)</td><td>1.34757E2<br>(1.21204E2)</td><td>0 (0)</td><td>-5.51965E2<br>(-5.41035E2)</td><td></td><td></td><td></td></tr><tr><td>0</td><td>0 (0)</td><td>0 (0)</td><td>7.80321E-7<br>(1.10949E-6)</td><td>6.74215E-8<br>(8.23192E-8)</td><td>1.0</td><td>0 (0)</td><td>3.94828E-7 (4.65746E-7)</td></tr><tr><td>100</td><td>0 (0)</td><td>1.86230E-5<br>(2.07786E-5)</td><td>0 (0)</td><td>-5.06754E-6<br>(-5.06469E-6)</td><td>0.1 1.0</td><td>7.21928E-6 (7.52315E-6)</td><td>0 (0)</td></tr><tr><td>0</td><td>100 (100)</td><td>3.19830E2<br>(3.19830E2)</td><td>4.08308E2<br>(4.29562E2)</td><td>-9.04538E1<br>(-8.04201E1)</td><td>0.1<sup>1</sup><br>1.0<sup>1</sup></td><td>1.01805E2 (1.01805E2)<br>1.01805E2 (1.01805E2)</td><td>1.44260E2 (1.43533E2)<br>1.48042E2 (1.48491E2)</td></tr><tr><td>0.5</td><td>3.67879E1<br>(3.67879E1)</td><td>3.59202E2<br>(3.54089E2)</td><td>4.29584E2<br>(4.47036E2)</td><td>2.75553E2<br>(2.51574E2)</td><td>0.1<sup>2</sup><br>1.0<sup>2</sup></td><td>1.01805E2 (1.01805E2)<br>1.01805E2 (1.01805E2)</td><td>1.28040E2 (1.28049E2)<br>1.48042E2 (1.48491E2)</td></tr><tr><td>1</td><td>1.35335E1<br>(1.35335E1)</td><td>3.02551E2<br>(3.01334E2)</td><td>5.94576E2<br>(5.94576E2)</td><td>7.29211E2<br>(7.16085E2)</td><td>0.1<sup>2</sup><br>1.0<sup>2</sup></td><td>1.39768E2 (1.32691E2)<br>1.06594E2 (1.06507E2)<br>1.08586E2 (1.06358E2)<br>1.06594E2 (1.06507E2)</td><td>9.93391E1 (1.03192E2)<br>1.58781E2 (1.58878E2)<br>9.11815E1 (9.42503E1)<br>1.58781E2 (1.58878E2)</td></tr><tr><td>0</td><td>100 (100)</td><td>3.19830E2<br>(3.19830E2)</td><td>3.12496E2<br>(3.12557E2)</td><td>-1.86266E2<br>(-1.68353E2)</td><td>0.1<sup>1</sup><br>1.0<sup>1</sup></td><td>1.01805E2 (1.01805E2)<br>1.01805E2 (1.01805E2)</td><td>1.38954E2 (1.37730E2)<br>9.50946E1 (9.57364E1)</td></tr><tr><td>0.5</td><td>3.67879E1<br>(3.67879E1)</td><td>3.53152E2<br>(3.50544E2)</td><td>2.62063E2<br>(2.68132E2)</td><td>1.01982E2<br>(1.01239E2)</td><td>0.1<sup>2</sup><br>1.0<sup>2</sup></td><td>1.01805E2 (1.01805E2)<br>1.01805E2 (1.01805E2)</td><td>1.22733E2 (1.22246E2)<br>9.50946E1 (9.57364E1)</td></tr><tr><td>1</td><td>1.35335E1<br>(1.35335E1)</td><td>2.88526E2<br>(2.92065E2)</td><td>3.02060E2<br>(3.05598E2)</td><td>4.22669E2<br>(4.09372E2)</td><td>0.1<sup>2</sup><br>1.0<sup>2</sup></td><td>1.34026E2 (1.27253E2)<br>1.06004E2 (1.06128E2)<br>1.02844E2 (1.00920E2)<br>1.06004E2 (1.06128E2)</td><td>8.74508E1 (9.06699E1)<br>8.80010E1 (8.87819E1)<br>7.92932E1 (8.17285E1)<br>8.80010E1 (8.87819E1)</td></tr><tr><td>0</td><td>100 (100)</td><td>3.19830E2<br>(3.19830E2)</td><td>3.88301E2<br>(4.04252E2)</td><td>-1.10462E2<br>(-9.98806E1)</td><td>0.1<sup>1</sup><br>1.0<sup>1</sup></td><td>1.01805E2 (1.01805E2)<br>1.01805E2 (1.01805E2)</td><td>1.43168E2 (1.42202E2)<br>1.37774E2 (1.38021E2)</td></tr><tr><td>0.5</td><td>3.67879E1<br>(3.67879E1)</td><td>3.57939E2<br>(3.53264E2)</td><td>3.94602E2<br>(4.07798E2)</td><td>2.39307E2<br>(2.17334E2)</td><td>0.1<sup>2</sup><br>1.0<sup>2</sup></td><td>1.01805E2 (1.01805E2)<br>1.01805E2 (1.01805E2)</td><td>1.26924E2 (1.26686E2)<br>1.37774E2 (1.38021E2)</td></tr><tr><td>1</td><td>1.35335E1<br>(1.35335E1)</td><td>2.99623E2<br>(2.99001E2)</td><td>5.33492E2<br>(5.29161E2)</td><td>6.65197E2<br>(6.38454E2)</td><td>0.1<sup>2</sup><br>1.0<sup>2</sup></td><td>1.38597E2 (1.31422E2)<br>1.06471E2 (1.06421E2)<br>1.07387E2 (1.05061E2)<br>1.06471E2 (1.06421E2)</td><td>9.67191E1 (1.00117E2)<br>1.45108E2 (1.45074E2)<br>8.84716E1 (9.10763E1)<br>1.45108E2 (1.45074E2)</td></tr></tbody></table>

and VZA (34°–55°). We perform TOA radiance calculations (representing satellite measurements) in two spectral regions: 12,585–13,285 cm⁻<fnref n="1" />, representing a region of solar beam forcing in the shortwave infrared, and 500–700 cm⁻<fnref n="1" />, representing a region of thermal emission. Calculations in these spectral bands are done at resolutions 0.0025 and 0.001 cm⁻<fnref n="1" />, respectively. There are, therefore, 280,000 and 200,000 wavelength points, respectively, for the RT calculations in the solar and thermal regions. Across the profiles, the surface albedo varies from 0.29 to 0.54 in the shortwave infrared and from 0.01 (emissivity = 0.99) to 0.04 (emissivity = 0.96) in the thermal infrared. The aerosol optical depth varies from 0.005 to 0.07 in the shortwave infrared and from 0.0005 to 0.01 in the thermal infrared. The RT calculations are performed using the most realistic treatment of scattering, with spherical treat-

ments for both the incoming solar beam and outgoing line of sight paths.

Fig. 12 shows an example of TOA radiances in the two spectral regions. Fig. 13 shows relative residual errors for the 2S-ESS model (compared to LIDORT results obtained with 16 half-space discrete ordinates). Errors are on the order of 0.1% in the shortwave infrared and 0.01% in the thermal infrared.

## 6. Comparisons with DISORT test problems

The DISORT code package distribution comes with a set of test problems that evaluate various aspects of the RTE solution. This set includes tests for (1) purely scattering slabs with successively more complicated phase functions, (2) purely absorbing scenar-

<footnote n="1">¹ $\varphi = 0^\circ$</footnote>

<footnote n="2">² $\varphi = 90^\circ$</footnote>
<caption>**Table 15**</caption>

DISORT Test Problem 8 (Absorbing/Isotropic Scattering Medium, Two Inhomogeneous Layers). The exact values (shown in parentheses) are based on LIDORT calculations using 32 half-space streams. The upwelling (downwelling) fluxes at the top (bottom) of the atmosphere calculated using LIDORT agree with the data in Table 1 of Ref. [40].

<table><thead><tr><th>τ</th><th>F<sup>+</sup><sub>dir</sub></th><th>F<sup>+</sup></th><th>F<sup>-</sup></th><th>df/dτ</th><th>μ</th><th>I<sup>+</sup></th><th>I<sup>-</sup></th></tr></thead><tbody><tr><td colspan="8" style="text-align: center; font-weight: bold;">Test Problem 8a</td></tr><tr><td>0</td><td>0 (0)</td><td>1 (1)</td><td>1.12557E-1<br>(9.30102E-2)</td><td>1.11256<br>(1.12461)</td><td>0.2</td><td>3.18310E-1 (3.18310E-1)</td><td>5.99538E-2 (5.62281E-2)</td></tr><tr><td>0.25</td><td>0 (0)</td><td>6.94657E-1<br>(7.22031E-1)</td><td>3.51390E-2<br>(2.79270E-2)</td><td>7.29796E-1<br>(6.54013E-1)</td><td>1.0</td><td>3.18310E-1 (3.18310E-1)</td><td>2.10385E-2 (1.94534E-2)</td></tr><tr><td>0.5</td><td>0 (0)</td><td>4.55128E-1<br>(5.12812E-1)</td><td>0 (0)</td><td>6.37180E-1<br>(5.62315E-1)</td><td>0.2</td><td>1.40773E-1 (1.37058E-1)</td><td>2.08488E-2 (1.85044E-2)</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">Test Problem 8b</td></tr><tr><td>0</td><td>0 (0)</td><td>1 (1)</td><td>2.61404E-1<br>(2.25247E-1)</td><td>5.04562E-1<br>(5.12905E-1)</td><td>0.2</td><td>3.18310E-1 (3.18310E-1)</td><td>1.26472E-1 (1.23664E-1)</td></tr><tr><td>0.25</td><td>0 (0)</td><td>7.75316E-1<br>(7.95235E-1)</td><td>1.45372E-1<br>(1.26410E-1)</td><td>3.68275E-1<br>(3.56968E-1)</td><td>1.0</td><td>3.18310E-1 (3.18310E-1)</td><td>5.08931E-2 (4.95641E-2)</td></tr><tr><td>0.5</td><td>0 (0)</td><td>6.10820E-1<br>(6.50263E-1)</td><td>0 (0)</td><td>6.10820E-2<br>(5.67329E-2)</td><td>0.2</td><td>1.86772E-1 (1.84035E-1)</td><td>8.59185E-2 (8.35507E-2)</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">Test Problem 8c</td></tr><tr><td>0</td><td>0 (0)</td><td>1 (1)</td><td>4.15103E-1<br>(3.78577E-1)</td><td>5.66041E-1<br>(5.65092E-1)</td><td>0.2</td><td>3.18310E-1 (3.18310E-1)</td><td>1.55446E-1 (1.49312E-1)</td></tr><tr><td>1</td><td>0 (0)</td><td>4.39030E-1<br>(4.86215E-1)</td><td>2.48746E-1<br>(2.43436E-1)</td><td>2.75111E-1<br>(2.76667E-1)</td><td>1.0</td><td>3.18310E-1 (3.18310E-1)</td><td>1.06662E-1 (1.04770E-1)</td></tr><tr><td>3</td><td>0 (0)</td><td>1.14965E-1<br>(1.59989E-1)</td><td>2.17992E-17 (0)</td><td>1.14965E-2<br>(1.35673E-2)</td><td>0.2</td><td>1.01160E-1 (9.87827E-2)</td><td>9.32271E-2 (9.65271E-2)</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">Test Problem 8d</td></tr><tr><td>0</td><td>0 (0)</td><td>1 (1)</td><td>4.15103E-1<br>(3.78577E-1)</td><td>5.66041E-1<br>(5.65092E-1)</td><td>0.2</td><td>3.18310E-1 (3.18310E-1)</td><td>1.55446E-1 (1.49312E-1)</td></tr><tr><td>1</td><td>0 (0)</td><td>4.39030E-1<br>(4.86215E-1)</td><td>2.48746E-1<br>(2.43436E-1)</td><td>2.75111E-1<br>(2.76667E-1)</td><td>1.0</td><td>3.18310E-1 (3.18310E-1)</td><td>1.06662E-1 (1.04770E-1)</td></tr><tr><td>3</td><td>0 (0)</td><td>1.14965E-1<br>(1.59989E-1)</td><td>2.17992E-17 (0)</td><td>1.14965E-2<br>(1.35673E-2)</td><td>0.2</td><td>1.01160E-1 (9.87827E-2)</td><td>9.32271E-2 (9.65271E-2)</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">Test Problem 8e</td></tr><tr><td>0</td><td>0 (0)</td><td>1 (1)</td><td>4.15103E-1<br>(3.78577E-1)</td><td>5.66041E-1<br>(5.65092E-1)</td><td>0.2</td><td>3.18310E-1 (3.18310E-1)</td><td>1.55446E-1 (1.49312E-1)</td></tr><tr><td>1</td><td>0 (0)</td><td>4.39030E-1<br>(4.86215E-1)</td><td>2.48746E-1<br>(2.43436E-1)</td><td>2.75111E-1<br>(2.76667E-1)</td><td>1.0</td><td>3.18310E-1 (3.18310E-1)</td><td>1.06662E-1 (1.04770E-1)</td></tr><tr><td>3</td><td>0 (0)</td><td>1.14965E-1<br>(1.59989E-1)</td><td>2.17992E-17 (0)</td><td>1.14965E-2<br>(1.35673E-2)</td><td>0.2</td><td>1.01160E-1 (9.87827E-2)</td><td>9.32271E-2 (9.65271E-2)</td></tr></tbody></table>

ios with increasing complexity of boundary and internal sources, (3) scattering/emitting/absorbing slabs with Lambertian and non-Lambertian surface reflectance models, (4) multi-layer media with an isotropic flux at the upper boundary, and (5) inhomogeneous media with multiple layers. Many of these problems also have accurate numerical results tabulated in peer-reviewed literature. A detailed description of the test problems can be found in Ref. [23]. In our work, Tables 8–16 compare 2S-ESS results with published results, if available, and otherwise with those obtained using LIDORT. Plane-parallel geometry is assumed, to be consistent with the prior work being compared against. Appropriate references are supplied in the tables, as well as the number of streams for those cases where LIDORT was employed for the comparisons. Physical considerations are used where appropriate; for example, upwelling fluxes are identically zero at the bottom of the atmosphere when bounded from below by a black surface. Therefore, these are set to zero even if LIDORT returns small values. The diffuse fluxes tabulated below have the following definitions:

$$F^{+}(\tau) = 2\pi \int_{0}^{1} \mu I^{0}(\tau, \mu) d\mu; \qquad (68a)$$

$$F^{-}(\tau) = 2\pi \int_{0}^{1} \mu I^{0}(\tau, -\mu) d\mu. \qquad (68b)$$

In the 2S approximation, these reduce to:

$$F^{+}(\tau) = 2\pi \bar{\mu} \mathbb{I}^{+}; \qquad (69a)$$

$$F^{-}(\tau) = 2\pi \bar{\mu} \mathbb{I}^{-}. \qquad (69b)$$

In addition, the direct downwelling flux can be evaluated as follows:

$$F_{\text{dir}}^{+}(\tau) = \mu_{0}F_{0}e^{-\tau/\mu_{0}}. \qquad (70)$$

The flux divergence can be calculated as [1]:

$$\frac{dF}{d\tau} = 4\pi (1 - \omega) [\bar{I}(\tau) - B(\tau)]. \qquad (71)$$

The mean intensity $\bar{I}$ is defined as:

$$\bar{I}(\tau) = \frac{1}{2} \int_{0}^{1} [I^{0}(\tau, \mu) + I^{0}(\tau, -\mu)] d\mu + \frac{1}{4\pi} F_{0} e^{-\tau/\mu_{0}}, \qquad (72)$$

where, the last term on the right-hand side of Eq. (72) refers to the direct beam contribution to the downward mean intensity.

In the 2S approximation, Eq. (72) reduces to:

$$\bar{I}(\tau) = \frac{1}{2}(\mathbb{I}^{+} + \mathbb{I}^{-}) + \frac{1}{4\pi} F_{0} e^{-\tau/\mu_{0}}. \qquad (73)$$

The results are broadly consistent with those in Sections 4 and 5.

## 7. Computational considerations

The high computational efficiency of the 2S-ESS model is a result of several factors. The ESS calculations do not involve matrix inversion, which is typically the most significant time sink in solving the RTE. The layer homogeneous and particular 2S solutions are obtained analytically. The boundary value problem is solved using penta-diagonal solvers, which are more efficient than matrix routines used for multi-stream scenarios. The computation time for discrete-ordinate solutions increases approximately as the cube of the number of streams used [42]. Table 17 shows the mean computational time for LIDORT, LIDORT (in 2S mode), and 2S-ESS calculations for the scenarios described in Section 5. The calculations were done on a single core of a Dell PowerEdge R620 server with 256 GB of RAM, running Redhat Enterprise Linux. Note that, while the calculation of optical properties such as gas absorption cross sections can be time-consuming, in our comparisons we consider only the times taken by the RT models when determining the radiation fields. The 2S-ESS model provides a mean computational speed-up of 2650 for the collection of solar beam scenarios, and an average speed-up of 2800 for thermal emission scenarios. In addition, 2S-ESS is about 8 and 20 times faster, respectively, than LIDORT in 2S mode for the same scenarios. This result demonstrates the impact of employing analytic formulas for 2S calculations as opposed to the matrix computations performed by LIDORT in 2S mode.

## 8. Differences between 2S-ESS and other established schemes

The 2S RTE is simple enough to be solved analytically, and here we adopt the same techniques as used by other authors, finding
<caption>**Table 16**</caption>

DISORT Test Problem 9 (Emitting/Absorbing/Scattering Medium, Multiple Inhomogeneous Layers). The exact values (shown in parentheses) are based on LIDORT calculations using 32 (problems 9a and 9b) and 64 (problem 9e) half-space streams, respectively. The upwelling (downwelling) fluxes at the top (bottom) of the atmosphere calculated using LIDORT for problems 9a and 9b agree with the data in Tables VI and VII of Ref. [41].

<table><thead><tr><th>τ</th><th>F<sup>+</sup><sub>dir</sub></th><th>F<sup>+</sup></th><th>F<sup>-</sup></th><th>dF/dτ</th><th>μ</th><th>I<sup>+</sup></th><th>I<sup>-</sup></th></tr></thead><tbody><tr><td colspan="8" style="text-align: center; font-weight: bold;">Test Problem 9a</td></tr><tr><td>0</td><td>0 (0)</td><td>1 (1)</td><td>2.60013E-1<br>(2.27952E-1)</td><td>8.82009E-1<br>(8.82116E-1)</td><td>0.2</td><td>3.18310E-1 (3.18310E-1)</td><td>1.05660E-1 (9.98685E-2)</td></tr><tr><td>1.05</td><td>0 (0)</td><td>2.92606E-1<br>(3.55246E-1)</td><td>8.56761E-2<br>(8.75404E-2)</td><td>2.26969E-1<br>(2.32396E-1)</td><td>0.2</td><td>3.18310E-1 (3.18310E-1)</td><td>6.11585E-2 (5.91379E-2)</td></tr><tr><td>2.1</td><td>0 (0)</td><td>9.27844E-2<br>(1.44255E-1)</td><td>2.76523E-2<br>(3.61809E-2)</td><td>7.22620E-2<br>(9.33650E-2)</td><td>0.2</td><td>1.55629E-1 (1.53497E-1)</td><td>2.02160E-2 (2.31986E-2)</td></tr><tr><td>6</td><td>0 (0)</td><td>1.77564E-3<br>(6.71536E-3)</td><td>6.78320E-4<br>(2.19334E-3)</td><td>1.22698E-3<br>(3.92868E-3)</td><td>0.2</td><td>1.71323E-2 (2.09208E-2)</td><td>1.10484E-2 (1.48554E-2)</td></tr><tr><td>21</td><td>0 (0)</td><td>1.81817E-8<br>(6.17135E-7)</td><td>-5.19733E-24 (0)</td><td>3.63634E-9<br>(1.02525E-7)</td><td>0.2</td><td>6.92215E-2 (7.06664E-2)</td><td>6.61159E-3 (9.72231E-3)</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">Test Problem 9b</td></tr><tr><td>0</td><td>0 (0)</td><td>1 (1)</td><td>1.47242E-1<br>(1.00045E-1)</td><td>8.03070E-1<br>(8.04568E-1)</td><td>0.2</td><td>3.18310E-1 (3.18310E-1)</td><td>7.36290E-2 (7.38886E-2)</td></tr><tr><td>1.05</td><td>0 (0)</td><td>3.78125E-1<br>(4.52455E-1)</td><td>6.43614E-2<br>(4.52170E-2)</td><td>2.65492E-1<br>(2.55317E-1)</td><td>0.2</td><td>3.18310E-1 (3.18310E-1)</td><td>2.11284E-2 (1.32785E-2)</td></tr><tr><td>2.1</td><td>0 (0)</td><td>1.56023E-1<br>(2.36472E-1)</td><td>2.74784E-2<br>(2.41999E-2)</td><td>1.10101E-1<br>(1.31017E-1)</td><td>0.2</td><td>6.00555E-2 (5.91905E-2)</td><td>3.17544E-2 (3.00285E-2)</td></tr><tr><td>6</td><td>0 (0)</td><td>7.75973E-3<br>(2.76496E-2)</td><td>1.87563E-3<br>(4.16071E-3)</td><td>4.81768E-3<br>(1.36242E-2)</td><td>0.2</td><td>1.92929E-1 (1.96606E-1)</td><td>9.53919E-3 (7.05806E-3)</td></tr><tr><td>21</td><td>0 (0)</td><td>2.33796E-6<br>(7.41883E-5)</td><td>-1.66315E-22 (0)</td><td>4.67591E-7<br>(1.22028E-5)</td><td>0.2</td><td>1.08508E-1 (1.15481E-1)</td><td>4.23567E-3 (4.06938E-3)</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">Test Problem 9c</td></tr><tr><td>0</td><td>1.57080<br>(1.57080)</td><td>6.09217<br>(6.09217)</td><td>4.77677<br>(4.68410)</td><td>3.47067<br>(3.49559)</td><td>0.2<sup>1</sup></td><td>1.93920 (1.93920)</td><td>1.62423 (1.61851)</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">Test Problem 9d</td></tr><tr><td>1.05</td><td>1.92354E-1<br>(1.92354E-1)</td><td>4.84435<br>(4.97288)</td><td>4.25390<br>(4.24382)</td><td>8.81645E-1<br>(8.81133E-1)</td><td>0.2<sup>1</sup></td><td>1.93920 (1.93920)</td><td>1.44392 (1.43873)</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">Test Problem 9e</td></tr><tr><td>2.1</td><td>2.35550E-2<br>(2.35550E-2)</td><td>4.33478<br>(4.46617)</td><td>4.14656<br>(4.16942)</td><td>3.04441E-1<br>(3.50114E-1)</td><td>0.2<sup>1</sup></td><td>1.93920 (1.93920)</td><td>1.58486 (1.57890)</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">Test Problem 9f</td></tr><tr><td>6</td><td>9.65131E-6<br>(9.65131E-6)</td><td>4.21639<br>(4.22731)</td><td>4.28367<br>(4.30667)</td><td>5.57724E-3<br>(1.93491E-2)</td><td>0.2<sup>1</sup></td><td>1.93920 (1.93920)</td><td>1.57163 (1.56554)</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">Test Problem 9g</td></tr><tr><td>21</td><td>9.03133E-19<br>(9.03133E-19)</td><td>4.78641<br>(4.73767)</td><td>5.13962<br>(5.11524)</td><td>7.81157E-2<br>(7.16719E-2)</td><td>0.2<sup>1</sup></td><td>1.93920 (1.93920)</td><td>1.44392 (1.43873)</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">Test Problem 9h</td></tr><tr><td>0</td><td>1.57080<br>(1.57080)</td><td>6.09217<br>(6.09217)</td><td>4.77677<br>(4.68410)</td><td>3.47067<br>(3.49559)</td><td>0.2<sup>1</sup></td><td>1.93920 (1.93920)</td><td>1.62423 (1.61851)</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">Test Problem 9i</td></tr><tr><td>1.05</td><td>1.92354E-1<br>(1.92354E-1)</td><td>4.84435<br>(4.97288)</td><td>4.25390<br>(4.24382)</td><td>8.81645E-1<br>(8.81133E-1)</td><td>0.2<sup>1</sup></td><td>1.93920 (1.93920)</td><td>1.44392 (1.43873)</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">Test Problem 9j</td></tr><tr><td>2.1</td><td>2.35550E-2<br>(2.35550E-2)</td><td>4.33478<br>(4.46617)</td><td>4.14656<br>(4.16942)</td><td>3.04441E-1<br>(3.50114E-1)</td><td>0.2<sup>1</sup></td><td>1.93920 (1.93920)</td><td>1.58486 (1.57890)</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">Test Problem 9k</td></tr><tr><td>6</td><td>9.65131E-6<br>(9.65131E-6)</td><td>4.21639<br>(4.22731)</td><td>4.28367<br>(4.30667)</td><td>5.57724E-3<br>(1.93491E-2)</td><td>0.2<sup>1</sup></td><td>1.93920 (1.93920)</td><td>1.57163 (1.56554)</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">Test Problem 9l</td></tr><tr><td>21</td><td>9.03133E-19<br>(9.03133E-19)</td><td>4.78641<br>(4.73767)</td><td>5.13962<br>(5.11524)</td><td>7.81157E-2<br>(7.16719E-2)</td><td>0.2<sup>1</sup></td><td>1.93920 (1.93920)</td><td>1.44392 (1.43873)</td></tr></tbody></table>

<caption>**Table 17**</caption>

Computation time comparisons.

<table><thead><tr><th rowspan="2">Scenario</th><th rowspan="2">Number of monochromatic RT calculations</th><th colspan="3">Mean computation time (s)</th><th colspan="2">2S-ESS speed-up compared with</th></tr><tr><th>2S-ESS</th><th>LIDORT</th><th>LIDORT (2S mode)</th><th>LIDORT</th><th>LIDORT (2S mode)</th></tr></thead><tbody><tr><td>Solar</td><td>280,000</td><td>43.510</td><td>115,074.808</td><td>336.173</td><td>2645</td><td>7.7</td></tr><tr><td>Thermal</td><td>200,000</td><td>8.635</td><td>24,069.026</td><td>170.156</td><td>2787</td><td>19.7</td></tr></tbody></table>

<footnote n="1">¹ $\varphi = 0^\circ$</footnote>

<footnote n="2">² $\varphi = 90^\circ$</footnote>

<footnote n="3">³ $\varphi = 180^\circ$</footnote>
particular integral solutions either by substitution with an exponential expression having dependence on optical depth mirroring that of the solar source term, or by substitution of a linear-in-optical-depth expression akin to the thermal emission linear parameterization. In the former case, the use of the "average secant" formulation for solar attenuation follows the practice adopted in the LIDORT and VLIDORT models with the pseudo-spherical approximation. To the best of our knowledge, the use of a Green's function solution in the 2S model is new; however, the solution presented here is based on the LIDORT Green's function treatment (first introduced in Ref. [17], and more recently developed for the VLIDORT model).

As far as the ESS model is concerned, the main new approach here is with the use of a quadrature scheme to solve the SS RTE in a fully spherical geometrical scenario; this scheme was first introduced into the LIDORT and VLIDORT models some years ago [43], and a similar but more general ESS model (including limb viewing) has recently been validated explicitly against a spherical Monte Carlo code [44].

2S calculates the MS radiation field in the two-stream discrete-ordinate approximation. First, 2S determines the complete (SS and MS) discrete-ordinate field (in the stream directions +/-) obtained from truncation of the phase function implicit in the solution method. However, for the post-processed radiation field in other polar directions (obtained by source-function integration), only the MS part of this field is evaluated. Indeed, the post-processed SS part of the 2S field is ignored since it is based on severe phase function truncation, and is replaced by a fully accurate independent result from the ESS model, which is based on the full phase function, and is not subject to the strictures of plane-parallel scattering. This procedure is described in detail in Ref. [12].

## 9. Conclusions

In this paper, we describe the capabilities and limitations of the generic, flexible and computationally efficient 2S-ESS RT model for vertically inhomogeneous media that accounts for spherical geometry effects. The theoretical bases of the 2S and ESS approximations are covered in rigorous detail. We perform detailed comparisons of the accuracy of the upwelling and downwelling radiances from the model for both slab and multi-layer scenarios. The results show that, for practical applications, the errors due to the 2S-ESS approximation are on the order of 0.1% for solar beam and 0.01% for thermal emission scenarios. Furthermore, the 2S-ESS model significantly outperforms the Full-2S model that truncates the forward peak of strongly anisotropic scattering, with negligible loss of computational efficiency. We also perform comparisons with the test cases included in the DISORT code package.

We remark here that the 2S-ESS RT code is fully linearized: in addition to generating the radiation field, the model can also produce analytically-derived weighting functions with respect to any atmospheric property (e.g., trace gas absorber number density profiles, or bulk properties such as total aerosol optical depth) or any surface quantity (e.g., albedo, or wind speed in sunglitt BRDF scenarios). This facility is useful for remote sensing retrieval applications. The Jacobian output can be compared directly against related output from the LIDORT model; the latter also has a full Jacobian facility. A second paper on the performance and accuracy of linearized output is currently in preparation.

The 2S-ESS code is publicly available on GitHub: https://github.com/vnatraj1/2S-ESS. All the data relevant to the tables and figures in the paper can be downloaded from zenodo: https://doi.org/10.5281/zenodo.6780812. We anticipate that the code will be of considerable interest for Earth, planetary and exoplanetary science applications.

## Declaration of Competing Interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

Vijay Natraj reports financial support and article publishing charges were provided by Jet Propulsion Laboratory.

## CRediT authorship contribution statement

**V. Natraj:** Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Writing – original draft, Supervision.

**R. Spurr:** Software, Writing – review & editing.

**A. Gao:** Validation.

**T. Le:** Formal analysis, Visualization.

**Z.C. Zeng:** Writing – review & editing.

**S. Fan:** Data curation.

**Y.L. Yung:** Conceptualization, Writing – review & editing.

**Data availability**

Code, documentation and sample driver programs are on GitHub. Codes to generate all the Tables and Figures in the manuscript have been uploaded to Zenodo. The links are provided in the manuscript.

## Acknowledgments

A portion of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). VN acknowledges support from the Tropospheric Ozone and its Precursors from Earth System Sounding (TROPESS) project. The authors gratefully acknowledge the insightful and constructive comments from the two anonymous reviewers, which improved the clarity and quality of the manuscript.

## References

[1] Stamnes K, Thomas GE, Stamnes JJ. Radiative transfer in the atmosphere and ocean. 2nd edition. Cambridge University Press; 2017. doi:10.1017/9781316148549.

[2] Shettle EP, Weinman JA. The transfer of solar irradiance through inhomogeneous turbid atmospheres evaluated by Eddington's approximation. J Atmos Sci 1970;27:1048–55.

[3] Liou KN. Analytic two-stream and four-stream solutions for radiative transfer. J Atmos Sci 1974;31:1473–5.

[4] Joseph JH, Wiscombe WJ, Weinman JA. The delta-Eddington approximation for radiative flux transfer. J Atmos Sci 1976;33:2452–9.

[5] Meador WE, Weaver WR. Two-stream approximations to radiative transfer in planetary atmospheres: a unified description of existing methods and a new improvement. J Atmos Sci 1980;37:630–43.

[6] King MD, Harshvardhan. Comparative accuracy of selected multiple scattering approximations. J Atmos Sci 1986;43(8):784–801.

[7] Toon O, McKay CP, Ackerman TP, Santhanam K. Rapid calculation of radiative heating rates and photodissociation rates in inhomogeneous multiple scattering atmospheres. J Geophys Res 1989;94:16287–301.

[8] Kylling A, Stamnes K, Tsay SC. A reliable and efficient two-stream algorithm for spherical radiative transfer: documentation of accuracy in realistic layered media. J Atmos Chem 1995;21:115–50.

[9] Gabriel P, Stephens GL, Wittmeyer IL. Adjoint perturbation and selection rule methods for solar broadband two-stream fluxes in multi-layer media. J Quant Spectrosc Radiat Transf 2000;65:693–728.

[10] Chen-Chen H, Perez-Hoyos S, Sanchez-Lavega A. Assessing multi-stream radiative transfer schemes for the calculation of aerosol radiative forcing in the Martian atmosphere. J Geophys Res 2021;126(10):e2021JE006889. doi:10.1029/2021JE006889.

[11] Min QL, Harrison LC. Retrieval of atmospheric optical depth profiles from downward-looking high-resolution O₂ A-band measurements: optically thin conditions. J Atmos Sci 2004;61:2469–77.

[12] Spurr R, Natraj V. A linearized two-stream radiative transfer code for fast approximation of multiple-scatter fields. J Quant Spectrosc Radiat Transf 2011;112:2630–7.

[13] Zhang Q, Natraj V, Li KF, Shia RL, Fu D, Pongetti TJ, et al. Accounting for aerosol scattering in the CLARS retrieval of column averaged CO₂ mixing ratios. J Geophys Res 2015;120:7205–18. doi:10.1002/2015JD023499.
[14] Zeng ZC, Zhang Q, Natraj V, Margolis JS, Shia RL, Newman S, et al. Aerosol scattering effects on water vapor retrievals over the Los Angeles Basin. Atmos Chem Phys 2017;17:2495–508. doi:10.5194/acp-17-2495-2017.

[15] Zeng ZC, Natraj V, Xu F, Pongetti TJ, Shia RL, Kort EA, et al. Constraining aerosol vertical profile in the boundary layer using hyperspectral measurements of oxygen absorption. Geophys Res Lett 2018;45:10772–80. doi:10.1029/2018GL079286.

[16] Zeng ZC, Xu F, Natraj V, Pongetti TJ, Shia RL, Zhang Q, et al. Remote sensing of angular-dependent scattering of aerosols in a North American megacity. Remote Sens Environ 2020;242:111760. doi:10.1016/j.rse.2020.111760.

[17] Spurr R. Simultaneous radiative transfer derivation of intensities and weighting functions in a general pseudo-spherical treatment. J Quant Spectrosc Radiat Transf 2002;75(4):129–75. doi:10.1016/S0022-4073(01)00245-X.

[18] Stamnes K, Tsay SC, Wiscombe W, Jayaweera K. Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting layered media. Appl Opt 1988;27:2502–9.

[19] Stamnes K, Conklin P. A new multi-layer discrete ordinate approach to radiative transfer in vertically inhomogeneous atmospheres. J Quant Spectrosc Radiat Transf 1984;31:273.

[20] Van Oss RF, Spurr RJD. Fast and accurate 4 and 6 stream linearized discrete ordinate radiative transfer models for ozone profile retrieval. J Quant Spectrosc Radiat Transf 2002;75:177–220.

[21] Siewert CE. A concise and accurate solution to Chandrasekhar's basic problem in radiative transfer. J Quant Spectrosc Radiat Transf 2000;64:109–30.

[22] Spurr R, Christi M, Kokhanovsky A. The LIDORT and VLIDORT linearized scalar and vector discrete ordinate radiative transfer models: updates in the last 10 years. Light scattering reviews, 12. Springer; 2019.

[23] Laszlo I, Stamnes K, Wiscombe WJ, Tsay SC, Kokhanovsky A. The discrete ordinate algorithm, DISORT for radiative transfer. Light scattering reviews, 11, Berlin, Heidelberg: Springer; 2016. Springer Praxis Books. doi:10.1007/978-3-662-49538-4_1.

[24] Natraj V, Jiang X, Shia RL, Huang X, Margolis JS, Yung YL. Application of principal component analysis to high spectral resolution radiative transfer: a case study of the O₂ A band. J Quant Spectrosc Radiat Transf 2005;95(4):539–56. doi:10.1016/j.jqsrt.2004.12.024.

[25] Natraj V, Shia RL, Yung YL. On the use of principal component analysis to speed up radiative transfer calculations. J Quant Spectrosc Radiat Transf 2010;111:810–16. doi:10.1016/j.jqsrt.2009.11.004.

[26] Spurr R, Natraj V, Lerot C, Van Roozendael M, Loyola D. Linearization of the principal component analysis method for radiative transfer acceleration: application to retrieval algorithms and sensitivity studies. J Quant Spectrosc Radiat Transf 2013;125:1–17. doi:10.1016/j.jqsrt.2013.04.002.

[27] Efremenko D, Doicu A, Loyola D, Trautmann T. Optical property dimensionality reduction techniques for accelerated radiative transfer performance: application to remote sensing total ozone retrievals. J Quant Spectrosc Radiat Transf 2014;133:128–35. doi:10.1016/j.jqsrt.2013.07.023.

[28] Efremenko DS, Loyola DG, Spurr R, Doicu A. Acceleration of radiative transfer model calculations for the retrieval of trace gases under cloudy conditions. J Quant Spectrosc Radiat Transf 2014;135:58–65. doi:10.1016/j.jqsrt.2013.11.014.

[29] Kopparla P, Natraj V, Spurr RJD, Shia RL, Yung YL, Crisp D. A fast and accurate PCA based radiative transfer model: extension to the broadband shortwave region. J Quant Spectrosc Radiat Transf 2016;173:65–71. doi:10.1016/j.jqsrt.2016.01.014.

[30] Kopparla P, Natraj V, Limpasuvan D, Spurr R, Crisp D, Shia RL, et al. PCA-based radiative transfer: improvements to aerosol scheme, vertical layering and spectral binning. J Quant Spectrosc Radiat Transf 2017;198:104–11. doi:10.1016/j.jqsrt.2017.05.005.

[31] Somkuti P, Bösch H, Natraj V, Kopparla P. Application of a PCA-based fast radiative transfer model to XCO₂ retrievals in the shortwave infrared. J Geophys Res 2017;122(19):10268–87. doi:10.1002/2017JD027013.

[32] Zeng ZC, Natraj V, Xu F, Chen S, Gong FY, Pongetti TJ, et al. GFIT3: a full physics retrieval algorithm for remote sensing of greenhouse gases in the presence of aerosols. Atmos Meas Technol 2021;14:6483–507. doi:10.5194/amt-14-6483-2021.

[33] Kokhanovsky AA, Kokhanovsky AA. Asymptotic radiative transfer. Light scattering reviews, Berlin, Heidelberg: Springer; 2006. Springer Praxis Books. doi:10.1007/3-540-37672-0_6.

[34] Grell GA, Peckham SE, Schmitz R, McKeen SA, Frost G, Skamarock WC, et al. Fully coupled "online" chemistry within the WRF model. Atmos Environ 2005;39(37):6957–75. doi:10.1016/j.atmosenv.2005.04.027.

[35] Fast JD, Gustafson Jr WI, Easter RC, Zaveri RA, Barnard JC, Chapman EG, et al. Evolution of ozone, particulates, and aerosol direct radiative forcing in the vicinity of Houston using a fully coupled meteorology-chemistry-aerosol model. J Geophys Res 2006;111(D21).

[36] Peckham S., Fast J., Schmitz R., Grell G., Gustafson W., McKeen S., et al. WRF/Chem version 3.3 users guide, 94 p. 2011.

[37] Van de Hulst HC. Multiple light scattering, tables, formulas and applications, 1. New York: Academic Press; 1980.

[38] Van de Hulst HC. Multiple light scattering, tables, formulas and applications, 2. New York: Academic Press; 1980.

[39] Garcia RDM, Siewert CE. Benchmark results in radiative transfer. Transport Theo Stat Phys 1985;14:437–83.

[40] Ozisik M, Shouman S. Source function expansion method for radiative transfer in a two-layer slab. J Quant Spectrosc Radiat Transf 1980;24:441–9.

[41] Devaux C, Grandjean P, Ishiguro Y, Siewert CE. On multi-region problems in radiative transfer. Astrophys Space Sci 1979;62:225–33.

[42] Stamnes K., Tsay S.C., Wiscombe W., Laszlo I. DISORT, a general-purpose Fortran program for discrete-ordinate-method radiative transfer in scattering and emitting layered media: documentation of methodology, DISORT report v1.1, 2000.

[43] Spurr RL, Kokhanovsky AA. VLIDORT: linearized pseudo-spherical scalar and vector discrete ordinate radiative transfer models for use in remote sensing retrieval problems. Light scattering reviews, 3, Berlin, Heidelberg: Springer; 2008. Springer Praxis Books. doi:10.1007/978-3-540-48456-9_7.

[44] Spurr R, Natraj V, Colosimo SF, Stutz J, Christi M, Korkin S. VLIDORT-QS: a quasi-spherical vector radiative transfer model. J Quant Spectrosc Radiat Transf 2022;291:108341. doi:10.1016/j.jqsrt.2022.108341.