Efficient radiative transfer model for thermal infrared brightness temperature simulation in cloudy atmospheres

Wenwen Li; Feng Zhang; Feng Zhang; Yi-Ning Shi; Hironobu Iwabuchi; Mingwei Zhu; Jiangnan Li; Wei Han; Husi Letu; Hiroshi Ishimoto

doi:10.1364/OE.400130

1. Introduction

As a fundamental ingredient of the earth-atmosphere system, clouds affect short-range weather processes, climate change, and atmospheric circulation by regulating the global radiation budget and influencing latent heating and cooling caused by phase transformation [1–3]. However, the interaction between cloud and radiation is complicated by the macro and micro properties of different types of clouds and their complex temporal and spatial variations. Satellite observation provides an effective way to monitor the earth-atmospheric system over regional and global scales in high spatiotemporal resolution. The cloud properties such as cloud top height, cloud effective radius, cloud optical thickness, and cloud water path can be retrieved from the visible to the infrared observation of satellite by various algorithms [4–9]. A forward radiative transfer model with high accuracy and computational efficiency is a critical component in the process of the cloud inversion, and it can also be used to calibrate satellite instrument and evaluate numerical model [10–15].

The solution of the radiative transfer equation is a crucial issue when handling the radiative processes in the radiative transfer forward model. Since the radiative transfer equation is an integro-differential equation, there is no exact solution in the realistic atmosphere. Many accurate radiative transfer schemes, such as the discrete ordinate radiative transfer (DISORT) method [16], adding-doubling algorithm [17,18], and Monte-Carlo technique [19,20] have been proposed and widely used, but these accurate schemes are limited in remote sensing applications due to the highly time-consuming process of single simulation and inefficient processing of high spectral and spatial satellite dataset. Since the scattering is much weaker in infrared radiative transfer than in shortwave radiative transfer, numerous rapid infrared radiative transfer schemes are developed by simplifying scattering or reducing streams [21–23]. Among these methods, an adding method of infrared four-stream discrete ordinates method (DOM) (4DDA) is designed for the flux calculation of infrared radiative transfer in the climate model [24]. In this method, the single-layer solutions for the infrared radiative transfer are solved by using the four-stream DOM. Then the radiative transfer in multiple layers with a vertically inhomogeneous atmosphere is settled by the analytical adding method based on the invariance principle. Infrared radiance measurements from satellites are the reliable source of earth-atmosphere observation on a near global scale. Since the brightness temperature (BT) is determined by radiative intensity instead of radiative flux and this physical quantity is essential to the satellite remote sensing. Therefore, in this study, 4DDA is extended to be able to calculate the thermal infrared radiative intensity varying with the zenith angle for remote sensing applications.

Furthermore, atmospheric gaseous absorption has a significant effect on thermal radiation. Absorption lines of gases are irregularly distributed across the spectrum, and their intensities and shapes strongly depend on pressure and temperature [25]. Due to the complexity of the gaseous absorption lines, a detailed line-by-line radiative transfer model (LBLRTM), is proposed to settle the inhomogeneous atmospheric radiation faultlessly [26]. However, LBLRTM has a far higher computational cost, that can not be afforded by remote sensing applications. Therefore, several approaches called channel-based forward models such as correlated k-distribution (CKD) [27], spectral sampling (OSS) [28], and radiance sampling method (RSM) [29] are proposed to reduce the computational time by using nonlinear functions of gas and temperature profiles, or reducing monochromatic calculations [30]. These approaches are beneficial for the hyperspectral satellite data processing and climate modeling community. J Li and Barker [31] proposed the alternate mapping correlated k-distribution (AMCKD) method, which is superior in handling the overlap of gaseous absorption lines. Most recently, AMCKD has been extended to the application of satellite simulation under clear sky condition and achieved high accuracy and efficiency [32]. The AMCKD method is employed in this study to improve computational efficiency for the satellite simulation under cloudy atmospheres.

This study develops an efficient radiative transfer model (ERTM) to simulate thermal infrared (TIR) brightness temperatures observed by the Advanced Himawari Imager (AHI) onboard the Himawari-8 satellites for cloudy atmospheres. In section 2, the ERTM is introduced. The AMCKD is applied to settle the gaseous absorption; a parameterization of cloud optical property is established for both ice and water clouds; the rapid radiative transfer scheme (4DDA) is extended to the satellite simulation, and incorporated into ERTM. Then the ERTM is integrated by the three parts with the requisite input meteorological dataset. section 3 validates the ERTM by comparing the simulated results with rigorous results from the benchmark model and AHI observations. Finally, a summary is given in section 4.

2. Efficient radiative transfer model (ERTM)

In this section, the ERTM that simulates thermal infrared brightness temperatures observed by AHI for cloudy atmospheres is constructed. AHI onboard the Himawari-8 satellites can provide high spectral, spatial, and temporal resolution data for multispectral bands, including six visible & near-infrared bands and ten infrared bands [33]. Table 1 shows the parameters of each AHI TIR channel, including center wavelength, bandwidth, and considered absorbers. For each channel, up to three kinds of gases are considered.

Table 1. Specification of the AHI TIR channels.

View Table | View all tables in this article

2.1 Gaseous absorption by AMCKD

AMCKD is employed to calculate the gaseous absorption for radiative transfer processes, and integrated into ERTM seamlessly. In AMCKD, the band-averaged gaseous transmittance in a homogenous layer of the atmosphere with a single gas is defined as

(1)$$T{r_\phi }(u) = \int_{\Delta v} {\phi (v){e^{ - uk(v)}}dv} \textrm{/}\int_{\Delta v} {\phi (v)dv}, $$

where $k(v)$ is the gaseous absorption coefficient at wavenumber v. $\phi (v)$ is the spectral response function of the measurement instrument, and u denotes the gas amount. The Eq. (1) can also be written in k space [34] as

(2)$$T{r_\phi }(u) = \int_\textrm{0}^\infty {{e^{ - uk}}{f_\phi }(k)dk}, $$

where ${f_\phi }(k)$ is the normalized probability distribution function for $k(v)$ considering the spectral response function.

By defining the cumulative probability function $\hat{g}(k) = \int_0^k {{f_\phi }} (k^{\prime})dk^{\prime}$, the integration over k in Eq. (2) can be replaced by an integration over $\hat{g}$. Dividing $\hat{g}$ space into ${N_{\hat{g}}}$ points ${\hat{g}_i}$ (i = 0, 1, 2, …, ${N_{\hat{g}}}$), with ${\hat{g}_0} = 0$, ${\hat{g}_{{N_{\hat{g}}}}} = 1$, and letting $\Delta {\hat{g}_i} = {\hat{g}_i} - {\hat{g}_{i - 1}}$, with $\sum\limits_{i = 1}^{{N_g}} {\Delta {{\hat{g}}_i}} = 1$, the integral in (2) becomes

(3)$$T{r_\phi }(u) = \int_0^1 {{e^{ - uk(\hat{g})}}d\hat{g}} \textrm{ = }\sum\limits_{i = 1}^{{N_{\hat{g}}}} {\int_{{{\hat{g}}_{i - 1}}}^{{{\hat{g}}_i}} {{e^{ - uk(\hat{g})}}} d\hat{g}}. $$

The transmission function in AMCKD is calculated in an approximate format as

(4)$$T{r_\phi }(u) = \sum\limits_{i = 1}^{{N_{\hat{g}}}} {{e^{ - u{{\bar{k}}_i}}}\Delta {{\hat{g}}_i}}, $$

where the absorption coefficient ${\overline k _i}{\kern 1pt}$in i-th subinterval can be obtained through fitting to the LBLRTM results of Eq. (3) in the same domain for a suitable range of u.

As the inequality relation is given by

(5)$$\int_{{{\hat{g}}_{i - 1}}}^{{{\hat{g}}_i}} {{e^{ - uk(\hat{g})}}} d\hat{g}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \ge {e^{ - \frac{u}{{\Delta {{\hat{g}}_i}}}\int_{{{\hat{g}}_{i - 1}}}^{{{\hat{g}}_i}} {k(\hat{g})} d\hat{g}}}\Delta {\hat{g}_i}{\kern 1pt} {\kern 1pt} {\kern 1pt}, $$

which indicates that ${\overline k _i}\textrm{ = }\frac{{{\alpha _i}}}{{\Delta {{\hat{g}}_i}}}\int_{{{\hat{g}}_{i - 1}}}^{{{\hat{g}}_i}} {k(\hat{g})} d\hat{g}{\kern 1pt}$, and ${\alpha _i} < 1$. The adjusted factor ${\alpha _i}$ is pre-calculated through the comparison of transmittance with LBLRTM results [31,32]. The calculation of ${\alpha _i}$ represents the lion’s share of the work of building the AMCKD. In the realistic inhomogeneous atmosphere at various temperature T and pressure ${P_{atm}}$, ${\bar{k}_i}$ becomes ${\bar{k}_i}({P_{atm}},T)$, which can be parameterized as a polynomial in temperature for certain reference pressures. Then ${\bar{k}_i}({P_{atm}},T)$ at any arbitrary pressure can be acquired by the linear interpolation between two neighboring reference pressures.

The alternate mapping method is used in the AMCKD to handle the overlap of gaseous absorption lines. Figure 1 illustrates the treatments of overlapping absorption lines for the TIR B10 (7.35 µm) channel. Three absorbing gases considered for the B10 channel are H₂O, N₂O, and CH₄. In AMCKD, a primary gas is chosen and sorted in each subinterval, and the same sorting rules are applied to the other gases in the same subinterval. Then the optical thickness of gaseous absorption for the i-th subinterval $\Delta {\hat{g}_i}$ at an arbitrary atmospheric layer is as follows:

(6)$${\tau _{gas\_i}} = {u_{{H_2}O}} \cdot {\bar{k}_{{H_2}O\_i}}({P_{atm}},T)\textrm{ + }{u_{{N_2}O}} \cdot {\bar{k}_{{N_2}O\_i}}({P_{atm}},T) + {u_{C{H_4}}} \cdot {\bar{k}_{C{H_4}\_i}}({P_{atm}},T). $$

A detailed description of AMCKD can be found in [31,32]. Besides, the accuracy and computational efficiency of AMCKD in applications for AHI are evaluated by [32], the absolute errors of BTs simulated by AMCKD for all TIR channels of AHI under the clear sky are bounded by 0.44 K, compared to the rigorous results of LBLRTM based on the standard atmospheric profile. Moreover, AMCKD has a comparable level of computational efficiency to Radiative Transfer for the Television Observation Satellite Operational Vertical Sounder (RTTOV). Therefore, AMCKD can be effectively applied for the satellite simulation.

Fig. 1. Example of the treatment of overlapping absorption lines based on AMCKD for the TIR B10 (7.35 µm) channel. (a) Gaseous absorption coefficients of H₂O, N₂O, and CH₄ as a function of wavenumber. (b) Gaseous absorption coefficients as a function of cumulative probability sorted by AMCKD, considering the spectral response function.

Download Full Size | PDF

2.2 Cloud optical property parameterization

The scattering and absorption of clouds need to be taken into consideration in the cloudy radiative transfer process. For the sake of convenience, a parameterization for cloud optical property is established. A highly irregular ice particle model called Voronoi aggregate developed by [35] and [36] is applied. The size distribution of ice particles is assumed to be the gamma distribution. A number of effective lengths ${L_{eff}}$ (given by Eq. (7)) ranging from 4 to 200 µm are calculated for 50 different size distributions,

(7)$${L_{eff}} = \frac{3}{2}\frac{{\int_{{L_{\min }}}^{{L_{\max }}} {V(L )n(L )dL} }}{{\int_{{L_{\min }}}^{{L_{\max }}} {S(L )n(L )dL} }}, $$

where L represents the maximum dimension of a single ice crystal; $V(L )$ and $S(L )$ are the volume and projected area of a single ice crystal, respectively. $n(L )$ is the size distribution of ice particles. Then the optical properties of ice clouds, including asymmetry factor, extinction coefficient, and single-scattering albedo, are calculated for the same 50 size distributions at the given 99 wavelengths in the Voronoi aggregate dataset [37].

Since the optical properties vary with the spectral wavelength and effective length of ice clouds, the optical properties are parameterized by the effective length for each wavelength. For a particular wavelength $\lambda$, the fitting polynomials of the optical properties are defined as

(8)$${g_{ice}}({\lambda ,{L_{eff}}} )= \sum\limits_{i = 0}^{10} {{a_{ice\_i}}(\lambda )L_{eff}^{5 - i}}, $$

(9)$${\beta _{e\_ice}}({\lambda ,{L_{eff}}} )= \sum\limits_{i = 0}^{10} {{b_{ice\_i}}(\lambda )L_{eff}^{5 - i}}, $$

(10)$$1 - {\omega _{ice}}({\lambda ,{L_{eff}}} )= \sum\limits_{i = 0}^{10} {{c_{ice\_i}}(\lambda )L_{eff}^{5 - i}}, $$

where ${g_{ice}}({\lambda ,{L_{eff}}} )$, ${\beta _{e\_ice}}({\lambda ,{L_{eff}}} )$, and ${\omega _{ice}}({\lambda ,{L_{eff}}} )$ represent the asymmetry factor, extinction coefficient (µm^-1), and single-scattering albedo of ice clouds, respectively. ${a_{ice\_i}}(\lambda )$, ${b_{ice\_i}}(\lambda )$, ${c_{ice\_i}}(\lambda )$ are the fitting coefficients at each wavelength.

Both spectral response function and Planck function vary with wavelength, and their variations should be considered for the simulation of the radiance at the top of the atmosphere (TOA). Therefore, the channel-averaged optical property (e.g., ${g_{ice\_\bar{\lambda }}}$) is calculated as

(11)$${g_{ice\_\bar{\lambda }}}({L_{eff}}) = \frac{{\int_{{\lambda _{\min }}}^{{\lambda _{\max }}} {{g_{ice}}(\lambda ,{L_{eff}})\phi (\lambda )B(\lambda ,{T_c})d\lambda } }}{{\int_{{\lambda _{\min }}}^{{\lambda _{\max }}} {\phi (\lambda )B(\lambda ,{T_c})d\lambda } }}\textrm{ = }\sum\limits_{i = 0}^{10} {{{\bar{a}}_{ice\_i}}L_{eff}^{5 - i}}, $$

where ${\lambda _{\max }}$ and ${\lambda _{\min }}$ are the maximum and minimum wavelength of the channel, respectively; $B({\lambda ,{T_c}} )$ refers to the Planck function; ${T_c}$ is equal to 287 K; ${\bar{a}_{ice\_i}}$ is equal to $\frac{{\int_{{\lambda _{\min }}}^{{\lambda _{\max }}} {{a_{ice\_i}}(\lambda )\phi (\lambda )B({\lambda ,{T_c}} )d\lambda } }}{{\int_{{\lambda _{\min }}}^{{\lambda _{\max }}} {\phi (\lambda )B({\lambda ,{T_c}} )d\lambda } }}$ (${\bar{b}_{ice\_i}}$ and ${\bar{c}_{ice\_i}}$ are obtained using a similar formula). The channel-averaged fitting coefficients for the optical properties of ice clouds at each TIR channel are available at the Harvard Dataverse (https://doi.org/10.7910/DVN/DEK4UW).

Figures 2(a)–2(c) present the channel-averaged optical properties of Voronoi aggregate at each TIR channel based on the established parameterization. Figure 2(a) shows that the asymmetry factor increases gradually with the effective length but ultimately stabilizes when the effective length is larger than 40 µm. In Fig. 2(b), the extinction coefficient drops sharply at almost all TIR channels when the effective length is less than 40 µm. It can be seen from Fig. 2(c) that single-scattering albedo converges to 0.5 as the effective length increases to 100 µm or larger. Figures 2(d)–2(f) display the relative errors (%) of the channel-averaged optical properties of ice clouds based on the parameterization as compared to the exact result. The absolute values of relative error of the asymmetry factor and single-scattering albedo are both less than 0.5%, while the maximum absolute value of the relative error of the extinction coefficient is 3.81%.

Fig. 2. Channel-averaged optical properties of ice cloud model (Voronoi aggregate): (a) asymmetry factor, (b) extinction coefficient, and (c) single-scattering albedo, based on the established parameterization, and their relative errors (%) (d-f) compared to the exact result, as a function of effective length ranging from 4 to 200 µm at each TIR channel.

Download Full Size | PDF

For water clouds, optical properties are obtained based on Lorenz-Mie theory [38], assuming the water droplets size distribution to be the gamma distribution. A number of effective radii ranging from 2 to 60 µm are calculated for 59 different size distributions, and the optical properties of water clouds assuming the same size distributions are obtained at 483 wavelengths ranging from 0.2 to 100 µm.

Then the parameterization for the optical properties of water clouds as a function of effective radius is established using similar fitting polynomials as ice clouds.

(12)$${g_{wat\_\bar{\lambda }}}({{R_{eff}}} )= \sum\limits_{i = 0}^{10} {{{\bar{a}}_{wat\_i}}R_{eff}^{5 - i}}, $$

(13)$${\beta _{e\_wat\_\bar{\lambda }}}({{R_{eff}}} )= \sum\limits_{i = 0}^{10} {{{\bar{b}}_{wat\_i}}R_{eff}^{5 - i}}, $$

(14)$$1 - {\omega _{wat\_\bar{\lambda }}}({{R_{eff}}} )= \sum\limits_{i = 0}^{10} {{{\bar{c}}_{wat\_i}}R_{eff}^{5 - i}}, $$

where ${R_{eff}}$ (µm) is the effective radius of water clouds; ${g_{wat\_\bar{\lambda }}}({{R_{eff}}} )$, ${\beta _{e\_wat\_\bar{\lambda }}}({{R_{eff}}} )$, and ${\omega _{wat\_\bar{\lambda }}}({{R_{eff}}} )$ represent the channel-averaged asymmetry factor, extinction coefficient (µm^-1), and single-scattering albedo of water clouds, respectively; ${\bar{a}_{wat\_i}}$, ${\bar{b}_{wat\_i}}$, ${\bar{c}_{wat\_i}}$ are the channel-averaged fitting coefficients for the optical properties of water clouds. The exact values of ${\bar{a}_{wat\_i}}$, ${\bar{b}_{wat\_i}}$, ${\bar{c}_{wat\_i}}$ are also shown at the aforementioned publicly-accessible website.

Figures 3(a)–3(c) show the channel-averaged optical properties of water clouds at each TIR channel based on the established parameterization. As shown in Fig. 3, the variation of optical properties of water clouds with effective radius shows a similar characteristic to that of ice clouds. Figures 3(d)–3(f) display the relative errors (%) of the channel-averaged optical properties of water clouds based on the established parameterization as compared to the exact result. The absolute value of the relative error of the asymmetry factor is close to 0%. The maximum absolute value of relative error is 0.40% for single-scattering albedo, while it is 1.23% for the extinction coefficient.

Fig. 3. Channel-averaged optical properties of water clouds: (a) asymmetry factor, (b) extinction coefficient, and (c) single-scattering albedo, based on the established parameterization, and their relative errors (%) (d-f) compared to the exact result, as a function of effective radius ranging from 2 to 60 µm at each TIR channel.

Download Full Size | PDF

2.3 Radiative transfer calculation

Clouds play a crucial role in the infrared radiative transfer, and the determination of the cloud top pressure ${P_{cld\_top}}$(hPa) and cloud base pressure ${P_{cld\_bas}}$(hPa) in the radiative transfer model is critical. The cloud with a larger optical depth tends to be geometrically thicker [39]. The geometric thickness H (m) is calculated by a function of cloud water path [9]. Then the cloud base pressure can be inferred from the geometric thickness and atmospheric profiles by the following formula

(15)$${P_{cld\_bas}} = {P_{cld\_top}}\exp \left[ {\frac{{M{g_{Earth}}}}{{({1 + 0.608q} )R\Gamma }}\ln \left( {\frac{{{T_{cld\_top}} + \Gamma H}}{{{T_{cld\_top}}}}} \right)} \right], $$

where ${T_{cld\_top}}$(K) is the cloud top temperature, and q (kg/kg) is the average specific humidity between cloud top and cloud base. The constants are molar mass of dry air M = 28.96 g/mol, gravitational acceleration ${g_{Earth}}$ = 9.8 m/s², gas constant R = 8.314 J/mol/K, and temperature lapse rate $\Gamma = \textrm{0}\textrm{.0065 K/m}$. For ice clouds, the ice water path of the l-th cloud layer $\Delta {W_{ice}}$ between the P_l₊₁ and P_l can be calculated from the vertical accumulated ice water path ${W_{ice}}$ by

(16)$$\Delta {W_{ice}} = \frac{{\ln {P_{l + 1}} - \ln {P_l}}}{{\ln {P_{cld\_bas}} - \ln {P_{cld\_top}}}} \times {W_{ice}}\begin{array}{c} {}\\ {} \end{array}({{P_{cld\_top}} \le {P_l} \le {P_{cld\_bas}}} ), $$

and the liquid water path of the l-th cloud layer $\Delta {W_{wat}}$ can be calculated using a similar formula.

Then the cloud extinction optical thickness of the l-th cloud layer ${\tau _{ice\_\bar{\lambda }}}$(${\tau _{wat\_\bar{\lambda }}}$) for ice clouds (water clouds) is calculated as

(17)$${\tau _{ice\_\bar{\lambda }}} = {\beta _{e\_ice\_\bar{\lambda }}}\Delta {W_{ice}}/{\rho _{ice}}, $$

(18)$${\tau _{wat\_\bar{\lambda }}} = {\beta _{e\_wat\_\bar{\lambda }}}\Delta {W_{wat}}/{\rho _{wat}}, $$

where ${\rho _{ice}}$ is the ice density, and ${\rho _{wat}}$ is the water density.

Finally, the total optical parameters of an arbitrary atmospheric layer with gaseous and cloud matter for $\Delta {\hat{g}_i}$ are

(19)$$\tau = {\tau _{gas\_i}}\textrm{ + }{\tau _{ice\_\bar{\lambda }}} + {\tau _{wat\_\bar{\lambda }}}, $$

(20)$$\omega = \frac{{{\omega _{ice\_\bar{\lambda }}}{\tau _{ice\_\bar{\lambda }}} + {\omega _{wat\_\bar{\lambda }}}{\tau _{wat\_\bar{\lambda }}}}}{{{\tau _{ice\_\bar{\lambda }}}\textrm{ + }{\tau _{wat\_\bar{\lambda }}}\textrm{ + }{\tau _{gas\_i}}}}, $$

(21)$$g = \frac{{{g_{ice\_\bar{\lambda }}}{\omega _{ice\_\bar{\lambda }}}{\tau _{ice\_\bar{\lambda }}} + {g_{wat\_\bar{\lambda }}}{\omega _{wat\_\bar{\lambda }}}{\tau _{wat\_\bar{\lambda }}}}}{{{\omega _{ice\_\bar{\lambda }}}{\tau _{ice\_\bar{\lambda }}} + {\omega _{wat\_\bar{\lambda }}}{\tau _{wat\_\bar{\lambda }}}}}. $$

The basic radiative transfer equation for thermal infrared radiation in a plane-parallel scattering atmosphere [40] is

(22)$$\mu \frac{{dI(\tau ,\mu )}}{{d\tau }} = I(\tau ,\mu ) - \frac{\omega }{2}\int_{ - 1}^1 {P(\mu ,\mu ^{\prime})} I(\tau ,\mu ^{\prime})d\mu ^{\prime} - (1 - \omega )B(\tau ), $$

where $\tau$, $\mu$and $\omega$ are the optical thickness, cosine of the zenith angle, and single-scattering albedo, respectively; $I({\tau ,\mu } )$ is the radiative intensity; P is the azimuthal independent phase function. The Planck function $B(\tau )$ is approximated exponentially in optical thickness as

(23)$$B(\tau ) = {B_1}{e^{ - {f_1}\tau }}, $$

where ${f_1}$ is equal to $- [{\ln ({{B_2}/{B_1}} )} ]/{\tau _1}$, and ${B_1}$ and ${B_2}$ are Planck functions at the top ($\tau = 0$) and bottom ($\tau = {\tau _1}$) of the layer.

The four-stream DOM solution which strikes the perfect balance between accuracy and efficiency is employed. Therefore, the radiative transfer equation is obtained by using two-node Gaussian quadrature, and Eq. (22) becomes

(24)$${\mu _i}\frac{{dI(\tau ,{\mu _i})}}{{d\tau }} = I(\tau ,{\mu _i}) - \frac{\omega }{2}\sum\limits_{l = 0}^3 {{\omega _l}{P_l}({{\mu_i}} )\sum\limits_{j ={-} 2}^2 {{a_j}{P_l}({{\mu_j}} )I(\tau ,{\mu _j})} } - (1 - \omega ){B_1}{e^{ - {f_1}\tau }}, $$

where i=±1, ±2, ${\mu _\textrm{1}} ={-} {\mu _{ - 1}} = 0.2113248$, ${\mu _\textrm{2}} ={-} {\mu _{ - \textrm{2}}} = 0.\textrm{7886752}$, and a_j = 0.5 (j=±1, ±2).${\omega _l}$ denotes the Legendre coefficients of the phase function.

Then the radiative transfer in a vertically inhomogeneous multiple-layer atmosphere is settled by the analytical adding method based on the infrared invariance principle. Finally, we can obtain the radiative intensity ${{\mathbf I}_n} = \left[ {\begin{array}{c} {I({\tau_{1,n}},{\mu_2})}\\ {I({\tau_{1,n}},{\mu_1})}\\ {I({\tau_{1,n}},{\mu_{ - 1}})}\\ {I({\tau_{1,n}},{\mu_{ - 2}})} \end{array}} \right]$ (n=1, 2, …, N, N+1) at the Gaussian quadrature points at each atmospheric level, where N is the total number of atmospheric layer, and ${\tau _{1,n}}$ represents the optical thickness between level 1 and level n (${\tau _{1,n}} = \sum\limits_{i = 1}^{n - 1} {{\tau _i}}$). A detailed description of the analytical adding method can be found in [24].

To calculate the intensity at arbitrary zenith angles in multiple layers, Eq. (24) for arbitrary zenith angles at layer n (from level n to level n+1) is written as

(25)$$\begin{array}{l} \mu \frac{{dI({\tau _{1,n}} + \tau ,\mu )}}{{d\tau }} = I({\tau _{1,n}} + \tau ,\mu ) - \frac{{{\omega _n}}}{2}\sum\limits_{l = 0}^3 {{\omega _{ln}}{P_l}(\mu )\sum\limits_{j ={-} 2}^2 {{a_j}{P_l}({{\mu_j}} )I({\tau _{1,n}} + \tau ,{\mu _j})} } \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} - (1 - {\omega _n}){B_n}{e^{ - {f_n}\tau }} \end{array}, $$

where ${f_n} ={-} [{\ln ({{B_{n + 1}}/{B_n}} )} ]/{\tau _n}$, and $0 \le \tau \le {\tau _n}$.

To solve the Eq. (25), the explicit expression of $I({\tau _{1,n}} + \tau ,{\mu _j})$ (j=±1, ±2) is needed. Therefore, $I({\tau _{1,n}} + \tau ,{\mu _j})$ with dependence on $\tau$ is expressed as

(26)$$\left[ {\begin{array}{c} {I({{\tau_{1,n}}\textrm{ + }\tau ,{\mu_2}} )}\\ {I({{\tau_{1,n}}\textrm{ + }\tau ,{\mu_1}} )}\\ {I({{\tau_{1,n}}\textrm{ + }\tau ,{\mu_{ - 1}}} )}\\ {I({{\tau_{1,n}}\textrm{ + }\tau ,{\mu_{ - 2}}} )} \end{array}} \right] = \left[ {\begin{array}{cccc} {\varphi_{2n}^ + {e_{2n}}}&{\varphi_{1n}^ + {e_{1n}}}&{\varphi_{1n}^ - {e_{3n}}}&{\varphi_{2n}^ - {e_{4n}}}\\ {\Phi _{2n}^ + {e_{2n}}}&{\Phi _{1n}^ + {e_{1n}}}&{\Phi _{1n}^ - {e_{3n}}}&{\Phi _{2n}^ - {e_{4n}}}\\ {\Phi _{2n}^ - {e_{2n}}}&{\Phi _{1n}^ - {e_{1n}}}&{\Phi _{1n}^ + {e_{3n}}}&{\Phi _{2n}^ + {e_{4n}}}\\ {\varphi_{2n}^ - {e_{2n}}}&{\varphi_{1n}^ - {e_{1n}}}&{\varphi_{1n}^ + {e_{3n}}}&{\varphi_{2n}^ + {e_{4n}}} \end{array}} \right]{{\mathbf G}_{1n}} + \left[ {\begin{array}{c} {Z_{2n}^ + }\\ {Z_{1n}^ + }\\ {Z_{1n}^ - }\\ {Z_{2n}^ - } \end{array}} \right]{e^{ - {f_n}\tau }}, $$

where ${e_1}_n = {e^{ - {k_{1n}}\tau }}$, ${e_{2n}} = {e^{ - {k_{2n}}\tau }}$, ${e_{3n}} = {e^{ - {k_{1n}}({{\tau_n} - \tau } )}}$, and ${e_{4n}} = {e^{ - {k_{2n}}({{\tau_n} - \tau } )}}$. The terms of $\varphi _{1n}^ \pm$, $\varphi _{2n}^ \pm$, $\Phi _{1n}^ \pm$, $\Phi _{2n}^ \pm$, $Z_{1n}^ \pm$, $Z_{2n}^ \pm$, ${k_{1n}}$ and ${k_{2n}}$ are defined the same as in [24], and listed in their appendix A. ${{\mathbf G}_{1n}} = {[{G_{2n}^1G_{1n}^1G_{ - 1n}^1G_{ - 2n}^1} ]^{\mathbf {\rm T}}}$, with the superscript ${\mathbf T}$ indicating a matrix transpose. ${{\mathbf G}_{1n}}$ can be obtained by

(27)$${{\mathbf G}_{1n}} = {\mathbf A}_n^{ - 1}({{{\mathbf I}_n} - {{\mathbf H}_n}} ), $$

where ${{\mathbf A}_n} = \left[ {\begin{array}{cccc} {\varphi_{2n}^ + }&{\varphi_{1n}^ + }&{\varphi_{1n}^ - {\psi_{1n}}}&{\varphi_{2n}^ - {\psi_{2n}}}\\ {\Phi _{2n}^ + }&{\Phi _{1n}^ + }&{\Phi _{1n}^ - {\psi_{1n}}}&{\Phi _{2n}^ - {\psi_{2n}}}\\ {\Phi _{2n}^ - }&{\Phi _{1n}^ - }&{\Phi _{1n}^ + {\psi_{1n}}}&{\Phi _{2n}^ + {\psi_{2n}}}\\ {\varphi_{2n}^ - }&{\varphi_{1n}^ - }&{\varphi_{1n}^ + {\psi_{1n}}}&{\varphi_{2n}^ + {\psi_{2n}}} \end{array}} \right]$, ${{\mathbf H}_n} = \left[ {\begin{array}{c} {Z_{2n}^ + }\\ {Z_{1n}^ + }\\ {Z_{1n}^ - }\\ {Z_{2n}^ - } \end{array}} \right]$, ${\psi _{1n}} = {e^{ - {k_{1n}}{\tau _n}}}$, and ${\psi _{2n}} = {e^{ - {k_{2n}}{\tau _n}}}$.

By substituting Eq. (26) to Eq. (25), the solution of Eq. (25) can be obtained. The intensity at arbitrary zenith angles at TOA $I({0,\mu } )$ is derived in an upward path from the level N+1 to level 1 as follows:

(28)$$\begin{array}{l} I({{\tau_{1,n}},\mu } )= I({{\tau_{1,n\textrm{ + }1}},\mu } ){e^{ - \frac{{{\tau _n}}}{\mu }}} - \frac{{{\omega _n}}}{{2\mu }}\sum\limits_{l = 0}^3 {{\omega _{ln}}{P_l}(\mu )\sum\limits_{j ={-} 2}^2 {{a_j}{P_l}({{\mu_j}} )} } {S_{jn}}\\ \;\;\;\;\;\;\;\;\;\;\;\;{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} + \frac{{(1 - {\omega _n}){B_n}}}{{\mu {f_n} + 1}}(1 - {e^{ - ({f_n} + \frac{1}{\mu }){\tau _n}}}){\kern 1pt} \end{array}, $$

where n = N, N-1, …, 2, 1. The term ${S_{jn}}$ (j=±1, ±2) is

(29)$$\begin{array}{l} {S_{ {\pm} 1n}} = \Phi _{2n}^ \pm G_{2n}^1({e^{ - ({k_{2n}} + 1/\mu ){\tau _n}}} - 1)/({k_{2n}} + 1/\mu ) + \Phi _{1n}^ \pm G_{1n}^1({e^{ - ({k_{1n}} + 1/\mu ){\tau _n}}} - 1)/({k_{1n}} + 1/\mu )\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} + \Phi _{1n}^ \mp G_{ - 1n}^1({e^{ - {k_{1n}}{\tau _n}}} - {e^{ - {\tau _n}/\mu }})/({k_{1n}} - 1/\mu ) + \Phi _{2n}^ \mp G_{ - 2n}^1({e^{ - {k_{2n}}{\tau _n}}} - {e^{ - {\tau _n}/\mu }})/({k_{2n}} - 1/\mu )\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} + Z_{1n}^ \pm ({e^{ - ({f_n} + 1/\mu ){\tau _n}}} - 1)/({f_n} + 1/\mu ) \end{array}, $$

(30)$$\begin{array}{l} {S_{ {\pm} 2n}} = \varphi _{2n}^ \pm G_{2n}^1({e^{ - ({k_{2n}} + 1/\mu ){\tau _n}}} - 1)/({k_{2n}} + 1/\mu ) + \varphi _{1n}^ \pm G_{1n}^1({e^{ - ({k_{1n}} + 1/\mu ){\tau _n}}} - 1)/({k_{1n}} + 1/\mu )\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} + \varphi _{1n}^ \mp G_{ - 1n}^1({e^{ - {k_{1n}}{\tau _n}}} - {e^{ - {\tau _n}/\mu }})/({k_{1n}} - 1/\mu ) + \varphi _{2n}^ \mp G_{ - 2n}^1({e^{ - {k_{2n}}{\tau _n}}} - {e^{ - {\tau _n}/\mu }})/({k_{2n}} - 1/\mu )\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} + Z_{2n}^ \pm ({e^{ - ({f_n} + 1/\mu ){\tau _n}}} - 1)/({f_n} + 1/\mu ) \end{array}, $$

and the surface boundary condition is considered as

(31)$$I({\tau _{1,N + 1}},\mu ) = 2(1 - {\varepsilon _s})\sum\limits_{j = 1}^2 {{a_j}{\mu _j}I({\tau _{1,N + 1}},{\mu _{ - j}})} + {\varepsilon _s}{B_s}, $$

where ${\varepsilon _s}$ and ${B_s}$ are the surface emissivity and Planck function evaluated at the surface temperature, respectively. Moreover, the δ-function adjustment is employed in the radiative scheme to enhance the accuracy [41], and hereafter the adding algorithm of the delta-four-stream DOM is denoted as δ-4DDA.

Based on the acquired optical parameters and the cosine of satellite zenith angle ${\mu _{sat}}$, the intensity at TOA denoted as ${I_{TOA}}$ can be obtained from δ-4DDA by ${I_{TOA}} = I({0,{\mu_{sat}}} )$. Following the AMCKD, each channel is divided into six subintervals, and the radiative transfer process needs to be dealt with in each sub-interval, separately. The channel-averaged radiance is obtained as

(32)$${\bar{I}_{TOA}} = \sum\limits_{i = 1}^{{N_{\hat{g}}}} {{I_{TOA\_i}}} \cdot \Delta {\hat{g}_i}, $$

where ${N_{\hat{g}}} = 6$, ${I_{TOA\_i}}$ is the intensity at TOA of the i-th subinterval. Then, the BT at TOA, ${T_B}$, can be obtained from the channel averaged radiance ${\bar{I}_{TOA}}$ [42] by

(33)$${T_B} = \frac{{hc}}{{K\ln (\frac{{2h{c^2}}}{{{{\bar{I}}_{TOA}}{{\hat{\lambda }}^5}}} + 1)\hat{\lambda }}}, $$

where h, c, K is the Planck’s constant, speed of light, and Boltzmann’s constant, respectively. $\hat{\lambda }$ is the center wavelength of the channel.

The δ-64 DISORT is inserted into LBLRTM with spectral resolution 0.1 cm^-1 as a benchmark model. The simulated BTs of ERTM are evaluated by comparing them with the results of the benchmark model. The brightness temperature difference (BTD) is defined as the BT bias between the ERTM and the benchmark model. Figure 4 illustrates the BTDs for various cloud conditions based on the mid-latitude summer atmospheric profile with a surface emissivity of 0.99 under two viewing zenith angles (20° and 50°). The details of the cloud conditions are listed in the figure. Here B10, B12, and B14 are chosen to discuss in particular, that is, two absorption channels (B10 and B12) and one window channel (B14). The absolute values of BTDs at all three channels are less than 0.3 K for any cloud condition assumption and zenith angle. Besides, the BTDs are more dependent on the zenith angle and less dependent on the cloud effective radius or cloud top pressure. Table 2 lists the root mean square error (RMSE) of simulated BTs by ERTM, AMCKD with δ-4 DISORT, and AMCKD with δ-32 DISORT, compared to benchmark results at each TIR channel under the cloud assumptions illustrated in Fig. 4. For ERTM, the maximum RMSE is only 0.21 K at the B16 channel, and the minimum RMSE is 0.04 K at the B13 channel. Also, ERTM generally has similar accuracy to AMCKD with δ-4 DISORT, but is slightly less accurate than AMCKD with δ-32 DISORT. Table 3 compares the runtime of ERTM, AMCKD with δ-4 DISORT, AMCKD with δ-32 DISORT, AMCKD with δ-64 DISORT, and benchmark model. The simulation is implemented on a computing cluster running the 64-bit Linux operating system (CentOS Linux release 7.2.1511). The cluster is composed of 1 login node (2 Intel Xeon Silver 4110 CPU; 16 cores; 61 GB of memory) and 3 compute nodes (2 Intel Xeon E5-2680 v4 CPU per node; 28 cores per node; 125 GB of memory per node). All the compute nodes in the cluster are used for the simulation, and it takes an average of 4.98 hours per core to accomplish the simulation based on the parallel computing technology. The computational efficiency of ERTM is approximately five orders of magnitude higher than that of the benchmark model, but two orders of magnitude higher than that of AMCKD with δ-64 DISORT (or δ-32 DISORT, δ-4 DISORT). Furthermore, compared with AMCKD with δ-4 DISORT, ERTM can gain equivalent accuracy but greatly reduce the computational cost.

Fig. 4. Brightness temperature differences as a function of ice water path (IWP) for ice clouds (left column) and liquid water path (LWP) for water clouds (right column) under two viewing zenith angles (20° and 50°) at the TIR B10 (7.35 µm), B12 (9.63 µm) and B14 (11.24 µm) channels.

Download Full Size | PDF

Table 2. RMSE of simulated brightness temperatures by ERTM, AMCKD with δ-4 DISORT, and AMCKD with δ-32 DISORT, compared to benchmark results at each TIR channel.

View Table | View all tables in this article

Table 3. Comparison of the runtime of ERTM, AMCKD with δ-4 DISORT, AMCKD with δ-32 DISORT, AMCKD with δ-64 DISORT, and benchmark model.

View Table | View all tables in this article

2.4 ERTM integration

The data used in ERTM are summarized in Table 4. Atmospheric profiles including the specific humidity, ozone mass mixing ratio, and temperature are from ERA-Interim product with a vertical resolution of 37 pressure levels, a temporal resolution of 6 h, and a spatial resolution of 0.75°×0.75°. ERA-Interim product is a global atmospheric reanalysis from the data assimilation system based on the European Centre for Medium-Range Weather Forecasts (ECMWF) Integrated Forecasting System (IFS) [43]. The land surface emissivity and surface temperature data in a 0.05° climate modeling grid are obtained from the MODIS land 8-day mean level 3 product (MYD11C2), retrieved through the day/night algorithm [44]. The land surface temperature products have been validated within ±1 K for the absolute errors under clear sky conditions by in-situ measurements. Emissivity products are highly consistent with the surface emissivity spectra measured by the sun-shadow method [45]. Sea surface temperature data are from the NOAA 1/4° daily Optimum Interpolation Sea Surface Temperature (daily OISST), which is an analysis built by uniting observations from various platforms including satellites, ships, and buoys using optimum interpolation (OI) [46]. RMSE of the daily OISST product relative to the buoys is about 0.3 K [47]. The sea surface emissivity is set as 0.99. The cloud properties, including cloud top pressure, cloud top temperature, cloud water path, cloud thermodynamic phase, cloud optical thickness, and cloud effective particle radius, are retrieved from Integrated Cloud Analysis System (ICAS) via an optimal estimation (OE) approach [48] using eight TIR bands of AHI [8]. The retrievals of cloud property from ICAS are highly consistent with collocated active remote sensing counterparts especially for single-layer clouds [8]. The measurements used in this paper are the Himawari-8 level 1 gridded data with a temporal resolution of 10 min and a spatial resolution of 0.05°×0.05°. All the meteorological field data and cloud property data are interpolated in space according to the measurement data.

Table 4. Summary of data used in ERTM.

View Table | View all tables in this article

A flow chart of Fig. 5 is used to illustrate the application of ERTM in the simulation of AHI observation. The TIR emission is generated from the ERA-interim temperature profiles and MODIS and OISST surface temperature. Gaseous absorption is obtained from ERA-interim atmospheric profile based on the AMCKD. Cloud optical properties are calculated by the pre-established cloud optical property parameterization. Then the BTs are simulated in δ-4DDA with these essential inputs. Finally, the simulated results are compared with rigorous results and AHI measurements to estimate the efficiency and accuracy of ERTM.

Fig. 5. Flow chart of the application of ERTM (AMCKD scheme, a cloud optical property parameterization and δ-4DDA) in the simulation of AHI observation.

Download Full Size | PDF

3. Results and discussion

Typhoon Mujigae on 3 October 2015 at 6:00 UTC is used as a practical case in this study to validate the ERTM. Typhoon Mujigae is the strongest typhoon in 2015, it led to severe weather events that caused heavy casualties and economic losses [49]. Figure 6 shows the cloud properties retrieved by ICAS for Typhoon Mujigae on 3 October 2015 at 6:00 UTC. Figure 6(c) indicates that the whole typhoon is covered by high clouds since the cloud top pressure values at 680hPa and 440hPa are used as a threshold to separate the cloud into the low, middle and high cloud [50]. Besides, the large cloud optical thickness (> 6) in the eyewall and spiral rainband of the typhoon demonstrates the existence of strong convection because the high clouds are further classified into cirrus, cirrostratus, and deep convective clouds using cloud optical thickness boundaries at 1 and 6 (Fig. 6(a)) . Figures 6(b) and 6(d) show that the typhoon is largely covered by ice clouds of large cloud effective radius only with little water clouds at the margin.

Fig. 6. Cloud properties for Typhoon Mujigae on 3 October 2015 at 6:00 UTC (a) cloud optical thickness, (b) cloud effective radius (µm), (c) cloud top pressure (hPa) and (d) cloud type retrieved from AHI observations.

Download Full Size | PDF

One thousand cloudy pixels are randomly selected from the simulated region to compare the simulated BTs by ERTM with the benchmark results in the practical case. The simulation is also implemented by the aforementioned computing cluster, and it takes an average of 8.93 hours per core to accomplish the simulation. The corresponding results are shown in Fig. 7, in which the simulations by ERTM have a high agreement with the standard results, and almost all pixels fall on the 1:1 line at each TIR channel. Moreover, the maximum RMSE of the simulated BTs is only 0.3665 at the TIR B10 channel. Then we systematically evaluate the δ-4DDA in ERTM by comparing the simulated result between ERTM and AMCKD with δ-64 DISORT in the whole study domain. These whole domain simulations take an average of 10.85 hours per core using the computing cluster. As shown in Fig. 8, the BTDs between ERTM and AMCKD with δ-64 DISORT are less than 0.3 K at all nine channels, which demonstrates the prominent accuracy of ERTM for the brightness temperature simulation in cloudy atmospheres.

Fig. 7. Comparison between brightness temperatures simulated by ERTM and results of the benchmark model for one thousand cloudy pixels at each TIR channel.

Download Full Size | PDF

Fig. 8. Comparison between brightness temperatures simulated by ERTM and results of the AMCKD with δ-64 DISORT at each TIR channel.

Download Full Size | PDF

To further assess the performance of the ERTM, simulated BTs are compared with observations from AHI. Figures 9–11 show the comparison between the observations and simulations of BTs at each TIR channels. The BTs of clear sky at B08-10 (water vapor absorption), B12 (ozone absorption), and B16 (carbon dioxide absorption) channels are smaller than those at the other channels (atmospheric window) where the gaseous absorption is weak. As shown in the Figs. 9–11, the BTDs of clear sky at B08 (strong water vapor absorption) channel are slightly larger than those at the other channels, which is probably due to the uncertainties of water vapor in reanalysis data. Generally, the distribution pattern of observed brightness temperatures is highly consistent with that of simulated brightness temperatures, demonstrating that the ERTM has excellent capability to simulate BTs for the clear and cloudy atmosphere. Then we quantitatively compare the observations with simulations over the simulated region. Figure 12 shows the violin plot of absolute bias between the AHI observations and simulated brightness temperatures for ice clouds and water clouds at each TIR channel. For all TIR channels, the mean absolute biases are less than 2, and 90% of pixels have an absolute bias of less than 3.

Fig. 9. Comparison between AHI observations (first row) and simulated brightness temperatures (second row) at the TIR B08 (6.24 µm), B09 (6.94 µm), and B10 (7.35 µm) channels. The differences between AHI observations and simulated brightness temperatures are shown in the third row.

Download Full Size | PDF

Fig. 10. Comparison between AHI observations (first row) and simulated brightness temperatures (second row) at the TIR B11 (8.59 µm), B12 (9.63 µm), and B13 (10.40 µm) channels. The differences between AHI observations and simulated brightness temperatures are shown in the third row.

Download Full Size | PDF

Fig. 11. Comparison between AHI observations (first row) and simulated brightness temperatures (second row) at the TIR B14 (11.24 µm), B15 (12.38 µm), and B16 (13.28 µm) channels. The differences between AHI observations and simulated brightness temperatures are shown in the third row.

Download Full Size | PDF

Fig. 12. The violin plot of absolute bias between the AHI observations and simulated brightness temperatures for ice clouds (a) and water clouds (b) at each TIR channel. Plots indicate the distribution (shaded), 10th and 90th percentiles (black cross), median (red horizontal dashed line) and mean (yellow horizontal line).

Download Full Size | PDF

The disparity between simulation and observation is indeed associated with many factors: (1) The uncertainties may come from ERTM itself including approximate rapid radiative transfer scheme and channel-averaged method for gaseous absorption, which have been proved negligible by comparing simulated results with benchmark results shown in the previous section. (2) Disparities may be caused by different ice crystal habits because ice particle is hypothesized to be 8-element column aggregate with severely rough surface in ICAS while Voronoi aggregate in ERTM. However, these disparities are generally under a few percent according to the previous work that TIR measurements are insensitive to the ice habit assumption [51,52]. (3) The uncertainties may be induced by approximate plane-parallel radiative transfer model without consideration of 3-D transfer effect [53] and cloud heterogeneity (e.g., horizontally inhomogeneity (sub-pixel variation), vertical variation of particle microphysical properties) (4) The deviations of external input data such as cloud properties retrieved by ICAS and meteorological field data (e.g., background surface properties, atmospheric profiles) are a major source of disparities. However, factor (3) and (4) are beyond the scope of this study.

4. Summary

In this study, an efficient radiative transfer model (ERTM) is developed to simulate thermal infrared brightness temperatures observed by the Advanced Himawari Imager (AHI). The ERTM contains an alternate mapping correlated k-distribution scheme, a parameterization for cloud optical property and an efficient infrared radiative transfer scheme. The AMCKD is employed to calculate the gaseous absorption in the inhomogeneous thermodynamic atmosphere. The optical properties of clouds are parameterized by the cloud effective length for ice clouds based on the Voronoi model, and by the cloud effective radius for water clouds based on the Lorenz-Mie theory. The δ-4DDA is exploited in ERTM to obtain the radiance at TOA. The δ-64 DISORT is inserted into LBLRTM as a benchmark model. A series of comparisons between ERTM and benchmark model under various cloud conditions are discussed to assess the performance of ERTM. Under the standard atmospheric profiles, the RMSEs of simulated BTs reach a maximum of 0.21K at the B16 channel, while the minimum RMSE is only 0.04K at the B13 channel. Moreover, the computational efficiency of ERTM is approximately five orders of magnitude higher than that of the benchmark model. Finally, Typhoon Mujigae on 3 October 2015 at 6:00 UTC is used as a practical case to validate the ERTM. The simulated brightness temperatures by ERTM are highly consistent with the rigorous results and observations from AHI, which verifies the excellent performance of ERTM for the realistic atmosphere. Given its great accuracy and high computational efficiency, ERTM is suitable for model assessment, data assimilation, and remote sensing retrieval. The ERTM is available from the corresponding author.

Funding

National Key Research and Development Program of China (2018YFC1507002); National Natural Science Foundation of China (41675003).

Acknowledgments

The authors appreciate four anonymous reviewers for their constructive comments and Prof. Pengwang Zhai for his editorial effort. The authors also thank Dr. Kun Wu for help in revising the manuscript. The ERA-Interim products used in this study are available from the European Centre for Medium-Range Weather Forecasts (ECMWF) forecast models and data assimilation systems (https://www.ecmwf.int/en/forecasts/datasets/archive-datasets/browse-reanalysis-datasets). MODIS land 8-day mean level 3 product (MYD11C2) are provided by the Land Processes Distributed Active Archive Center (https://lpdaac.usgs.gov/). Optimum Interpolation Sea Surface Temperature (daily OISST) is obtained from National Oceanic and Atmospheric Administration (https://www.ncdc.noaa.gov/oisst). Himawari L1 Gridded data is provided by the Japan Meteorological Agency (http://www.jma.go.jp/jma/indexe.html).

WENWEN LI: Writing - original draft, Software, Validation, Visualization. FENG ZHANG: Writing - original draft, Writing - review & editing, Conceptualization, Methodology. YI-NING SHI: Review & editing. HIRONOBU IWABUCHI: Review & editing, Providing the cloud property data retrieved from Integrated Cloud Analysis System (OCAS). MINGWEI ZHU: Review & editing. JIANGNAN LI: Review & editing. WEI HAN: Review & editing. HUSI LETU: Review & editing. HIROSHI ISHIMOTO: Review & editing, Providing the ice particle optical property database of Voronoi aggregate model.

Disclosures

The authors declare no conflicts of interest.

References

1. V. Ramanathan, R. Cess, E. Harrison, P. Minnis, B. Barkstrom, E. Ahmad, and D. Hartmann, “Cloud-radiative forcing and climate: Results from the Earth Radiation Budget Experiment,” Science 243(4887), 57–63 (1989). [CrossRef]

2. A. Arking, “The radiative effects of clouds and their impact on climate,” Bull. Am. Meteorol. Soc. 72(6), 795–813 (1991). [CrossRef]

3. Z. Li, H. W. Barker, and L. Moreau, “The variable effect of clouds on atmospheric absorption of solar radiation,” Nature 376(6540), 486–490 (1995). [CrossRef]

4. T. Nakajima and M. D. King, “Determination of the optical thickness and effective particle radius of clouds from reflected solar radiation measurements. Part I: Theory,” J. Atmos. Sci. 47(15), 1878–1893 (1990). [CrossRef]

5. S. Platnick, M. D. King, S. A. Ackerman, W. P. Menzel, B. A. Baum, J. C. Riédi, and R. A. Frey, “The MODIS cloud products: Algorithms and examples from Terra,” IEEE Trans. Geosci. Remote Sensing 41(2), 459–473 (2003). [CrossRef]

6. C. Wang, P. Yang, B. A. Baum, S. Platnick, A. K. Heidinger, Y. Hu, and R. E. Holz, “Retrieval of ice cloud optical thickness and effective particle size using a fast infrared radiative transfer model,” J. Appl. Meteorol. Climatol. 50(11), 2283–2297 (2011). [CrossRef]

7. P. Minnis, S. Sun-Mack, D. F. Young, P. W. Heck, D. P. Garber, Y. Chen, D. A. Spangenberg, R. F. Arduini, Q. Z. Trepte, and W. L. Smith, “CERES Edition-2 cloud property retrievals using TRMM VIRS and Terra and Aqua MODIS data—Part I: Algorithms,” IEEE Trans. Geosci. Remote Sensing 49(11), 4374–4400 (2011). [CrossRef]

8. H. Iwabuchi, N. S. Putri, M. Saito, Y. Tokoro, M. Sekiguchi, P. Yang, and B. A. Baum, “Cloud Property Retrieval from Multiband Infrared Measurements by Himawari-8,” J. Meteorol. Soc. Japan 96B, 27–42 (2018). [CrossRef]

9. H. Iwabuchi, M. Saito, Y. Tokoro, N. S. Putri, and M. Sekiguchi, “Retrieval of radiative and microphysical properties of clouds from multispectral infrared measurements,” Prog. in Earth and Planet. Sci. 3(1), 32 (2016). [CrossRef]

10. C. Liu, P. Yang, S. L. Nasiri, S. Platnick, K. G. Meyer, C. Wang, and S. Ding, “A fast Visible Infrared Imaging Radiometer Suite simulator for cloudy atmospheres,” J. Geophys. Res.: Atmos. 120(1), 240–255 (2015). [CrossRef]

11. B. Yao, C. Liu, Y. Yin, P. Zhang, M. Min, and W. Han, “Radiance-Based Evaluation of WRF Cloud Properties Over East Asia: Direct Comparison With FY-2E Observations,” J. Geophys. Res.: Atmos. 123(9), 4613–4629 (2018). [CrossRef]

12. J.-P. Chaboureau, O. Nuissier, and C. Claud, “Verification of ensemble forecasts of Mediterranean high-impact weather events against satellite observations,” Nat. Hazards Earth Syst. Sci. 12(8), 2449–2462 (2012). [CrossRef]

13. M. Han, S. A. Braun, T. Matsui, and C. R. Williams, “Evaluation of cloud microphysics schemes in simulations of a winter storm using radar and radiometer measurements,” J. Geophys. Res.: Atmos. 118(3), 1401–1419 (2013). [CrossRef]

14. K. L. Chan, A. Hartl, Y. F. Lam, P. H. Xie, W. Q. Liu, H. M. Cheung, J. Lampel, D. Poehler, A. Li, and J. Xu, “Observations of tropospheric NO2 using ground based MAX-DOAS and OMIn measurements during the Shanghai World Expo 2010,” Atmos. Environ. 119, 45–58 (2015). [CrossRef]

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

16. K. Stamnes, S.-C. Tsay, W. Wiscombe, and K. Jayaweera, “Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting layered media,” Appl. Opt. 27(12), 2502–2509 (1988). [CrossRef]

17. S. Twomey, H. Jacobowitz, and H. Howell, “Matrix methods for multiple-scattering problems,” J. Atmos. Sci. 23(3), 289–298 (1966). [CrossRef]

18. J. Hovenier, “Symmetry relationships for scattering of polarized light in a slab of randomly oriented particles,” J. Atmos. Sci. 26(3), 488–499 (1969). [CrossRef]

19. J. R. Howell, “The Monte Carlo method in radiative heat transfer,” J. Heat Transfer 120(3), 547–560 (1998). [CrossRef]

20. N. Gentile, “Implicit Monte Carlo diffusion—an acceleration method for Monte Carlo time-dependent radiative transfer simulations,” J. Comput. Phys. 172(2), 543–571 (2001). [CrossRef]

21. J. Li and Q. Fu, “Absorption approximation with scattering effect for infrared radiation,” J. Atmos. Sci. 57(17), 2905–2914 (2000). [CrossRef]

22. J. Li, “Accounting for unresolved clouds in a 1D infrared radiative transfer model. Part I: Solution for radiative transfer, including cloud scattering and overlap,” J. Atmos. Sci. 59(23), 3302–3320 (2002). [CrossRef]

23. J. Li and V. Ramaswamy, “Four-stream spherical harmonic expansion approximation for solar radiative transfer,” J. Atmos. Sci. 53(8), 1174–1186 (1996). [CrossRef]

24. F. Zhang, K. Wu, J. Li, Q. Yang, J.-Q. Zhao, and J. Li, “Analytical infrared delta-four-stream adding method from invariance principle,” J. Atmos. Sci. 73(10), 4171–4188 (2016). [CrossRef]

25. E. J. Mlawer, M. J. Iacono, R. Pincus, H. W. Barker, L. Oreopoulos, and D. L. Mitchell, “Contributions of the ARM program to radiative transfer modeling for climate and weather applications,” Meteor. Monogr. 57, 15.1–15.19 (2016). [CrossRef]

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

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

28. J.-L. Moncet, G. Uymin, and H. E. Snell, “Atmospheric radiance modeling using the optimal spectral sampling (OSS) method,” in Algorithms and Technologies for Multispectral, Hyperspectral, and Ultraspectral Imagery X, (International Society for Optics and Photonics, 2004), 368–374.

29. S. A. Tjemkes and J. Schmetz, “Synthetic satellite radiances using the radiance sampling method,” J. Geophys. Res.: Atmos. 102(D2), 1807–1818 (1997). [CrossRef]

30. X. Liu, W. L. Smith, D. K. Zhou, and A. Larar, “Principal component-based radiative transfer model for hyperspectral sensors: theoretical concept,” Appl. Opt. 45(1), 201–209 (2006). [CrossRef]

31. J. Li and H. Barker, “A radiation algorithm with correlated-k distribution. Part I: Local thermal equilibrium,” J. Atmos. Sci. 62(2), 286–309 (2005). [CrossRef]

32. F. Zhang, M. Zhu, J. Li, W. Li, D. Di, Y. Shi, and K. Wu, “Alternate Mapping Correlated k-Distribution Method for Infrared Radiative Transfer Forward Simulation,” Remote Sens. 11(9), 994 (2019). [CrossRef]

33. K. Bessho, K. Date, M. Hayashi, A. Ikeda, T. Imai, H. Inoue, Y. Kumagai, T. Miyakawa, H. Murata, and T. Ohno, “An introduction to Himawari-8/9—Japan’s new-generation geostationary meteorological satellites,” J. Meteorol. Soc. Japan 94(2), 151–183 (2016). [CrossRef]

34. R. M. Goody and Y. L. Yung, Atmospheric radiation: theoretical basis (Oxford university press, 1995).

35. H. Ishimoto, K. Masuda, Y. Mano, N. Orikasa, and A. Uchiyama, “Irregularly shaped ice aggregates in optical modeling of convectively generated ice clouds,” J. Quant. Spectrosc. Radiat. Transfer 113(8), 632–643 (2012). [CrossRef]

36. H. Letu, H. Ishimoto, J. Riedi, T. Y. Nakajima, L. C-Labonnote, A. J. Baran, T. M. Nagao, and M. Sekiguchi, “Investigation of ice particle habits to be used for ice cloud remote sensing for the GCOM-C satellite mission,” Atmos. Chem. Phys. 16(18), 12287–12303 (2016). [CrossRef]

37. H. Zhang, Q. Chen, and B. Xie, “A new parameterization for ice cloud optical properties used in BCC-RAD and its radiative impact,” J. Quant. Spectrosc. Radiat. Transfer 150, 76–86 (2015). [CrossRef]

38. G. Mie, “Beiträge zur Optik trüber Medien, speziell kolloidaler Metallösungen,” Ann. Phys. 330(3), 377–445 (1908). [CrossRef]

39. P. Veglio and T. Maestri, “Statistics of vertical backscatter profiles of cirrus clouds,” Atmos. Chem. Phys. 11(24), 12925–12943 (2011). [CrossRef]

40. K. Liou and Y. Yung, “An introduction to atmospheric radiation,” Q. J. R. Meteorol. Soc. 129(590), 1741 (2003). [CrossRef]

41. W. Wiscombe, “The Delta–M method: Rapid yet accurate radiative flux calculations for strongly asymmetric phase functions,” J. Atmos. Sci. 34(9), 1408–1422 (1977). [CrossRef]

42. M. Williams and R. A. Houze, “Satellite-Observed Characteristics of Winter Monsoon Cloud Clusters,” Mon. Wea. Rev. 115(2), 505–519 (1987). [CrossRef]

43. D. P. Dee, S. Uppala, A. Simmons, P. Berrisford, P. Poli, S. Kobayashi, U. Andrae, M. Balmaseda, G. Balsamo, and D. P. Bauer, “The ERA-Interim reanalysis: Configuration and performance of the data assimilation system,” Q. J. R. Meteorol. Soc. 137(656), 553–597 (2011). [CrossRef]

44. Z. Wan and Z.-L. Li, “A physics-based algorithm for retrieving land-surface emissivity and temperature from EOS/MODIS data,” IEEE Trans. Geosci. Remote Sensing 35(4), 980–996 (1997). [CrossRef]

45. Z. Wan, “New refinements and validation of the MODIS land-surface temperature/emissivity products,” Remote Sens. Environ. 112(1), 59–74 (2008). [CrossRef]

46. V. Banzon, T. M. Smith, T. M. Chin, C. Liu, and W. Hankins, “A long-term record of blended satellite and in situ sea-surface temperature for climate monitoring, modeling and environmental studies,” Earth Syst. Sci. Data 8(1), 165–176 (2016). [CrossRef]

47. R. W. Reynolds and D. B. Chelton, “Comparisons of daily sea surface temperature analyses for 2007–08,” J. Clim. 23(13), 3545–3562 (2010). [CrossRef]

48. C. D. Rodgers, Inverse methods for atmospheric sounding: theory and practice (World scientific, 2000), Vol. 2.

49. B. Wang, M. Wei, W. Hua, Y. Zhang, X. Wen, J. Zheng, N. Li, H. Li, Y. Wu, and J. Zhu, “Characteristics and possible formation mechanisms of severe storms in the outer rainbands of Typhoon Mujiga (1522),” J. Meteorol. Res. 31(3), 612–624 (2017). [CrossRef]

50. T. Inoue and H. Kamahori, “Statistical relationship between ISCCP cloud type and vertical relative humidity profile,” J. Meteorol. Soc. Japan 79(6), 1243–1256 (2001). [CrossRef]

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

52. T. S. L’Ecuyer, P. Gabriel, K. Leesman, S. J. Cooper, and G. L. Stephens, “Objective assessment of the information content of visible and infrared radiance measurements for cloud microphysical property retrievals over the global oceans. Part I: Liquid clouds,” J. Appl. Meteorol. Climatol. 45(1), 20–41 (2006). [CrossRef]

53. D. Kimes and J. Kirchner, “Radiative transfer model for heterogeneous 3-D scenes,” Appl. Opt. 21(22), 4119–4129 (1982). [CrossRef]

TIR Channels	Center wavelength (µm)	Bandwidth (µm)	Absorbing gas
B08	6.2383	0.8219	H₂O
B09	6.9395	0.4019	H₂O
B10	7.3471	0.1871	H₂O, N₂O, CH₄
B11	8.5905	0.3727	H₂O, N₂O, CH₄
B12	9.6347	0.3779	H₂O, CO₂, O₃
B13	10.4029	0.4189	H₂O, CO₂
B14	11.2432	0.6678	H₂O, CO₂
B15	12.3828	0.9656	H₂O, CO₂, O₃
B16	13.2844	0.5638	H₂O, CO₂, O₃

TIR Channels	B08	B09	B10	B11	B12	B13	B14	B15	B16
AMCKD + δ-4DDA	0.1254	0.1042	0.0595	0.0951	0.1118	0.0413	0.0533	0.0772	0.2090
AMCKD + δ-4 DISORT	0.1051	0.1017	0.0712	0.0917	0.1292	0.0400	0.0510	0.0768	0.2133
AMCKD + δ-32 DISORT	0.1045	0.0977	0.0688	0.0444	0.1321	0.0144	0.0276	0.0795	0.2216

Data name/source	Quantity	Spatial resolution	Temporal resolution
ERA-Interim	Atmospheric profiles	0.75°	6h
MODIS level 3 (MYD11C2)	Land surface temperature/ emissivity	0.05°	8-day mean
OISST	Sea surface temperature	0.25°	Daily
ICAS	Cloud property	2km	10min
Himawari-8 level 1	Brightness temperature	0.05°	10min

TIR Channels	Center wavelength (µm)	Bandwidth (µm)	Absorbing gas
B08	6.2383	0.8219	H₂O
B09	6.9395	0.4019	H₂O
B10	7.3471	0.1871	H₂O, N₂O, CH₄
B11	8.5905	0.3727	H₂O, N₂O, CH₄
B12	9.6347	0.3779	H₂O, CO₂, O₃
B13	10.4029	0.4189	H₂O, CO₂
B14	11.2432	0.6678	H₂O, CO₂
B15	12.3828	0.9656	H₂O, CO₂, O₃
B16	13.2844	0.5638	H₂O, CO₂, O₃

TIR Channels	B08	B09	B10	B11	B12	B13	B14	B15	B16
AMCKD + δ-4DDA	0.1254	0.1042	0.0595	0.0951	0.1118	0.0413	0.0533	0.0772	0.2090
AMCKD + δ-4 DISORT	0.1051	0.1017	0.0712	0.0917	0.1292	0.0400	0.0510	0.0768	0.2133
AMCKD + δ-32 DISORT	0.1045	0.0977	0.0688	0.0444	0.1321	0.0144	0.0276	0.0795	0.2216

Efficient radiative transfer model for thermal infrared brightness temperature simulation in cloudy atmospheres

Abstract

1. Introduction

2. Efficient radiative transfer model (ERTM)

2.1 Gaseous absorption by AMCKD

2.2 Cloud optical property parameterization

2.3 Radiative transfer calculation

2.4 ERTM integration

3. Results and discussion

4. Summary

Funding

Acknowledgments

Disclosures

References

Cited By

Figures (12)

Tables (4)

Equations (33)

Optics Express