Scalar thermal radiation using the adding-doubling method

Bingqiang Sun; Bingqiang Sun; Chenxu Gao; Chenxu Gao; Robert Spurr

doi:10.1364/OE.462580

1. Introduction

Atmospheric radiation is crucial for remote sensing, weather, and climate studies [1–3]. The theoretical foundation for the calculation of atmospheric radiation is the radiative transfer equation (RTE) [4]. The RTE can be obtained under the incoherence and far-field approximations from Maxwell’s equations [5,6]. In general, the RTE is a three-dimensional (3D) vector integro-differential equation with no analytic solutions [3,4]. The 3D RTE has been extensively discussed in the literature; see for example [7]. However, the calculation of 3D RTE is time-consuming and not considered in this study. The plane-parallel assumption is employed in this research after ignoring the impact of terrestrial curvature. For real calculations of the RTE, the atmosphere is divided into a stratification of layers, for which the single scattering albedos and scattering phase matrices are constant within each layer. In each layer, physical processes include molecular absorption by gas species, Rayleigh scattering by air molecules, and scattering and absorption of particulate matters (aerosols and clouds) [1–3].

The molecular absorption coefficient can be derived either with the line-by-line integration [8], the correlated K-distribution method [9–11], or some other band model methods [1], with spectroscopic data from the HITRAN database [12]. Rayleigh scattering coefficient is proportional to the reciprocal of the fourth power of wavelength and can be straightforwardly obtained using the number density of air molecules, the index of refraction and the depolarization ratio of air [2]. Monochromatic specification or elastic scattering is usually required for the scattering and absorption of particulate matters [13]. Single scattering by particulate matters is actually a bulk property determined by averaging single particle scattering in random orientations over particle size distribution [14]. Single particle scattering properties can be found by solving Maxwell’s equations [13]. Scattering by spherical particles can be obtained analytically using the Lorenz-Mie theory [13,15,16]. Scattering by non-spherical particles has to be determined using suitable numerically accurate approaches, such as the finite-difference time-domain method [17–19], the discrete-dipole approximation method [20,21], the T-matrix method [22–26], and the geometric optics methods [27–35]. In narrowband models, monochromatic scattering of particulate matters and Rayleigh scattering have to be spectrally integrated for the coupling with molecular absorption using the correlated K-distribution method [9,14]. Accordingly, the total single scattering albedo, scattering phase matrix, and optical depth in spectral average can be provided in each layer. Furthermore, the radiation effect of air molecules and particulate matters is handled as multiple scattering process using the RTE.

Multiple scattering in the RTE can be solved numerically by the invariant-imbedding method [4], the successive order of scattering method [36–38], the discrete ordinate method [39,40], the matrix operator method [41,42], the Monte Carlo method [43,44], adding-doubling methods [45–47], and so on [48]. The adding-doubling (A-D) algorithm is easy to implement and can be cast in general form for scalar and vector scenarios [49–53]. With the involvement of viewing zenith angles, the repeated reflections between two layers in the adding or doubling process are dealt with the summation of multiple reflections using the Taylor expansion in the classical A-D algorithm literature [45]. In thermal radiation, the Planck function is usually assumed to vary within a layer in terms of variations in the vertical temperature profile.

In this work, the repeated reflections are directly computed in terms of matrix inversions with the involvement of arbitrary viewing zenith angles so the radiance output in any viewing zenith angle can be quickly obtained. Moreover, the adding-doubling technique is derived in detail for thermal radiation source, and the doubling algorithm is specifically treated in the linear and exponential parameterizations. A state-of-the-art MATLAB program implements the algorithm numerically. The program is fast and a general scenario can be completed in a matter of seconds after optimization. Here we have implemented just the scalar scenario in the MATLAB program since polarization is weak in the thermal IR regime and the flux is invariant either for scalar or vector cases. The algorithm can be generalized in a straightforward manner to the vector RTE.

This study is organized as follows: the adding-doubling algorithm and the doubling techniques for linear and exponential parameterizations are shown in detail in section 2; Verification of the model and discussions on the results are presented in section 3; Concluding remarks are given in section 4.

2. Method

The plane-parallel assumption is employed in this study and a schematic figure for thermal radiation is shown in Fig. 1. The atmosphere is divided into n layers or (n+1) levels, where the scattering phase function ${p_i}$ and single scattering albedo ${\omega _i}$ are constant within a layer i ($i = 1,\ldots ,n$). The total optical depth is $\tau _t^i$ for layer i ($i = 1,\ldots ,n$), and the temperature profile ${T_i}$ is given at the level i ($i = 1,\ldots ,n + 1$). For the boundary conditions, the top of atmosphere (TOA) has a 2.725K (T_CMB) cosmic microwave background and the ground surface is characterized by a surface albedo ${\alpha _S}$ either with a Lambertian case or with a Bidirectional Reflectance Distribution Function (BRDF) R_s, and in addition there is surface thermal emission for surface temperature ${T_S}$. The scalar radiative transfer equation in the thermal radiation regime can then be written as:

(1)$$u\frac{{\partial I(\tau ,u)}}{{\partial \tau }} = I(\tau ,u) - \frac{\omega }{2}\int_{ - 1}^1 {p(u,u\prime )I(\tau ,u\prime )du\prime } - (1 - \omega )B(\tau ),$$

and the boundary conditions are

(2)$$I(0, - \mu ) = B({T_{\textrm{CMB}}}),$$

(3)$$I({\tau _S},\mu ) = {\alpha _S}\int_0^1 {{R_S}(\mu , - \mu \prime )I({\tau _S}, - \mu \prime )d\mu \prime } ,$$

where the layer or level number is suppressed for clarity. Also, ${\tau _s} = \sum\limits_{i = 1}^n {\tau _t^i} $ and $\int_0^1 {{R_S}(\mu , - \mu \prime )d\mu \prime } = 1$. Zenith angle $\theta $ is defined relative to the z-axis as shown in Fig. 1 and the convention $u = \cos \theta \textrm{ and }\mu = |u|$ are used throughout this paper. The phase function $p(u,u\prime )$ is azimuth-independent because the Planck function is isotropic; the wavenumber (or wavelength) related quantities $I,\;\tau ,\;p,\;\omega ,\;B,{\alpha _S},\textrm{ and }{R_S}$ are averaged over a narrow wavenumber (or wavelength) band. All quantities are represented in wavenumber space in the study.

Fig. 1. Schematic figure for scalar thermal radiation in the plane-parallel assumption.

Download Full Size | PDF

The phase function can be expanded in terms of Legendre polynomials as:

(4)$$P(u,u\prime ,\varphi - \varphi \prime ) = P(x) = \sum\limits_{n = 0}^N {(2n + 1){\alpha _n}{P_n}(x)} ,(x = \cos \Theta )$$

(5)$$\cos \Theta = \mu \mu \prime + \sqrt {1 - {\mu ^2}} \sqrt {1 - \mu {\prime ^2}} \cos (\varphi - \varphi \prime ),$$

(6)$${\alpha _n} = \frac{1}{2}\int_{ - 1}^1 {P(x){P_n}(x)dx} ,$$

where ${P_n}(u)$ denotes the Legendre polynomial and the azimuthal-independent phase function is related to the total phase function through

(7)$$p(u,u\prime ) = \frac{1}{{2\pi }}\int_0^{2\pi } {P(u,u\prime ,\varphi - \varphi \prime )d({\varphi - \varphi \prime } )} .$$

(8)$$p(u,u\prime ) = \sum\limits_{n = 0}^N {(2n + 1){\alpha _n}{P_n}(u){P_n}(u\prime )} .$$

The Planck function per wavenumber is given by

(9)$${B_\nu }(T) = \frac{{2h{\nu ^3}{c^2}}}{{{e^{h\nu c/{k_B}T}} - 1}},$$

where the constants (in SI units) are the Planck constant $h = 6.62607015\ast {10^{ - 34}}J \cdot s$, the light speed $c = 2.99792458\ast {10^8}m \cdot {s^{ - 1}}$, and the Boltzmann constant ${k_B} = 1.380648\ast {10^{ - 23}}J \cdot {K^{ - 1}}$. The Planck function can be integrated over any wavenumber band as follows

(10)$$\begin{aligned} {B_{\Delta \nu }}(T) = \int_{{\nu _1}}^{{\nu _2}} {{B_\nu }(T)d\nu } &= {c_1}\left[ {\frac{{6({\textrm{L}{\textrm{i}_4}({e^{ - {c_2}{\nu_1}}}) - \textrm{L}{\textrm{i}_4}({e^{ - {c_2}{\nu_2}}})} )}}{{c_2^4}} + \frac{{\textrm{3}({\nu_1^2\textrm{L}{\textrm{i}_2}({e^{ - {c_2}{\nu_1}}}) - \nu_2^2\textrm{L}{\textrm{i}_2}({e^{ - {c_2}{\nu_2}}})} )}}{{c_2^2}}} \right] + \\& {c_1}\left[ {\frac{{6({{\nu_1}\textrm{L}{\textrm{i}_3}({e^{ - {c_2}{\nu_1}}}) - {\nu_2}\textrm{L}{\textrm{i}_3}({e^{ - {c_2}{\nu_2}}})} )}}{{c_2^3}} + \frac{{\nu_1^3\textrm{L}{\textrm{i}_1}({e^{ - {c_2}{\nu_1}}}) - \nu_2^3\textrm{L}{\textrm{i}_1}({e^{ - {c_2}{\nu_2}}})}}{{{c_2}}}} \right], \end{aligned}$$

where the units for wavenumber and temperature in Eq. (10) are cm^-1 and K, and ${c_1} = 2h{c^2} = 1.19104297\ast {10^{ - 8}}W \cdot {m^{ - 2}} \cdot {(c{m^{ - 1}})^{ - 1}} \cdot {(c{m^{ - 3}})^{ - 1}}$ and ${c_2} = hc/({k_B}T) = 1.4387768775/T\;cm$ are two constants. The poly-logarithm function Li is defined as

(11)$$\textrm{L}{\textrm{i}_n}(x) = \sum\limits_{s = 1}^\infty {\frac{{{x^s}}}{{{s^n}}}} .$$

The explicit expression for n = 1 is $\textrm{L}{\textrm{i}_1}(x) ={-} \ln (1 - x)$. The source function $B(\tau )$ in Eq. (1) is given in terms of ${B_{\Delta \nu }}(T)$ in Eq. (10). The linear and exponential parameterizations for blackbody radiation can be respectively given as

(12)$$B(\tau ) = a + b\tau ,\;a = {B_0},\;b = \frac{{{B_1} - {B_0}}}{{{\tau _\textrm{t}}}},$$

(13)$$B(\tau ) = a{e^{b\tau }},\;a = {B_0},\;b = \ln ({B_1}/{B_0})/{\tau _\textrm{t}},$$

where ${\tau _\textrm{t}}$ is the optical thickness of the layer; B₀ and B₁ are the top and bottom blackbody radiation of a layer, respectively. The constant parameterization can be straightforwardly reached by setting b = 0 and a = constant in Eq. (12) or (13).

The flux F and the actinic flux F_a at any optical depth are defined by

(14)$${F^ \pm }(\tau ) = 2\pi \int_0^1 {{I^ \pm }(\tau , \pm \mu )\mu d\mu } ,$$

(15)$$F_a^ \pm (\tau ) = \frac{1}{2}\int_0^1 {{I^ \pm }(\tau , \pm \mu )d\mu } .$$

The total fluxes can be defined as follows:

(16)$$F = 2\pi \int_{ - 1}^1 {I(\tau ,u)udu} = {F^ + } - {F^ - },$$

(17)$${F_a} = \frac{1}{2}\int_{ - 1}^1 {I(\tau ,u)du} = F_a^ +{+} F_a^ - .$$

Integrating both side of Eq. (1) as $\frac{1}{2}\int_{ - 1}^1 {du}$ and using the orthogonality of the Legendre polynomial, one can get the Gershun’s law in thermal radiation as

(18)$$\frac{{\partial F}}{{\partial \tau }} ={-} (1 - \omega )[{4\pi (B - {F_a})} ],$$

which states not only the physical meaning of the actinic flux, but also the energy conservation in thermal radiation.

Solutions for the upwelling and downwelling radiances associated with internal blackbody radiation sources can be formally expressed as:

(19)$$I = {\Delta _{ds}} + {\Delta _{ss}} + {\Delta _{ms}},$$

where “$\Delta $” represents the radiation from the blackbody radiation. The first term represents the direct thermal radiation, and can be generated explicitly. Upwelling radiances (denoted by superscript “+”) for the direct radiation at any optical depth and at the top of the layer are

(20)$$\begin{array}{l} \Delta _{ds}^ + (\tau ,\mu ) = (1 - \omega )d_0^ + (\tau ,\mu ),\\ \Delta _{ds}^ + (0,\mu ) = (1 - \omega ){r_0}(\mu ), \end{array}$$

where for linear parameterization,

(21)$$\begin{array}{l} {d^ + }(\tau ,\mu ) = \left\{ {a\left( {1 - {e^{ - \frac{{{\tau_\textrm{t}} - \tau }}{\mu }}}} \right) + b\left[ {({\tau + \mu } )- ({{\tau_\textrm{t}} + \mu } ){e^{ - \frac{{{\tau_\textrm{t}} - \tau }}{\mu }}}} \right]} \right\},\\ {r_0}(\mu ) = \left\{ {a\left( {1 - {e^{ - \frac{{{\tau_\textrm{t}}}}{\mu }}}} \right) + b\left[ {\mu - ({{\tau_\textrm{t}} + \mu } ){e^{ - \frac{{{\tau_\textrm{t}}}}{\mu }}}} \right]} \right\}, \end{array}$$

and for exponential parameterization,

(22)$$\begin{aligned} d_0^ + (\tau ,\mu ) &= \frac{a}{{1 - b\mu }}\left( {{e^{b\tau }} - {e^{b{\tau_\textrm{t}}}}{e^{ - \frac{{{\tau_\textrm{t}} - \tau }}{\mu }}}} \right),\\ {r_0}(\mu ) &= \left\{ \begin{array}{l} \frac{a}{{1 - b\mu }}[{1 - {e^{ - (1 - b\mu ){\tau_\textrm{t}}/\mu }}} ],\;\;1 - b\mu \ne 0\\ a\frac{{{\tau_\textrm{t}}}}{\mu },\qquad\qquad\qquad\qquad1 - b\mu = 0 \end{array} \right.. \end{aligned}$$

Similarly, downwelling radiances (denoted by superscript “-” for the direct radiation at any optical depth and at the bottom of the layer are

(23)$$\begin{array}{l} \Delta _{ds}^ - (\tau ,\mu ) = (1 - \omega )d_0^ - (\tau ,\mu ),\\ \Delta _{ds}^ - ({\tau _\textrm{t}}, - \mu ) = (1 - \omega ){t_0}(\mu ), \end{array}$$

where now the coefficients for linear parameterization are given by

(24)$$\begin{array}{l} d_0^ - (\tau ,\mu ) = \left\{ {a\left( {1 - {e^{ - \frac{\tau }{\mu }}}} \right) + b\left[ {({\tau - \mu } )+ \mu {e^{ - \frac{\tau }{\mu }}}} \right]} \right\},\\ {t_0}(\mu ) = \left\{ {a\left( {1 - {e^{ - \frac{{{\tau_\textrm{t}}}}{\mu }}}} \right) + b\left[ {({{\tau_\textrm{t}} - \mu } )+ \mu {e^{ - \frac{{{\tau_\textrm{t}}}}{\mu }}}} \right]} \right\}, \end{array}$$

and the coefficients for exponential parameterization are

(25)$$\begin{array}{l} d_0^ - (\tau ,\mu ) = \frac{a}{{1 + b\mu }}\left( {{e^{b\tau }} - {e^{ - \frac{\tau }{\mu }}}} \right),\\ {t_0}(\mu ) = \frac{a}{{1 + b\mu }}\left( {{e^{b{\tau_\textrm{t}}}} - {e^{ - \frac{{{\tau_\textrm{t}}}}{\mu }}}} \right). \end{array}$$

Single scattering for blackbody radiation can be expressed formally or calculated numerically. Upwelling and downwelling radiances for single scattering at layer top and bottom are then:

(26)$$\begin{array}{l} \Delta _{ss}^ + (0,\mu ) = \frac{{\omega (1 - \omega )}}{2}\int_0^1 {[{p(\mu ,\mu \prime ){r_1}(\mu ,\mu \prime ) + p(\mu , - \mu \prime ){r_2}(\mu ,\mu \prime )} ]d\mu \prime } ,\\ \Delta _{ss}^ - ({\tau _\textrm{t}}, - \mu ) = \frac{{\omega (1 - \omega )}}{2}\int_0^1 {[{p( - \mu ,\mu \prime ){t_1}(\mu ,\mu \prime ) + p( - \mu , - \mu \prime ){t_2}(\mu ,\mu \prime )} ]d\mu \prime } , \end{array}$$

where for the linear parameterization, the coefficients are

(27)$$\begin{aligned} \mu &\ne \mu \prime ,\qquad{r_1}(\mu ,\mu \prime ) = \;\;a\left[ {\left( {1 - {e^{ - \frac{{{\tau_\textrm{t}}}}{\mu }}}} \right) - \frac{{\mu \prime }}{{\mu - \mu \prime }}\left( {{e^{ - \frac{{{\tau_\textrm{t}}}}{\mu }}} - {e^{ - \frac{{{\tau_\textrm{t}}}}{{\mu \prime }}}}} \right)} \right]\\ \ &+ b\left[ {\mu + \mu \prime - ({{\tau_\textrm{t}} + \mu + \mu \prime } ){e^{ - \frac{{{\tau_\textrm{t}}}}{\mu }}} - ({{\tau_\textrm{t}} + \mu \prime } )\frac{{\mu \prime }}{{\mu - \mu \prime }}\left( {{e^{ - \frac{{{\tau_\textrm{t}}}}{\mu }}} - {e^{ - \frac{{{\tau_\textrm{t}}}}{{\mu \prime }}}}} \right)} \right], \end{aligned}$$

(28)$$\begin{aligned} \mu = \mu \prime ,\quad{r_1}&(\mu ,\mu \prime ) = \;\;a\left[ {\left( {1 - {e^{ - \frac{{{\tau_\textrm{t}}}}{\mu }}}} \right) - {e^{ - \frac{{{\tau_\textrm{t}}}}{\mu }}}\frac{{{\tau_\textrm{t}}}}{\mu }} \right]\\ &+ b\left[ {2\mu - \left( {\frac{{\tau_\textrm{t}^2}}{\mu } + 2{\tau_\textrm{t}} + 2\mu } \right){e^{ - \frac{{{\tau_\textrm{t}}}}{\mu }}}} \right], \end{aligned}$$

(29)$$\begin{array}{l} {r_2}(\mu ,\mu \prime ) = \;\;a\left[ {\left( {1 - {e^{ - \frac{{{\tau_\textrm{t}}}}{\mu }}}} \right) - \frac{{\mu \prime }}{{\mu + \mu \prime }}\left( {1 - {e^{ - \frac{{{\tau_\textrm{t}}}}{\mu } - \frac{{{\tau_\textrm{t}}}}{{\mu \prime }}}}} \right)} \right]\\ - b\left[ {({\mu - \mu \prime } )- ({{\tau_\textrm{t}} + \mu - \mu \prime } ){e^{ - \frac{{{\tau_\textrm{t}}}}{\mu }}} + \mu \prime \frac{{\mu \prime }}{{\mu + \mu \prime }}\left( {1 - {e^{ - \frac{{{\tau_\textrm{t}}}}{\mu } - \frac{{{\tau_\textrm{t}}}}{{\mu \prime }}}}} \right)} \right], \end{array}$$

(30)$$\begin{array}{l} {t_1}(\mu ,\mu \prime ) = \;\;a\left[ {\left( {1 - {e^{ - \frac{{{\tau_\textrm{t}}}}{\mu }}}} \right) - \frac{{\mu \prime }}{{\mu + \mu \prime }}\left( {1 - {e^{ - \frac{{{\tau_\textrm{t}}}}{\mu } - \frac{{{\tau_\textrm{t}}}}{{\mu \prime }}}}} \right)} \right]\\ + b\left[ {({\mu - \mu \prime } ){e^{ - \frac{{{\tau_\textrm{t}}}}{\mu }}} + {\tau_\textrm{t}} - ({\mu - \mu \prime } )- ({{\tau_\textrm{t}} + \mu \prime } )\frac{{\mu \prime }}{{\mu + \mu \prime }}\left( {1 - {e^{ - \frac{{{\tau_\textrm{t}}}}{\mu } - \frac{{{\tau_\textrm{t}}}}{{\mu \prime }}}}} \right)} \right], \end{array}$$

(31)$$\begin{aligned} \mu \ne \mu \prime \quad\quad{t_2}(\mu ,\mu \prime ) = \;\;a\left[ {\left( {1 - {e^{ - \frac{{{\tau_\textrm{t}}}}{\mu }}}} \right) - \frac{{\mu \prime }}{{\mu - \mu \prime }}\left( {{e^{ - \frac{{{\tau_\textrm{t}}}}{\mu }}} - {e^{ - \frac{{{\tau_\textrm{t}}}}{{\mu \prime }}}}} \right)} \right]\\ \;\;\;\; - b\left[ {({\mu + \mu \prime } ){e^{ - \frac{{{\tau_\textrm{t}}}}{\mu }}} + {\tau_\textrm{t}} - ({\mu + \mu \prime } )+ \mu \prime \frac{{\mu \prime }}{{\mu - \mu \prime }}\left( {{e^{ - \frac{{{\tau_\textrm{t}}}}{\mu }}} - {e^{ - \frac{{{\tau_\textrm{t}}}}{{\mu \prime }}}}} \right)} \right], \end{aligned}$$

(32)$$\begin{aligned} \mu = \mu \prime \quad{t_2}(\mu ,\mu \prime ) = \;\;a\left[ {\left( {1 - {e^{ - \frac{{{\tau_\textrm{t}}}}{\mu }}}} \right) - {e^{ - \frac{{{\tau_\textrm{t}}}}{\mu }}}\frac{{{\tau_\textrm{t}}}}{\mu }} \right]\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; - b\left[ {2\mu {e^{ - \frac{{{\tau_\textrm{t}}}}{\mu }}} + {\tau_\textrm{t}} - 2\mu + {e^{ - \frac{{{\tau_\textrm{t}}}}{\mu }}}{\tau_\textrm{t}}} \right], \end{aligned}$$

and for the exponential parameterization, the coefficients are:

(33)$$\begin{array}{c} {1 - b\mu \ne 0}\\ {1 - b\mu \prime \ne 0}\\ {\mu \ne \mu \prime } \end{array}\;,\;\;{r_1}(\mu ,\mu \prime ) = \;\;\frac{a}{{1 - b\mu \prime }}\left[ {\frac{{1 - {e^{ - (1 - b\mu )\frac{{{\tau_\textrm{t}}}}{\mu }}}}}{{1 - b\mu }} - {e^{b{\tau_\textrm{t}}}}\frac{{\mu \prime }}{{\mu - \mu \prime }}\left( {{e^{ - \frac{{{\tau_\textrm{t}}}}{\mu }}} - {e^{ - \frac{{{\tau_\textrm{t}}}}{{\mu \prime }}}}} \right)} \right],$$

(34)$$\begin{array}{c} {1 - b\mu = 0}\\ {1 - b\mu \prime \ne 0}\\ {\mu \ne \mu \prime } \end{array}\;,\;\;{r_1}(\mu ,\mu \prime ) = \;\;\frac{a}{{1 - \mu \prime /\mu }}\left[ {\frac{{{\tau_\textrm{t}}}}{\mu } - \frac{{\mu \prime }}{{\mu - \mu \prime }}\left( {1 - {e^{\frac{{{\tau_\textrm{t}}}}{\mu } - \frac{{{\tau_\textrm{t}}}}{{\mu \prime }}}}} \right)} \right],$$

(35)$$\begin{array}{c} {1 - b\mu \ne 0}\\ {1 - b\mu \prime = 0}\\ {\mu \ne \mu \prime } \end{array}\;,\;\;{r_1}(\mu ,\mu \prime ) = \;\;\frac{a}{{1 - \mu /\mu \prime }}\left[ {\frac{{{\tau_\textrm{t}}}}{{\mu \prime }} - \frac{\mu }{{\mu \prime - \mu }}\left( {1 - {e^{ - \frac{{{\tau_\textrm{t}}}}{\mu } + \frac{{{\tau_\textrm{t}}}}{{\mu \prime }}}}} \right)} \right],$$

(36)$$\begin{array}{c} {1 - b\mu \ne 0}\\ {1 - b\mu \prime \ne 0}\\ {\mu = \mu \prime } \end{array},\;\;\;{r_1}(\mu ,\mu \prime ) = \;\;\frac{a}{{1 - b\mu }}\left[ {\frac{{1 - {e^{ - (1 - b\mu )\frac{{{\tau_\textrm{t}}}}{\mu }}}}}{{1 - b\mu }} - {e^{ - (1 - b\mu )\frac{{{\tau_\textrm{t}}}}{\mu }}}\frac{{{\tau_\textrm{t}}}}{\mu }} \right],$$

(37)$$\begin{array}{c} {1 - b\mu = 0}\\ {1 - b\mu \prime = 0}\\ {\mu = \mu \prime } \end{array},\;\;\;{r_1}(\mu ,\mu \prime ) = \;\;a\frac{{\tau _\textrm{t}^2}}{{2{\mu ^2}}},$$

(38)$$1 - b\mu \ne 0,\;\;\;{r_2}(\mu ,\mu \prime ) = \;\;\frac{a}{{1 + b\mu \prime }}\left[ {\frac{{1 - {e^{ - (1 - b\mu )\frac{{{\tau_\textrm{t}}}}{\mu }}}}}{{1 - b\mu }} - \frac{{\mu \prime }}{{\mu + \mu \prime }}\left( {1 - {e^{ - \frac{{{\tau_\textrm{t}}}}{\mu } - \frac{{{\tau_\textrm{t}}}}{{\mu \prime }}}}} \right)} \right],$$

(39)$$1 - b\mu = 0,\;\;\;{r_2}(\mu ,\mu \prime ) = \;\;\frac{a}{{1 + \mu \prime /\mu }}\left[ {\frac{{{\tau_\textrm{t}}}}{\mu } - \frac{{\mu \prime }}{{\mu + \mu \prime }}\left( {1 - {e^{ - \frac{{{\tau_\textrm{t}}}}{\mu } - \frac{{{\tau_\textrm{t}}}}{{\mu \prime }}}}} \right)} \right],$$

(40)$$1 - b\mu \prime \ne 0,\;\;\;{t_1}(\mu ,\mu \prime ) = \;\;\frac{a}{{1 - b\mu \prime }}\left[ {{e^{b{\tau_\textrm{t}}}}\frac{{1 - {e^{ - (1 + b\mu )\frac{{{\tau_\textrm{t}}}}{\mu }}}}}{{1 + b\mu }} - {e^{b{\tau_\textrm{t}}}}\frac{{\mu \prime }}{{\mu + \mu \prime }}\left( {1 - {e^{ - \frac{{{\tau_\textrm{t}}}}{\mu } - \frac{{{\tau_\textrm{t}}}}{{\mu \prime }}}}} \right)} \right],$$

(41)$$1 - b\mu \prime = 0,\;\;\;{t_1}(\mu ,\mu \prime ) = \;\;\frac{a}{{1 + \mu /\mu \prime }}\left[ {\frac{\mu }{{\mu + \mu \prime }}\left( {{e^{\frac{{{\tau_\textrm{t}}}}{{\mu \prime }}}} - {e^{ - \frac{{{\tau_\textrm{t}}}}{\mu }}}} \right) - \frac{{{\tau_\textrm{t}}}}{{\mu \prime }}{e^{ - \frac{{{\tau_\textrm{t}}}}{\mu }}}} \right],$$

(42)$$\mu \ne \mu \prime \;\;\;{t_2}(\mu ,\mu \prime ) = \;\;\frac{a}{{1 + b\mu \prime }}\left[ {{e^{b{\tau_\textrm{t}}}}\frac{{1 - {e^{ - (1 + b\mu )\frac{{{\tau_\textrm{t}}}}{\mu }}}}}{{1 + b\mu }} - \frac{{\mu \prime }}{{\mu - \mu \prime }}\left( {{e^{ - \frac{{{\tau_\textrm{t}}}}{\mu }}} - {e^{ - \frac{{{\tau_\textrm{t}}}}{{\mu \prime }}}}} \right)} \right],$$

(43)$$\mu = \mu \prime \;\;\;{t_2}(\mu ,\mu \prime ) = \;\;\frac{a}{{1 + b\mu \prime }}\left[ {{e^{b{\tau_\textrm{t}}}}\frac{{1 - {e^{ - (1 + b\mu )\frac{{{\tau_\textrm{t}}}}{\mu }}}}}{{1 + b\mu }} - {e^{ - \frac{{{\tau_\textrm{t}}}}{\mu }}}\frac{{{\tau_\textrm{t}}}}{\mu }} \right].$$

The direct transmission of incident radiation is given through

(44)$$T_0^\ast (\mu ,\mu \prime ) = {T_0}(\mu ,\mu \prime ) = {e ^{ - \frac{{{\tau _\textrm{t}}}}{{\mu \prime }}}}\delta (\mu - \mu \prime ).$$

The reflection and transmission of single scattering associated with the incident radiation are given as

(45)$${R_1}(\mu ,\mu \prime ) = \frac{\omega }{{4({\mu + \mu \prime } )}}p(\mu , - \mu \prime )\left( {1 - {e^{ - \frac{{{\tau_\textrm{t}}}}{\mu } - \frac{{{\tau_\textrm{t}}}}{{\mu \prime }}}}} \right),$$

(46)$${T_1}(\mu ,\mu \prime ) = \frac{\omega }{{4({\mu - \mu \prime } )}}p( - \mu , - {\mu _0})\left( {{e^{ - \frac{{{\tau_\textrm{t}}}}{\mu }}} - {e^{ - \frac{{{\tau_\textrm{t}}}}{{\mu \prime }}}}} \right).$$

(47)$$R_1^\ast (\mu ,\mu \prime ) = {R_1}(\mu ,\mu \prime ),\;T_1^\ast (\mu ,\mu \prime ) = T(\mu ,\mu \prime ),$$

where ‘*’ denotes the incident direction is from bottom. The direct and single scattering given in Eqs. (20)–(47) can be used as the starting parameters of the doubling process.

The upwelling and downwelling radiances for a homogeneous layer are shown schematically in Fig. 2. The upwelling and downwelling radiances for the thermal layer can be expressed as:

(48)$$\left( {\begin{array}{c} {{I^ + }({\tau_1})}\\ {{I^ - }({\tau_2})} \end{array}} \right) = \left( {\begin{array}{cc} R&{{T^\ast }}\\ T&{{R^\ast }} \end{array}} \right)\left( {\begin{array}{c} {{I^ - }({\tau_1})}\\ {{I^ + }({\tau_2})} \end{array}} \right) + \left( {\begin{array}{c} {\Delta _\textrm{t}^ + }\\ {\Delta _\textrm{t}^ - } \end{array}} \right).$$

where ${I^ - }({\tau _1})$ and ${I^ + }({\tau _2})$ are the incoming radiances from the layer top and bottom, ${I^ + }({\tau _1})$ and ${I^ - }({\tau _2})$ are the outgoing radiances from the layer top and bottom, $\Delta _t^ + $ and $\Delta _t^ - $ are the internal radiation at the layer top and bottom, and all the zenith angle argument is suppressed. Radiance is represented by a vector, and the reflection and transmission are represented by matrices. The matrix-matrix or matrix-vector multiplication represent the discretization of a convolution as

(49)$$Z(\mu ,\mu \prime ) = 2\int_0^1 {X(\mu ,\mu \prime \prime )Y(\mu \prime \prime ,\mu \prime )\mu \prime \prime d\mu \prime \prime } ,$$

where X is a matrix, Y is a matrix or vector, and Z has the same dimension as Y.

Fig. 2. Schematic figure of upwelling and downwelling radiances in a homogeneous thermal layer. The linear parameterization inside the layer is assumed for denoting the blackbody radiation.

Download Full Size | PDF

The adding principle for the thermal radiation is shown in Fig. 3, where the quantities in the top layer, the bottom layer, and the added layer are denoted by subscript “a”, “b”, “a + b”, respectively. The top, middle, and bottom levels are at optical depth ${\tau _1}$, ${\tau _2}$, and ${\tau _3}$, where ${\tau _2}$ is shown both in top and bottom layers for clarity in Fig. 3. According to the adding principle, the upwelling and downwelling radiances for two thermal layers are derived as [1]:

(50)$$\left( {\begin{array}{c} {{I^ + }({\tau_1})}\\ {{I^ - }({\tau_3})} \end{array}} \right) = \left( {\begin{array}{cc} {{R_{a + b}}}&{T_{a + b}^\ast }\\ {{T_{a + b}}}&{R_{a + b}^\ast } \end{array}} \right)\left( {\begin{array}{c} {{I^ - }({\tau_1})}\\ {{I^ + }({\tau_3})} \end{array}} \right) + \left( {\begin{array}{c} {\Delta _{a + b}^ + }\\ {\Delta _{a + b}^ - } \end{array}} \right),$$

where the incoming radiances are from top and bottom levels, and

(51)$$\left( {\begin{array}{cc} {{R_{a + b}}}&{T_{a + b}^\ast }\\ {{T_{a + b}}}&{R_{a + b}^\ast } \end{array}} \right) = \left( {\begin{array}{cc} {{R_a} + T_a^\ast {{(I - {R_b}R_a^\ast )}^{ - 1}}{R_b}{T_a}}&{T_a^\ast {{(I - {R_b}R_a^\ast )}^{ - 1}}T_b^\ast }\\ {{T_b}{{(I - R_a^\ast {R_b})}^{ - 1}}{T_a}}&{R_b^\ast{+} {T_b}{{(I - R_a^\ast {R_b})}^{ - 1}}R_a^\ast T_b^\ast } \end{array}} \right),$$

(52)$$\left( {\begin{array}{c} {\Delta _{a + b}^ + }\\ {\Delta _{a + b}^ - } \end{array}} \right) = \left( {\begin{array}{c} {\Delta _a^ + }\\ {\Delta _b^ - } \end{array}} \right) + \left( {\begin{array}{cc} {T_a^\ast {{(I - {R_b}R_a^\ast )}^{ - 1}}{R_b}}&{T_a^\ast {{(I - {R_b}R_a^\ast )}^{ - 1}}}\\ {{T_b}{{(I - R_a^\ast {R_b})}^{ - 1}}}&{{T_b}{{(I - R_a^\ast {R_b})}^{ - 1}}R_a^\ast } \end{array}} \right)\left( {\begin{array}{c} {\Delta _a^ - }\\ {\Delta _b^ + } \end{array}} \right).$$

Fig. 3. Schematic figure of the adding principle for two homogeneous layers.

Download Full Size | PDF

When the incoming radiances are the same from top and bottom level, the upwelling and downwelling radiances at the middle layer are:

(53)$$\left( {\begin{array}{c} {{I^ + }({\tau_2})}\\ {{I^ - }({\tau_2})} \end{array}} \right) = \left( {\begin{array}{cc} {{R_{{\tau_2}}}}&{T_{{\tau_2}}^\ast }\\ {{T_{{\tau_2}}}}&{R_{{\tau_2}}^\ast } \end{array}} \right)\left( {\begin{array}{c} {{I^ - }({\tau_1})}\\ {{I^ + }({\tau_3})} \end{array}} \right) + \left( {\begin{array}{c} {\Delta _{{\tau_2}}^ + }\\ {\Delta _{{\tau_2}}^ - } \end{array}} \right),$$

(54)$$\left( {\begin{array}{cc} {{R_{{\tau_2}}}}&{T_{{\tau_2}}^\ast }\\ {{T_{{\tau_2}}}}&{R_{{\tau_2}}^\ast } \end{array}} \right) = \left( {\begin{array}{cc} {{{(I - {R_b}R_a^\ast )}^{ - 1}}{R_b}{T_a}}&{{{(I - {R_b}R_a^\ast )}^{ - 1}}T_b^\ast }\\ {{{(I - R_a^\ast {R_b})}^{ - 1}}{T_a}}&{{{(I - R_a^\ast {R_b})}^{ - 1}}R_a^\ast T_b^\ast } \end{array}} \right),$$

(55)$$\left( {\begin{array}{c} {\Delta _{{\tau_2}}^ + }\\ {\Delta _{{\tau_2}}^ - } \end{array}} \right) = \left( {\begin{array}{cc} {{{(I - {R_b}R_a^\ast )}^{ - 1}}{R_b}}&{{{(I - {R_b}R_a^\ast )}^{ - 1}}}\\ {{{(I - R_a^\ast {R_b})}^{ - 1}}}&{{{(I - R_a^\ast {R_b})}^{ - 1}}R_a^\ast } \end{array}} \right)\left( {\begin{array}{c} {\Delta _a^ - }\\ {\Delta _b^ + } \end{array}} \right).$$

When the bottom layer is a surface, the schematic figure for the adding principle is shown in Fig. 4. The surface can only provide reflection R_b and the surface thermal radiation $\Delta _b^ + $. The adding steps in Eqs. (50)–(55) can be further simplified. For the adding of a layer and a surface boundary, the upwelling radiance can be expressed as

(56)$$\begin{aligned} {I^ + }({\tau _1}) &= \quad[{{R_a} + T_a^\ast {{(I - {R_b}R_a^\ast )}^{ - 1}}{R_b}{T_a}} ]{I^ - }({\tau _1})\\ & + \Delta _a^ +{+} [{T_a^\ast {{(I - {R_b}R_a^\ast )}^{ - 1}}{R_b}} ]\Delta _a^ -{+} [{T_a^\ast {{(I - {R_b}R_a^\ast )}^{ - 1}}} ]\Delta _b^ + . \end{aligned}$$

Fig. 4. Schematic figure of the adding principle for a layer and a surface boundary.

Download Full Size | PDF

The upwelling and downwelling radiances at the surface boundary are

(57)$$ \begin{aligned} I^{+}\left(\tau_{2}\right)=& {\left[\left(I-R_{b} R_{a}^{*}\right)^{-1} R_{b} T_{a}\right] I^{-}\left(\tau_{1}\right) } \\ &+\left[\left(I-R_{b} R_{a}^{*}\right)^{-1} R_{b}\right] \Delta_{a}^{-}+\left[\left(I-R_{b} R_{a}^{*}\right)^{-1}\right] \Delta_{b}^{+}, \end{aligned} $$

(58)$$ \begin{aligned} I^{-}\left(\tau_{2}\right)=& {\left[\left(I-R_{a}^{*} R_{b}\right)^{-1} T_{a}\right] I^{-}\left(\tau_{1}\right) } \\ &+\left[\left(I-R_{a}^{*} R_{b}\right)^{-1}\right] \Delta_{a}^{-}+\left[\left(I-R_{a}^{*} R_{b}\right)^{-1} R_{a}^{*}\right] \Delta_{b}^{+} \end{aligned} $$

The reflection and transmission convolution in Eq. (49) can be calculated using Gaussian or Radau-Gaussian quadrature and then be transferred into super-matrix multiplication. The matrix element for the super-matrix can be expressed as [45]:

(59)$${X_{i,j}} = {c_i}{X_{ij}}({\mu _i},{\mu _j}){c_j},\;\;{c_i} = \left\{ {\begin{array}{c} {\sqrt {2{\omega_i}{\mu_i}} ,\;}\\ {1,\quad} \end{array}} \right.\begin{array}{c} {1 \le i \le n}\\ {(n + 1) \le i \le (n + {n_0})} \end{array},$$

where ${\omega _i}$ is the quadrature weight, n is the quadrature number, and n₀ is the number of viewing zenith angles. For convenience of calculation, such a matrix can be divided into four parts:

(60)$$X = \left( {\begin{array}{cc} {{X_{11}}(n \times n)}&{{X_{12}}(n \times {n_0})}\\ {{X_{21}}({n_0} \times n)}&{{X_{22}}({n_0} \times {n_0})} \end{array}} \right).$$

Given this notation super-matrix multiplication using the super-matrix can then be written as

(61)$$XY = \left( {\begin{array}{cc} {{X_{11}}{Y_{11}}}&{{X_{11}}{Y_{12}}}\\ {{X_{21}}{Y_{11}}}&{{X_{21}}{Y_{12}}} \end{array}} \right).$$

Matrix inversions can be directly expressed as

(62)$${(I - X)^{ - 1}} = \left( {\begin{array}{cc} {{{(I - {X_{11}})}^{ - 1}}}&{{{(I - {X_{11}})}^{ - 1}}{X_{12}}}\\ {{X_{21}}{{(I - {X_{11}})}^{ - 1}}}&{I + {X_{22}} + {X_{21}}{{(I - {X_{11}})}^{ - 1}}{X_{12}}} \end{array}} \right).$$

The doubling process is actually to repeat the adding process when the top and bottom layers exactly have the same quantities. The doubling principle is illustrated in Fig. 5. Consequently, the reflection and transmission matrices after doubling are then given by:

(63)$$\left( {\begin{array}{cc} {{R_{n + 1}}}&{T_{n + 1}^\ast }\\ {{T_{n + 1}}}&{R_{n + 1}^\ast } \end{array}} \right) = \left( {\begin{array}{cc} {{R_n} + T_n^\ast {{(I - {R_n}R_n^\ast )}^{ - 1}}{R_n}{T_n}}&{T_n^\ast {{(I - {R_n}R_n^\ast )}^{ - 1}}T_n^\ast }\\ {{T_n}{{(I - R_n^\ast {R_n})}^{ - 1}}{T_n}}&{R_n^\ast{+} {T_n}{{(I - R_n^\ast {R_n})}^{ - 1}}R_n^\ast T_n^\ast } \end{array}} \right),$$

Fig. 5. Schematic figure of the doubling principle, where the numbers are the multiple of the starting layer and can be denoted as ${2^n}$(n is the doubling power).

Download Full Size | PDF

With internal blackbody radiation, the doubling principle cannot be directly applied because the Planck function is assumed to vary in linear or exponential parameterizations with optical depth. The doubling realization is associated with different parameterizations. For linear parameterization, the doubling process is illustrated in Fig. 6, where the Planck function decomposition is shown in the layers. The corresponding doubling process can be formulated as

(64)$$\begin{array}{l} \left( {\begin{array}{c} {\Gamma _{n + 1}^ + }\\ {\Gamma _{n + 1}^ - } \end{array}} \right) = \left( {\begin{array}{c} {\Gamma _n^ + }\\ {\Gamma _n^ - } \end{array}} \right) + \left( {\begin{array}{cc} {T_n^\ast {{(I - {R_n}R_n^\ast )}^{ - 1}}{R_n}}&{T_n^\ast {{(I - {R_n}R_n^\ast )}^{ - 1}}}\\ {{T_n}{{(I - R_n^\ast {R_n})}^{ - 1}}}&{{T_n}{{(I - R_n^\ast {R_n})}^{ - 1}}R_n^\ast } \end{array}} \right)\left( {\begin{array}{c} {\Gamma _n^ - }\\ {\Gamma _n^ + } \end{array}} \right),\\ \left( {\begin{array}{c} {\Delta _{n + 1}^ + }\\ {\Delta _{n + 1}^ - } \end{array}} \right) = \left( {\begin{array}{c} {\Delta _n^ + }\\ {\Delta _n^ -{+} \Gamma _n^ - } \end{array}} \right) + \left( {\begin{array}{cc} {T_n^\ast {{(I - {R_n}R_n^\ast )}^{ - 1}}{R_n}}&{T_n^\ast {{(I - {R_n}R_n^\ast )}^{ - 1}}}\\ {{T_n}{{(I - R_n^\ast {R_n})}^{ - 1}}}&{{T_n}{{(I - R_n^\ast {R_n})}^{ - 1}}R_n^\ast } \end{array}} \right)\left( {\begin{array}{c} {\Delta _n^ - }\\ {\Delta _n^ +{+} \Gamma _n^ + } \end{array}} \right), \end{array}$$

where ${\Gamma ^ \pm }$ is the doubling process for $b\Delta $ and the doubling process for linear parameterization is the superposition of two doubling processes. The doubling process and the Planck function decomposition for exponential parameterization is illustrated in Fig. 7. For the doubling process, we have instead:

(65)$$\left( {\begin{array}{c} {\Delta _{n + 1}^ + }\\ {\Delta _{n + 1}^ - } \end{array}} \right) = \left( {\begin{array}{c} {\Delta _n^ + }\\ {{e^{{2^n}b\Delta \tau }} \cdot \Delta _n^ - } \end{array}} \right) + \left( {\begin{array}{cc} {T_n^\ast {{(I - {R_n}R_n^\ast )}^{ - 1}}{R_n}}&{T_n^\ast {{(I - {R_n}R_n^\ast )}^{ - 1}}}\\ {{T_n}{{(I - R_n^\ast {R_n})}^{ - 1}}}&{{T_n}{{(I - R_n^\ast {R_n})}^{ - 1}}R_n^\ast } \end{array}} \right)\left( {\begin{array}{c} {\Delta _n^ - }\\ {{e^{{2^n}b\Delta \tau }} \cdot \Delta _n^ + } \end{array}} \right).$$

Fig. 6. Doubling realization for linear parameterization of the Planck function.

Download Full Size | PDF

Fig. 7. Doubling realization for exponential parameterization of the Planck function.

Download Full Size | PDF

Eventually, Eq. (1) in the scenario shown in Fig. 1 can numerically solved by combining the doubling and adding processes described in this section.

3. Results and discussions

3.1 MATLAB program

The algorithm described in section 2 is implemented in a state-of-the-art MATLAB program. The equations and figures in section 2 are realized by the MATLAB functions as shown in Table 1.

Table 1. MATLAB functions and their corresponding equations and figures

View Table

The program is composed of three processes: the input, adding-doubling, and output as indicated in Fig. 8. The schematic figure is shown in Fig. 1. The input parameters include the wavenumber limits ${\nu _1}\textrm{ and }{\nu _2}$, the number of layers n, the profiles of temperature, single scattering albedo, optical thickness, and scattering phase function, the surface albedo, and the surface temperature. In practice, it is the Legendre polynomial expansion coefficients $(2n + 1){\alpha _n}$ in Eq. (4) are used as input, instead of the actual scattering phase functions. There is no scattering when the single scattering albedo is equal to zero in certain layer. The Lambertian surface is employed in this study. The radiative transfer in a single layer is calculated by direct computation or with the doubling process, depending on the single scattering albedo as zero or not. The radiative transfer in multiple layers is derived with the adding process. The radiances at user-defined directions and at layer boundaries or partial levels within a layer can be directly computed and output using the adding process. The flux or actinic flux can be extracted easily once the radiances are obtained. The radiance and flux profiles can be directly plotted in the MATLAB program. The scalar thermal radiation using adding-doubling algorithm in MATLAB has the acronym “STRADM”.

Fig. 8. Flowchart for the adding-doubling MATLAB program.

Download Full Size | PDF

3.2 Verification and discussions

As an example application, the atmosphere is divided into 23 layers and the input temperature profile and two scattering phase functions used in cases 1 and 2 are shown in Fig. 9. The scattering phase functions denoted as phasef1 and phasef2 in Fig. 9(b) represent weak and strong scattering, respectively. The wavenumber limits are 2499.5 cm^-1 and 2500.5 cm^-1. An Intel Core i7-8700CPU @ 3.20GHZ Dell PC computer is used for the computation. The expansion coefficients in Eq. (6) associated with phasef1 and phasef2 in Fig. 9(b) are precomputed and saved. Verification comparisons are carried out between the STRADM and the discrete ordinate radiative transfer code (denoted as DORT in this study) (Spurr, 2006). The linear Planck function parameterization is used in cases 1 and 2. The comparisons between linear and exponential parameterizations are given in case 3 for a single layer.

a) case 1: the single scattering albedo is 0.402856 for all layers, the optical depths are 2i (i, the layer number enumerating from top as shown in Fig. 1), and the strong scattering phase function denoted in phasef2 in Fig. 9(b) is used for the scattering phase function in all layers. The surface albedo and temperature are 0.5 and 300K, respectively. The surface emission is present. The Gaussian quadrature number n is 200 and the viewing zenith angle number n₀ is 40. The initial optical depth for doubling is 1.0e-6. The computational time for case 1 is 2.56s.
The flux, actinic flux, and the relative differences between the STRADM and the DORT are shown in Fig. 10. The upper panels show the flux and actinic flux comparisons, where the legends flux_up and flux_dn denote upwelling and downwelling fluxes, respectively, the legends aflux_up and aflux_dn denote upwelling and downwelling actinic fluxes, respectively, and the legends _D and _M denote DORT and STRADM simulations, respectively. The lower panels shown the flux relative difference between DORT and STRADM results, which is defined as (flux(DORT)-flux(STRADM))/flux(DORT)*1000, so does the actinic flux relative difference. All differences are less than 0.01%.
Figures 11 and 12 show the radiances and their relative differences for output at the 24 level boundaries. The legends D and M denote the DORT and the STRADM simulations too. The radiance relative difference definition is similar as the flux relative difference definition in Fig. 10. The radiance comparisons between the two models are shown all on the left y-axes and the relative differences are shown all on the right y-axes. The unit for the differences is 0.01% and once again, all the differences are less than 0.01%.
b) case 2: the single scattering albedos in layers 16 and 21 are 0.6 and 0.4, respectively, and in all other layers these are zero. Optical depths in layer 16 and 21 are 15 and 5 respectively, and in all other layers are 1; phase functions phasef1 and phasef2 in Fig. 9(b) are used for layers 21 and 16, respectively. The surface emission is setting on. The surface albedo and temperature are 0.5 and 300K, respectively. The Gaussian quadrature number n is 200 and the viewing zenith angle number n₀ is 40. The initial optical depth for doubling is around 1.0e-6. The computational time for this case is 0.42s.
The legends in Fig. 13 are the same as the ones in Fig. 10. The relative differences are still less than 0.01%. The upwelling, downwelling, and net radiance profiles are shown in Fig. 14, where net radiance is defined as the upwelling radiance minus the downwelling radiance. The left and right color bars are applied to the upwelling/downwelling radiances and net radiance profiles, respectively. There is a strong radiation region around the surface for the upwelling radiance profile because the surface emission is on and the surface temperature is high. At top level, there is strong net radiance because the downwelling radiance here is zero.
For verification purposes, case 1 is an extreme scenario in which all layers are assigned strong scattering. The computational time for the extreme situations is around several seconds. Case 2 is a normal scenario, where a general case can actually be finished within 1 second using STRADM.
c) Case 3: Only one layer is used to compare the linear and exponential parameterizations of the Planck function. The bottom temperature is constant for the following figures, where the horizontal axes in Figs. 15 and 16 show the boundary temperature difference (bottom minus top temperature) and the vertical axes indicate the layer optical depth. All the differences are defined as (flux(exponential)-flux(linear))/flux(linear)*100.

Fig. 9. Temperature profile and two scattering phase functions used in cases 1 and 2.

Download Full Size | PDF

Fig. 10. Flux, actinic flux, and relative differences between the STRADM and DORT.

Download Full Size | PDF

Fig. 11. Radiances and relative differences between the STRADM and the DORT for the upper 12 levels enumerated from top.

Download Full Size | PDF

Fig. 12. Radiances and relative differences between the STRADM and the DORT for the lower 12 levels enumerated from top.

Download Full Size | PDF

Fig. 13. Flux, actinic flux, and relative differences between the DORT and the STRADM.

Download Full Size | PDF

Fig. 14. Upwelling, downwelling, and net radiance profiles, where the net radiance profile is defined as. $I(\tau ,\mu ) - I(\tau , - \mu )$

Download Full Size | PDF

Fig. 15. Relative differences between the linear and exponential parameterizations with the bottom temperature as 288K. Panels (a)-(c) are for no scattering, scattering phasef1 with single scattering albedo 0.6, and scattering phasef2 with single scattering albedo 0.6, respectively.

Download Full Size | PDF

Fig. 16. All the same as Fig. 15 except that the bottom temperature as 270 K.

Download Full Size | PDF

Figures 15 and 16 show the differences between the exponential and linear parameterizations for bottom temperature as 288K and 270K with no scattering for panel (a), scattering with phasef1 and single scattering albedo 0.6 for panel (b), and scattering with phasef2 and single scattering albedo 0.6 for panel (b). The differences are mainly influenced by the temperature difference, where the difference becomes large with the temperature difference and the difference can reach around 6% with 15 K difference. For large optical depth, the differences are converged to a small value, which can be obtained from Eqs. (21) and (22). The incorporation of scattering makes the convergence slower as shown in panels (b) and (c) than no scattering. Even if the temperature difference is the same, the lower bottom temperature case has larger differences by comparing Figs. 15 and 16.

Figure 17 shows the relative differences of two parameterizations for fixed top and bottom temperatures with different optical depths, single scattering albedos, and scattering phase functions. The differences of two parameterizations are mainly affected by the temperature difference. The differences with optical depth increase are convergent but the convergent rate is reduced with the enhanced scattering.

Fig. 17. Relative differences between linear and exponential parameterizations with bottom temperature as 270 K and top temperature as 255 K (a) for phasef1 and (b) for phasef2.

Download Full Size | PDF

4. Conclusion

Thermal radiation using the adding-doubling algorithm in scalar radiative transfer is derived in detail in this study. The Planck function is integrated over the wavenumber range and expressed analytically in term of logarithmic polynomial functions. The Planck function within a layer can be parameterized either with linear, exponential, and constant approximations. The doubling algorithms for these different parameterizations are derived in detail. The radiances at any zenith angles can be directly obtained by matrix inversion. The algorithm is implemented into a state-of-the-art MATLAB program, the scalar thermal radiation using adding-doubling algorithm in MATLAB (STRADM). The STRADM is easy to use and is public available, obtained either from the website at https://faculty.fudan.edu.cn/bqsun/zh_CN/yjgk/411343/list/index.htm or by contacting the corresponding author. The program is validated against discrete ordinate radiative transfer code and differences are always less than 0.01%. The STRADM program is highly efficient, with a general case calculated within 1 second. Differences between the linear and the exponential black-body parameterizations mainly depend on the temperature difference. Larger temperature differences can lead to bigger differences between results from the two parameterizations, and it can reach 6% when the temperature difference is in 15K. The difference is convergent to a small value with the increase of optical depth and the convergent rate is reduced with enhanced scattering.

Funding

Natural Science Foundation of Shanghai (19ZR1404100); National Natural Science Foundation of China (41975021); Shanghai Science and Technology Development Foundation (20dz1200700).

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. R. M. Goody and Y. L. Yung, Atmospheric radiation: theoretical basis (Oxford University, 1995).

2. K.-N. Liou, An introduction to atmospheric radiation (Elsevier, 2002).

3. H. C. van de Hulst, Multiple light scattering: tables, formulas, and applications (Academic Press, 1980), Vol. 1.

4. S. Chandrasekhar, Radiative transfer (Dover Publications, 1960).

5. L. Tsang, J. A. Kong, and R. T. Shin, Theory of microwave remote sensing (Wiley, 1985).

6. M. I. Mishchenko, “Vector radiative transfer equation for arbitrarily shaped and arbitrarily oriented particles: a microphysical derivation from statistical electromagnetics,” Appl. Opt. 41(33), 7114–7134 (2002). [CrossRef]

7. A. Marshak and A. Davis, 3D radiative transfer in cloudy atmospheres (Springer, 2005).

8. S. A. Clough, M. J. Iacono, and J. L. Moncet, “Line-by-line calculations of atmospheric fluxes and cooling rates: Application to water vapor,” J. Geophys. Res. 97(D14), 15761–15785 (1992). [CrossRef]

9. A. A. Lacis and V. Oinas, “A description of the correlated k distribution method for modeling nongray gaseous absorption, thermal emission, and multiple scattering in vertically inhomogeneous atmospheres,” J. Geophys. Res. 96(D5), 9027–9063 (1991). [CrossRef]

10. Q. Fu and K. Liou, “On the correlated k-distribution method for radiative transfer in nonhomogeneous atmospheres,” J. Atmos. Sci. 49(22), 2139–2156 (1992). [CrossRef]

11. H. Zhang, T. Nakajima, G. Shi, T. Suzuki, and R. Imasu, “An optimal approach to overlapping bands with correlated k distribution method and its application to radiative calculations,” J. Geophys. Res. Atmos. 108 (2003).

12. L. S. Rothman, I. E. Gordon, A. Barbe, D. C. Benner, P. F. Bernath, M. Birk, V. Boudon, L. R. Brown, A. Campargue, and J.-P. Champion, “The HITRAN 2008 molecular spectroscopic database,” J. Quant. Spectrosc. Radiat. Transfer 110(9-10), 533–572 (2009). [CrossRef]

13. H. C. van de Hulst, Light scattering by small particles (John Wiley and Sons, 1957).

14. B. A. Baum, P. Yang, A. J. Heymsfield, S. Platnick, M. D. King, Y. Hu, and S. T. Bedka, “Bulk scattering properties for the remote sensing of ice clouds. Part II: Narrowband models,” J. Appl. Meteorol. 44(12), 1896–1911 (2005). [CrossRef]

15. G. Mie, “Sättigungsstrom und Stromkurve einer schlecht leitenden Flüssigkeit,” Ann. Phys. 331(8), 597–614 (1908). [CrossRef]

16. M. I. Mishchenko, L. D. Travis, and A. A. Lacis, Scattering, absorption, and emission of light by small particles (Cambridge University, 2002).

17. A. Taflove and S. C. Hagness, Computational electrodynamics: the finite-difference time-domain method (Artech House, 2005).

18. K. Yee, “Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media,” IEEE Trans. Antennas Propag. 14(3), 302–307 (1966). [CrossRef]

19. P. Yang and K. Liou, “Finite-difference time domain method for light scattering by small ice crystals in three-dimensional space,” J. Opt. Soc. Am. A 13(10), 2072–2085 (1996). [CrossRef]

20. B. T. Draine and P. J. Flatau, “Discrete-dipole approximation for scattering calculations,” J. Opt. Soc. Am. A 11(4), 1491–1499 (1994). [CrossRef]

21. M. A. Yurkin and A. G. Hoekstra, “The discrete dipole approximation: an overview and recent developments,” J. Quant. Spectrosc. Radiat. Transfer 106(1-3), 558–589 (2007). [CrossRef]

22. P. Waterman, “Matrix formulation of electromagnetic scattering,” Proc. IEEE 53(8), 805–812 (1965). [CrossRef]

23. P. C. Waterman, “Symmetry, unitarity, and geometry in electromagnetic scattering,” Phys. Rev. D 3(4), 825–839 (1971). [CrossRef]

24. M. I. Mishchenko, L. D. Travis, and D. W. Mackowski, “T-matrix method and its applications to electromagnetic scattering by particles: A current perspective,” J. Quant. Spectrosc. Radiat. Transfer 111(11), 1700–1703 (2010). [CrossRef]

25. L. Bi, P. Yang, G. W. Kattawar, and M. I. Mishchenko, “Efficient implementation of the invariant imbedding T-matrix method and the separation of variables method applied to large nonspherical inhomogeneous particles,” J. Quant. Spectrosc. Radiat. Transfer 116, 169–183 (2013). [CrossRef]

26. B. Sun, L. Bi, P. Yang, M. Kahnert, and G. Kattawar, Invariant Imbedding T-matrix method for light scattering by nonspherical and inhomogeneous particles (Elsevier, 2019).

27. A. Macke, M. I. Mishchenko, and B. Cairns, “The influence of inclusions on light scattering by large ice particles,” J. Geophys. Res. 101(D18), 23311–23316 (1996). [CrossRef]

28. K. Muinonen, “Scattering of light by crystals: a modified Kirchhoff approximation,” Appl. Opt. 28(15), 3044–3050 (1989). [CrossRef]

29. P. Yang and K. Liou, “Geometric-optics–integral-equation method for light scattering by nonspherical ice crystals,” Appl. Opt. 35(33), 6568–6584 (1996). [CrossRef]

30. P. Yang and K. Liou, “Light scattering by hexagonal ice crystals: solutions by a ray-by-ray integration algorithm,” J. Opt. Soc. Am. A 14(9), 2278–2289 (1997). [CrossRef]

31. L. Bi, P. Yang, G. W. Kattawar, Y. Hu, and B. A. Baum, “Scattering and absorption of light by ice particles: solution by a new physical-geometric optics hybrid method,” J. Quant. Spectrosc. Radiat. Transfer 112(9), 1492–1508 (2011). [CrossRef]

32. B. Sun, P. Yang, G. W. Kattawar, and X. Zhang, “Physical-geometric optics method for large size faceted particles,” Opt. Express 25(20), 24044–24060 (2017). [CrossRef]

33. K. F. Ren, F. Onofri, C. Rozé, and T. Girasole, “Vectorial complex ray model and application to two-dimensional scattering of plane wave by a spheroidal particle,” Opt. Lett. 36(3), 370–372 (2011). [CrossRef]

34. B. Sun, G. W. Kattawar, P. Yang, and K. F. Ren, “Rigorous 3-D vectorial complex ray model applied to light scattering by an arbitrary spheroid,” J. Quant. Spectrosc. Radiat. Transfer 179, 1–10 (2016). [CrossRef]

35. Q. Duan, X. e. Han, S. Idlahcen, and K. F. Ren, “Three-dimensional light scattering by a real liquid jet: VCRM simulation and experimental validation,” J. Quant. Spectrosc. Radiat. Transfer 239, 106677 (2019). [CrossRef]

36. J. Lenoble, M. Herman, J. Deuzé, B. Lafrance, R. Santer, and D. Tanré, “A successive order of scattering code for solving the vector equation of transfer in the earth's atmosphere with aerosols,” J. Quant. Spectrosc. Radiat. Transfer 107(3), 479–507 (2007). [CrossRef]

37. M. Duan, Q. Min, and D. Lü, “A polarized radiative transfer model based on successive order of scattering,” Adv. Atmos. Sci. 27(4), 891–900 (2010). [CrossRef]

38. P.-W. Zhai, Y. Hu, C. R. Trepte, and P. L. Lucker, “A vector radiative transfer model for coupled atmosphere and ocean systems based on successive order of scattering method,” Opt. Express 17(4), 2057–2079 (2009). [CrossRef]

39. K. Stamnes, S.-C. Tsay, W. Wiscombe, and I. Laszlo, “DISORT, a general-purpose Fortran program for discrete-ordinate-method radiative transfer in scattering and emitting layered media: documentation of methodology,” (2000).

40. R. J. Spurr, “VLIDORT: A linearized pseudo-spherical vector discrete ordinate radiative transfer code for forward model and retrieval studies in multilayer multiple scattering media,” J. Quant. Spectrosc. Radiat. Transfer 102(2), 316–342 (2006). [CrossRef]

41. G. N. Plass, G. W. Kattawar, and F. E. Catchings, “Matrix operator theory of radiative transfer. 1: Rayleigh scattering,” Appl. Opt. 12(2), 314–329 (1973). [CrossRef]

42. G. W. Kattawar, G. N. Plass, and F. E. Catchings, “Matrix operator theory of radiative transfer. 2: Scattering from maritime haze,” Appl. Opt. 12(5), 1071–1084 (1973). [CrossRef]

43. G. N. Plass and G. W. Kattawar, “Monte carlo calculations of radiative transfer in the earth's atmosphere-ocean system: I. flux in the atmosphere and ocean,” J. Phys. Oceanogr. 2(2), 139–145 (1972). [CrossRef]

44. P.-W. Zhai, G. W. Kattawar, and P. Yang, “Impulse response solution to the three-dimensional vector radiative transfer equation in atmosphere-ocean systems. I. Monte Carlo method,” Appl. Opt. 47(8), 1037–1047 (2008). [CrossRef]

45. J. F. de Haan, P. Bosma, and J. Hovenier, “The adding method for multiple scattering calculations of polarized light,” Astron. Astrophys. 183, 371–391 (1987).

46. J. Hovenier, “Multiple scattering of polarized light in planetary atmospheres,” Astron. Astrophys. 13, 7 (1971).

47. W. Wiscombe, “Extension of the doubling method to inhomogeneous sources,” J. Quant. Spectrosc. Radiat. Transfer 16(6), 477–489 (1976). [CrossRef]

48. J. Lenoble, Radiative transfer in scattering and absorbing atmospheres: standard computational procedures (A. Deepak Publishing, 1985), Vol. 300.

49. Q. Liu and F. Weng, “Advanced doubling–adding method for radiative transfer in planetary atmospheres,” J. Atmos. Sci. 63(12), 3459–3465 (2006). [CrossRef]

50. Z. Zhang, P. Yang, G. Kattawar, H.-L. A. Huang, T. Greenwald, J. Li, B. A. Baum, D. K. Zhou, and Y. Hu, “A fast infrared radiative transfer model based on the adding–doubling method for hyperspectral remote-sensing applications,” J. Quant. Spectrosc. Radiat. Transfer 105(2), 243–263 (2007). [CrossRef]

51. W. Bai, P. Zhang, W. Zhang, G. Ma, and C. Qi, “A model for accurately calculating hyper-spectral, middle-shortwave infrared radiative transfer for remote sensing,” Sci. China Earth Sci. 61(3), 317–326 (2018). [CrossRef]

52. J. Ding, P. Yang, M. D. King, S. Platnick, X. Liu, K. G. Meyer, and C. Wang, “A fast vector radiative transfer model for the atmosphere-ocean coupled system,” J. Quant. Spectrosc. Radiat. Transfer 239, 106667 (2019). [CrossRef]

53. W. Li, F. Zhang, Y.-N. Shi, H. Iwabuchi, M. Zhu, J. Li, W. Han, H. Letu, and H. Ishimoto, “Efficient radiative transfer model for thermal infrared brightness temperature simulation in cloudy atmospheres,” Opt. Express 28(18), 25730–25749 (2020). [CrossRef]

zenphasef.m	Eq. (8)
planckf_band.m	Eq. (10)
noscattering_lin.m	Eqs. (20), (21), (23), and (24)
noscattering_exp.m	Eqs. (20), (22), (23), and (25)
coefsingle_lin.m	Eqs. (27)-(32)
coefsingle_exp.m	Eqs. (33)-(43)
doubling_lin.m	Eqs. (44)–(47), (64), Figs. 5 and 6
doubling_exp.m	Eqs. (44)–(47), (65), Figs. 5 and 7
adding.m	Eqs. (50)-(52), Fig. 3
addingm.m	Eqs. (53)–(55), Fig. 3
addingb.m	Eq. (56), Fig. 4
addingbm.m	Eqs. (57)–(58), Fig. 4

Scalar thermal radiation using the adding-doubling method

Abstract

1. Introduction

2. Method

3. Results and discussions

3.1 MATLAB program

3.2 Verification and discussions

4. Conclusion

Funding

Disclosures

Data availability

References

Data availability

Cited By

Figures (17)

Tables (1)

Equations (65)

Optics Express