Multi-layer solar radiative transfer considering the vertical variation of inherent microphysical properties of clouds

Yi-Ning Shi; Yi-Ning Shi; Feng Zhang; Feng Zhang; Ka Lok Chan; Thomas Trautmann; Jiangnan Li

doi:10.1364/OE.27.0A1569

1. Introduction

Clouds as well as aerosols play a key role in radiative transfer (RT) and climate simulations [1–7]. In current RT models used for climate simulations, e.g., RRTMG [8], the Fu-Liao RT model [9, 10], the Beijing Climate Center RT model [11–13], the Canadian Climate Center RT model [14, 15], clouds are vertically divided into a few homogeneous layers and the Eddington approximation [16–20] is used to handle solar RT in each layer. In order to reduce computational demand, the vertical resolution of many climate models is rather coarse, e.g., CSIRO-MK3.6.0 [21] as well as CESM [22] have only 18 levels for total atmosphere, CanAM4 [23] have 32 levels, MPI-ESM [24] have 47 levels.

On the other hand, many observations (e.g., [25, 26]) show that the inherent microphysical properties of cloud, e.g., liquid water content and droplet radius, vary with height. It leads to vertically inhomogeneous optical properties within clouds. In addition, optical properties within aerosol layers are also vertically inhomogeneous according to observations [27–32]. Therefore, the vertical variation of inherent microphysical properties is an essential character of clouds and aerosols.

There are many ways to obtain the vertical profiles of liquid water content (LWC) as well as droplet radius which can be used to acquire optical properties within clouds (e.g., [33–36]) in reality. However, due to the coarse vertical resolution, current climate models usually ignore the vertical variation of microphysical properties within clouds. This assumption has been reported to contribute >10 % error for radiances and several percents for irradiances in the RT calculation [37]. In order to take the vertical variation into account, Li et al. used a Monte Carlo RT model to simulate the reflectance/transmittance of a vertically inhomogeneous layer [38]. The results show that the error related to vertical inhomogeneity of inherent microphysical properties is up to 15 % in the reflectance. Considering the low computational efficiency of the Monte Carlo RT model, this technique is difficult to be applied to climate simulations [19]. In addition, vertical inhomogeneity of inherent microphysical properties can also be considered by increasing the vertical resolution of climate models. However, it is not only computationally expensive but also difficult to be implemented as a climate model which involves many other processes. Therefore, providing a feasible RT solution to consider the effect corresponding to the vertical variability of inherent microphysical properties in climate models is an interesting research topic. It could be an essential step forward for the use of RT models in current climate models.

Zhang et al. proposed a new single-layer solar RT solution in which exponential expressions are used to represent the optical properties which are obtained from inhomogeneous microphysical properties and the perturbation method coupled with the Eddington approximation is applied to solve RT equations [4]. In order to consider the radiation effect related to the vertical variation of microphysical properties of clouds in climate simulations, it is necessary to have a multi-layer RT scheme based on the new inhomogeneous solution. Therefore, we proposed the multi-layer RT scheme by including a modified adding method based on Chandrasekhar’s invariance principle. We also compare the proposed scheme and the conventional Eddington approximation to the Eddington approximation scheme with the high vertical resolution to show the superiority of the proposed scheme.

Besides, it is well-known that delta scaling [39] improves the accuracy of reflectance/transmittance result by separating the contribution to the forward scattered energy from the optical properties. The delta scaling has been applied to many RT schemes (e.g. [15, 40–47]). In the proposed scheme, the optical properties of each layer are also adjusted following the delta scaling. The basic RT equation, coupled with the delta scaling, is derived in the following section. A perturbation method coupled with the Eddington approximation is used to solve the RT equation (the detailed derivation is shown in the Appendix). In section 3, we present the modified adding method to resolve multi-layer RT. Investigation of the accuracy and the computational efficiency of the proposed scheme for two double-layer cases as well as for a realistic atmospheric profile are presented in section 4. A summary is given in section 5.

2. Single-layer RT solution

2.1. Basic RT equation

The RT in the plane-parallel atmosphere can be described by Eq. (1).

μ \frac{d I (τ, μ)}{d τ} = I (τ, μ) - \frac{ω (τ)}{2} \int_{- 1}^{1} I (τ, μ^{'}) P (τ, μ, μ^{'}) d μ^{'}

where

I (τ, μ)

is the radiative intensity, τ

represents the optical depth, μ is the cosine of the zenith angle [ $μ > 0$ ( $μ < 0$ ) denotes the upward (downward) light]. $ω (τ)$ refers to the single scattering albedo and $P (τ, μ, μ^{'})$ refers to the azimuthally averaged scattering phase function.

According to the theory of delta scaling [39], the normalized azimuthally averaged scattering phase function can be approximated as

P (τ, μ, μ^{'}) = 2 f δ (μ - μ^{'}) + (1 - f) [1 + 3 g^{*} (τ) μ μ^{'}]

where f is the truncated fraction and

g^{*} (τ)

is the asymmetry factor after delta scaling,

g^{*} (τ) = \frac{g (τ) - f}{1 - f}

By substituting Eq. (2) to Eq. (1), we obtain

μ \frac{d I (τ, μ)}{d τ} = I (τ, μ) [1 - ω (τ) f] - \frac{ω (τ) (1 - f)}{2} \int_{- 1}^{1} I (τ, μ^{'}) [1 + 3 g^{*} (τ) μ μ^{'}] d μ^{'}

In current climate simulations, almost all RT models employ fixed $ω (τ)$ and $g (τ)$ for cloudy layers. Thus the vertical variations of optical properties caused by inhomogeneous microphysical properties within clouds are often ignored in the RTcalculation. In order to account for the vertical variation, Zhang et al. use exponential expressions to represent the vertical variation of the single scattering albedo and the asymmetry factor [4]

ω (τ) = \hat{ω} + ε_{ω} (e^{- a_{1} τ} - e^{- a_{1} τ_{0} / 2})

g (τ) = \hat{g} + ε_{g} (e^{- a_{2} τ} - e^{- a_{2} τ_{0} / 2})

where ε_ω, a₁, ε_g and a₂ are parameters indicating the vertical variation of

ω (τ)

and

g (τ)

. A schematic illustration is shown in Fig. 1. Usually, optical depth (τ) is used to represent the height (z) in the RT calculation. We assume τ = 0 at the top and

τ = τ_{0}

at the bottom of a cloudy layer. There is a one-to-one match between each height (z) and each optical depth (τ). By substituting Eq. (6) to Eq. (3), we rewrite the expression of the asymmetry factor as

g^{*} (τ) = {\hat{g}}^{*} + ε_{g}^{*} (e^{- a_{2} τ} - e^{- a_{2} τ_{0} / 2})

where

{\hat{g}}^{*} = (\hat{g} - f) / (1 - f)

and

ε_{g}^{*} = ε_{g} / (1 - f)

. According to [39], we choose

f = {\hat{g}}^{2}

.

Fig. 1 The schematic illustration of solar radiative transfer in a cloudy layer.

Download Full Size | PDF

The intensity is divided into two parts: I^dir and I^sca. I^dir is related to the direct sunlight and I^sca is related to the scattered sunlight. Correspondingly, we divided the RT equation shown in Eq. (4) also into two parts.

2.2. RT solution of direct sunlight

The RT equation of direct sunlight is

μ \frac{d I^{dir} (τ, μ)}{d τ} = I^{dir} (τ, μ) [1 - ω (τ) f]

I^{dir} (0, - μ) = δ (μ - μ_{0}) \frac{F_{0}}{2 π}

I^{dir} (τ_{0}, μ) = 0

where μ₀ refers to the cosine of solar zenith angle,

F_{0}

denotes the incident solar flux at the top of the atmosphere (TOA). By substituting Eq. (5) to Eq. (8), the solution can be rewritten as Eqs. (11) - (13)

I^{dir} (τ, μ) = 0

I^{dir} (τ, - μ) = δ (μ - μ_{0}) \frac{F_{0}}{2 π} e^{(- n_{1} τ / μ)} J (τ, μ)

with, J (τ, μ) = exp {[f ε_{ω} (\frac{1}{a_{1}} - \frac{1}{a_{1}} e^{- a_{1} τ} - e^{- a_{1} τ_{0} / 2} τ)] / μ}

where

n_{1} = 1 - \hat{ω} f

. The expression for radiative fluxes is

F_{\pm} (τ) = 2 π \int_{0}^{\pm 1} I (τ, μ) μ d μ

where ± refer to upward/downward direction. Therefore, the upward/downward fluxes of direct sunlight are

F_{+}^{dir} (0) = 0

F_{-}^{dir} (τ_{0}) = μ_{0} F_{0} e^{- n_{1} τ_{0} / μ_{0}} J (τ_{0}, μ_{0})

The transmittance of direct sunlight is

t^{dir} (μ_{0}) = e^{- n_{1} τ_{0} / μ_{0}} J (τ_{0}, μ_{0})

2.3. RT solution of scattered sunlight

There are two types of RT mechanisms for the scattered sunlight: one (direct incident radiation) is associated with the incident direct sunlight, and the other (diffuse incident radiation) is associated with the incident diffuse sunlight. The diffuse incident radiation plays critical roles in the adding method which is used to solve the multi-layer RT. The relationship between the direct incident radiation and the diffuse incident radiation has been discussed in [48].

2.3.1. Direct incident radiation

The RT equation associated with direct incident radiation is

\begin{matrix} μ \frac{d I^{s c a} (τ, μ)}{d τ} = I^{sca} (τ, μ) [1 - ω (τ) f] - \frac{ω (τ) (1 - f)}{2} \\ \times \int_{- 1}^{1} [I^{dir} (τ, μ) + I^{sca} (τ, μ)] [1 + 3 g^{*} (τ) μ μ^{'}] d μ^{'} \end{matrix}

I^{sca} (0, - μ) = 0, I^{sca} (τ_{0}, μ) = 0

We use the Eddington approximation to solve the Eqs.. According to the Eddington approximation, the intensity can be separated by using Legendre polynomials

I^{sca} (τ, μ) = \sum_{l = 0}^{1} I_{l}^{sca} (τ) P_{l} (μ)

where

P_{l} (μ)

is the l th-order Legendre polynomial and

I_{l}^{sca} (τ)

is components of intensity on each Legendre polynomials. By substituting Eqs. (11) - (13) and (20) to Eq. (18) and then using Eq. (14) to convert intensity to flux, the resulting Eqs. can be rewritten as

\frac{d F_{+}^{sca} (τ)}{d τ} = γ_{1} (τ) F_{+}^{sca} (τ) - γ_{2} (τ) F_{-}^{sca} (τ) - γ_{3}^{+} (τ) e^{(- n_{1} τ / μ_{0})} J (τ, μ_{0})

\frac{d F_{-}^{sca} (τ)}{d τ} = γ_{2} (τ) F_{+}^{sca} (τ) - γ_{1} (τ) F_{-}^{sca} (τ) + γ_{3}^{-} (τ) e^{(- n_{1} τ / μ_{0})} J (τ, μ_{0})

F_{-}^{sca} (0) = 0, F_{+}^{sca} (τ_{0}) = 0

where

γ_{1} (τ) = 1.75 - (0.75 f + 1) ω (τ) - 0.75 (1 - f) ω (τ) g^{*} (τ)

γ_{2} (τ) = - 0.25 - (0.75 f - 1) ω (τ) - 0.75 (1 - f) ω (τ) g^{*} (τ)

γ_{3}^{\pm} (τ) = (1 - f) ω (τ) [0.5 \mp 0.75 g^{*} (τ) μ_{0}] F_{0}

The perturbation method is applied to solve the Eqs. The solution of Eqs. (21) - (23) as well as the detailed derivation are presented in the Appendix A.

The reflectance and transmittance of direct incident radiation are

r (μ_{0}) = 2 \int_{0}^{1} R (μ, μ_{0}) μ d μ = \frac{F_{+}^{sca} (0)}{μ_{0} F_{0}}

t (μ_{0}) = 2 \int_{0}^{1} T (μ, μ_{0}) μ d μ = \frac{F_{-}^{sca} (τ_{0})}{μ_{0} F_{0}}

where

R (μ, μ_{0}) = \frac{π I^{s c a} (τ = 0, μ)}{μ_{0} F_{0}}

and

T (μ, μ_{0}) = \frac{π I^{s c a} (τ = τ_{0}, - μ)}{μ_{0} F_{0}}

are the zenith angle dependent reflectance and transmittance of direct incident radiation, respectively.

2.3.2. Diffuse incident radiation

According to [49, 50], the incoming radiance is assumed to be isotropic (the incoming radiance from each angle is the same) and both terms $γ_{3}^{\pm} (τ) e^{(- n_{1} τ / μ_{0})} J (τ, μ_{0})$ associated with direct incident radiation are equal to 0 in the RT equation. For a vertically homogeneous layer, the reflectance/transmittance of RT when the light source at the top are the same as when the light source at the bottom. However, this is not valid for a vertically inhomogeneous layer. We have to calculate the reflectance/transmittance of these two cases separately.

The RT equation when the light source at the top is

\frac{d {\bar{F}}_{+}^{sca} (τ)}{d τ} = γ_{1} (τ) {\bar{F}}_{+}^{sca} (τ) - γ_{2} (τ) {\bar{F}}_{-}^{sca} (τ)

\frac{d {\bar{F}}_{-}^{sca} (τ)}{d τ} = γ_{2} (τ) {\bar{F}}_{+}^{sca} (τ) - γ_{1} (τ) {\bar{F}}_{-}^{sca} (τ)

{\bar{F}}_{-}^{sca} (0) = 1, {\bar{F}}_{+}^{sca} (τ_{0}) = 0

The RT equation when the light source at the bottom is

\frac{d {\bar{F}}_{+}^{sca, *} (τ)}{d τ} = γ_{1} (τ) {\bar{F}}_{+}^{sca, *} (τ) - γ_{2} (τ) {\bar{F}}_{-}^{sca, *} (τ)

\frac{d {\bar{F}}_{-}^{sca, *} (τ)}{d τ} = γ_{2} (τ) {\bar{F}}_{+}^{sca, *} (τ) - γ_{1} (τ) {\bar{F}}_{-}^{sca, *} (τ)

{\bar{F}}_{-}^{sca, *} (0) = 0, {\bar{F}}_{+}^{sca, *} (τ_{0}) = 1

The solution of Eqs. (29) - (31) and Eqs. (32) - (34) as well as the detailed derivation are presented in the Appendix B. The reflectance and transmittance of diffuse incident radiation are

\bar{r} = 2 \int_{0}^{1} \bar{R} (μ) μ d μ = {\bar{F}}_{+}^{sca} (0), \bar{t} = 2 \int_{0}^{1} \bar{T} (μ) μ d μ = {\bar{F}}_{-}^{sca} (τ_{0})

{\bar{r}}^{*} = 2 \int_{0}^{1} {\bar{R}}^{*} (μ) μ d μ = {\bar{F}}_{-}^{sca, *} (τ_{0}), {\bar{t}}^{*} = 2 \int_{0}^{1} {\bar{T}}^{*} (μ) μ d μ = {\bar{F}}_{+}^{sca, *} (0)

where

\bar{R} (μ) = π {\bar{I}}^{sca} (τ = 0, μ)

,

\bar{T} (μ) = π {\bar{I}}^{sca} (τ = τ_{0}, - μ)

are the zenith angle dependent reflectance and transmittance of RT when the light source at the top, respectively.

{\bar{R}}^{*} (μ) = π {\bar{I}}^{sca, *} (τ = τ_{0}, - μ)

and

{\bar{T}}^{*} (μ) = π {\bar{I}}^{sca, *} (τ = 0, μ)

are reflectance and transmittance of RT when the light source at the bottom.

3. Multi-layer RT calculation

In this section, we present the adding method to solve multi-layer RT. The technique is first applied to calculate reflectance/transmittance of the combination of two layers and then further extended to multi-layer RT cases.

3.1. Combination of two layers

The schematic diagram of double-layer case has been shown in [50]. $R_{1, 2}$ ( $T_{1, 2}$ ) is the zenith angle-dependent reflectance (transmittance) of the combination. U and D are the upward and downward dimensionless intensities at the interface between two layers. The meanings of $R_{1, 2}^{*}$ , $T_{1, 2}^{*}$ , $U^{*}$ , $D^{*}$ are same as $R_{1, 2}$ , $T_{1, 2}$ , U, D but for the case with the incident light coming from the bottom.

3.1.1. Direct incident radiation

In RT associated with the direct incident radiation, the expressions of $u (μ_{0})$ , $d (μ_{0})$ are

u (μ_{0}) = 2 \int_{0}^{1} U (μ, μ_{0}) μ d μ, d (μ_{0}) = 2 \int_{0}^{1} D (μ, μ_{0}) μ d μ

where

U (μ, μ_{0}) = \frac{π I^{sca} (τ = τ_{1}, μ)}{μ_{0} F_{0}}

and

D (μ, μ_{0}) = \frac{π I^{sca} (τ = τ_{1}, - μ)}{μ_{0} F_{0}}

Following [19, 51], the four principles of invariance are

U (μ, μ_{0}) = R_{2} (μ, μ_{0}) t_{1}^{dir} (μ_{0}) + 2 \int_{0}^{1} {\bar{R}}_{2} (μ, μ^{'}) D (μ^{'}, μ_{0}) μ^{'} d μ^{'}

D (μ, μ_{0}) = T_{1} (μ, μ_{0}) + 2 \int_{0}^{1} {\bar{R}}_{1}^{*} (μ, μ^{'}) U (μ^{'}, μ_{0}) μ^{'} d μ^{'}

R_{1, 2} (μ, μ_{0}) = R_{1} (μ, μ_{0}) + 2 \int_{0}^{1} {\bar{T}}_{1}^{*} (μ, μ^{'}) U (μ^{'}, μ_{0}) μ^{'} d μ^{'}

T_{1, 2} (μ, μ_{0}) = T_{2} (μ, μ_{0}) t_{1}^{dir} (μ_{0}) + 2 \int_{0}^{1} {\bar{T}}_{2} (μ, μ^{'}) D (μ^{'}, μ_{0}) μ^{'} d μ^{'}

The subscripts 1 and 2 refer to the first layer and the second layer, respectively. By multiplying $2 μ$ in both sides of Eqs. (38) - (41) and integrating over μ from 0 to 1, we can obtain the following Eqs. from Eqs. (27), (28), (35), (36) and (37):

u (μ_{0}) = r_{2} (μ_{0}) t_{1}^{dir} (μ_{0}) + {\bar{r}}_{2} d (μ_{0})

d (μ_{0}) = t_{1} (μ_{0}) + {\bar{r}}_{1}^{*} u (μ_{0})

r_{1, 2} (μ_{0}) = r_{1} (μ_{0}) + {\bar{t}}_{1}^{*} u (μ_{0})

t_{1, 2} (μ_{0}) = t_{2} (μ_{0}) t_{1}^{dir} (μ_{0}) + {\bar{t}}_{2} d (μ_{0})

From Eqs. (42) - (45), we obtain

u (μ_{0}) = \frac{r_{2} (μ_{0}) t_{1}^{dir} (μ_{0}) + {\bar{r}}_{2} t_{1} (μ_{0})}{1 - {\bar{r}}_{2} {\bar{r}}_{1}^{*}}

d (μ_{0}) = t_{1} (μ_{0}) + {\bar{r}}_{1}^{*} \frac{r_{2} (μ_{0}) t_{1}^{dir} (μ_{0}) + {\bar{r}}_{2} t_{1} (μ_{0})}{1 - {\bar{r}}_{2} {\bar{r}}_{1}^{*}}

r_{1, 2} (μ_{0}) = r_{1} (μ_{0}) + \frac{{\bar{t}}_{1}^{*} [r_{2} (μ_{0}) t_{1}^{dir} (μ_{0}) + {\bar{r}}_{2} t_{1} (μ_{0})]}{1 - {\bar{r}}_{2} {\bar{r}}_{1}^{*}}

t_{1, 2} (μ_{0}) = t_{2} (μ_{0}) t_{1}^{dir} (μ_{0}) + {\bar{t}}_{2} t_{1} (μ_{0}) + \frac{{\bar{t}}_{2} {\bar{r}}_{1}^{*} [r_{2} (μ_{0}) t_{1}^{dir} (μ_{0}) + {\bar{r}}_{2} t_{1} (μ_{0})]}{1 - {\bar{r}}_{2} {\bar{r}}_{1}^{*}}

3.1.2. Diffuse incident radiation

The four principles of invariance are also applied to RT associated with the diffuse incident radiation. For downwelling RT,

\bar{U} (μ) = 2 \int_{0}^{1} {\bar{R}}_{2} (μ, μ^{'}) \bar{D} (μ^{'}) μ^{'} d μ^{'}

\bar{D} (μ) = {\bar{T}}_{1} (μ) + 2 \int_{0}^{1} {\bar{R}}_{1}^{*} (μ, μ^{'}) \bar{U} (μ^{'}) μ^{'} d μ^{'}

{\bar{R}}_{1, 2} (μ) = {\bar{R}}_{1} (μ) + 2 \int_{0}^{1} {\bar{T}}_{1}^{*} (μ, μ^{'}) \bar{U} (μ^{'}) μ^{'} d μ^{'}

{\bar{T}}_{1, 2} (μ) = 2 \int_{0}^{1} {\bar{T}}_{2} (μ, μ^{'}) \bar{D} (μ^{'}) μ^{'} d μ^{'}

By applying the similar procedure as Eqs. (38) - (41) to Eqs. (42) - (45), we obtain

\bar{u} = {\bar{r}}_{2} \bar{d}

\bar{d} = {\bar{t}}_{1} + {\bar{r}}_{1}^{*} \bar{u}

{\bar{r}}_{1, 2} = {\bar{r}}_{1} + {\bar{t}}_{1}^{*} \bar{u}

{\bar{t}}_{1, 2} = {\bar{t}}_{2} \bar{d}

From Eqs. (54) - (57), we can derive the reflectance/transmittance of the combination of downwelling RT

{\bar{r}}_{1, 2} = {\bar{r}}_{1} + \frac{{\bar{t}}_{1}^{*} {\bar{r}}_{2} {\bar{t}}_{1}}{1 - {\bar{r}}_{2} {\bar{r}}_{1}^{*}}

{\bar{t}}_{1, 2} = \frac{{\bar{t}}_{1} {\bar{t}}_{2}}{1 - {\bar{r}}_{2} {\bar{r}}_{1}^{*}}

Similar to downwelling RT, the reflectance/transmittance of the combination of upwelling RT are

{\bar{r}}_{1, 2}^{*} = {\bar{r}}_{2}^{*} + \frac{{\bar{t}}_{2}^{*} {\bar{r}}_{1} {\bar{t}}_{2}}{1 - {\bar{r}}_{2} {\bar{r}}_{1}^{*}}

{\bar{t}}_{1, 2}^{*} = \frac{{\bar{t}}_{1}^{*} {\bar{t}}_{2}^{*}}{1 - {\bar{r}}_{2} {\bar{r}}_{1}^{*}}

3.2. Multiple layers

To solve the multi-layer RT, Eqs. (49) and (60) are extended from 2 layers to k layers ( $k > 2$ ) following a path from TOA (layer 1) to the ground (layer N).

\begin{matrix} t_{1, k} (μ_{0}) = t_{k} (μ_{0}) t_{1, k - 1}^{dir} (μ_{0}) + {\bar{t}}_{k} t_{1, k - 1} (μ_{0}) \\ + \frac{{\bar{t}}_{k} {\bar{r}}_{1, k - 1}^{*} [r_{k} (μ_{0}) t_{1, k - 1}^{d i r} (μ_{0}) + {\bar{r}}_{k} t_{1, k - 1} (μ_{0})]}{1 - {\bar{r}}_{k} {\bar{r}}_{1, k - 1}^{*}} \end{matrix}

{\bar{r}}_{1, k}^{*} = {\bar{r}}_{k}^{*} + \frac{{\bar{t}}_{k}^{*} {\bar{r}}_{1, k - 1} {\bar{t}}_{k}}{1 - {\bar{r}}_{k} {\bar{r}}_{1, k - 1}^{*}}

where

t_{1, k - 1}^{dir} (μ_{0}) = \prod_{i = 1}^{k - 1} t_{i}^{dir} (μ_{0})

. We also extend Eqs. (48) and (58) from 2 layers to k layers following a path from the ground (layer N) to TOA (layer 1).

r_{k, N} (μ_{0}) = r_{k} (μ_{0}) + \frac{{\bar{t}}_{k}^{*} [r_{k + 1, N} (μ_{0}) t_{k}^{dir} (μ_{0}) + {\bar{r}}_{k + 1, N} t_{k} (μ_{0})]}{1 - {\bar{r}}_{k + 1, N} {\bar{r}}_{k}^{*}}

{\bar{r}}_{k, N} = {\bar{r}}_{k} + \frac{{\bar{t}}_{k}^{*} {\bar{r}}_{k + 1, N} {\bar{t}}_{k}}{1 - {\bar{r}}_{k + 1, N} {\bar{r}}_{k}^{*}}

where

r_{N} (μ_{0})

and

{\bar{r}}_{N}

are the direct and diffuse surface albedo of the ground, respectively.

u (μ_{0})

and

d (μ_{0})

at level

k + 1

(the bottom of layer k) can be obtained from Eqs. (46) and (47).

u_{k + 1} (μ_{0}) = \frac{r_{k + 1, N} (μ_{0}) t_{1, k}^{dir} (μ_{0}) + {\bar{r}}_{k + 1, N} t_{1, k} (μ_{0})}{1 - {\bar{r}}_{k + 1, N} {\bar{r}}_{1, k}^{*}}

d_{k + 1} (μ_{0}) = t_{1, k} (μ_{0}) + {\bar{r}}_{1, k}^{*} \frac{r_{k + 1, N} (μ_{0}) t_{1, k}^{dir} (μ_{0}) + {\bar{r}}_{k + 1, N} t_{1, k} (μ_{0})}{1 - {\bar{r}}_{k + 1, N} {\bar{r}}_{1, k}^{*}}

The upward and downward radiative fluxes at level $k + 1$ are

F_{+, k + 1} = μ_{0} F_{0} u_{k + 1} (μ_{0})

F_{-, k + 1} = μ_{0} F_{0} d_{k + 1} (μ_{0}) + μ_{0} F_{0} t_{1, k}^{dir} (μ_{0})

The upward and downward radiative fluxes at TOA (level 1) are

F_{+, 1} = μ_{0} F_{0} r_{1, N} (μ_{0})

F_{-, 1} = μ_{0} F_{0}

4. Comparisons and discussions

4.1. Combinations of two layers

In this subsection, the accuracy of the proposed scheme (inhomogeneous scheme) is investigated by comparing it to the conventional Eddington scheme (homogeneous scheme) for two double-layer cases and each layer with vertically inhomogeneous optical properties. According to [37], we divide each layer into 100 homogeneous sub-layers, and the optical properties vary with each sub-layer to represent the vertical inhomogeneity. The Eddington approximation is applied to handle RT in the benchmark.

The reason for using the Eddington approximation instead of a standard RT algorithm (e.g. DIScrete Ordinates Radiative Transfer (DISORT)) in benchmark is that the Eddington approximation has been widely used in current climate models. Therefore, the differences between the Eddington approximation and DISORT should be negligible. Using the Eddington method as reference also helps to identify the discrepancies due to the vertical inhomogeneity of cloud optical properties. Otherwise, it is difficult to separate sources of discrepancies. Besides, the accuracy of the Eddington approximation already has been investigated for idealized cases and real cases in [39, 52–54] by comparing to double adding method, DISORT and 1D/3DMonte Carlo RT model. The results show that the Eddington approximation follows exact solution well, especially for large solar zenith angle cases.

Fig. 2 Relative differences between the homogeneous solution and the benchmark via optical depth of combination τ_total and solar zenith angle μ₀ for (a) reflectance and (c) absorptance. (b) and (d) show the relative differences between the inhomogeneous solution and the benchmark for reflectance and absorptance, respectively.

Download Full Size | PDF

Fig. 3 Same as Fig. 2 but for the combination with ω1 (τ) = 0.95, ω₂ (τ) = 0.93, g1 (τ) = 0.85 + 0.04(e^−0.1^τ − e^−0.1^τ₀^/2)) and g₂ (τ) = 0.8+0.02(e^−0.05^τ − e^−0.05^τ₀^/2).

Download Full Size | PDF

We consider a combination with

ω_{1} (τ) = 0.98 + 0.02 (e^{- 0.2 τ} - e^{- 0.2 τ_{1} / 2})

ω_{2} (τ) = 0.97 + 0.01 (e^{- 0.1 τ} - e^{- 0.1 τ_{2} / 2})

g_{1} (τ) = 0.85, g_{2} (τ) = 0.8

We choose following single scattering albedo and asymmetry factor for the homogeneous solutions.

ω_{1} (τ) = 0.98, ω_{2} (τ) = 0.97

g_{1} (τ) = 0.85, g_{2} (τ) = 0.8

The subscripts 1 and 2 refer to the first layer and the second layer, respectively. The optical depth of the first layer and the second layer are $τ_{1} = 0.4 τ_{total}$ and $τ_{2} = 0.6 τ_{total}$ . The reflectance and absorptance are usedto evaluate the accuracy of the proposed scheme. The definitions of reflectance and absorptance are

Reflectance = F_{+, 1} / μ_{0} F_{0}

Absorptance = 1 - (F_{+, 1} + F_{-, 3}) / μ_{0} F_{0}

where

F_{+, 1}

refers to upward flux at level 1 and

F_{-, 3}

refers to downward flux at level 3 (note that there are three levels in the two layer case). The results of reflectance and absorptance are shown in Fig. 2. We found that the relative errors of reflectance of both schemes increase with

τ_{total}

when the latter increases from 0 to 50. The relative differences of reflectance between the homogeneous scheme and the benchmark range from 0 % to -13.8 % (Fig. 2(a)) while the maximum relative difference of the inhomogeneous scheme is only -0.8 % (Fig. 2(b)). It indicates that the proposed inhomogeneous scheme reduces the error related to the vertical inhomogeneity of optical properties. Similar results are shown in Figs. 2(c) and 2(d). The maximum difference of absorptance between the homogeneous scheme and the benchmark (Fig. 2(c)) is 29.2 %, while the maximum difference of the inhomogeneous solution (Fig. 2(d)) is only 1.7 %.

We consider another case with the following optical properties

ω_{1} (τ) = 0.95, ω_{2} (τ) = 0.93

g_{1} (τ) = 0.85 + 0.04 (e^{- 0.1 τ} - e^{- 0.1 τ_{1} / 2})

g_{2} (τ) = 0.8 + 0.02 (e^{- 0.05 τ} - e^{- 0.05 τ_{2} / 2})

We choose the following single scattering albedo and asymmetry factor for the homogeneous solutions.

ω_{1} (τ) = 0.95, ω_{2} (τ) = 0.93

g_{1} (τ) = 0.85, g_{2} (τ) = 0.8

The optical depth of the first layer and the second layer are $τ_{1} = 0.5 τ_{total}$ and $τ_{2} = 0.5 τ_{total}$ in this case. Ignoring the vertical variation of asymmetry factor leads to a large bias of $9.3 %$ for the reflectance (Fig. 3(a)) and of $- 5.2 %$ for the absorptance (Fig. 3(c)). Again the proposed inhomogeneous scheme is more accurate than the homogeneous scheme (Figs. 3(b) and 3(d)).

Table. 1. The Experiments Setup Description.

View Table | View all tables in this article

4.2. Multilayer atmosphere

The Canadian Climate Center RT model [14, 15] which handle RT in CanAM4 [23] is used for the calculation in this subsection. The RT model uses the correlated-k distribution method for gaseous transmission calculations. Several gases are considered in the model, including H₂O, O₃, CH₄, CO₂, and NO₂. The calculations are performed in the following four solar wavenumber bands 2500-4200, 4200-8400, 8400-14500 and 14500-50000 cm $^{- 1}$ (2.38-4.00, 1.19-2.38, 0.69-1.19 and 0.2-0.69 μm). The parameterizations of [56] are adopted to calculate the optical properties of water clouds. The superiority of this radiative transfer model has been investigated by comparing to other radiative transfer models [52]. The results from the Coupled Model Inter comparison Project Phase 5 (CMIP5) show that this RT model matches well with the NASA Clouds and the Earth’s Radiant Energy System (CERES) satellite observations and the model simulation is mostly close to the averaged results of all models [57].

Fig. 4 The schematic diagram of liquid water content (a) and droplet radius (b) in the inhomogeneous/benchmark schemes and the homogeneous scheme for case 1. The red line and the blue line represent the values of liquid water content and effective radius in the inhomogeneous/benchmark schemes, respectively. The shaded areas represent the liquid water content and droplet radius in the homogeneous scheme.

Download Full Size | PDF

According to [25, 26], the liquid water content (LWC) as well as the droplet radius (r_e) increase with height above the cloud base due to the water vapor condensation. In order to take the vertical variation of the microphysical properties of cloud into account, we assume LWC ( $g m^{- 3}$ ) increases with distance from the cloud base. Two boundary layer marine stratocumulus profiles are chosen. The test setup is summarized in Table 1, where z is the height from the cloud base to the cloud top in meter. LWC ranges from 0.25 $g m^{- 3}$ to 0.28 $g m^{- 3}$ in both case 1 and case 2. These behaviors are consistent with observations [25, 26]. The liquid water path (LWP) which is defined as $\bar{m} L W P = \int_{Cloud_base}^{Cloud_top} LWC d z$ is 260 $g m^{- 2}$ in both cases. It represents the low cloud case [40]. The mean values of LWC are chosen for each layer in the homogeneous scheme.

Under the adiabatic assumption, the expression of $r_{e} (z)$ (m) can be written as [58]

r_{e} (z) = {[\frac{1}{N} \frac{LWC (z)}{\frac{4}{3} π ρ_{w}}]}^{\frac{1}{3}}

where

ρ_{w} = 10^{6} g m^{- 3}

is the density of water and N is the cloud droplet number concentration. According to [35], N is set to

7.5 \times 10^{7} m^{- 3}

and

2.8 \times 10^{8} m^{- 3}

in case 1 and case 2, respectively. In order to make the test setup more clear, a schematic diagram of LWC and r_e for the inhomogeneous/benchmark scheme and the homogeneous scheme in case 1 is shown in Fig. 4. Similar results have also been shown for case 2 (not shown in this paper).

Fig. 5 The benchmark profiles of downward flux (a), upward flux (c) and heating rate (e) for case 1. The differences between homogeneous/inhomogeneous scheme and benchmark for downward flux (b), upward flux (d) and heating rate (f).

Download Full Size | PDF

Fig. 6 Same as Fig. 5 but for case 2.

Download Full Size | PDF

Before running the proposed inhomogeneous scheme, a simple regression code is used in the RT model to obtain the parameter in Eqs. (5) and (6) from bulk optical properties calculated by the parameterizations of [56]. The correlation coefficients for each band are close to 1, which evaluates the validity of Eqs. (5) and (6) when fitting the vertical variation of inherent optical properties.

Figure 5 shows the benchmark of downward flux (Fig. 5(a)), upward flux (Fig. 5(c)), heating rate (Fig. 5(e)) as well as the differences between the homogeneous/inhomogeneous scheme and the benchmark (Figs. 5(b), 5(d) and 5(f)) for case 1. Ignoring the vertical variation of microphysical properties leads to a significant bias of 5.77 W/m² for the downward flux (Fig. 5(b)) and 8.97 W/m² for the upward flux (Fig. 5(d)) while the proposed inhomogeneous scheme reduced the bias to 0.05 W/m² for the downward flux and 0.22 W/m² for the upward flux. The homogeneous scheme caused a <-0.27 K/day bias at the bottom of the cloud while the maximum bias of inhomogeneous scheme is only -0.05 K/day (Fig. 5(f)).

The results for case 2 are shown in Fig. 6. The maximum difference between the homogeneous scheme and the benchmark is 6.04 W/m² for the downward flux (Fig. 6(b)) and -8.15 W/m² for the upward flux (Fig. 6(d)) while the maximum difference between the inhomogeneous scheme and the benchmark is -0.20 W/m² for the downward flux and -0.24 W/m² for the upward flux. For the heating rate calculation (Fig. 6(f)), the inhomogeneous scheme is also more accurate than the homogeneous scheme. All results indicate that the inhomogeneous scheme minimized the error related to the vertical variation of microphysical properties.

Table 2. Computational Time of the Homogeneous/Inhomogeneous radiative transfer scheme and the Benchmark (normalized to the Computational Time of the Homogeneous scheme).

View Table | View all tables in this article

The computational efficiency is essential for any RT scheme which is used in climate models. The computational time for the homogeneous scheme, the inhomogeneous scheme, and the benchmark is listed in Table 2. The results are computed by Dell G7 7588 laptop computer with 8 Intel(R) Core(TM) i5-8300H CPUs, a 64-bit operating system with 8 GB memory. The Intel Fortran (version: 19.0.0.117) which is installed on Ubuntu (version: 18.04.1 LTS) is used for calculations. In Table 2, “Algorithm only” means that we only consider pure radiative transfer solution without gaseous transmission calculation and parameterization of cloud optical properties. “Radiative transfer model” means that we consider the complete radiative transfer process which includes gaseous transmission calculation, parameterization of cloud optical properties, and radiative transfer solution. The inhomogeneous scheme only takes triple the CPU time of homogeneous scheme when we only consider the algorithm itself and twice the CPU time when we consider the complete radiative transfer process. Compared to the benchmark, the proposed inhomogeneous scheme is more efficient to be applied to climate models to consider the radiative effect related to the vertical variation of the microphysical properties. In addition, we note that the increase of vertical resolution can complicate other physical processes in climate models. Therefore, the superiority of our proposed scheme is not confined to increase computing efficiency.

5. Summary and conclusions

The present study aims to propose a multi-layer RT scheme which is able to take into account the vertical variation of the microphysical properties of clouds for climate simulations. The scheme is based on the inhomogeneous solution proposed in [4] in which exponential expressions are used to represent the vertical variation of optical properties caused by inhomogeneous microphysical properties, and a perturbation method is applied to solve the RT equation. The delta scaling [39] is introduced to adjust the optical properties to improve the accuracy of single-layer RT calculation. As the reflectance/transmittance of upwelling RT are different from downwelling RT in a vertically inhomogeneous layer, we modified the adding method to solve the multi-layer RT which takes the vertical inhomogeneity of microphysical properties into account.

Two cases of the combination of two vertically inhomogeneous layers are used to investigate the accuracy of the proposed scheme with a wide range of solar zenith angles and optical depths. The results show that the bias related to the vertical variation of optical properties reaches 13.8 % for reflectance and 29.2 % for absorptance while the bias of the proposed scheme is only -0.8 % for reflectance and 1.7 % for absorptance.

We also apply the proposed scheme to the Canadian Climate Center RT model [14, 15] which handle RT in CanAM4 [23] and compare the results to the ones calculated by the conventional Eddington approximation. Two boundary marine layer stratocumulus profiles are used. Due to ignoring vertical variation of inherent microphysical properties, the conventional Eddington approximation leads to a significant bias of 5.77 W/m² for downward flux, 8.97 W/m² for upward flux and -0.27 K/day for heating rate under case 1 condition. The bias of downward flux, upward flux, and heating rate are -0.05 W/m², -0.22 W/m² and 0.05 K/day, respectively, for calculations with the proposed scheme, which shows the proposed scheme is more accurate than the conventional Eddington approximation. Similar results are also obtained from case 2. All results indicate that the proposed scheme reduced the error associated with the vertical variation of inherent microphysical properties. In terms of computational efficiency, the proposed scheme is approximately three times more computationally expensive compared to the Eddington approximation considering the algorithm itself. The computational time for the proposed scheme is doubled as compared to the Eddington approximation when applying it to the full RT model. Compared to the benchmark, the proposed inhomogeneous scheme is more efficient to be applied to climate models to consider the radiative effect related to the vertical variation of inherent microphysical properties.

There are still some limitations of the present work. The vertical variability of cloud inherent microphysical properties can be highly variable in the realistic atmosphere, and it is not strictly of the form explored in the cases shown in the paper. For example, the vertical profile of cumulus congestus clouds can differ significantly from this type of profile once precipitation processes take over. We will propose more flexible expressions to represent the optical properties of various clouds in the future. In addition, we note that the Eddington approximation makes a great balance between accuracy and computational efficiency in a global simulation scenario. The four-stream inhomogeneous RT scheme which has much higher accuracy can also be proposed by a similar technique. However, according to [49, 50], the computational time of four-stream solver is about four to five times longer than the Eddington approximation under a homogeneous layer assumption. Therefore, the four-stream inhomogeneous RT scheme cannot be applied to global climate simulations in the current state of the art.

Appendix

A. Solution of Equations (21) - (23) in the manuscript

Equations (21) - (23) in the manuscript are

\frac{d F_{+}^{sca} (τ)}{d τ} = γ_{1} (τ) F_{+}^{sca} (τ) - γ_{2} (τ) F_{-}^{sca} (τ) - γ_{3}^{+} (τ) e^{(- n_{1} τ / μ_{0})} J (τ, μ_{0})

\frac{d F_{-}^{sca} (τ)}{d τ} = γ_{2} (τ) F_{+}^{sca} (τ) - γ_{1} (τ) F_{-}^{sca} (τ) + γ_{3}^{-} (τ) e^{(- n_{1} τ / μ_{0})} J (τ, μ_{0})

F_{-}^{sca} (0) = 0, F_{+}^{sca} (τ_{0}) = 0

where

γ_{1} (τ) = 1.75 - (0.75 f + 1) ω (τ) - 0.75 (1 - f) ω (τ) g^{*} (τ)

,

γ_{2} (τ) = - 0.25 - (0.75 f - 1) ω (τ) - 0.75 (1 - f) ω (τ) g^{*} (τ)

,

γ_{3}^{\pm} (τ) = (1 - f) ω (τ) [0.5 \mp 0.75 g^{*} (τ) μ_{0}] F_{0}

. Following the perturbation method proposed by [59], the upward/downward flux are expanded as

F_{+}^{sca} (τ) = F_{+}^{sca, 0} (τ) + ε_{ω} F_{+}^{sca, 1} (τ) + ε_{g}^{*} F_{+}^{sca, 2} (τ)

F_{-}^{sca} (τ) = F_{-}^{sca, 0} (τ) + ε_{ω} F_{-}^{sca, 1} (τ) + ε_{g}^{*} F_{-}^{sca, 2} (τ)

After substituting Eqs. (5) and (6) to $γ_{1} (τ)$ , $γ_{2} (τ)$ , $γ_{3}^{\pm} (τ)$ and then removing terms containing second and higher order terms of ε_ω, ε_g, we obtain

γ_{1} (τ) = γ_{10} + γ_{11} ε_{ω} (e^{- a_{1} τ} - e^{- a_{1} τ_{0} / 2}) + γ_{12} ε_{g}^{*} (e^{- a_{2} τ} - e^{- a_{2} τ_{0} / 2})

γ_{2} (τ) = γ_{20} + γ_{21} ε_{ω} (e^{- a_{1} τ} - e^{- a_{1} τ_{0} / 2}) + γ_{22} ε_{g}^{*} (e^{- a_{2} τ} - e^{- a_{2} τ_{0} / 2})

γ_{3}^{\pm} (τ) = γ_{30}^{\pm} + γ_{31}^{\pm} ε_{ω} (e^{- a_{1} τ} - e^{- a_{1} τ_{0} / 2}) + γ_{32}^{\pm} ε_{g}^{*} (e^{- a_{2} τ} - e^{- a_{2} τ_{0} / 2})

where

γ_{10} = 1.75 - (0.75 f + 1) \hat{ω} - 0.75 (1 - f) \hat{ω} {\hat{g}}^{*}

,

γ_{20} = - 0.25 - (0.75 f - 1) \hat{ω} - 0.75 (1 - f) \hat{ω} {\hat{g}}^{*}

,

γ_{30}^{\pm} = (1 - f) (0.5 \mp 0.75 {\hat{g}}^{*} μ_{0}) \hat{ω} F_{0}

,

γ_{11} = - (0.75 f + 1) - 0.75 (1 - f) {\hat{g}}^{*}

,

γ_{21} = - (0.75 f - 1) - 0.75 (1 - f) {\hat{g}}^{*}

,

γ_{31}^{\pm} = γ_{30}^{\pm} / \hat{ω}

,

γ_{12} = γ_{22} = - 0.75 (1 - f) \hat{ω}

and

γ_{32}^{\pm} = \mp 0.75 (1 - f) μ_{0} F_{0} \hat{ω}

. In addition, we use the Taylor polynomial to expand the term

J (τ, μ_{0})

in Eqs. (85) and (86), and we remove the terms containing second and higher order terms of ε_ω.

J (τ, μ_{0}) = 1 + \frac{f ε_{ω}}{μ_{0}} (\frac{1}{a_{1}} - \frac{1}{a_{1}} e^{- a_{1} τ} - e^{- a_{1} τ_{0} / 2} τ)

By substituting Eqs. (88) - (93) to Eqs. (85) and (86), we obtain

\begin{matrix} \frac{d F_{+}^{s c a} (τ)}{d τ} = [γ_{10} + γ_{11} ε_{ω} (e^{- a_{1} τ} - e^{- a_{1} τ_{0} / 2}) + γ_{12} ε_{g}^{*} (e^{- a_{2} τ} - e^{- a_{2} τ_{0} / 2})] \\ \times [F_{+}^{sca, 0} (τ) + ε_{ω} F_{+}^{sca, 1} (τ) + ε_{g}^{*} F_{+}^{sca, 2} (τ)] - [γ_{20} + γ_{21} ε_{ω} (e^{- a_{1} τ} - e^{- a_{1} τ_{0} / 2}) \\ + γ_{22} ε_{g}^{*} (e^{- a_{2} τ} - e^{- a_{2} τ_{0} / 2})] [F_{-}^{sca, 0} (τ) + ε_{ω} F_{-}^{sca, 1} (τ) + ε_{g}^{*} F_{-}^{sca, 2} (τ)] \\ - [γ_{30}^{+} + γ_{31}^{+} ε_{ω} (e^{- a_{1} τ} - e^{- a_{1} τ_{0} / 2}) + γ_{32}^{+} ε_{g}^{*} (e^{- a_{2} τ} - e^{- a_{2} τ_{0} / 2})] \\ \times e^{- n_{1} τ_{0} / μ_{0}} [1 + \frac{f ε_{ω}}{μ_{0}} (\frac{1}{a_{1}} - \frac{1}{a_{1}} e^{- a_{1} τ} - e^{- a_{1} τ_{0} / 2} τ)] \end{matrix}

\begin{matrix} \frac{d F_{-}^{s c a} (τ)}{d τ} = [γ_{20} + γ_{21} ε_{ω} (e^{- a_{1} τ} - e^{- a_{1} τ_{0} / 2}) + γ_{22} ε_{g}^{*} (e^{- a_{2} τ} - e^{- a_{2} τ_{0} / 2})] \\ \times [F_{+}^{sca, 0} (τ) + ε_{ω} F_{+}^{sca, 1} (τ) + ε_{g}^{*} F_{+}^{sca, 2} (τ)] - [γ_{10} + γ_{11} ε_{ω} (e^{- a_{1} τ} - e^{- a_{1} τ_{0} / 2}) \\ + γ_{12} ε_{g}^{*} (e^{- a_{2} τ} - e^{- a_{2} τ_{0} / 2})] [F_{-}^{sca, 0} (τ) + ε_{ω} F_{-}^{sca, 1} (τ) + ε_{g}^{*} F_{-}^{sca, 2} (τ)] \\ + [γ_{30}^{-} + γ_{31}^{-} ε_{ω} (e^{- a_{1} τ} - e^{- a_{1} τ_{0} / 2}) + γ_{32}^{-} ε_{g}^{*} (e^{- a_{2} τ} - e^{- a_{2} τ_{0} / 2})] \\ \times e^{- n_{1} τ_{0} / μ_{0}} [1 + \frac{f ε_{ω}}{μ_{0}} (\frac{1}{a_{1}} - \frac{1}{a_{1}} e^{- a_{1} τ} - e^{- a_{1} τ_{0} / 2} τ)] \end{matrix}

After removing terms containing second or higher order terms of ε_ω, ε_g, we separate Eqs. (94) and (95) into three Eqs. involving $F_{\pm}^{sca, 0} (τ)$ , $F_{\pm}^{sca, 1} (τ)$ and $F_{\pm}^{sca, 2} (τ)$ .

The resulting Eqs. for $F_{\pm}^{sca, 0} (τ)$ are

\frac{d F_{+}^{sca, 0} (τ)}{d τ} = γ_{10} F_{+}^{sca, 0} (τ) - γ_{20} F_{-}^{sca, 0} (τ) - γ_{30}^{+} e^{- \frac{n_{1} τ}{μ_{0}}}

\frac{d F_{-}^{sca, 0} (τ)}{d τ} = γ_{20} F_{+}^{sca, 0} (τ) - γ_{10} F_{-}^{sca, 0} (τ) + γ_{30}^{-} e^{- \frac{n_{1} τ}{μ_{0}}}

F_{-}^{sca, 0} (0) = 0, F_{+}^{sca, 0} (τ_{0}) = 0

Equations (96) - (98) are the RT equation based on the Eddington approximation in a vertically homogeneous layer. According to [41], the solutions of Eqs. (96) and (97) are

F_{+}^{sca, 0} (τ) = α^{+} C_{10} e^{k (τ - τ_{0})} + α^{-} C_{20} e^{- k τ} + G_{1}^{+} e^{- \frac{n_{1} τ}{μ_{0}}}

F_{-}^{sca, 0} (τ) = α^{-} C_{10} e^{k (τ - τ_{0})} + α^{+} C_{20} e^{- k τ} + G_{1}^{-} e^{- \frac{n_{1} τ}{μ_{0}}}

where

k^{2} = (γ_{10} - γ_{20}) (γ_{10} + γ_{20})

,

α^{\pm} = 0.5 (1 \pm \frac{γ_{10} - γ_{20}}{k})

,

G_{1}^{\pm} = \pm \frac{μ_{0}^{2}}{n_{1}^{2} - μ_{0}^{2} k^{2}} (\frac{n_{1}}{μ_{0}} γ_{30}^{\pm} \mp γ_{10} γ_{30}^{\pm} \mp γ_{20} γ_{30}^{\mp})

. C₁₀ and C₂₀ are determined by the boundary conditions. By substituting Eqs. (99) and (100) to Eq. (98), we get

C_{10} = \frac{α^{-} G_{1}^{-} e^{- k τ_{0}} - α^{+} G_{1}^{+} e^{- \frac{n_{1} τ_{0}}{μ_{0}}}}{{(α^{+})}^{2} - {(α^{-})}^{2} e^{- 2 k τ_{0}}}

C_{20} = \frac{- G_{1}^{-} - α^{-} C_{10} e^{- k τ_{0}}}{α^{+}}

Likewise, the Eqs. for $F_{\pm}^{sca, 1} (τ)$ and $F_{\pm}^{sca, 2} (τ)$ are

\begin{matrix} \frac{d F_{+}^{s c a, i} (τ)}{d τ} = γ_{10} F_{+}^{sca, i} (τ) - γ_{20} F_{-}^{sca, i} (τ) + (γ_{1 i} F_{+}^{sca, 0} (τ) - γ_{2 i} F_{-}^{sca, 0} (τ) - γ_{3 i}^{+} e^{- \frac{n_{1} τ}{μ_{0}}}) \\ \times (e^{- a_{i} τ} - e^{- a_{i} τ_{0} / 2}) - γ_{4 i}^{+} e^{- \frac{n_{1} τ}{μ_{0}}} (\frac{1}{a_{1}} - \frac{1}{a_{1}} e^{- a_{1} τ} - e^{- a_{1} τ_{0} / 2} τ) \end{matrix}

\begin{matrix} \frac{d F_{-}^{s c a, i} (τ)}{d τ} = γ_{20} F_{+}^{sca, i} (τ) - γ_{10} F_{-}^{sca, i} (τ) + (γ_{2 i} F_{+}^{sca, 0} (τ) - γ_{1 i} F_{-}^{sca, 0} (τ) + γ_{3 i}^{-} e^{- \frac{n_{1} τ}{μ_{0}}}) \\ \times (e^{- a_{i} τ} - e^{- a_{i} τ_{0} / 2}) + γ_{4 i}^{-} e^{- \frac{n_{1} τ}{μ_{0}}} (\frac{1}{a_{1}} - \frac{1}{a_{1}} e^{- a_{1} τ} - e^{- a_{1} τ_{0} / 2} τ) \end{matrix}

F_{-}^{sca, i} (0) = 0, F_{+}^{sca, i} (τ_{0}) = 0

where

γ_{41}^{\pm} = \frac{f}{μ_{0}} γ_{30}^{\pm}

,

γ_{42}^{\pm} = 0

and

i = 1, 2

. If we define

M_{i} (τ) = F_{+}^{sca, i} (τ) + F_{-}^{sca, i} (τ)

,

N_{i} (τ) = F_{+}^{sca, i} (τ) - F_{-}^{sca, i} (τ)

then Eqs. (103) and (104) yield

\begin{matrix} \frac{d M_{i} (τ)}{d τ} = (γ_{10} + γ_{20}) N_{i} (τ) + χ_{1 i}^{+} e^{k (τ - τ_{0})} + χ_{2 i}^{+} e^{- k τ} + χ_{3 i}^{+} e^{(k - a_{i}) τ - k τ_{0}} + χ_{4 i}^{+} e^{- (k + a_{i}) τ} \\ + χ_{5 i}^{+} e^{- \frac{n_{1}}{μ_{0}} τ} + χ_{6 i}^{+} e^{- (\frac{n_{1}}{μ_{0}} + a_{i}) τ} + χ_{7 i}^{+} τ e^{- \frac{n_{1}}{μ_{0}} τ} \end{matrix}

\begin{matrix} \frac{d N_{i} (τ)}{d τ} = (γ_{10} - γ_{20}) M_{i} (τ) + χ_{1 i}^{-} e^{k (τ - τ_{0})} + χ_{2 i}^{-} e^{- k τ} + χ_{3 i}^{-} e^{(k - a_{i}) τ - k τ_{0}} + χ_{4 i}^{-} e^{- (k + a_{i}) τ} \\ + χ_{5 i}^{-} e^{- \frac{n_{1}}{μ_{0}} τ} + χ_{6 i}^{-} e^{- (\frac{n_{1}}{μ_{0}} + a_{i}) τ} + χ_{7 i}^{-} τ e^{- \frac{n_{1}}{μ_{0}} τ} \end{matrix}

where

χ_{1 i}^{\pm} = - (γ_{1 i} \pm γ_{2 i}) (α^{+} \mp α^{-}) C_{10} e^{- a_{i} τ_{0} / 2}

,

χ_{2 i}^{\pm} = - (γ_{1 i} \pm γ_{2 i}) (α^{-} \mp α^{+}) C_{20} e^{- a_{i} τ_{0} / 2}

,

χ_{3 i}^{\pm} = (γ_{1 i} \pm γ_{2 i}) (α^{+} \mp α^{-}) C_{10}

,

χ_{4 i}^{\pm} = (γ_{1 i} \pm γ_{2 i}) (α^{-} \mp α^{+}) C_{20}

,

χ_{5 i}^{\pm} = - [(γ_{1 i} \pm γ_{2 i}) (G_{1}^{+} \mp G_{1}^{-}) e^{- a_{i} τ_{0} / 2} - (γ_{3 i}^{+} \mp γ_{3 i}^{-}) e^{- a_{i} τ_{0} / 2} + (γ_{4 i} \mp γ_{4 i}) / a_{i}]

,

χ_{6 i}^{\pm} = (γ_{1 i} \pm γ_{2 i}) (G_{1}^{+} \mp G_{1}^{-}) - (γ_{3 i}^{+} \mp γ_{3 i}^{-}) + (γ_{4 l}^{+} \mp γ_{4 l}^{-}) / a_{i}

and

χ_{7 i}^{\pm} = (γ_{4 i}^{+} \mp γ_{4 i}^{-}) e^{- a_{i} τ_{0} / 2}

. From Eqs. (106) and (107), we can get

\begin{matrix} \frac{d^{2} M_{i} (τ)}{d τ^{2}} = k^{2} M_{i} (τ) + ϕ_{1 i}^{+} e^{k (τ - τ_{0})} + ϕ_{2 i}^{+} e^{- k τ} + ϕ_{3 i}^{+} e^{(k - a_{i}) τ - k τ_{0}} + ϕ_{4 i}^{+} e^{- (k + a_{i}) τ} \\ + ϕ_{5 i}^{+} e^{- \frac{n_{1}}{μ_{0}} τ} + ϕ_{6 i}^{+} e^{- (\frac{n_{1}}{μ_{0}} + a_{i}) τ} + ϕ_{7 i}^{+} τ e^{- \frac{n_{1}}{μ_{0}} τ} \end{matrix}

\begin{matrix} \frac{d^{2} N_{i} (τ)}{d τ^{2}} = k^{2} N_{i} (τ) + ϕ_{1 i}^{-} e^{k (τ - τ_{0})} + ϕ_{2 i}^{-} e^{- k τ} + ϕ_{3 i}^{-} e^{(k - a_{i}) τ - k τ_{0}} + ϕ_{4 i}^{-} e^{- (k + a_{i}) τ} \\ + ϕ_{5 i}^{-} e^{- \frac{n_{1}}{μ_{0}} τ} + ϕ_{6 i}^{-} e^{- (\frac{n_{1}}{μ_{0}} + a_{i}) τ} + ϕ_{7 i}^{-} τ e^{- \frac{n_{1}}{μ_{0}} τ} \end{matrix}

where

ϕ_{1 i}^{\pm} = (γ_{10} \pm γ_{20}) χ_{1 i}^{\mp} + k χ_{1 i}^{\pm}

,

ϕ_{2 i}^{\pm} = (γ_{10} \pm γ_{20}) χ_{2 i}^{\mp} - k χ_{2 i}^{\pm}

,

ϕ_{3 i}^{\pm} = (γ_{10} \pm γ_{20}) χ_{3 i}^{\mp} + (k - a_{i}) χ_{3 i}^{\pm}

,

ϕ_{4 i}^{\pm} = (γ_{10} \pm γ_{20}) χ_{4 i}^{\mp} - (k + a_{i}) χ_{4 i}^{\pm}

,

ϕ_{5 i}^{\pm} = (γ_{10} \pm γ_{20}) χ_{5 i}^{\mp} - \frac{n_{1}}{μ_{0}} χ_{5 i}^{\pm} + χ_{7 i}^{\pm}

,

ϕ_{6 i}^{\pm} = (γ_{10} \pm γ_{20}) χ_{6 i}^{\mp} - \frac{n_{1} + a_{1} μ_{0}}{μ_{0}} χ_{6 i}^{\pm}

, and

ϕ_{7 i}^{\pm} = (γ_{10} \pm γ_{20}) χ_{7 i}^{\mp} - \frac{n_{1}}{μ_{0}} χ_{7 i}^{\pm}

. The solutions of Eqs. (108) and (109) are

\begin{matrix} M_{i} (τ) = C_{1 i} e^{k (τ - τ_{0})} + C_{2 i} e^{- k τ} + P_{1 i}^{+} τ e^{k (τ - τ_{0})} + P_{2 i}^{+} τ e^{- k τ} + P_{3 i}^{+} e^{(k - a_{i}) τ - k τ_{0}} \\ + P_{4 i}^{+} e^{- (k + a_{i}) τ} + P_{5 i}^{+} e^{- \frac{n_{1}}{μ_{0}} τ} + P_{6 i}^{+} e^{- (\frac{n_{1}}{μ_{0}} + a_{i}) τ} + P_{7 i}^{+} τ e^{- \frac{n_{1}}{μ_{0}} τ} \end{matrix}

\begin{matrix} N_{i} (τ) = C_{3 i} e^{k (τ - τ_{0})} + C_{4 i} e^{- k τ} + P_{1 i}^{-} τ e^{k (τ - τ_{0})} + P_{2 i}^{-} τ e^{- k τ} + P_{3 i}^{-} e^{(k - a_{i}) τ - k τ_{0}} \\ + P_{4 i}^{-} e^{- (k + a_{i}) τ} + P_{5 i}^{-} e^{- \frac{n_{1}}{μ_{0}} τ} + P_{6 i}^{-} e^{- (\frac{n_{1}}{μ_{0}} + a_{i}) τ} + P_{7 i}^{-} τ e^{- \frac{n_{1}}{μ_{0}} τ} \end{matrix}

where

P_{1 i}^{\pm} = \frac{ϕ_{1 i}^{\pm}}{2 k}

,

P_{2 i}^{\pm} = - \frac{ϕ_{2 i}^{\pm}}{2 k}

,

P_{3 i}^{\pm} = \frac{ϕ_{3 i}^{\pm}}{{(k - a_{i})}^{2} - k^{2}}

,

P_{4 i}^{\pm} = \frac{ϕ_{4 i}^{\pm}}{{(k + a_{i})}^{2} - k^{2}}

,

P_{7 i}^{\pm} = \frac{μ_{0}^{2} ϕ_{7 i}^{\pm}}{n_{1}^{2} - μ_{0}^{2} k^{2}}

,

P_{5 i}^{\pm} = \frac{μ_{0}^{2} ϕ_{5 i}^{\pm} + 2 μ_{0} n_{1} P_{7 i}^{\pm}}{n_{1}^{2} - μ_{0}^{2} k^{2}}

and

P_{6 i}^{\pm} = \frac{μ_{0}^{2} ϕ_{6 i}^{\pm}}{{(n_{1} + a_{1} μ_{0})}^{2} - μ_{0}^{2} k^{2}}

. By using the definition of

M_{i} (τ)

and

N_{i} (τ)

, we arrive at the expressions for

F_{\pm}^{sca, i} (τ)

.

\begin{matrix} F_{+}^{sca, i} (τ) = D_{1 i}^{+} e^{k (τ - τ_{0})} + D_{2 i}^{+} e^{- k τ} + σ_{1 i}^{+} τ e^{k (τ - τ_{0})} + σ_{2 i}^{+} τ e^{- k τ} + σ_{3 i}^{+} e^{(k - a_{i}) τ - k τ_{0}} \\ + σ_{4 i}^{+} e^{- (k + a_{i}) τ} + σ_{5 i}^{+} e^{- \frac{n_{1}}{μ_{0}} τ} + σ_{6 i}^{+} e^{- (\frac{n_{1}}{μ_{0}} + a_{i}) τ} + σ_{7 i}^{+} τ e^{- \frac{n_{1}}{μ_{0}} τ} \end{matrix}

\begin{matrix} F_{-}^{sca, i} (τ) = D_{1 i}^{-} e^{k (τ - τ_{0})} + D_{2 i}^{-} e^{- k τ} + σ_{1 i}^{-} τ e^{k (τ - τ_{0})} + σ_{2 i}^{-} τ e^{- k τ} + σ_{3 i}^{-} e^{(k - a_{i}) τ - k τ_{0}} \\ + σ_{4 i}^{-} e^{- (k + a_{i}) τ} + σ_{5 i}^{-} e^{- \frac{n_{1}}{μ_{0}} τ} + σ_{6 i}^{-} e^{- (\frac{n_{1}}{μ_{0}} + a_{i}) τ} + σ_{7 i}^{-} τ e^{- \frac{n_{1}}{μ_{0}} τ} \end{matrix}

where

X_{i} = \frac{χ_{1 i}^{-} - P_{1 i}^{-}}{2 k}

,

Y_{i} = \frac{χ_{2 i}^{-} - P_{2 i}^{-}}{2 k}

,

D_{1 i}^{\pm} = α^{\pm} C_{1 i} \pm X_{i}

,

D_{2 i}^{\pm} = α^{\mp} C_{2 i} \mp Y_{i}

and

σ_{j i}^{\pm} = \frac{1}{2} (P_{j i}^{+} \pm P_{j i}^{-}) (j = 1, 2, 3, 4, 5, 6, 7)

. By substituting Eqs. (112) and (113) to the Eq. (105) which refers to the boundary conditions, we get C_1i and C_2i

\begin{matrix} C_{1 i} = \frac{- α^{-} (X_{i} e^{- k τ_{0}} - Y_{i} - σ_{3 i}^{-} e^{- k τ_{0}} - σ_{4 i}^{-} - σ_{5 i}^{-} - σ_{6 i}^{-}) e^{- k τ_{0}}}{{(α^{+})}^{2} - {(α^{-})}^{2} e^{- 2 k τ_{0}}} \\ + \frac{α^{+} (- X_{i} + Y_{i} e^{- k τ_{0}} - σ_{1 i}^{+} τ_{0} - σ_{2 i}^{+} τ_{0} e^{- k τ_{0}} - σ_{3 i}^{+} e^{- a_{i} τ_{0}})}{{(α^{+})}^{2} - {(α^{-})}^{2} e^{- 2 k τ_{0}}} \\ - \frac{α^{+} (σ_{4 i}^{+} e^{- (k + a_{i}) τ_{0}} + σ_{5 i}^{+} e^{- \frac{n_{1}}{μ_{0}} τ_{0}} + σ_{6 i}^{+} e^{- (\frac{n_{1}}{μ_{0}} + a_{i}) τ_{0}} + σ_{7 i}^{+} τ_{0} e^{- \frac{n_{1}}{μ_{0}} τ_{0}})}{{(α^{+})}^{2} - {(α^{-})}^{2} e^{- 2 k τ_{0}}} \end{matrix}

C_{2 i} = \frac{X_{i} e^{- k τ_{0}} - Y_{i} - σ_{3 i}^{-} e^{- k τ_{0}} - σ_{4 i}^{-} - σ_{5 i}^{-} - σ_{6 i}^{-} - α^{-} C_{1 i} e^{- k τ_{0}}}{α^{+}}

Finally we obtain the solution of Eqs. (21) - (23) in the manuscript by using the Eqs. (88) and (89). In addition, we also propose an inhomogeneous solution in which linear expressions are used to represent the vertical variation of the single scattering albedo and the asymmetry factor. The detailed process has been shown in [60].

B. Solution of Equations (29) - (31) and (32) - (34) in the manuscript

Equations (29) - (31) in the manuscript are

\frac{d {\bar{F}}_{+}^{sca} (τ)}{d τ} = γ_{1} (τ) {\bar{F}}_{+}^{sca} (τ) - γ_{2} (τ) {\bar{F}}_{-}^{sca} (τ)

\frac{d {\bar{F}}_{-}^{sca} (τ)}{d τ} = γ_{2} (τ) {\bar{F}}_{+}^{sca} (τ) - γ_{1} (τ) {\bar{F}}_{-}^{sca} (τ)

{\bar{F}}_{-}^{sca} (0) = 1, {\bar{F}}_{+}^{sca} (τ_{0}) = 0

By applying the similar procedure as Eqs. (85) and (86) to Eqs. (94) and (95), Eqs. (116) and (117) yields

\begin{matrix} \frac{d {\bar{F}}_{+}^{s c a} (τ)}{d τ} = [γ_{10} + γ_{11} ε_{ω} (e^{- a_{1} τ} - e^{- a_{1} τ_{0} / 2}) + γ_{12} ε_{g}^{*} (e^{- a_{2} τ} - e^{- a_{2} τ_{0} / 2})] \\ \times [{\bar{F}}_{+}^{sca, 0} (τ) + ε_{ω} {\bar{F}}_{+}^{sca, 1} (τ) + ε_{g}^{*} {\bar{F}}_{+}^{sca, 2} (τ)] - [γ_{20} + γ_{21} ε_{ω} (e^{- a_{1} τ} - e^{- a_{1} τ_{0} / 2}) \\ + γ_{22} ε_{g}^{*} (e^{- a_{2} τ} - e^{- a_{2} τ_{0} / 2})] [{\bar{F}}_{-}^{sca, 0} (τ) + ε_{ω} {\bar{F}}_{-}^{sca, 1} (τ) + ε_{g}^{*} {\bar{F}}_{-}^{sca, 2} (τ)] \end{matrix}

\begin{matrix} \frac{d {\bar{F}}_{-}^{s c a} (τ)}{d τ} = [γ_{20} + γ_{21} ε_{ω} (e^{- a_{1} τ} - e^{- a_{1} τ_{0} / 2}) + γ_{22} ε_{g}^{*} (e^{- a_{2} τ} - e^{- a_{2} τ_{0} / 2})] \\ \times [{\bar{F}}_{+}^{sca, 0} (τ) + ε_{ω} {\bar{F}}_{+}^{sca, 1} (τ) + ε_{g}^{*} {\bar{F}}_{+}^{sca, 2} (τ)] - [γ_{10} + γ_{11} ε_{ω} (e^{- a_{1} τ} - e^{- a_{1} τ_{0} / 2}) \\ + γ_{12} ε_{g}^{*} (e^{- a_{2} τ} - e^{- a_{2} τ_{0} / 2})] [{\bar{F}}_{-}^{sca, 0} (τ) + ε_{ω} {\bar{F}}_{-}^{sca, 1} (τ) + ε_{g}^{*} {\bar{F}}_{-}^{sca, 2} (τ)] \end{matrix}

By ignoring terms with second or higher order of ε_ω, ε_g, we separate Eqs. (119) and (120) to three Eqs. involving ${\bar{F}}_{\pm}^{sca, 0} (τ)$ , ${\bar{F}}_{\pm}^{sca, 1} (τ)$ and ${\bar{F}}_{\pm}^{sca, 2} (τ)$ .

The Eqs. for ${\bar{F}}_{\pm}^{sca, 0} (τ)$ are

\frac{d {\bar{F}}_{+}^{sca, 0} (τ)}{d τ} = γ_{10} {\bar{F}}_{+}^{sca, 0} (τ) - γ_{20} {\bar{F}}_{-}^{sca, 0} (τ)

\frac{d {\bar{F}}_{-}^{sca, 0} (τ)}{d τ} = γ_{20} {\bar{F}}_{+}^{sca, 0} (τ) - γ_{10} {\bar{F}}_{-}^{sca, 0} (τ)

{\bar{F}}_{-}^{sca, 0} (0) = 1, {\bar{F}}_{+}^{sca, 0} (τ_{0}) = 0

Similar to the solution of Eqs. (96) - (98), the solutions of Eqs. (121) and (122) are

{\bar{F}}_{+}^{sca, 0} (τ) = α^{+} {\bar{C}}_{10} e^{k (τ - τ_{0})} + α^{-} {\bar{C}}_{20} e^{- k τ}

{\bar{F}}_{-}^{sca, 0} (τ) = α^{-} {\bar{C}}_{10} e^{k (τ - τ_{0})} + α^{+} {\bar{C}}_{20} e^{- k τ}

After substituting Eqs. (124) and (125) to Eq. (123) which refers to the boundary conditions, we obtain

{\bar{C}}_{10} = \frac{- α^{-} e^{- k τ_{0}}}{{(α^{+})}^{2} - {(α^{-})}^{2} e^{- 2 k τ_{0}}}

{\bar{C}}_{20} = \frac{1 - α^{-} {\bar{C}}_{10} e^{- k τ_{0}}}{α^{+}}

The Eqs. of ${\bar{F}}_{\pm}^{sca, i} (τ)$ ( $i = 1, 2$ ) are same as Eqs. (103) - (105) but without terms containing $e^{- \frac{n_{1}}{μ_{0}} τ}$ . Therefore, the expression of ${\bar{F}}_{\pm}^{sca, i} (τ)$ is same as Eqs. (112) and (113) but: 1. The expressions of ${\bar{F}}_{\pm}^{sca, i} (τ)$ are exclusive of terms containing $e^{- \frac{n_{1}}{μ_{0}} τ}$ ; 2. All parameters (e.g. χ, ϕ, P) in the expressions of ${\bar{F}}_{\pm}^{sca, i} (τ)$ are based on ${\bar{C}}_{10}$ and ${\bar{C}}_{20}$ instead of C₁₀ and C₂₀.

The solutions of Eqs. (32) - (34) in the manuscript are same as the solution of Eqs. (29) - (31) in the manuscript but all parameters are based on ${\bar{C}}_{10}^{*}$ and ${\bar{C}}_{20}^{*}$ . The expressions of ${\bar{C}}_{10}^{*}$ and ${\bar{C}}_{20}^{*}$ are

{\bar{C}}_{10}^{*} = \frac{α^{+}}{{(α^{+})}^{2} - {(α^{-})}^{2} e^{- 2 k τ_{0}}}

{\bar{C}}_{20}^{*} = \frac{- α^{-} {\bar{C}}_{10}^{*} e^{- k τ_{0}}}{α^{+}}

Funding

National Key R & D Program of China (2018YFC1507002); National Natural Science Foundation of China (41675003, 41675056); Postgraduate Research & Practice innovation Program of Jiangsu Province (KYCX18_1001).

Acknowledgments

We appreciate 3 anonymous reviewers very much for their constructive comments.

References

1. X. Guo, D. Fu, X. Li, Z. Hu, H. Lei, H. Xiao, and Y. Hong, “Advances in cloud physics and weather modification in China,” Adv. Atmospheric Sci. 32, 230–249 (2015). [CrossRef]

2. J. Li, W.-C. Wang, X. Dong, and J. Mao, “Cloud-radiation-precipitation associations over the Asian monsoon region: an observational analysis,” Clim. Dyn. 49, 3237–3255 (2017). [CrossRef]

3. Y. Yan, Y. Liu, and J. Lu, “Cloud vertical structure, precipitation, and cloud radiative effects over Tibetan Plateau and its neighboring regions,” J. Geophys. Res. Atmospheres 121, 5864–5877 (2016). [CrossRef]

4. F. Zhang, J.-R. Yan, J. Li, K. Wu, H. Iwabuchi, and Y.-N. Shi, “A new radiative transfer method for solar radiation in a vertically internally inhomogeneous medium,” J. Atmospheric Sci. 75, 41–55 (2018). [CrossRef]

5. K. Chan, M. Wiegner, M. Wenig, and D. Pöhler, “Observations of tropospheric aerosols and NO₂ in Hong Kong over 5 years using ground based MAX-DOAS,” Sci. The Total. Environ. 619-620, 1545–1556 (2018). [CrossRef]

6. K. L. Chan, M. Wiegner, H. Flentje, I. Mattis, F. Wagner, J. Gasteiger, and A. Geiß, “Evaluation of ECMWF-IFS (version 41R1) operational model forecasts of aerosol transport by using ceilometer network measurements,” Geosci. Model. Dev. 11, 3807–3831 (2018). [CrossRef]

7. Y. Sun, C. Liu, K. Chan, W. Wang, C. Shan, Q. Hu, and J. Liu, “The influence of instrumental line shape degradation on the partial columns of O₃, CO, CH₄ and N₂O derived from high-resolution FTIR spectrometry,” Remote. Sens. 10, 2041 (2018). [CrossRef]

8. M. J. Iacono, J. S. Delamere, E. J. Mlawer, M. W. Shephard, S. A. Clough, and W. D. Collins, “Radiative forcing by long-lived greenhouse gases: calculations with the AER radiative transfer models,” J. Geophys. Res. Atmospheres 113, D13 (2008). [CrossRef]

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

10. Q. Fu and K. N. Liou, “Parameterization of the radiative properties of cirrus clouds,” J. Atmospheric Sci. 50, 2008–2025 (1993). [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. Atmospheres 108, D20 (2003). [CrossRef]

12. H. Zhang, G. Shi, T. Nakajima, and T. Suzuki, “The effects of the choice of the k-interval number on radiative calculations,” J. Quant. Spectrosc. Radiat. Transf. 98, 31–43 (2006). [CrossRef]

13. H. Zhang, T. Suzuki, T. Nakajima, G. Shi, X. Zhang, and Y. Liu, “Effects of band division on radiative calculations,” Opt. Eng. 45, 016002 (2006). [CrossRef]

14. J. Li and H. W. Barker, “A radiation algorithm with correlated-k distribution. part i: local thermal equilibrium,” J. Atmospheric Sci. 62, 286–309 (2005). [CrossRef]

15. J. Li, S. Dobbie, P. Räisänen, and Q. Min, “Accounting for unresolved clouds in a 1-D solar radiative-transfer model,” Q. J. Royal Meteorol. Soc. 131, 1607–1629 (2005). [CrossRef]

16. V. Sobolev, “Chapter 7 - linear integral equations for the reflection and transmission coefficients,” in Light Scattering in Planetary Atmospheres, 76 (Pergamon, 1975), pp. 126 – 152. [CrossRef]

17. J. Lenoble, Atmospheric radiative transfer (A Deepak Publishing, 1993).

18. K. Y. Kondratyev, V. V. Kozoderov, and O. I. Smokty, “Theoretical models of the non-polarized optical radiation in the atmosphere-surface system,” in Remote Sensing of the Earth from Space: Atmospheric Correction, (Springer, 1992), pp. 69–140.

19. K. Liou, Radiation and Cloud Processes in the Atmosphere: Theory, Observation and Modeling (Oxford University, 1992).

20. R. Guzzi, “Radiative transfer, solution techniques,” in Encyclopedia of Remote Sensing, E. G. Njoku, ed. (Springer, 2014). [CrossRef]

21. M. Collier, S. Jeffrey, L. Rotstayn, K.-H. Wong, S. Dravitzki, C. Moseneder, C. Hamalainen, J. Syktus, R. Suppiah, J. Antony, A. ElZein, and M. Atif, “The CSIRO-Mk3. 6.0 Atmosphere-Ocean GCM: participation in CMIP5 and data publication,”in International Congress on Modelling and Simulation–MODSIM, (2011).

22. R. B. Neale, A. Gettelman, S. Park, A. J. Conley, D. Kinnison, D. Marsh, A. K. Smith, F. Vitt, H. Morrison, P. Cameron-smith, W. D. Collins, M. J. Iacono, R. C. Easter, X. Liu, M. A. Taylor, C. chieh Chen, P. H. Lauritzen, D. L. Williamson, R. Garcia, J. francois Lamarque, M. Mills, S. Tilmes, S. J. Ghan, P. J. Rasch, and M. Meteorology, “Description of the NCAR Community Atmosphere Model (CAM 5.0),” Tech. rep., NCAR Tech. (2010).

23. K. von Salzen, J. F. Scinocca, N. A. McFarlane, J. Li, J. N. S. Cole, D. Plummer, D. Verseghy, M. C. Reader, X. Ma, M. Lazare, and L. Solheim, “The Canadian Fourth Generation Atmospheric Global Climate Model (CanAM4). part i: representation of physical processes,” Atmosphere-Ocean 51, 104–125 (2013). [CrossRef]

24. V. Brovkin, L. Boysen, T. Raddatz, V. Gayler, A. Loew, and M. Claussen, “Evaluation of vegetation cover and land-surface albedo in MPI-ESM CMIP5 simulations,” J. Adv. Model. Earth Syst. 5, 48–57 (2013). [CrossRef]

25. R. Chen, R. Wood, Z. Li, R. Ferraro, and F.-L. Chang, “Studying the vertical variation of cloud droplet effective radius using ship and space-borne remote sensing data,” J. Geophys. Res. Atmospheres 113, D8 (2008). [CrossRef]

26. V. R. Noonkester, “Droplet spectra observed in marine stratus cloud layers,” J. Atmospheric Sci. 41, 829–845 (1984). [CrossRef]

27. K. Chan and K. Chan, “Aerosol optical depths and their contributing sources in Taiwan,” Atmospheric Environ. 148, 364 – 375 (2017). [CrossRef]

28. C. Xing, C. Liu, S. Wang, K. L. Chan, Y. Gao, X. Huang, W. Su, C. Zhang, Y. Dong, G. Fan, T. Zhang, Z. Chen, Q. Hu, H. Su, Z. Xie, and J. Liu, “Observations of the vertical distributions of summertime atmospheric pollutants and the corresponding ozone production in Shanghai, China,” Atmospheric Chem. Phys. 17, 14275–14289 (2017). [CrossRef]

29. Y. Zhu, K. L. Chan, Y. F. Lam, M. Horbanski, D. Pöhler, J. Boll, I. Lipkowitsch, S. Ye, and M. Wenig, “Analysis of spatial and temporal patterns of on-road NO₂ concentrations in Hong Kong,” Atmospheric Meas. Tech. 11, 6719–6734 (2018). [CrossRef]

30. M. Saito and H. Iwabuchi, “A new method of measuring aerosol optical properties from digital twilight photographs,” Atmospheric Meas. Tech. 8, 4295–4311 (2015). [CrossRef]

31. M. Saito and H. Iwabuchi, “Cloud discrimination from sky images using a clear-sky index,” J. Atmospheric Ocean. Technol. 33, 1583–1595 (2016). [CrossRef]

32. T. Hayasaka, S. Satake, A. Shimizu, N. Sugimoto, I. Matsui, K. Aoki, and Y. Muraji, “Vertical distribution and optical properties of aerosols observed over Japan during the Atmospheric Brown Clouds-East Asia Regional Experiment 2005,” J. Geophys. Res. Atmospheres 112, D22 (2007). [CrossRef]

33. G. M. Martin, D. W. Johnson, and A. Spice, “The measurement and parameterization of effective radius of droplets in warm stratocumulus clouds,” J. Atmospheric Sci. 51, 1823–1842 (1994). [CrossRef]

34. R. Wood, “Parametrization of the effect of drizzle upon the droplet effective radius in stratocumulus clouds,” Q. J. Royal Meteorol. Soc. 126, 3309–3324 (2000). [CrossRef]

35. N. L. Miles, J. Verlinde, and E. E. Clothiaux, “Cloud droplet size distributions in low-level stratiform clouds,” J. Atmospheric Sci. 57, 295–311 (2000). [CrossRef]

36. F.-L. Chang and Z. Li, “Estimating the vertical variation of cloud droplet effective radius using multispectral near-infrared satellite measurements,” J. Geophys. Res. Atmospheres 107, D15 (2002). [CrossRef]

37. M. Duan, Q. Min, and K. Stamnes, “Impact of vertical stratification of inherent optical properties on radiative transfer in a plane-parallel turbid medium,” Opt. Express 18, 5629–5638 (2010). [CrossRef] [PubMed]

38. J. Li, D. J. W. Geldart, and P. Chỳlek, “Solar radiative transfer in clouds with vertical internal inhomogeneity,” J. Atmospheric Sci. 51, 2542–2552 (1994). [CrossRef]

39. J. H. Joseph, W. J. Wiscombe, and J. A. Weinman, “The Delta-Eddington approximation for radiative flux transfer,” J. Atmospheric Sci. 33, 2452–2459 (1976). [CrossRef]

40. Q. Fu, K. N. Liou, M. C. Cribb, T. P. Charlock, and A. Grossman, “Multiple scattering parameterization in thermal infrared radiative transfer,” J. Atmospheric Sci. 54, 2799–2812 (1997). [CrossRef]

41. W. E. Meador and W. R. Weaver, “Two-stream approximations to radiative transfer in planetary atmospheres: a unified description of existing methods and a new improvement,” J. Atmospheric Sci. 37, 630–643 (1980). [CrossRef]

42. K. Stamnes, G. E. Thomas, and J. J. Stamnes, Radiative Transfer in the Atmosphere and Ocean (Cambridge University, 2017). [CrossRef]

43. P. Räisänen, “Two-stream approximations revisited: A new improvement and tests with GCM data,” Q. J. Royal Meteorol. Soc. 128, 2397–2416 (2002). [CrossRef]

44. B. Ritter and J.-F. Geleyn, “A comprehensive radiation scheme for numerical weather prediction models with potential applications in climate simulations,” Mon. Weather. Rev. 120, 303–325 (1992). [CrossRef]

45. G. L. Stephens, S.-C. Tsay, P. W. Stackhouse, and P. J. Flatau, “The relevance of the microphysical and radiative properties of cirrus clouds to climate and climatic feedback,” J. Atmospheric Sci. 47, 1742–1754 (1990). [CrossRef]

46. H. Zhang, F. Zhang, Q. Fu, Z. Shen, and P. Lu, “Two- and four-stream combination approximations for computation of diffuse actinic fluxes,” J. Atmospheric Sci. 67, 3238–3252 (2010). [CrossRef]

47. F. Zhang, Y.-N. Shi, J. Li, K. Wu, and H. Iwabuchi, “Variational iteration method for infrared radiative transfer in a scattering medium,” J. Atmospheric Sci. 74, 419–430 (2017). [CrossRef]

48. F. Zhang, J. Li, H. Zhang, L. Ma, and X. Zhou, “On the relationship between direct and diffuse radiation,” J. Quant. Spectrosc. Radiat. Transf. 115, 60 – 65 (2013). [CrossRef]

49. F. Zhang, Z. Shen, J. Li, X. Zhou, and L. Ma, “Analytical delta-four-stream doubling-adding method for radiative transfer parameterizations,” J. Atmospheric Sci. 70, 794–808 (2013). [CrossRef]

50. F. Zhang and J. Li, “Doubling-adding method for delta-four-stream spherical harmonic expansion approximation in radiative transfer parameterization,” J. Atmospheric Sci. 70, 3084–3101 (2013). [CrossRef]

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

52. C. A. Randles, S. Kinne, G. Myhre, M. Schulz, P. Stier, J. Fischer, L. Doppler, E. Highwood, C. Ryder, B. Harris, J. Huttunen, Y. Ma, R. T. Pinker, B. Mayer, D. Neubauer, R. Hitzenberger, L. Oreopoulos, D. Lee, G. Pitari, G. Di Genova, J. Quaas, F. G. Rose, S. Kato, S. T. Rumbold, I. Vardavas, N. Hatzianastassiou, C. Matsoukas, H. Yu, F. Zhang, H. Zhang, and P. Lu, “Intercomparison of shortwave radiative transfer schemes in global aerosol modeling: results from the AeroCom Radiative Transfer Experiment,” Atmospheric Chem. Phys. 13, 2347–2379 (2013). [CrossRef]

53. M. D. King and Harshvardhan, “Comparative accuracy of selected multiple scattering approximations,” J. Atmospheric Sci. 43, 784–801 (1986). [CrossRef]

54. H. W. Barker, J. N. S. Cole, J. Li, B. Yi, and P. Yang, “Estimation of errors in two-stream approximations of the solar radiative transfer equation for cloudy-sky conditions,” J. Atmospheric Sci. 72, 4053–4074 (2015). [CrossRef]

55. R. McClatchey, R. Fenn, J. A. Selby, F. Volz, and J. Garing, Optical Properties of the Atmosphere (3 Edition) (Air Force Cambridge Research Laboratories, 1972).

56. J. S. Dobbie, J. Li, and P. Chỳlek, “Two- and four-stream optical properties for water clouds and solar wavelengths,” J. Geophys. Res. Atmospheres 104, 2067–2079 (1999). [CrossRef]

57. E. K. Dolinar, X. Dong, B. Xi, J. H. Jiang, and H. Su, “Evaluation of CMIP5 simulated clouds and TOA radiation budgets using NASA satellite observations,” Clim. Dyn. 44, 2229–2247 (2015). [CrossRef]

58. R. Bennartz, “Global assessment of marine boundary layer cloud droplet number concentration from satellite,” J. Geophys. Res. Atmospheres 112, D2 (2007).

59. T. Kato, Perturbation Theory for Linear Operators (Springer, 1995). [CrossRef]

60. Y.-N. Shi, F. Zhang, K. L. Chan, W. Li, H. Lin, and M. Duan, “An improved Eddington approximation method for irradiance calculation in a vertical inhomogeneous medium,” J. Quant. Spectrosc. Radiat. Transf. 226, 40 – 50 (2019). [CrossRef]

	Homogeneous	Inhomogeneous	Benchmark
Atmospheric Profile	The midlatitude winter atmosphere [55] is used with a subdivision of 400 layers each of which having a geometrical thickness of 0.25 km.
Ground albedo	0.1
Solar zenith angle	1.0
Cloud location	1 - 2 km
The vertical resolution of the cloudy area	4 layers with geometrical thickness of 0.25 km		100 layers with geometrical thickness of 0.01 km
LWC within the cloud ( $g m^{- 3}$ )	Case 1: 0.276,0.269, 0.261,0.254 Case 2: 0.276,0.266, 0.258,0.252	Case 1: $0.25 + 0.00003 z$ Case 2: $0.25 + 0.000001 z^{1.5}$

	Homogeneous	Benchmark	Inhomogeneous
Algorithm only	1.00	5.89	2.89
Radiative transfer model	1.00	4.25	2.01

	Homogeneous	Inhomogeneous	Benchmark
Atmospheric Profile	The midlatitude winter atmosphere [55] is used with a subdivision of 400 layers each of which having a geometrical thickness of 0.25 km.
Ground albedo	0.1
Solar zenith angle	1.0
Cloud location	1 - 2 km
The vertical resolution of the cloudy area	4 layers with geometrical thickness of 0.25 km		100 layers with geometrical thickness of 0.01 km
LWC within the cloud ( $g m^{- 3}$ )	Case 1: 0.276,0.269, 0.261,0.254 Case 2: 0.276,0.266, 0.258,0.252	Case 1: $0.25 + 0.00003 z$ Case 2: $0.25 + 0.000001 z^{1.5}$

	Homogeneous	Benchmark	Inhomogeneous
Algorithm only	1.00	5.89	2.89
Radiative transfer model	1.00	4.25	2.01

Multi-layer solar radiative transfer considering the vertical variation of inherent microphysical properties of clouds

Abstract

1. Introduction

2. Single-layer RT solution

2.1. Basic RT equation

2.2. RT solution of direct sunlight

2.3. RT solution of scattered sunlight

2.3.1. Direct incident radiation

2.3.2. Diffuse incident radiation

3. Multi-layer RT calculation

3.1. Combination of two layers

3.1.1. Direct incident radiation

3.1.2. Diffuse incident radiation

3.2. Multiple layers

4. Comparisons and discussions

4.1. Combinations of two layers

4.2. Multilayer atmosphere

5. Summary and conclusions

Appendix

A. Solution of Equations (21) - (23) in the manuscript

B. Solution of Equations (29) - (31) and (32) - (34) in the manuscript

Funding

Acknowledgments

References

Cited By

Figures (6)

Tables (2)

Equations (129)

Optics Express