Perturbation Method for Solar/Infrared Radiative Transfer in a Scattering Medium with Vertical Inhomogeneity in Internal Optical Properties

Abstract

A new scheme based on perturbation method is presented to solve the problem of solar/infrared radiative transfer (SRT/IRT) in a scattering medium, in which the inherent optical properties (IOPs) are vertically inhomogeneous. The Eddington approximation for SRT and the two-stream approximation for IRT are used as the zeroth-order solution, and multiple-scattering effect of inhomogeneous IOPs is included in the first-order solution. Observations show that the stratocumulus clouds are vertically inhomogeneous, and the accuracy of SRT/IRT for stratocumulus clouds by different solutions is evaluated. In the spectral band of 0.25–0.69 μm, the relative error in absorption with inhomogeneous SRT solution is 1.4% at most, but with the homogeneous SRT solution, it can be up to 7.4%. In the spectral band of 5–8 μm, the maximum relative error of downward emissivity can reach −11% for the homogeneous IRT solution but only −2% for the inhomogeneous IRT solution.

1. Introduction

Solving the radiative transfer equation (RTE) is a key issue in radiation scheme for climate model and remote sensing. In most numerical radiative transfer algorithms, the atmosphere is divided into many homogeneous layers. The inherent optical properties (IOPs) are then fixed within each layer and the variations of IOPs inside each layer are ignored, effectively regarding each layer as internally homogeneous. The standard solar/infrared radiative transfer (SRT/IRT) solutions are based on this assumption of internal homogeneity [1, 2, 3, 4], which cannot resolve within-layer vertical inhomogeneity.

It has been well established by observation that cumulus and stratocumulus clouds (hereinafter, collectively referred to as cumulus clouds) are inhomogeneous, both horizontally and vertically [5, 6, 7, 8, 9]. Inside a cumulus cloud, the liquid water content (LWC) and the cloud droplet size distribution vary with height, and so the IOPs of cloud droplets depend on vertical height.

How to deal with vertical internal inhomogeneity in SRT/IRT models is an interesting topic for researchers. Li developed a Monte Carlo cloud model that can be used to investigate photon transport in inhomogeneous clouds by considering an internal variation of the optical properties [10]. Their model showed that when overcast clouds become broken clouds, the difference in reflectance at large solar zenith angles between vertically inhomogeneous clouds and their plane-parallel counterparts can be as much as 10%.

However, the Monte Carlo method is very expensive in computing and not applicable to climate models or remote sensing [11]. The albedo of inhomogeneous mixed-phase clouds at visible wavelengths could be obtained by using a Monte Carlo method to compare such clouds with plane-parallel homogeneous clouds [12].

In principle, the vertical inhomogeneity problem of the SRT/IRT process can be solved by increasing the number of layers of the climate model. However, it is time-consuming to increase the vertical resolution of a climate model. Typically, there are only 30–100 layers in a climate model [13], which is not high enough to resolve the cloud vertical inhomogeneity. To completely address the problem of vertical inhomogeneity by using a limited number of layers in a climate model, the standard SRT method must be extended to deal with the vertical inhomogeneity inside each model layer. The primary purpose of this study is to introduce a new inhomogeneous SRT/IRT solution presented by Zhang and Shi. This solution follows a perturbation method: the zeroth-order solution is the standard Eddington approximation for SRT and two-stream approximation for IRT, with a first-order perturbation to account for the inhomogeneity effect. In Section 2, the basic theory of SRT/IRT is introduced, and the new inhomogeneous SRT/IRT solution is presented. In Section 3, the inhomogeneous SRT/IRT solution is applied to cloud as realistic examples to demonstrate the practicality of this new method. A summary is given in Section 4.

2. SRT/IRT solution for an inhomogeneous layer

2.1. SRT solution

The azimuthally averaged solar radiative transfer equation [1, 2, 3, 4, 10, 11, 12] is

$$\mu \frac{{\mathit{dI}}_S\left(\tau ,\mu \right)}{\mathit{d\tau }}=I_S\left(\tau ,\mu \right)-\frac{\omega \left(\tau \right)}{2}{\int }_{-1}^1{I}_S\left(\tau ,\mu \right)P\left(\tau ,\mu ,{\mu }^{\prime }\right)d{\mu }^{\prime }-\frac{\omega \left(\tau \right)}{4\pi }{F}_{0}P\left(\tau ,\mu ,-{\mu }_0\right)e^{-\frac{\tau }{{\mu }_0}}$$ E1

where μ is the cosine of the zenith angle (μ > 0 and μ < 0 refer to upward and downward radiation, respectively), $P\left(\tau ,\mu ,{\mu }^{\prime }\right)$ is the scattering phase function, $\tau$ is the optical depth ( $\tau =0$ and $\tau ={\tau }_0$ refer to the top and bottom of the medium, respectively), $\omega \left(\tau \right)$ is the single-scattering albedo, and $F_0$ is the incoming solar flux. For the Eddington approximation, $P\left(\tau ,\mu ,{\mu }^{\prime }\right)=1+3g\left(\tau \right)\mu {\mu }^{\prime }$ (−1 < μ <1) and $g\left(\tau \right)$ are the asymmetry factors. For the scattering atmosphere, the irradiance fluxes in the upward and downward directions can be written as

$$F_S^{\pm }\left(\tau \right)=2\pi {\int }_0^{\pm 1}{I}_S\left(\tau ,\mu \right)\mathit{\mu d\mu }$$ E2

To simulate a realistic medium such as cloud or snow, we consider $\omega \left(\tau \right)$ and $g\left(\tau \right)$ to vary with $\tau$ , and we use exponential expressions here to simplify the process. The single-scattering albedo and asymmetry factor are written as

$$\omega \left(\tau \right)=\widehat{\omega }+{\epsilon }_{\omega }\left(e^{-a_1\tau }-e^{-a_1{\tau }_0/2}\right)$$ E3a

$$g\left(\tau \right)=\widehat{g}+{\epsilon }_g\left(e^{-a_2\tau }-e^{-a_2{\tau }_0/2}\right)$$ E3b

where ${\tau }_0$ is the optical depth of the layer, $\widehat{\omega }$ is the single-scattering albedo at ${\tau }_0/2$ , and $\widehat{g}$ is the asymmetry factor at the same place. Both ${\epsilon }_g$ and ${\epsilon }_{\omega }$ are small parameters that are far less than $\widehat{g}$ and $\widehat{\omega }$ , respectively, in a realistic medium.

According to the Eddington approximation, the radiative intensity $I_S\left(\tau ,\mu \right)$ can be written as

$$I_S\left(\tau ,\mu \right)=I_{S0}\left(\tau \right)+I_{S1}\left(\tau \right)\mu$$ E4

Using Eqs. (1), (2), and (4), we obtain

$$\frac{{\mathit{dF}}_S^{+}\left(\tau \right)}{\mathit{d\tau }}={\gamma }_1\left(\tau \right)F_S^{+}\left(\tau \right)-{\gamma }_2\left(\tau \right)F_S^{-}\left(\tau \right)-{\gamma }_3\left(\tau \right)\omega \left(\tau \right)F_0{e}^{-\frac{\tau }{{\mu }_0}}$$ E5a

$$\frac{{\mathit{dF}}_S^{-}\left(\tau \right)}{\mathit{d\tau }}={\gamma }_2\left(\tau \right)F_S^{+}\left(\tau \right)-{\gamma }_1\left(\tau \right)F_S^{-}\left(\tau \right)+\left[1-{\gamma }_3\left(\tau \right)\right]\omega \left(\tau \right)F_0{e}^{-\frac{\tau }{{\mu }_0}}$$ E5b

$$F_S^{-}\left(0\right)=0,F_S^{+}\left({\tau }_0\right)=R_{\mathit{dif}}{F}_S^{-}\left({\tau }_0\right)+R_{\mathit{dir}}{\mu }_0{F}_0{e}^{-\frac{{\tau }_0}{{\mu }_0}}$$ E5c

where ${\gamma }_1\left(\tau \right)=\frac{1}{4}\left\{7-\left[4+3g\left(\tau \right)\right]\omega \left(\tau \right)\right\}$ , ${\gamma }_2\left(\tau \right)=\frac{-1}{4}\left\{1-\left[4-3g\left(\tau \right)\right]\omega \left(\tau \right)\right\}$ , and ${\gamma }_3\left(\tau \right)=\frac{1}{4}\left[2-3g\left(\tau \right){\mu }_0\right]$ ; ${\tau }_0$ is the optical depth of the single layer; and $R_{\mathit{dif}}$ ( $R_{\mathit{dir}}$ ) is the diffuse (resp., direct) reflection from the layer below or the diffuse (direct) surface albedo. Substituting ${\gamma }_1\left(\tau \right)$ , ${\gamma }_2\left(\tau \right)$ , and ${\gamma }_3\left(\tau \right)$ into Eq. (3) and ignoring the small second-order parameters ${{\epsilon }_{\omega }}^2$ , ${{\epsilon }_g}^2$ , and ${\epsilon }_{\omega }{\epsilon }_g$ , we get

$${\gamma }_1\left(\tau \right)={\gamma }_1^0+{\gamma }_1^1{\epsilon }_{\omega }\left(e^{-a_1\tau }-e^{-a_1{\tau }_0/2}\right)+{\gamma }_1^2{\epsilon }_g\left(e^{-a_2\tau }-e^{-a_2{\tau }_0/2}\right)$$ E6a

$${\gamma }_2\left(\tau \right)={\gamma }_2^0+{\gamma }_2^1{\epsilon }_{\omega }\left(e^{-a_1\tau }-e^{-a_1{\tau }_0/2}\right)+{\gamma }_2^2{\epsilon }_g\left(e^{-a_2\tau }-e^{-a_2{\tau }_0/2}\right)$$ E6b

$${\gamma }_3\left(\tau \right)={\gamma }_3^0+{\gamma }_3^2{\epsilon }_g\left(e^{-a_2\tau }-e^{-a_2{\tau }_0/2}\right)$$ E6c

where ${\gamma }_1^0=\frac{1}{4}\left[7-\left(4+3\widehat{g}\right)\widehat{\omega }\right]$ , ${\gamma }_2^0=\frac{-1}{4}\left[1-\left(4-3\widehat{g}\right)\widehat{\omega }\right]$ , ${\gamma }_3^0=\frac{1}{4}\left(2-3\widehat{g}{\mu }_0\right)$ , ${\gamma }_1^1=\frac{-1}{4}\left(4+3\widehat{g}\right)$ , ${\gamma }_1^2=-\frac{3}{4}\widehat{\omega }$ , ${\gamma }_2^1=\frac{1}{4}\left(4-3\widehat{g}\right)$ , ${\gamma }_2^2=-\frac{3}{4}\widehat{\omega }$ , and ${\gamma }_3^2=-\frac{3}{4}{\mu }_0$ .

By perturbation theory [14], the corresponding flux can also be expanded by using the perturbation coefficients ${\epsilon }_{\omega }$ and ${\epsilon }_g$ :

$$F_S^{+}=F_{S0}^{+}+{\epsilon }_{\omega }{F}_{S1}^{+}+{\epsilon }_g{F}_{S2}^{+}$$ E7a

$$F_S^{-}=F_{S0}^{-}+{\epsilon }_{\omega }{F}_{S1}^{-}+{\epsilon }_g{F}_{S2}^{-}$$ E7b

Substituting Eqs. (6) and (7) into Eq. (5) yields

$$\begin{array}{c}\frac{{\mathit{dF}}_S^{+}}{\mathit{d\tau }}=\left[{\gamma }_1^0+{\gamma }_1^1{\epsilon }_{\omega }\left(e^{-a_1\tau }-e^{-a_1{\tau }_0/2}\right)+{\gamma }_1^2{\epsilon }_g\left(e^{-a_2\tau }-e^{-a_2{\tau }_0/2}\right)\right]\left(F_{S0}^{+}+{\epsilon }_{\omega }{F}_{S1}^{+}+{\epsilon }_g{F}_{S2}^{+}\right)\\ -\left[{\gamma }_2^0+{\gamma }_2^1{\epsilon }_{\omega }\left(e^{-a_1\tau }-e^{-a_1{\tau }_0/2}\right)+{\gamma }_2^2{\epsilon }_g\left(e^{-a_2\tau }-e^{-a_2{\tau }_0/2}\right)\right]\left(F_{S0}^{-}+{\epsilon }_{\omega }{F}_{S1}^{-}+{\epsilon }_g{F}_{S2}^{-}\right)\\ -\left[\widehat{\omega }{\gamma }_3^0+{\gamma }_3^0{\epsilon }_{\omega }\left(e^{-a_1\tau }-e^{-a_1{\tau }_0/2}\right)+\widehat{\omega }{\gamma }_3^2{\epsilon }_g\left(e^{-a_2\tau }-e^{-a_2{\tau }_0/2}\right)\right]F_0{e}^{-\frac{\tau }{{\mu }_0}}\end{array}$$ E8a

$$\begin{array}{c}\frac{{\mathit{dF}}_S^{-}}{\mathit{d\tau }}=\left[{\gamma }_2^0+{\gamma }_2^1{\epsilon }_{\omega }\left(e^{-a_1\tau }-e^{-a_1{\tau }_0/2}\right)+{\gamma }_2^2{\epsilon }_g\left(e^{-a_2\tau }-e^{-a_2{\tau }_0/2}\right)\right]\left(F_{S0}^{+}+{\epsilon }_{\omega }{F}_{S1}^{+}+{\epsilon }_g{F}_{S2}^{+}\right)\\ -\left[{\gamma }_1^0+{\gamma }_1^1{\epsilon }_{\omega }\left(e^{-a_1\tau }-e^{-a_1{\tau }_0/2}\right)+{\gamma }_1^2{\epsilon }_g\left(e^{-a_2\tau }-e^{-a_2{\tau }_0/2}\right)\right]\left(F_{S0}^{-}+{\epsilon }_{\omega }{F}_{S1}^{-}+{\epsilon }_g{F}_{S2}^{-}\right)\\ +\left[\widehat{\omega }{\gamma }_4^0+{\gamma }_4^0{\epsilon }_{\omega }\left(e^{-a_1\tau }-e^{-a_1{\tau }_0/2}\right)-\widehat{\omega }{\gamma }_3^2{\epsilon }_g\left(e^{-a_2\tau }-e^{-a_2{\tau }_0/2}\right)\right]F_0{e}^{-\frac{\tau }{{\mu }_0}}\end{array}$$ E8b

where ${\gamma }_4^0=1-{\gamma }_3^0$ . And, Eq. (8) can be rewritten as separate equations for $F_{S0}^{\pm }$ , $F_{S1}^{\pm }$ , and $F_{S2}^{\pm }$ . We obtain the following equations for the scattered flux $F_{S0}^{\pm }$ :

$$\frac{{\mathit{dF}}_{S0}^{+}}{\mathit{d\tau }}={\gamma }_1^0{F}_{S0}^{+}-{\gamma }_2^0{F}_{S0}^{-}-{\gamma }_3^0\widehat{\omega }{F}_0{e}^{-\frac{\tau }{{\mu }_0}}$$ E9a

$$\frac{{\mathit{dF}}_{S0}^{-}}{\mathit{d\tau }}={\gamma }_2^0{F}_{S0}^{+}-{\gamma }_1^0{F}_{S0}^{-}+{\gamma }_4^0\widehat{\omega }{F}_0{e}^{-\frac{\tau }{{\mu }_0}}$$ E9b

$$F_{S0}^{-}\left(0\right)=0,F_{S0}^{+}\left({\tau }_0\right)=R_{\mathit{dif}}{F}_{S0}^{-}\left({\tau }_0\right)+R_{\mathit{dir}}{\mu }_0{F}_0{e}^{-\frac{{\tau }_0}{{\mu }_0}}$$ E9c

Eq. (9) is the standard SRT equation for a homogeneous layer [15] and has the following solution:

$$F_{S0}^{+}=K_1{e}^{\mathit{k\tau }}+\mathrm{\Gamma }{K}_2{e}^{-\mathit{k\tau }}+G_1{e}^{-\frac{\tau }{{\mu }_0}}$$ E10a

$$F_{S0}^{-}=\mathrm{\Gamma }{K}_1{e}^{\mathit{k\tau }}+K_2{e}^{-\mathit{k\tau }}+G_2{e}^{-\frac{\tau }{{\mu }_0}}$$ E10b

where $K_1=\frac{\left(\mathrm{\Gamma }-R_{\mathit{dif}}\right)\mathrm{\Gamma }{G}_2{e}^{-k{\tau }_0}-\left(G_1-R_{\mathit{dir}}{G}_2-R_{\mathit{dir}}{\mu }_0{F}_0\right)e^{-\frac{{\tau }_0}{{\mu }_0}}}{\left(1-R_{\mathit{dif}}\mathrm{\Gamma }\right)e^{k{\tau }_0}-\left(\mathrm{\Gamma }-R_{\mathit{dif}}\right)\mathrm{\Gamma }{e}^{-k{\tau }_0}}$ , $K_2=-\mathrm{\Gamma }{K}_1-G_2$ , $G_1=\left[{\gamma }_3^0\left(\frac{1}{{\mu }_0}-{\gamma }_1^0\right)-{\gamma }_2^0{\gamma }_4^0\right]\frac{{{\mu }_0}^2\widehat{\omega }{F}_0}{1-{{\mu }_0}^2{k}^2}$ , $G_2=-\left[{\gamma }_4^0\left(\frac{1}{{\mu }_0}+{\gamma }_1^0\right)+{\gamma }_2^0{\gamma }_3^0\right]\frac{{{\mu }_0}^2\widehat{\omega }{F}_0}{1-{{\mu }_0}^2{k}^2}$ , $\mathrm{\Gamma }=1-\frac{2k}{{\gamma }_1^0+{\gamma }_2^0+k}$ , and $k^2=\left({\gamma }_1^0+{\gamma }_2^0\right)\left({\gamma }_1^0-{\gamma }_2^0\right)$ . And, the equations for the perturbation terms $F_{\mathit{Si}}^{\pm }$ (i = 1, 2) are

$$\frac{{\mathit{dF}}_{\mathit{Si}}^{+}}{\mathit{d\tau }}={\gamma }_1^0{F}_{\mathit{Si}}^{+}-{\gamma }_2^0{F}_{\mathit{Si}}^{-}+\left(e^{-a_i\tau }-e^{-a_i{\tau }_0/2}\right)\left({\gamma }_1^i{F}_{S0}^{+}-{\gamma }_2^i{F}_{S0}^{-}\right)-{\gamma }_3^{i-1}{F}_0\left(e^{-a_i\tau }-e^{-a_i{\tau }_0/2}\right)e^{-\frac{\tau }{{\mu }_0}}$$ E11a

$$\frac{{\mathit{dF}}_{\mathit{Si}}^{-}}{\mathit{d\tau }}={\gamma }_2^0{F}_{\mathit{Si}}^{+}-{\gamma }_1^0{F}_{\mathit{Si}}^{-}+\left(e^{-a_i\tau }-e^{-a_i{\tau }_0/2}\right)\left({\gamma }_2^i{F}_{S0}^{+}-{\gamma }_1^i{F}_{S0}^{-}\right)-{\gamma }_4^{i-1}{F}_0\left(e^{-a_i\tau }-e^{-a_i{\tau }_0/2}\right)e^{-\frac{\tau }{{\mu }_0}}$$ E11b

$$F_{\mathit{Si}}^{-}\left(0\right)=0,F_{\mathit{Si}}^{+}\left({\tau }_0\right)=R_{\mathit{dif}}{F}_{\mathit{Si}}^{-}\left({\tau }_0\right)$$ E11c

where ${\gamma }_3^1=-{\gamma }_4^1=\widehat{\omega }{\gamma }_3^2$ . Letting $M_i=F_{\mathit{Si}}^{+}+F_{\mathit{Si}}^{-}$ and $N_i=F_{\mathit{Si}}^{+}-F_{\mathit{Si}}^{-}$ , Eq. (11a) and (11b) yields

$$\begin{array}{c}\frac{{\mathit{dM}}_i}{\mathit{d\tau }}=\left({\gamma }_1^0+{\gamma }_2^0\right)N_i+\left({\psi }_i^{+}+{\psi }_i^{-}\right)e^{-\left(k+a_i\right)\tau }+\left({\zeta }_i^{+}+{\zeta }_i^{-}\right)e^{\left(k-a_i\right)\tau }+\left({\chi }_i^{+}+{\chi }_i^{-}\right)e^{-\left(a_i+\frac{1}{{\mu }_0}\right)\tau }\\ -e^{-a_i{\tau }_0/2}\left[\left({\psi }_i^{+}+{\psi }_i^{-}\right)e^{-\mathit{k\tau }}+\left({\zeta }_i^{+}+{\zeta }_i^{-}\right)e^{\mathit{k\tau }}+\left({\chi }_i^{+}+{\chi }_i^{-}\right)e^{-\frac{\tau }{{\mu }_0}}\right]\end{array}$$ E12a

$$\begin{array}{c}\frac{{\mathit{dN}}_i}{\mathit{d\tau }}=\left({\gamma }_1^0-{\gamma }_2^0\right)M_i+\left({\psi }_i^{+}-{\psi }_i^{-}\right)e^{-\left(k+a_i\right)\tau }+\left({\zeta }_i^{+}-{\zeta }_i^{-}\right)e^{\left(k-a_i\right)\tau }+\left({\chi }_i^{+}-{\chi }_i^{-}\right)e^{-\left(a_i+\frac{1}{{\mu }_0}\right)\tau }\\ -e^{-a_i{\tau }_0/2}\left[\left({\psi }_i^{+}-{\psi }_i^{-}\right)e^{-\mathit{k\tau }}+\left({\zeta }_i^{+}-{\zeta }_i^{-}\right)e^{\mathit{k\tau }}+\left({\chi }_i^{+}-{\chi }_i^{-}\right)e^{-\frac{\tau }{{\mu }_0}}\right]\end{array}$$ E12b

where ${\mathrm{\Psi }}_i^{+}=K_2\left({\gamma }_1^i\mathrm{\Gamma }-{\gamma }_2^i\right)$ , ${\mathrm{\Psi }}_i^{-}=K_2\left({\gamma }_2^i\mathrm{\Gamma }-{\gamma }_1^i\right)$ , ${\zeta }_i^{+}=K_1\left({\gamma }_1^i-{\gamma }_2^i\mathrm{\Gamma }\right)$ , ${\zeta }_i^{-}=K_1\left({\gamma }_2^i-{\gamma }_1^i\mathrm{\Gamma }\right)$ , ${\chi }_i^{+}={\gamma }_1^i{G}_1-{\gamma }_2^i{G}_2-{\gamma }_3^{i-1}{F}_0$ , and ${\chi }_i^{-}={\gamma }_2^i{G}_1-{\gamma }_1^i{G}_2+{\gamma }_4^{i-1}{F}_0$ .

From Eq. (12), we obtain

$$\frac{d^2{M}_i}{d{\tau }^2}=k^2{M}_i+{\eta }_{1i}^{+}{e}^{-\left(k+a_i\right)\tau }+{\eta }_{2i}^{+}{e}^{\left(k-a_i\right)\tau }+{\eta }_{3i}^{+}{e}^{-\left(a_i+\frac{1}{{\mu }_0}\right)\tau }+{\eta }_{4i}^{+}{e}^{-\mathit{k\tau }}+{\eta }_{5i}^{+}{e}^{\mathit{k\tau }}+{\eta }_{6i}^{+}{e}^{-\frac{\tau }{{\mu }_0}}$$ E13a

$$\frac{d^2{N}_i}{d{\tau }^2}=k^2{N}_i+{\eta }_{1i}^{-}{e}^{-\left(k+a_i\right)\tau }+{\eta }_{2i}^{-}{e}^{\left(k-a_i\right)\tau }+{\eta }_{3i}^{-}{e}^{-\left(a_i+\frac{1}{{\mu }_0}\right)\tau }+{\eta }_{4i}^{-}{e}^{-\mathit{k\tau }}+{\eta }_{5i}^{-}{e}^{\mathit{k\tau }}+{\eta }_{6i}^{-}{e}^{-\frac{\tau }{{\mu }_0}}$$ E13b

where ${\eta }_{1i}^{\pm }=\left({\gamma }_1^0\pm {\gamma }_2^0\right)\left({\psi }_i^{+}\mp {\psi }_i^{-}\right)-\left(k+a_i\right)\left({\psi }_i^{+}\pm {\psi }_i^{-}\right)$ , ${\eta }_{2i}^{\pm }=\left(k-a_i\right)\left({\zeta }_i^{+}\pm {\zeta }_i^{-}\right)+\left({\gamma }_1^0\pm {\gamma }_2^0\right)\left({\zeta }_i^{+}\mp {\zeta }_i^{-}\right)$ , ${\eta }_{3i}^{\pm }=\left({\chi }_i^{+}\mp {\chi }_i^{-}\right)\left({\gamma }_1^0\pm {\gamma }_2^0\right)-\left(a_i+\frac{1}{{\mu }_0}\right)\left({\chi }_i^{+}\mp {\chi }_i^{-}\right)$ , ${\eta }_{4i}^{\pm }=-e^{-a_i{\tau }_0/2}\left[\left({\gamma }_1^0\pm {\gamma }_2^0\right)\left({\psi }_i^{+}\mp {\psi }_i^{-}\right)-k\left({\psi }_i^{+}\pm {\psi }_i^{-}\right)\right]$ , ${\eta }_{5i}^{\pm }=-e^{-a_i{\tau }_0/2}\left[k\left({\zeta }_i^{+}\pm {\zeta }_i^{-}\right)+\left({\gamma }_1^0\pm {\gamma }_2^0\right)\left({\zeta }_i^{+}\mp {\zeta }_i^{-}\right)\right]$ , and ${\eta }_{6i}^{\pm }=-e^{-a_i{\tau }_0/2}\left[\left({\chi }_i^{+}\mp {\chi }_i^{-}\right)\left({\gamma }_1^0\pm {\gamma }_2^0\right)-\frac{1}{{\mu }_0}\left({\chi }_i^{+}\mp {\chi }_i^{-}\right)\right]$ .

The solutions of Eq. (13) are

$$M_i=A_i^{+}{e}^{-\mathit{k\tau }}+B_i^{+}{e}^{\mathit{k\tau }}+P_i^{+}{e}^{-\left(k+a_i\right)\tau }+Q_i^{+}{e}^{\left(k-a_i\right)\tau }+R_i^{+}{e}^{-\left(a_i+\frac{1}{{\mu }_0}\right)\tau }-\frac{{\eta }_{4i}^{+}}{2k}{e}^{-\mathit{k\tau }}+\frac{{\eta }_{5i}^{+}}{2k}{e}^{\mathit{k\tau }}+\frac{{\eta }_{6i}^{+}{\mu }_0^2}{1-{\mu }_0^2{k}^2}{e}^{-\frac{\tau }{{\mu }_0}}$$ E14a

$$N_i=A_i^{-}{e}^{-\mathit{k\tau }}+B_i^{-}{e}^{\mathit{k\tau }}+P_i^{-}{e}^{-\left(k+a_i\right)\tau }+Q_i^{-}{e}^{\left(k-a_i\right)\tau }+R_i^{-}{e}^{-\left(a_i+\frac{1}{{\mu }_0}\right)\tau }-\frac{{\eta }_{4i}^{-}}{2k}{e}^{-\mathit{k\tau }}+\frac{{\eta }_{5i}^{-}}{2k}{e}^{\mathit{k\tau }}+\frac{{\eta }_{6i}^{-}{\mu }_0^2}{1-{\mu }_0^2{k}^2}{e}^{-\frac{\tau }{{\mu }_0}}$$ E14b

where $P_i^{\pm }=\frac{{\eta }_{1i}^{\pm }}{{\left(k+a_i\right)}^2-k^2}$ , $Q_i^{\pm }=\frac{{\eta }_{2i}^{\pm }}{{\left(k-a_i\right)}^2-k^2}$ , and $R_i^{\pm }=\frac{{\eta }_{3i}^{\pm }}{{\left(a_i+\frac{1}{{\mu }_0}\right)}^2-k^2}{e}^{-\frac{\tau }{{\mu }_0}}$ . Finally, we can obtain $F_{\mathit{Si}}^{-}$ and $F_{\mathit{Si}}^{+}$ as

$$F_{\mathit{Si}}^{+}=D_{1i}^{+}{e}^{-\mathit{k\tau }}+D_{2i}^{+}{e}^{\mathit{k\tau }}+{\phi }_{1i}^{+}{e}^{-\left(k+a_i\right)\tau }+{\phi }_{2i}^{+}{e}^{\left(k-a_i\right)\tau }+{\phi }_{3i}^{+}{e}^{-\left(a_i+\frac{1}{{\mu }_0}\right)\tau }+{\phi }_{4i}^{+}\tau e^{-k_1\tau }+{\phi }_{5i}^{+}\tau e^{\mathit{k\tau }}+{\phi }_{6i}^{+}{e}^{-\frac{\tau }{{\mu }_0}}$$ E15a

$$F_{\mathit{Si}}^{-}=D_{1i}^{-}{e}^{-\mathit{k\tau }}+D_{2i}^{-}{e}^{\mathit{k\tau }}+{\phi }_{1i}^{-}{e}^{-\left(k+a_i\right)\tau }+{\phi }_{2i}^{-}{e}^{\left(k-a_i\right)\tau }+{\phi }_{3i}^{-}{e}^{-\left(a_i+\frac{1}{{\mu }_0}\right)\tau }+{\phi }_{4i}^{-}\tau e^{-k_1\tau }+{\phi }_{5i}^{-}\tau e^{\mathit{k\tau }}+{\phi }_{6i}^{-}{e}^{-\frac{\tau }{{\mu }_0}}$$ E15b

where $D_{1i}^{\pm }=A_i^{+}{\alpha }^{\mp }\pm X_i$ , $D_{2i}^{\pm }=B_i^{+}{\alpha }^{\pm }\pm Y$ , ${\alpha }^{\pm }=\frac{1}{2}\left(1\pm \frac{k}{{\gamma }_1^0+{\gamma }_2^0}\right)$ , $X_i=\frac{e^{-a_i{\tau }_0/2\left({\psi }_i^{+}+{\psi }_i^{-}\right)}}{2\left({\gamma }_1^0+{\gamma }_2^0\right)}-\frac{{\eta }_{4i}^{+}}{4k\left({\gamma }_1^0+{\gamma }_2^0\right)}$ ,

$Y_i=\frac{e^{-a_i{\tau }_0/2\left({\zeta }_i^{+}+{\zeta }_i^{-}\right)}}{2\left({\gamma }_1^0+{\gamma }_2^0\right)}+\frac{{\eta }_{5i}^{+}}{4k\left({\gamma }_1^0+{\gamma }_2^0\right)}$ , ${\varphi }_{1i}^{\pm }=\frac{1}{2}\left(P_i^{+}\pm P_i^{-}\right)$ , ${\varphi }_{2i}^{\pm }=\frac{1}{2}\left(Q_i^{+}\pm Q_i^{-}\right)$ , ${\varphi }_{3i}^{\pm }=\frac{1}{2}\left(R_i^{+}\pm R_i^{-}\right)$ , ${\varphi }_{4i}^{\pm }=-\frac{{\eta }_{4i}^{+}\pm {\eta }_{4i}^{-}}{4k}$ , ${\varphi }_{5i}^{\pm }=\frac{{\eta }_{5i}^{+}\pm {\eta }_{5i}^{-}}{4k}$ , and ${\varphi }_{6i}^{\pm }=\frac{\left({\eta }_{6i}^{+}\pm {\eta }_{6i}^{-}\right){\mu }_0^2}{1-{\mu }_0^2{k}^2}$ . $B_i$ and $A_i$ are determined by the boundary conditions as

$$\begin{array}{l}{B}_i^{+}\\ =-\frac{{\varphi }_{1i}^{-}+{\varphi }_{2i}^{-}+{\varphi }_{3i}^{-}+{\varphi }_{6i}^{-})\left({\alpha }^{-}-R_{\mathit{dif}}{\alpha }^{+}\right)e^{-k{\tau }_0}+{\alpha }^{+}\left[\left({\varphi }_{1i}^{+}-R_{\mathit{dif}}{\varphi }_{1i}^{-}\right)e^{-\left(k+a_i\right){\tau }_0}+\left({\varphi }_{2i}^{+}-R_{\mathit{dif}}{\varphi }_{2i}^{-}\right)e^{\left(k-a_i\right){\tau }_0}\right]}{{\alpha }^{-}\left({\alpha }^{-}-R_{\mathit{dif}}{\alpha }^{+}\right)e^{-k{\tau }_0}-{\alpha }^{+}\left({\alpha }^{+}-R_{\mathit{dif}}{\alpha }^{-}\right)e^{k{\tau }_0}}\hfill \\ +\frac{{\alpha }^{+}\left[{\varphi }_{3i}^{+}-R_{\mathit{dif}}{\varphi }_{3i}^{-}\right)e^{-\left(a_i+\frac{1}{{\mu }_0}\right){\tau }_0}+{\varphi }_{4i}^{+}-R_{\mathit{dif}}{\varphi }_{4i}^{-}){\tau }_0{e}^{-k{\tau }_0}+\left({\varphi }_{5i}^{+}-R_{\mathit{dif}}{\varphi }_{5i}^{-}\right){\tau }_0{e}^{k{\tau }_0}+\left({\varphi }_{6i}^{+}-R_{\mathit{dif}}{\varphi }_{6i}^{-}\right)e^{-\frac{{\tau }_0}{{\mu }_0}}]}{{\alpha }^{-}\left({\alpha }^{-}-R_{\mathit{dif}}{\alpha }^{+}\right)e^{-k{\tau }_0}-{\alpha }^{+}\left({\alpha }^{+}-R_{\mathit{dif}}{\alpha }^{-}\right)e^{k{\tau }_0}}\hfill \\ +\frac{\left(X_i+Y_i\right)\left({\alpha }^{-}-R_{\mathit{dif}}{\alpha }^{+}\right)e^{-k{\tau }_0}+\left(X_i+R_{\mathit{dif}}{X}_i\right){\alpha }^{+}{e}^{-k{\tau }_0}+\left(Y_i+R_{\mathit{dif}}{Y}_i\right){\alpha }^{+}{e}^{k{\tau }_0}}{{\alpha }^{-}\left({\alpha }^{-}-R_{\mathit{dif}}{\alpha }^{+}\right)e^{-k{\tau }_0}-{\alpha }^{+}\left({\alpha }^{+}-R_{\mathit{dif}}{\alpha }^{-}\right)e^{k{\tau }_0}}\hfill \end{array}$$ E16a

$$A_i^{+}=\frac{1}{{\alpha }^{+}}\left(X_i+Y_i-B_i{\alpha }^{-}-{\varphi }_{1i}^{-}-{\varphi }_{2i}^{-}-{\varphi }_{3i}^{-}-{\varphi }_{6i}^{-}\right)$$ E16b

All detailed calculation about solar radiation can be found at [16].

2.2. IRT solution

The azimuthally averaged infrared radiative transfer equation for intensity $I_{\mathrm{I}}\left(\tau ,\mu \right)$ is [1, 2, 3, 4, 10, 11, 12]

$$\mu \frac{{\mathit{dI}}_{\mathrm{I}}\left(\tau ,\mu \right)}{\mathit{d\tau }}=I_I\left(\tau ,\mu \right)-\frac{\omega \left(\tau \right)}{2}{\int }_{-1}^1{I}_{\mathrm{I}}\left(\tau ,\mu \right)P\left(\tau ,\mu ,{\mu }^{\text{'}}\right)d{\mu }^{\text{'}}-\left[1-\omega \left(\tau \right)\right]B\left(T\right)$$ E17

where $\mu$ , $\tau$ , $P\left(\tau ,\mu ,{\mu }^{\text{'}}\right)$ , and $\omega \left(\tau \right)$ are same as in Eq. (1). $B\left(T\right)$ is the Planck function at temperature $T$ , which represents the internal infrared emission of the medium.

The Planck function is approximated lineally as a function of optical depth [2] as

$$B\left[T\left(\tau \right)\right]=B_0+\mathit{\beta \tau }$$ E18

where $\beta =\left(B_1-B_0\right)/{\tau }_0$ and ${\tau }_0$ are the total optical depth of the medium. The Planck functions $B_0$ and $B_1$ are evaluated by using the temperature of the top ( $\tau =0$ ) and the bottom ( $\tau ={\tau }_0$ ) of the medium.

According to the two-stream approximation, the intensities can be written as $I_{\mathrm{I}}\left(\tau ,{\mu }_1\right)=I_{\mathrm{I}}^{+}\left(\tau \right)$ and $I_{\mathrm{I}}\left(\tau ,{\mu }_{-1}\right)=I_{\mathrm{I}}^{-}\left(\tau \right)$ , respectively, where ${\mu }_1=-{\mu }_{-1}=1/1.66$ is a diffuse factor that converts radiative intensity to flux [17]. ${\int }_{-1}^1{I}_I\left(\tau ,\mu \right)P\left(\tau ,\mu ,{\mu }^{\prime }\right)d{\mu }^{\prime }$ can be written as

$${\int }_{-1}^1{I}_I\left(\tau ,\mu \right)P\left(\tau ,\mu ,{\mu }^{\prime }\right)d{\mu }^{\prime }=\left[1+3g\left(\tau \right){\mathit{\mu \mu }}_1\right]I_I^{+}\left(\tau \right)+\left[1+3g\left(\tau \right){\mathit{\mu \mu }}_{-1}\right]I_I^{-}\left(\tau \right)$$ E19

where $g\left(\tau \right)$ is the asymmetry factor.

Using Eqs. (17) and (19), we can obtain

$$\frac{{\mathit{dI}}_{\mathrm{I}}^{+}\left(\tau \right)}{\mathit{d\tau }}={\gamma }_1\left(\tau \right)I_{\mathrm{I}}^{+}\left(\tau \right)-{\gamma }_2\left(\tau \right)I_{\mathrm{I}}^{-}\left(\tau \right)-{\gamma }_3\left(\tau \right)B\left(\tau \right)$$ E20a

$$\frac{{\mathit{dI}}_{\mathrm{I}}^{+}\left(\tau \right)}{\mathit{d\tau }}={\gamma }_2\left(\tau \right)I_{\mathrm{I}}^{+}\left(\tau \right)-{\gamma }_1\left(\tau \right)I_{\mathrm{I}}^{-}\left(\tau \right)+{\gamma }_3\left(\tau \right)B\left(\tau \right)$$ E20b

where ${\gamma }_1\left(\tau \right)=\frac{1‐\omega \left(\tau \right)\left(1+g\left(\tau \right)\right)/2}{{\mu }_1}$ , ${\gamma }_2\left(\tau \right)=\frac{\omega \left(\tau \right)\left[1-g\left(\tau \right)\right]}{2{\mu }_1}$ , and ${\gamma }_3\left(\tau \right)=\frac{1-\omega \left(\tau \right)}{{\mu }_1}$ .

For IRT, we also use Eq. (3) to represent an inhomogeneous medium such as cloud or snow, in which $\omega \left(\tau \right)$ and $g\left(\tau \right)$ vary with $\tau$ . By substituting Eq. (3) into ${\gamma }_1\left(\tau \right)$ , ${\gamma }_2\left(\tau \right)$ , and ${\gamma }_3\left(\tau \right)$ and by ignoring the second order of the small parameters of ${\epsilon }_{\omega }^2$ , ${\epsilon }_g^2$ , and ${\epsilon }_{\omega }{\epsilon }_g$ , we can obtain

$${\gamma }_1\left(\tau \right)={\gamma }_1^0+{\gamma }_1^1{\epsilon }_{\omega }\left(e^{-a_1\tau }-e^{-a_1{\tau }_0/2}\right)+{\gamma }_1^2{\epsilon }_g\left(e^{-a_2\tau }-e^{-a_2{\tau }_0/2}\right)$$ E21a

$${\gamma }_2\left(\tau \right)={\gamma }_2^0+{\gamma }_2^1{\epsilon }_{\omega }\left(e^{-a_1\tau }-e^{-a_1{\tau }_0/2}\right)+{\gamma }_2^2{\epsilon }_g\left(e^{-a_2\tau }-e^{-a_2{\tau }_0/2}\right)$$ E21b

$${\gamma }_3\left(\tau \right)={\gamma }_3^0+{\gamma }_3^1{\epsilon }_{\omega }\left(e^{-a_1\tau }-e^{-a_1{\tau }_0/2}\right)$$ E21c

In the above formula, ${\gamma }_i^0$ , ${\gamma }_i^1$ , and ${\gamma }_i^2$ (i = 1, 2, 3) are the known factors of $\widehat{\omega }$ and $\widehat{g}$ . These known factors are introduced for simplifying original expressions, in which ${\gamma }_1^0=\frac{1-\widehat{\omega }\left(1+\widehat{g}\right)/2}{{\mu }_1}$ , ${\gamma }_2^0=\frac{\widehat{\omega }\left(1-\widehat{g}\right)}{2{\mu }_1}$ , ${\gamma }_3^0=\frac{1-\widehat{\omega }}{{\mu }_1}$ , ${\gamma }_1^1=-\frac{1+\widehat{g}}{2{\mu }_1}$ , ${\gamma }_2^1=\frac{1-\widehat{g}}{2{\mu }_1}$ , ${\gamma }_3^1=-\frac{1}{{\mu }_1}$ , ${\gamma }_1^2={\gamma }_2^2=-\frac{\widehat{\omega }}{2{\mu }_1}$ , and ${\gamma }_3^2=0$ .

Same as in Eq. (7), the upward and downward intensity can be written as

$$I_{\mathrm{I}}^{+}=I_{\mathrm{I}0}^{+}+{\epsilon }_{\omega }{I}_{\mathrm{I}1}^{+}+{\epsilon }_g{I}_{\mathrm{I}1}^{+}$$ E22a

$$I_{\mathrm{I}}^{-}=I_{\mathrm{I}0}^{-}+{\epsilon }_{\omega }{I}_{\mathrm{I}1}^{-}+{\epsilon }_g{I}_{\mathrm{I}1}^{-}$$ E22b

By substituting Eqs. (21)–(22) into Eq. (20), we obtain

$$\begin{array}{c}\frac{{\mathit{dI}}_{\mathrm{I}}^{+}}{\mathit{d\tau }}=\left[{\gamma }_1^0+{\gamma }_1^1{\epsilon }_{\omega }\left(e^{-a_1\tau }-e^{-a_1{\tau }_0/2}\right)+{\gamma }_1^2{\epsilon }_g\left(e^{-a_2\tau }-e^{-a_2{\tau }_0/2}\right)\right]\left(I_{\mathrm{I}0}^{+}+{\epsilon }_{\omega }{I}_{\mathrm{I}1}^{+}+{\epsilon }_g{I}_{\mathrm{I}2}^{+}\right)\\ -\left[{\gamma }_2^0+{\gamma }_2^1{\epsilon }_{\omega }\left(e^{-a_1\tau }-e^{-a_1{\tau }_0/2}\right)+{\gamma }_2^2{\epsilon }_g\left(e^{-a_2\tau }-e^{-a_2{\tau }_0/2}\right)\right]\left(I_{\mathrm{I}0}^{-}+{\epsilon }_{\omega }{I}_{\mathrm{I}1}^{-}+{\epsilon }_g{I}_{\mathrm{I}2}^{-}\right)\\ -\left[{\gamma }_3^0+{\gamma }_3^1{\epsilon }_{\omega }\left(e^{-a_1\tau }-e^{-a_1{\tau }_0/2}\right)\right]B\left(\tau \right)\end{array}$$ E23a

$$\begin{array}{c}\frac{{\mathit{dI}}_{\mathrm{I}}^{-}}{d_{\tau }}=\left[{\gamma }_2^0+{\gamma }_2^1{\epsilon }_{\omega }\left(e^{-a_1\tau }-e^{-a_1{\tau }_0/2}\right)+{\gamma }_2^2{\epsilon }_g\left(e^{-a_2\tau }-e^{-a_2{\tau }_0/2}\right)\right]\left(I_{\mathrm{I}0}^{+}+{\epsilon }_{\omega }{I}_{\mathrm{I}1}^{+}+{\epsilon }_g{I}_{\mathrm{I}2}^{+}\right)\\ -\left[{\gamma }_1^0+{\gamma }_1^1{\epsilon }_{\omega }\left(e^{-a_1\tau }-e^{-a_1{\tau }_0/2}\right)+{\gamma }_1^2{\epsilon }_g\left(e^{-a_2\tau }-e^{-a_2{\tau }_0/2}\right)\right]\left(I_{\mathrm{I}0}^{-}+{\epsilon }_{\omega }{I}_{\mathrm{I}1}^{-}+{\epsilon }_g{I}_{\mathrm{I}2}^{-}\right)\\ +\left[{\gamma }_3^0+{\gamma }_3^1{\epsilon }_{\omega }\left(e^{-a_1\tau }-e^{-a_1{\tau }_0/2}\right)\right]B\left(\tau \right)\end{array}$$ E23b

By removing the second-order and higher-order perturbation terms, we can also separate Eq. (23) into three equations of $I_{\mathrm{I}i}^{\pm }$ (i = 0, 1, 2). The equations of $I_{\mathrm{I}0}^{\pm }$ can be written as

$$\frac{{\mathit{dI}}_{\mathrm{I}0}^{+}}{\mathit{d\tau }}={\gamma }_1^0{\mathit{dI}}_{\mathrm{I}0}^{+}-{\gamma }_2^0{\mathit{dI}}_{\mathrm{I}0}^{-}-{\gamma }_3^{0}B\left(\tau \right)$$ E24a

$$\frac{{\mathit{dI}}_{\mathrm{I}0}^{-}}{\mathit{d\tau }}={\gamma }_2^0{I}_{\mathrm{I}0}^{+}-{\gamma }_1^0{I}_{\mathrm{I}0}^{-}+{\gamma }_3^{0}B\left(\tau \right)$$ E24b

$$I_{\mathrm{I}0}^{-}\left(0\right)=0,I_{\mathrm{I}0}^{+}\left({\tau }_0\right)=\left(1-{\epsilon }_s\right)I_{\mathrm{I}0}^{-}\left({\tau }_0\right)+{\epsilon }_{s}B\left(T_s\right)$$ E24c

where $T_s$ and ${\epsilon }_s$ are surface temperature and surface emissivity, respectively. Eq. (24) is the standard homogeneous two-stream infrared radiative transfer equation [3, 15] with solutions

$$I_{\mathrm{I}0}^{+}={\alpha }^{+}{K}_0{e}^{-k\left({\tau }_0-\tau \right)}+{\alpha }^{-}{H}_0{e}^{-\mathit{k\tau }}+G_1\tau +G_2^{+}$$ E25a

$$I_{\mathrm{I}0}^{-}={\alpha }^{-}{K}_0{e}^{-k\left({\tau }_0-\tau \right)}+{\alpha }^{+}{H}_0{e}^{-\mathit{k\tau }}+G_1\tau +G_2^{-}$$ E25b

where $k^2=\left({\gamma }_1^0+{\gamma }_2^0\right)\left({\gamma }_1^0-{\gamma }_2^0\right)$ , ${\alpha }^{\pm }=\frac{1}{2}\left(1\pm \frac{k}{{\gamma }_1^0+{\gamma }_2^0}\right)$ , $G_1=\frac{{\gamma }_3^0}{{\gamma }_1^0-{\gamma }_2^0}\beta$ , $G_2^{\pm }=\frac{{\gamma }_3^0}{{\gamma }_1^0-{\gamma }_2^0}{B}_0\pm \frac{{\mathit{\beta \gamma }}_3^0}{k^2}$ , $H_0=\frac{{\alpha }^{-}{e}^{-k{\tau }_0}\left[\left(G_2^{+}-{\mathit{RG}}_2^{-}\right)+\left(1-R\right)G_1{\tau }_0-\left(1-R\right)B\left(T_s\right)\right]-\left({\alpha }^{+}-R{\alpha }^{-}\right)G_2^{-}}{{\alpha }^{+}\left({\alpha }^{+}-R{\alpha }^{-}\right)-{\alpha }^{-}\left({\alpha }^{-}-R_{\mathit{dif}}{\alpha }^{+}\right)e^{-2k{\tau }_0}}$ , $K_0=-\frac{{\alpha }^{+}{H}_0+G_2^{-}}{{\alpha }^{-}{e}^{-k{\tau }_0}}$ , and $R=1-{\epsilon }_s$ .

The equations for $I_{\mathrm{I}i}^{\pm }$ (i = 1, 2) are

$$\frac{{\mathit{dI}}_{\mathrm{I}i}^{+}}{\mathit{d\tau }}={\gamma }_1^0{I}_{\mathrm{I}i}^{+}-{\gamma }_2^0{I}_{\mathrm{I}i}^{-}+\left(e^{-a_i\tau }-e^{-a_i{\tau }_0/2}\right)\left[{\gamma }_1^i{I}_{\mathrm{I}0}^{+}-{\gamma }_2^i{I}_{\mathrm{I}0}^{-}-{\gamma }_3^{i}B\left(\tau \right)\right]$$ E26a

$$\frac{{\mathit{dI}}_{\mathrm{I}i}^{-}}{\mathit{d\tau }}={\gamma }_2^0{I}_{\mathrm{I}i}^{+}-{\gamma }_1^0{I}_{\mathrm{I}i}^{-}+\left(e^{-a_i\tau }-e^{-a_i{\tau }_0/2}\right)\left[{\gamma }_2^i{I}_{\mathrm{I}0}^{+}-{\gamma }_1^i{I}_{\mathrm{I}0}^{-}+{\gamma }_3^{i}B\left(\tau \right)\right]$$ E26b

$$I_{\mathrm{I}i}^{-}\left(0\right)=0,I_{\mathrm{I}i}^{+}\left({\tau }_0\right)=\left(1-{\epsilon }_s\right)I_{\mathrm{I}i}^{-}\left({\tau }_0\right)$$ E26c

Let $M_i=I_{\mathrm{I}i}^{+}+I_{\mathrm{I}i}^{-}$ and $N_i=I_{\mathrm{I}i}^{+}-I_{\mathrm{I}i}^{-}$ . Eq. (26a) and (26b) yields

$$\frac{{\mathit{dM}}_i}{\mathit{d\tau }}=\left({\gamma }_1^0+{\gamma }_2^0\right)N_i+{\chi }_{1i}^{+}{e}^{-k{\tau }_0+\left(k-a_i\right)\tau }+{\chi }_{2i}^{+}{e}^{-\left(k+a_i\right)\tau }+{\chi }_{3i}^{+}{e}^{-k\left({\tau }_0-\tau \right)}+{\chi }_{4i}^{+}{e}^{-\mathit{k\tau }}+{\chi }_{5i}^{+}+{\chi }_{6i}^{+}{e}^{-a_i\tau }$$ E27a

$$\begin{array}{c}\frac{{\mathit{dN}}_i}{\mathit{d\tau }}=\left({\gamma }_1^0-{\gamma }_2^0\right)M_i+{\chi }_{1i}^{-}{e}^{-k{\tau }_0+\left(k-a_i\right)\tau }+{\chi }_{2i}^{-}{e}^{-\left(k+a_i\right)\tau }+{\chi }_{3i}^{-}{e}^{-k\left({\tau }_0-\tau \right)}+{\chi }_{4i}^{-}{e}^{-\mathit{k\tau }}+{\chi }_{5i}^{-}+{\chi }_{6i}^{-}{e}^{-a_i\tau }\\ +{\chi }_{7i}^{-}\tau +{\chi }_{8i}^{-}\tau e^{-a_i\tau }\end{array}$$ E27b

where ${\chi }_{1i}^{\pm }=K_0\left({\alpha }^{+}\mp {\alpha }^{-}\right)\left({\gamma }_1^i\pm {\gamma }_2^i\right)$ , ${\chi }_{2i}^{\pm }=\mp H_0\left({\alpha }^{+}\mp {\alpha }^{-}\right)\left({\gamma }_1^i\pm {\gamma }_2^i\right)$ , ${\chi }_{3i}^{\pm }=-K_0\left({\alpha }^{+}\mp {\alpha }^{-}\right)\left({\gamma }_1^i\pm {\gamma }_2^i\right)e^{-a_i{\tau }_0/2}$ , ${\chi }_{4i}^{\pm }=\pm H_0\left({\alpha }^{+}\mp {\alpha }^{-}\right)\left({\gamma }_1^i\pm {\gamma }_2^i\right)e^{-a_i{\tau }_0/2}$ , ${\chi }_{5i}^{+}=-\left(G_2^{+}-G_2^{-}\right)\left({\gamma }_1^i+{\gamma }_2^i\right)e^{-a_i{\tau }_0/2}$ , ${\chi }_{5i}^{-}=-\left(G_2^{+}+G_2^{-}\right)\left({\gamma }_1^i-{\gamma }_2^i\right)e^{-a_i{\tau }_0/2}+2B_0{\gamma }_3^i{e}^{-a_i{\tau }_0/2}$ , ${\chi }_{6i}^{+}=\left(G_2^{+}-G_2^{-}\right)\left({\gamma }_1^i+{\gamma }_2^i\right)$ , ${\chi }_{6i}^{-}=\left(G_2^{+}+G_2^{-}\right)\left({\gamma }_1^i-{\gamma }_2^i\right)-2B_0{\gamma }_3^i$ , ${\chi }_{7i}^{-}=\left[-2G_1\left({\gamma }_1^i-{\gamma }_2^i\right)+2{\mathit{\beta \gamma }}_3^i\right]e^{-a_i{\tau }_0/2}$ , and ${\chi }_{8i}^{-}=2G_1\left({\gamma }_1^i-{\gamma }_2^i\right)-2{\mathit{\beta \gamma }}_3^i$ .

From Eq. (27), we can obtain

$$\begin{array}{c}\frac{d^2{M}_i}{d{\tau }^2}=k^2{M}_i+{\varphi }_{1i}^{+}{e}^{-k{\tau }_0+\left(k-a_i\right)\tau }+{\varphi }_{2i}^{+}{e}^{-\left(k+a_i\right)\tau }+{\varphi }_{3i}^{+}{e}^{-k\left({\tau }_0-\tau \right)}+{\varphi }_{4i}^{+}{e}^{-\mathit{k\tau }}+{\varphi }_{5i}^{+}+{\varphi }_{6i}^{+}{e}^{-a_i\tau }\\ +{\varphi }_{7i}^{+}\tau +{\varphi }_{8i}^{+}\tau e^{-a_i\tau }\end{array}$$ E28a

$$\begin{array}{c}\frac{d^2{N}_i}{d{\tau }^2}=k^2{N}_i+{\varphi }_{1i}^{-}{e}^{-k{\tau }_0+\left(k-a_i\right)\tau }+{\varphi }_{2i}^{-}{e}^{-\left(k+a_i\right)\tau }+{\varphi }_{3i}^{-}{e}^{-k\left({\tau }_0-\tau \right)}+{\varphi }_{4i}^{-}{e}^{-\mathit{k\tau }}+{\varphi }_{5i}^{-}+{\varphi }_{6i}^{-}{e}^{-a_i\tau }\\ +{\varphi }_{8i}^{-}\tau e^{-a_i\tau }\end{array}$$ E28b

where ${\varphi }_{1i}^{\pm }=\left({\gamma }_1^0\pm {\gamma }_2^0\right){\chi }_{1i}^{\mp }+\left(k-a_i\right){\chi }_{1i}^{\pm }$ , ${\varphi }_{2i}^{\pm }=\left({\gamma }_1^0\pm {\gamma }_2^0\right){\chi }_{2i}^{\mp }-\left(k+a_i\right){\chi }_{2i}^{\pm }$ , ${\varphi }_{3i}^{\pm }=\left({\gamma }_1^0\pm {\gamma }_2^0\right){\chi }_{3i}^{\mp }+k{\chi }_{3i}^{\pm }$ , ${\varphi }_{4i}^{\pm }=\left({\gamma }_1^0\pm {\gamma }_2^0\right){\chi }_{4i}^{\mp }-k{\chi }_{4i}^{\pm }$ , ${\varphi }_{5i}^{+}=\left({\gamma }_1^0+{\gamma }_2^0\right){\chi }_{5i}^{-}$ , ${\varphi }_{5i}^{-}=\left({\gamma }_1^0-{\gamma }_2^0\right){\chi }_{5i}^{+}+{\chi }_{7i}^{-}$ , ${\varphi }_{6i}^{+}=\left({\gamma }_1^0+{\gamma }_2^0\right){\chi }_{6i}^{-}-a_i{\chi }_{6i}^{+}$ , ${\varphi }_{6i}^{-}=\left({\gamma }_1^0-{\gamma }_2^0\right){\chi }_{6i}^{+}-a_i{\chi }_{6i}^{-}+{\chi }_{8i}^{-}$ , ${\varphi }_{7i}^{+}=\left({\gamma }_1^0+{\gamma }_2^0\right){\chi }_{7i}^{-}$ , ${\varphi }_{8i}^{+}=\left({\gamma }_1^0+{\gamma }_2^0\right){\chi }_{8i}^{-}$ , and ${\varphi }_{8i}^{-}=-a_i{\chi }_{8i}^{-}$ . Thus, the solutions are

$$\begin{array}{c}{M}_i=K_{1i}{e}^{-k\left({\tau }_0-\tau \right)}+H_{1i}{e}^{-\mathit{k\tau }}+P_{1i}^{+}{e}^{-k{\tau }_0+\left(k-a_i\right)\tau }+P_{2i}^{+}{e}^{-\left(k+a_i\right)\tau }+P_{3i}^{+}\tau e^{-k\left({\tau }_0-\tau \right)}+P_{4i}^{+}\tau e^{-\mathit{k\tau }}\\ +P_{5i}^{+}+P_{6i}^{+}{e}^{-a_i\tau }+P_{7i}^{+}\tau +P_{8i}^{+}\tau e^{-a_i\tau }\end{array}$$ E29a

$$\begin{array}{c}{N}_i=K_{2i}{e}^{-k\left({\tau }_0-\tau \right)}+H_{2i}{e}^{-\mathit{k\tau }}+P_{1i}^{-}{e}^{-k{\tau }_0+\left(k-a_i\right)\tau }+P_{2i}^{-}{e}^{-\left(k+a_i\right)\tau }+P_{3i}^{-}\tau e^{-k\left({\tau }_0-\tau \right)}+P_{4i}^{-}\tau e^{-\mathit{k\tau }}\\ +P_{5i}^{-}+P_{6i}^{-}{e}^{-a_i\tau }+P_{8i}^{-}\tau e^{-a_i\tau }\end{array}$$ E29b

where $P_{1i}^{\pm }=\frac{{\varphi }_{1i}^{\pm }}{{\left(k-a_i\right)}^2-k^2}$ , $P_{2i}^{\pm }=\frac{{\varphi }_{2i}^{\pm }}{{\left(k+a_i\right)}^2-k^2}$ , $P_{3i}^{\pm }=\frac{{\varphi }_{3i}^{\pm }}{2k}$ , $P_{4i}^{\pm }=-\frac{{\varphi }_{4i}^{\pm }}{2k}$ , $P_{5i}^{\pm }=-\frac{{\varphi }_{5i}^{\pm }}{k^2}$ , $P_{6i}^{\pm }=\frac{{\varphi }_{6i}^{\pm }\left(a_i^2-k^2\right)+2a_i{\varphi }_{8i}^{\pm }}{{\left(a_i^2-k^2\right)}^2}$ , $P_{7i}^{+}=-\frac{{\varphi }_{7i}^{+}}{k^2}$ , and $P_{8i}^{+}=\frac{{\varphi }_{8i}^{+}}{a_i^2-k^2}$ .

The expressions of $I_{\mathrm{I}i}^{\pm }$ are

$$\begin{array}{c}{I}_{\mathrm{I}i}^{+}=D_{1i}^{+}{e}^{-k\left({\tau }_0-\tau \right)}+D_{2i}^{+}{e}^{-\mathit{k\tau }}+{\sigma }_{1i}^{+}{e}^{-k{\tau }_0+\left(k-a_i\right)\tau }+{\sigma }_{2i}^{+}{e}^{-\left(k+a_i\right)\tau }+{\sigma }_{3i}^{+}\tau e^{-k\left({\tau }_0-\tau \right)}+{\sigma }_{4i}^{+}\tau e^{-\mathit{k\tau }}\\ +{\sigma }_{5i}^{+}+{\sigma }_{6i}^{+}{e}^{-a_i\tau }+{\sigma }_{7i}\tau +{\sigma }_{8i}^{+}\tau e^{-a_i\tau }\end{array}$$ E30a

$$\begin{array}{c}{I}_{\mathrm{I}i}^{-}=D_{1i}^{-}{e}^{-k\left({\tau }_0-\tau \right)}+D_{2i}^{-}{e}^{-\mathit{k\tau }}+{\sigma }_{1i}^{-}{e}^{-k{\tau }_0+\left(k-a_i\right)\tau }+{\sigma }_{2i}^{-}{e}^{-\left(k+a_i\right)\tau }+{\sigma }_{3i}^{-}\tau e^{-k\left({\tau }_0-\tau \right)}+{\sigma }_{4i}^{-}\tau e^{-\mathit{k\tau }}\\ +{\sigma }_{5i}^{-}+{\sigma }_{6i}^{-}{e}^{-a_i\tau }+{\sigma }_{7i}\tau +{\sigma }_{8i}^{-}\tau e^{-a_i\tau }\end{array}$$ E30b

where $D_{1i}^{\pm }=K_{1i}{\alpha }^{\pm }\pm X_i$ , $D_{2i}^{\pm }=H_{1i}{\alpha }^{\mp }\pm Y_i$ , $X_i=\frac{P_{3i}^{+}-{\chi }_{3i}^{+}}{2\left({\gamma }_1^0+{\gamma }_2^0\right)}$ , $Y_i=\frac{P_{4i}^{+}-{\chi }_{4i}^{+}}{2\left({\gamma }_1^0+{\gamma }_2^0\right)}$ , ${\sigma }_{\mathit{ji}}^{\pm }=\frac{1}{2}\left(P_{\mathit{ji}}^{+}\pm P_{\mathit{ji}}^{-}\right)$ . (j = 1, 2, 3, 4, 5, 6, 8), ${\sigma }_{7i}=\frac{1}{2}{P}_{7i}^{+}$ , and $K_{1i}$ and $H_{1i}$ are determined by boundary conditions. By substituting Eq. (30) into the boundary conditions of Eq. (26c), we can obtain

$$\begin{array}{l}{H}_{1i}=\frac{\left({\alpha }^{+}-R{\alpha }^{-}\right)\left(X_i{e}^{-k{\tau }_0}+Y_i-{\sigma }_{1i}^{-}{e}^{-k{\tau }_0}-{\sigma }_{2i}^{-}-{\sigma }_{5i}^{-}-{\sigma }_{6i}^{-}\right)+{\alpha }^{-}{e}^{-k{\tau }_0}\left[\left(R+1\right)X_i+\left(R+1\right)Y_i{e}^{-k{\tau }_0}\right]}{{\alpha }^{+}\left({\alpha }^{+}-R_{\mathit{dif}}{\alpha }^{-}\right)-{\alpha }^{-}\left({\alpha }^{-}-R_{\mathit{dif}}{\alpha }^{-}\right)e^{-2k{\tau }_0}}\\ -\frac{{\alpha }^{-}{e}^{-k{\tau }_0}\left[\left(R{\sigma }_{1i}^{-}-{\sigma }_{1i}^{+}\right)e^{-a_i{\tau }_0}+\left(R{\sigma }_{2i}^{-}-{\sigma }_{2i}^{+}\right)e^{-\left(k+a_i\right){\tau }_0}+\left(R{\sigma }_{3i}^{-}-{\sigma }_{3i}^{+}\right){\tau }_0+\left(R{\sigma }_{4i}^{-}-{\sigma }_{4i}^{+}\right){\tau }_0{e}^{-k{\tau }_0}\right]}{{\alpha }^{+}\left({\alpha }^{+}-R_{\mathit{dif}}{\alpha }^{-}\right)-{\alpha }^{-}\left({\alpha }^{-}-R_{\mathit{dif}}{\alpha }^{-}\right)e^{-2k{\tau }_0}}\\ -\frac{{\alpha }^{-}{e}^{-k{\tau }_0}\left[\left(R{\sigma }_{5i}^{-}-{\sigma }_{5i}^{+}\right)+\left(R{\sigma }_{6i}^{-}-{\sigma }_{6i}^{+}\right)e^{-a_i{\tau }_0}+\left(R-1\right){\sigma }_{7i}{\tau }_0+\left(R{\sigma }_{8i}^{-}-{\sigma }_{8i}^{+}\right){\tau }_0{e}^{-a_i{\tau }_0}\right]}{{\alpha }^{+}\left({\alpha }^{+}-R_{\mathit{dif}}{\alpha }^{-}\right)-{\alpha }^{-}\left({\alpha }^{-}-R_{\mathit{dif}}{\alpha }^{-}\right)e^{-2k{\tau }_0}}\end{array}$$ E31a

$$K_{1i}=\frac{1}{{\alpha }^{-}{e}^{-k{\tau }_0}}\left(X_i{e}^{-k{\tau }_0}+Y_i-{\alpha }^{+}{H}_{1i}-{\sigma }_{1i}^{-}{e}^{-k{\tau }_0}-{\sigma }_{2i}^{-}-{\sigma }_{5i}^{-}-{\sigma }_{6i}^{-}\right)$$ E31b

Finally, the upward and downward fluxes are obtained by

$$F_{\mathrm{I}}^{+}\left(0\right)=\pi I_{\mathrm{I}}^{+}\left(0\right)$$ E32a

$$F_{\mathrm{I}}^{-}\left({\tau }_0\right)=\pi I_{\mathrm{I}}^{-}\left({\tau }_0\right)$$ E32b

All detailed calculation about solar radiation can be found at [18].

3. Results and discussion

We apply the two schemes to idealized medium to investigate its accuracy, and the result has been shown on [16] and [18].

For true cloud medium, because ice clouds' optical properties strongly depend on the complex particle habits [19, 20, 21]. Therefore, we limit our discussion here to water cloud only. According to the observation, the internal LWC (g m⁻³) and droplet radius of the cloud tend to increase with height [22]. To take this feature into account, LWC and droplet cross-sectional area (DCA; cm⁻², m⁻³) should increase linearly from the cloud base to the position near the top of the cloud:

$$\mathrm{LWC}=0.22+0.00008z$$ E33a

$$\mathrm{DCA}=100+z$$ E33b

where 0 < z < z₀. The terms z and z₀ denote the height from the cloud base and the height of the cloud top, respectively. From Eq. (33a) to (33b), the cloud effective radius ( $r_e$ ; μm) and liquid water path (LWP; g m⁻²) can be obtained:

$$r_e\left(z\right)=\frac{3}{4\rho }\frac{\mathrm{LWC}}{\mathrm{DCA}}{10}^{10}$$ E34a

$$\mathrm{LWP}={\int }_0^{z_0}\mathrm{LWC}\mathit{dz}$$ E34b

where ρ (g m⁻³) is the liquid water density. In this case, LWC varies from 0.22 to 0.30 g m⁻³, and $r_e$ varies from 2.06 to 16.50 μm, in which both ranges are consistent with observation [23]. According to [24], we choose LWP = 260 (g m⁻²) to represent low cloud. In the benchmark calculations, z₀ is divided into 100 internal homogeneous sub-layers, although other numbers can be chosen (e.g., 200). In principle, more internal sub-layers should result in more accurate results. We use 100 internal sub-layers throughout this study because having any more makes little difference to the calculated results. Using 100 sub-layers are sufficiently accurate to resolve the vertical internal inhomogeneity of the medium. We use the optical properties of a water cloud in the solar spectral band of 0.25–0.69 μm and at 0.94 μm and in the infrared spectral band of 5–8 μm and 11 μm.

In Figure 1a and b, the benchmark values of the inhomogeneous IOPs and the parameterized results for the spectral band of 0.25–0.69 μm are shown. The parameterized inhomogeneous IOPs are

$$1-\omega \left(\tau \right)=3.979\times {10}^{-7}-1.897\times {10}^{-6}\left(e^{-0.1539\tau }-e^{-0.1539{\tau }_0/2}\right)$$ E35a

$$g\left(\tau \right)=0.8359+0.0289\left(e^{-0.1539\tau }-e^{-0.1539{\tau }_0/2}\right)$$ E35b

where ${\tau }_0=110.84$ . The corresponding results for reflection and absorption are shown in Figure 1c–f. For reflection, the relative error with the homogeneous solution increases from 0.25 to 0.71% as μ₀ increases from 0.01 to 1, whereas the relative error with the inhomogeneous solution increases from 0.05 to 0.14%. For absorption, the relative error is not sensitive to μ₀; it is around 7.4% with the homogeneous solution but around only 1.4% with the inhomogeneous solution.

Figure 1.

For the band of 0.25–0.69 μm, (a-b) show cloud asymmetry factor/single-scattering albedo versus cloud optical depth (a for asymmetry factor; b for single-scattering albedo), (c-d) show the reflectance/absorptance versus solar zenith angle (c for reflectance; d for absorptance) and (e-f) show the relative errors of the homogeneous and inhomogeneous solutions (e for reflectance error, f for absorptance error).

In Figure 2a and b, the benchmark values of the inhomogeneous IOPs and the parameterized results for the wavelength 0.94 μm are shown. The parameterized inhomogeneous IOPs are

$$1-\omega \left(\tau \right)=1.936\times {10}^{-4}-5.263\times {10}^{-4}\left(e^{-0.0357\tau }-e^{-0.0357{\tau }_0/2}\right)$$ E36a

$$g\left(\tau \right)=0.8321-0.0403\left(e^{0.0218\tau }-e^{0.0218{\tau }_0/2}\right)$$ E36b

where ${\tau }_0=54.46$ . Figure 2c–f shows the corresponding results for reflection and absorption. For reflection, the relative error with the homogeneous solution increases from 1.1 to 3.0% as μ₀ increases from 0.01 to 1, whereas the relative error with the inhomogeneous solution increases from 0.7 to 2.0%. For absorption, the relative error is not sensitive to μ₀; it is around 10% with the homogeneous solution but around only 5.7% with the inhomogeneous solution.

Figure 2.

Same as Figure 1 but for the wavelength 0.94 μm.

The benchmark values of IOPs and parameterized results for the band of 5–8 μm are shown in Figure 3a and b. Here, we assume

$$\omega \left(\tau \right)=0.6757-0.3697\left(e^{-0.0142\tau }-e^{-0.0142{\tau }_0/2}\right)$$ E37a

$$g\left(\tau \right)=0.8644+0.1023\left(e^{-0.0155\tau }-e^{0.0155{\tau }_0/2}\right)$$ E37b

where ${\tau }_0=55.85$ . For upward emissivity (Figure 3c and d), the relative errors of both solutions are not sensitive to $F_{\mathrm{I}}^{+}\left({\tau }_0\right)$ ; the errors are around −3% for homogeneous solution and around 1% for inhomogeneous solution. For downward emissivity (Figure 3e and f), the relative error of homogeneous solution is 4% when $F_{\mathrm{I}}^{+}\left({\tau }_0\right)=0$ , while the error of inhomogeneous solution is only 1%. With $F_{\mathrm{I}}^{+}\left({\tau }_0\right)$ increasing from 0 to $5\mathit{\pi B}\left(T\right)$ , the error of homogeneous solution decreases to 0 firstly but then negatively increases to around −10%. The error of inhomogeneous solution shows a similar decreasing-increasing pattern, but the negative increase only reaches about −2%.

Figure 3.

For the band of 5-8 μm, (a-b) show the cloud single-scattering albedo and asymmetry factor versus cloud optical depth, black dots represent the exact values and the blue lines is the fitting results (a for single-scattering albedo; b for asymmetry factor); (c-d) show the upward/downward emissivity versus the ratio of the radiation incident from the bottom to the internal infrared emission of the medium (c for upward emissivity; d for downward emissivity) and (e-f) show the relative errors of the homogeneous and inhomogeneous solutions (e for upward emissivity; f for downward emissivity).

The benchmark values of IOPs and parameterized results for the band of 11 μm are shown in Figure 4a and b. In this case, we assume

$$\omega \left(\tau \right)=0.4623-0.2155\left(e^{-0.1018\tau }-e^{-0.1018{\tau }_0/2}\right)$$ E38a

$$g\left(\tau \right)=0.9118-0.0083\left(e^{0.1087\tau }-e^{0.1087{\tau }_0/2}\right)$$ E38b

where ${\tau }_0=28.23$ . For upward emissivity (Figure 4c and d), the relative error of homogeneous solution is −1.2%, while the error of inhomogeneous solution is less than 0.5%. For downward emissivity (Figure 4e and f), with $F_{\mathrm{I}}^{+}\left({\tau }_0\right)$ increasing from 0 to $5\mathit{\pi B}\left(T\right)$ , the error of homogeneous (inhomogeneous) solution varies from 3 to −11% (from 0 to −1%).

Figure 4.

Same as Figure 3 but for the wavelength of 11 μm.

4. Summary and conclusions

In the above, we have considered the vertically inhomogeneous structures of only cloud and snow, whereas all physical quantities in the atmosphere are vertically inhomogeneous (e.g., the concentrations of all types of gases and aerosols). In current climate models, the vertical layer resolution is far from that required to resolve such vertical inhomogeneity. In this study, we have proposed a new inhomogeneous SRT/IRT solution to address the vertical inhomogeneity by introducing an internal variation of IOPs inside each model layer. This scheme is based on standard perturbation theory and allows us to use the standard solar Eddington solution and standard infrared two-stream solution for homogeneous layers to identify a zeroth-order equation and a first-order equation that includes the inhomogeneous effect. The new SRT/IRT solution can accurately express the inhomogeneous effect in each model layer, and it reduces to the standard solution when the medium is homogeneous.

The new inhomogeneous SRT/IRT solution is a good way to resolve cloud vertical inhomogeneity. In the spectral band of 0.25–0.69 μm, the relative error in the inhomogeneous SRT solution is no more than 1.4%, whereas the error with the homogeneous SRT solution can be up to 7.4%. At the specific wavelength of 0.94 μm, the relative error with the inhomogeneous solution is not more than 5.7% but can be up to 10% with the homogeneous SRT solution. In the band of 5–8 μm, the homogeneous IRT solution is not sensitive to $F_{\mathrm{I}}^{+}\left({\tau }_0\right)$ , and its relative error may reach −3.2% for upward emissivity, whereas the error of inhomogeneous IRT solution is only 1%. With $F_{\mathrm{I}}^{+}\left({\tau }_0\right)$ increasing from 0 to $5\mathit{\pi B}\left(T\right)$ , the error of downward emissivity for homogeneous solution varies from 4 to −10%, while the error ranges from 1 to −2% for inhomogeneous IRT solution. In the band of 11 μm, the relative error of homogeneous IRT solution is around −1.2% for upward emissivity, and the error of inhomogeneous IRT solution is only less than 0.5%. For downward emissivity, the maximum error of homogeneous IRT solution can be up to −11%, and the maximum error of inhomogeneous IRT solution is only around −1% when $F_{\mathrm{I}}^{+}\left({\tau }_0\right)=5\mathit{\pi B}\left(T\right)$ .

In specific spectral bands or at particular wavelengths, the vertical variations in IOPs can typically be fitted easily into Eq. (3) to obtain the required parameters. A simple fitting program can be easily incorporated into a climate model to produce the inhomogeneous IOPs of stratocumulus clouds. If no such cloud inhomogeneity information is available in the current climate models, the vertical variation rates of cloud LWC and DCA can be derived empirically from observations, which show that the vertical variation rates of LWC and DCA in stratocumulus clouds are not very different [5, 7, 8].

In this study, we presented only a single-layer inhomogeneous SRT/IRT solution. To implement the new solution in a climate model, the adding process for layer-to-layer connections has to be solved. Under the homogeneous condition, the single-layer result in reflection and transmission is the same for an upward path and a downward path, but this is not true for an inhomogeneous layer. Therefore, the adding process has to be modified. We will present an algorithm for this multilayer adding process in our next study, in which the climatic impact of inhomogeneous clouds and inhomogeneous snows will be explored. The code base for the inhomogeneous SRT/IRT solution is available from the authors upon request.

Acknowledgments

The work is supported by National Natural Science Foundation of China (41675003) and the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD).