Abstract
Satellite data assimilation requires a computationally fast and accurate radiative transfer model. Currently, three fast models are commonly used in the Numerical Weather Prediction models (NWP) for satellite data assimilation, including Radiative Transfer for TIROS Operational Vertical Sounder (RTTOV), Community Radiative Transfer Model (CRTM), and Advanced Radiative transfer Modeling System (ARMS). ARMS was initiated in 2018 and is now becoming the third pillar supporting many users in NWP and remote sensing fields. Its radiative transfer solvers (e.g. Doubling Adding method) is inherited from CRTM. In this study, we propose a Discrete Ordinate Adding Method (DOAM) to solve the radiative transfer equation including both solar and thermal source terms. In order to accelerate the DOAM computation, the single scattering approximation is used in the layer with an optical depth less than 10−8 or a single scattering albedo less than 10−10. From principles of invariance, the adding method is then applied to link the radiances between the layers. The accuracy of DOAM is evaluated through four benchmark cases. It is shown that the difference between DOAM and DIScrete Ordinate Radiative Transfer (DISORT) decreases with an increase of stream number. The relative bias of the 4-stream DOAM ranges from -5.03 % to 5.92 % in the triple layers of a visible wavelength case, while the maximum bias of the 8-stream DOAM is only about 1 %. The biases can be significantly reduced by the single scattering correction. Comparing to the visible case, the accuracy of the 4-stream DOAM is much higher in the thermal case with a maximum bias -1.69 %. Similar results are also shown in two multiple-layer cases. In the MacBook Pro (15-inch, 2018) laptop, the 2-stream DOAM only takes 1.68 seconds for calculating azimuthally independent radiance of 3000 profiles in the hyper-spectral oxygen A-band (wavelength ranges from 0.757 µm to 0.775 µm), while the 4-stream DOAM takes 4.06 seconds and the 16-stream DOAM takes 45.93 seconds. The time of the 2-, 4- and 16- stream DOAM are 0.86 seconds, 1.09 seconds and 4.34 seconds for calculating azimuthally averaged radiance. DISORT with 16 streams takes 1521.56 seconds and 127.64 seconds under the same condition. As a new solver, DOAM has been integrated into ARMS and is used to simulate the brightness temperatures at MicroWave Humidity Sounder (MWHS) as well as MicroWave Radiation Imager (MWRI) frequencies. The simulations by DOAM are compared to those by Doubling Adding method and accuracy of both solvers shows a general agreement. All the results show that the DOAM is accurate and computational efficient for applications in NWP data assimilation and satellite remote sensing.
© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
Radiative transfer models play a key role in the satellite remote sensing and data assimilation [1–7]. With the new generation of satellite instruments, the data volume is unprecedentedly huge and uses of the data require computationally fast and accurate observational operators [8–11]. Currently, many fast radiative transfer models, including Radiative Transfer for TIROS Operational Vertical Sounder (RTTOV) [12], Community Radiative Transfer Model (CRTM) [13,14], libRadtran [15,16], and Advanced Radiative transfer Modeling System (ARMS) [8,17], have been developed for radiance simulations. ARMS was initiated in 2018 and requested from user communities in China. Different from other fast radiative transfer models, ARMS includes some new modules which takes the scattering from spherical and non-spherical particles hydrometeors as well as the emission and scattering from rough ocean/land surfaces into account. ARMS also support all instruments onboard FengYun satellites. More details can be found in the earlier studies [8,17].
A fast radiative transfer model generally consists of a gaseous absorption module, an aerosol and cloud scattering module, a surface emission model and database, a radiative transfer solver as well as tangent linear and adjoints. In the gaseous absorption modules, the band gaseous absorption for a particular instrument is derived as the parameterized formulas with the coefficients which are generated from the regression against the Line-By-Line (LBL) radiative transfer computations [18]. Gaseous optical depth at each atmospheric layer is explicitly expressed as a function of temperature, pressure and gas content of the layer [19–21]. The optical properties of clouds are pre-calculated within a wide range of liquid/ice water content and effective radius. A look-up table is established based on the results [4,22,23]. Because of its complexity, the radiative transfer solver is difficult to be replaced by a regression module and its accuracy as well as computational efficiency in various conditions therefore determine the values of the observational operators.
Both DIScrete Ordinate Radiative Transfer (DISORT) [24] and Doubling Adding method (DA) [25–27] are popular radiative transfer solvers. In DISORT, the discrete ordinate method [24,28–31] is used to solve single-layer radiative transfer equation but without boundary conditions. The boundary conditions as well as continuity conditions at each layer interfaces are utilized to establish a system of linear equations. After using the numerical method to solve the system of linear equations, the intensity or the Stokes vector at each atmospheric level can be obtained. This treatment is conceptual complex and computational expensive, especially for a radiative transfer model with a high vertical resolution. In DA, atmosphere is divided into many sub-layers with very small optical depth (usually around 10$^{-8}$). The reflectance and transmittance of each sub-layer are calculated following a simple expression. The doubling and adding process is used to handle layer connection. While the DA concept is much more simple than DISORT, it is also computational expensive due to the large iteration loop in its algorithm. Zhang et al. (2013) [32] and (2016) [33] proposed a new four-stream solver, referred as 4DDA, to simulate irradiance in the visible and infrared spectrum. In the 4DDA, the four-stream discrete ordinate method is used to solve the single-layer radiative transfer equation. The adding method from four principles of invariance is applied to handle multiple layer radiative transfer. Therefore, the method overcomes the shortcoming of both DISORT and DA. Due to its high accuracy and efficiency, 4DDA has been used in climate simulations [34–36]. However, 4DDA only solves azimuthally averaged radiative transfer equation and it cannot be applied for radiance simulations.
In this study, the Discrete Ordinate Adding Method (DOAM) is developed as part of ARMS for Numerical Weather Prediction models (NWP) data assimilation and satellite remote sensing. In section 2, the basic azimuthally dependent radiative transfer equation is introduced and discretized. The multi-layer solar radiative transfer process is presented in section 3. The differences between solar and thermal radiative transfer are illustrated in section 4. Assessments of DOAM performance for two idealized triple-layer cases and two multiple-layer cases are made in section 5. Finally, DOAM is plugged into ARMS to simulate brightness temperatures at MicroWave Humidity Sounder (MWHS) as well as MicroWave Radiation Imager (MWRI) frequencies. The simulations between DOAM and DA are compared. A summary is presented in section 6.
2. Basic radiative transfer equation
The radiative transfer in a plane-parallel homogeneous atmospheric layer can be described by the following equation [37,38]:
According to the literature [37,38], the radiance can be divided into two parts as following:
where $I_\textrm {{dir}}(\tau ,\mu ,\phi )$ is the radiance of direct solar radiation which is refered as direct radiation thereafter and $I_\textrm {{dif}}(\tau ,\mu ,\phi )$ is the radiance of multiple scattering radiation, or diffuse radiation. The expression of direct radiation $I_\textrm {{dir}}$ isThe solar source is related to the solar flux from TOA while the thermal source of atmosphere can be derived from Planck function. Due to the difference of these two energy sources, the physical process of solar radiative transfer is also different from ones of thermal radiative transfer. Therefore, we should solve these two kinds of radiative transfer separately and then add radiance of two parts together [41].
3. Solar radiative transfer solution
The single-layer solar radiative transfer equation is the same as Eq. (5) but without the terms related to Planck function $B(T)$. We neglect superscript "$m$" in the following discussions, as the calculation process of each $I_\textrm {{dif}}^{m}(\tau ,\mu )$ is the same.
In the adding method, there are two kinds of diffuse radiation: one is related to the incident direct radiation, another is associated with the incident diffuse radiation [32,33,42]. The relationship between two diffuse radiation sources has already been analyzed in [43].
3.1 single layer solution
3.1.1 incident direct radiation
The single-layer solar radiative transfer related to incident direct radiation can be described by the following equations:
Then we can calculate the upward and downward output radiance matrix of incident direct radiation as follows
In order to accelerate the computational scheme, the single scattering approximation [44] is used to replace the discrete ordinate method in the layer with an optical depth less than $10^{-8}$ or a small single scattering albedo less than $10^{-10}$. The computational speed of the single scattering approximation is about 6 times faster than that of the discrete ordinate method. Following the approximation, multiple scattering is neglected and radiative transfer equation becomes
The solution of Eqs. (22)–(24) are
We note that there is a singularity in Eq. (27) when $\mu _{0}$ is very close to a Gaussian angle $\mu _{k}$. The L’Hopital’s rule is applied to handle it and the expression of $I_\textrm {{dif}}(\tau _{0},\mu _{-k})$ can be written as
3.1.2 incident diffuse radiation
According to the literature [32,33,42], the single-layer solar radiative transfer related to incident diffuse radiation can be described by the following equation:
In the $2N$ stream case, there are $N$ kinds of boundary conditions which refer to the incident diffuse radiation from $N$ Gaussian angles. The $j$-th boundary condition can be written as
Similar to the process of handling incident direct radiation, Eqs. (29) and (30) can be rewritten in matrix form
In the single scattering approximation, Eq. (37) can be simplified as
Following [32,33], we can define reflectance and transmittance matrix of incident diffuse radiation as
Before proceeding above single-layer solution, the Delta-scaling technology [45] is applied in order to improve the results of the solution.
3.2 Adding method from four principles of invariance
Adding method from four principles of invariance is used to handle the layer connection in this study. The method can provide an analytical solution of multiple-layer radiative transfer, which is different from the numerical solution in DISORT [24]. Compared to the numerical solution, the analytical solution generally more computational efficient. In order to detailly describe the process of the adding method, Let’s first consider a double-layer case. According to the literature [32,37,38,42], the four principles of invariance are
All variables have been defined in [32,37,38,42] and the subscripts "1", "2" and "1,2" represent the first layer, the second layer and the combination of two layers, respectively. By applying $2N$ Gaussian quadrature to handle the integration and then converting reflectance/transmittance to radiance, Eqs. (42)–(45) yields
The four principles of invariance can also be applied into radiative transfer related to incident diffuse radiation. After applying a similar procedure from Eqs. (42)–(45) to Eqs. (52)–(55), we obtain the reflectance and transmittance matrix of the combination of two layers [32,33,42]
Equations (55) and (58) can be extended from a double layer case to a multiple layer case through a path from TOA ($1$-th layer) to $k$-th layer [32,42]
In addition, if we applied DOAM to satellite remote sensing, only the upward radiance at TOA (e.g. $\boldsymbol {I}_{1,N}^{+}$) is needed and we just process Eqs. (62) and (63).
3.3 Radiance at viewing zenith angle
Calculating radiance at arbitrary viewing zenith angle is necessary for a radiative transfer scheme which is required in many applications. In order to consider the radiance at the viewing zenith angle $I_\textrm {{dif}}(\tau ,\mu _\textrm {sim})$ without affecting the results at other angles, the following expression is used to handle the integration in Eq. (8)
4. Difference between solar and thermal radiative transfer
The procedure of handling thermal radiative transfer is similar to ones of handling solar radiative transfer. Differences between these two procedures are described in the following subsections.
4.1 Difference in solving single-layer radiative transfer equation
Comparing to the solar radiative transfer equation, the thermal radiative transfer equation consists of the term related to the Planck function $B(T)$ but without term related to the incident solar flux $F_{0}$. In DOAM, we use a linear expression to describe the vertical variation of Planck function within a layer [46,47]
where $\beta =(B_{1}-B_{0})/\tau _{0}$. $B_{0}$ and $B_{1}$ is the Planck function at the top and bottom of the layer. Therefore, the solution of thermal radiative transfer equation is that [24]In the single scattering approximation, the solution of thermal radiative transfer equation is
As the thermal emission is azimuthally isotropic, we just need to solve the equation of $I_\textrm {{dif}}^{0}(\tau ,\mu )$ when handling thermal radiative transfer.
4.2 Difference in adding method
The four principles of invariance in thermal radiative transfer can be written as [33]
5. Comparisons and discussions
In this section, the accuracy and computational efficiency of DOAM are investigated in two idealized triple-layer cases and two multiple-layer cases. Results from DISORT [24] with 64 streams under the same optical condition are used as a benchmark. Then we apply DOAM into ARMS [17] and compare simulated brightness temperatures to those from DA [25,26] .
5.1 Triple layer cases
We first consider a triple layer case with the ground reflection of zero (blackbody) in the visible spectrum. The case consists of two Rayleigh scattering layers and a cloud layer in the middle. The Cloud C1 phase function [48] is used to describe the scattering in the cloud layer. The single scattering albedo and optical depth of both Rayleigh scattering layers are set to 0.9 and 0.3, respectively. The thermal emission is ignored. As DOAM is designed for both satellite data assimilation and remote sensing applications, we only consider the reflectance at TOA in the investigation
Figure 1 shows bias of reflectance at TOA calculated by the 4-stream (Fig. 1(a) and 1(b)) and the 8-stream (Fig. 1(c) and 1(d)) DOAM at different viewing zenith angles. The optical depth of cloud is set to 10 in Figs. 1(a), 1(c) and varies from 0.1 to 100 in Figs. 1(b), 1(d). The single scattering albedo of cloud varies from 0.9 to 0.99 in Figs. 1(a), 1(c) and is fixed at 0.93 in Figs. 1(b), 1(d). The cosine of solar zenith angle is 0.5. The solar and viewing azimuth angle are $0^{\circ }$ and $60^{\circ }$, respectively. The bias increases with the decrease of single scattering albedo in both the 4-stream and the 8-stream DOAM. The relative bias between the 4-stream DOAM and the benchmark can reach -4.85 % (Fig. 1(a)), while the maximum bias is only -1.12 % after applying the 8-stream DOAM (Fig. 1(c)). The bias also increases with the optical depth. The maximum bias of the 4-stream DOAM is -4.49 % (Fig. 1(b)) while the 8-stream DOAM reduce this bias to -1.10 % (Fig. 1(d)). After applying the single scattering correction [44], the maximum bias of the 4-stream DOAM reduce to -2.29 % in the Fig. 1(a) case, -2.42 % in the Fig. 1(b) case. The maximum bias of the 8-stream DOAM reduce to about 0.1 % in the both Fig. 1(c) and Fig. 1(d) cases.

Fig. 1. The relative bias of reflectance at TOA calculated by the 4-stream (a-b) and the 8-stream (c-d) DOAM. The horizontal coordinate represent the cosine of viewing zenith angle, while the vertical coordinate refers to optical depth of the cloud layer $\tau _{0}$ (b, d) or single scattering albedo of the cloud layer $\omega$ (a, c). The optical depth of cloud is set to 10 in (a), (c) and varies from 0.1 to 100 in (b), (d). The single scattering albedo of cloud varies from 0.9 to 0.99 in (a), (c) and is fixed at 0.93 in (b), (d). The cosine of solar zenith angle is 0.5. The solar and viewing azimuth angle are $0^{\circ }$ and $60^{\circ }$, respectively. Results from DISORT with 64 streams are used as benchmark.
Figure 2 shows the bias at different solar zenith angles. The cosine of viewing zenith angle, the solar azimuth angle and viewing azimuth angle are set to $0.65$, $0^{\circ }$ and $0^{\circ }$, respectively. Compared to the viewing zenith angle, the solar zenith angle has a greater impact on the result of reflectance at TOA. With the variation of single scattering albedo and solar zenith angle, the relative bias between the 4-stream DOAM and the benchmark ranges from -5.03 % to 5.92 % (Fig. 2(a)), while the bias only ranges from -1.34 % to 1.60 % (Fig. 2(c)) after applying the 8-stream DOAM. Simliar results are shown in the Fig. 2(b) and 2(d). The bias of the 4-stream DOAM ranges from -3.08 % to 4.44 % (Fig. 2(b)). Again, the 8-stream DOAM is much more accurate than the 4-stream DOAM (Fig. 2(d)). After applying the single scattering correction [44], the relative bias between the 4-stream DOAM and the benchmark ranges from -2.50 % to 3.21 % and -2.44 % to 2.52 % in the Fig. 2(a) and Fig. 2(c) case, respectively. The maximum bias of the 8-stream DOAM is less than 0.3 % in the both Fig. 2(b) and Fig. 2(d) cases.

Fig. 2. Same as Fig. 1 but the horizontal coordinate represent the cosine of solar zenith angle. The cosine of viewing zenith angle, the solar azimuth angle and viewing azimuth angle are set to $0.65$, $0^{\circ }$ and $0^{\circ }$, respectively.
Let’s consider another triple layer case which refers to infrared radiative transfer. In the case, the thermal emission is included and the solar energy is ignored. The case consists of two pure absorption layers and a cloud layer in the middle. The Henyey-Greenstein phase function [49] with an asymmetry factor $g=0.9$ is introduced to describe the cloud scattering. The optical depth of pure absorption layer is set to $1$. The effective emissivity is used to replace the reflectance in the investigation
The incident thermal emission from the surface is set to $0$ with $B(T_s)$=0, $\varepsilon _\textrm {surface}(\mu )=1$.Figure 3 shows bias of effective emissivity at TOA calculated by the 2-stream (a-b) and the 4-stream (c-d) DOAM at different viewing zenith angle. The optical depth of cloud is set to 1 in Figs. 3(a), 3(c) and varies from 0.1 to 100 in Figs. 3(b), 3(d). The single scattering albedo of cloud varies from 0. to 0.7 in Figs. 3(a), 3(c) and is fixed at 0.7 in Figs. 3(b), 3(d). The bias of the 2-stream DOAM increases with single scattering albedo and can be up to -3.88 % when $\omega =0.7$. The 4-stream DOAM can reduce this bias to -1.06 %. Similar results are shown in the Fig. 3(b) and 3(d). The maximum bias of the 2-stream DOAM is -4.04 % (Fig. 3(b)) while the maximum bias of the 4-stream DOAM is -1.69 % (Fig. 3(d)).

Fig. 3. The relative bias of effective emissivity calculated by the 2-stream (Fig. 3(a) and 3(b)) and the 4-stream (Fig. 3(c) and 3(d)) DOAM. The meanings of the horizontal/vertical coordinates are same as ones in Fig. 1. The optical depth of cloud is set to 1 in (a), (c) and varies from 0.1 to 100 in (b), (d). The single scattering albedo of cloud varies from 0 to 0.7 in (a), (c) and is fixed at 0.7 in (b), (d). Results from DISORT with 64 stream are used as benchmark.
5.2 Multiple layer cases
A Rayleigh scattering case and an aerosol case are used in this subsection. Optical properties of the Rayleigh scattering case is calculated in the hyper-spectral oxygen A-band. It contains 3000 spectrum from 0.757 $\mu$m to 0.775 $\mu$m. U.S. standard atmosphere profile [50] which contains 35 atmospheric layers is applied and we only consider Rayleigh scattering in this case. We use blackbody lower boundary surface and ignore the thermal emission. The cosine of solar zenith angle and viewing zenith angle are 0.51 and 0.4, respectively. Both solar azimuth angle and viewing azimuth angle are set to 0$^{\circ }$. Figure 4(a) show the benchmark of reflectance at TOA calculated by DISORT with 64 stream and Fig. 4(b) show the relative error between the 2-, 4-, 16- DOAM and benchmark. The bias of DOAM ranges from -30.17 % to -27.28 % in the 2-stream case while the maximum bias of the 4-stream DOAM is -2.54 %. The 16-stream DOAM can reduce this bias to 0.21 %.

Fig. 4. The benchmark of reflectance at TOA calculated by DISORT with 64 stream in the Rayleigh scattering case (a). The relative error between 2-stream, 4-stream, 16-stream DOAM and benchmark (b).
The optical properties of the aerosol case are calculated by Optical Properties of Aerosols and Clouds (OPAC) [51]. We choose the continental polluted aerosol with 50 % relative humidity. Henyey-Greenstein phase function [49] and blackbody lower boundary surface are used here for illustrative purposes. We also ignore the thermal emission in this case. The solar and viewing angles are same as ones in the Rayleigh scattering case. Results are shown in the Fig. 5. Figure 5(a) show the benchmark of reflectance at TOA calculated by DISORT with 64 stream and Fig. 5(b) show the Root Mean Square Error (RMSE) of the 2-16 stream DOAM. Similar to the Rayleigh scattering case, the RMSE of DOAM is gradually reduced with increasing stream numbers. The RMSE of the 2-stream DOAM is 0.0118 while the RMSE of the 4-stream DOAM is 0.0067. The RMSE of the 16-stream DOAM is only 0.00007.

Fig. 5. The benchmark of reflectance at TOA calculated by DISORT with 64 stream in the aerosol case (a). The RMSE of DOAM in terms of number of streams (b).

Fig. 6. The computational time of DOAM in terms of number of streams. The green line represent the case for calculating azimuthally dependent radiance while the red line refers to the case for calculating azimuthally averaged radiance.
Computational efficiency is another key issue for any radiative transfer scheme which is applied into satellite data assimilation and remote sensing. We investigate the computational time of DOAM for azimuthally dependent radiance and azimuthally averaged radiance calculation. The azimuthally dependent radiance is necessary for visable spectrum, while infrared/microwave spectrum only need azimuthally averaged radiance. Figure 6 shows two types of computational time in terms of number of stream. Above Rayleigh scattering case is also used here. The calculation is performed on MacBook Pro (15-inch, 2018) with the 2.9 GHz 6-Core Intel Core i9 processor and the 32GB 2400 MHz DDR4 memory. Intel Fortran (version: 19.0.0.117) is applied for programming. In calculating the azimuthally dependent radiance, the 2-stream DOAM only takes 1.68 seconds when calculating the reflectance of 3000 profiles while the 4-stream DOAM takes 4.06 seconds and the 16-stream DOAM takes 45.93 seconds. In the azimuthally averaged radiance calculation, the time of the 2-, 4- and 16- stream DOAM are 0.86 seconds, 1.09 seconds and 4.34 seconds, respectively. DISORT with 16 streams takes 1521.56 seconds and 127.64 seconds under the same condition. A higher computational efficiency is attributed to the following three factors: 1. Adding method is an analytical method which is generally much faster than the numerical approach used in DISORT [24] and other similar solvers. 2. DOAM uses the single scattering approximation to replace the discrete ordinate method in the layer with an optical depth less than $10^{-8}$ or a single scattering albedo less than $10^{-10}$. The single scattering approximation is about 6 times faster than the discrete ordinate method in above MacBook Pro laptop. 3. Also, the DOAM code is coded with a standard FORTRAN-90. The code performs some matrix operations which is much faster than the traditional looping calculation.
All the results show that DOAM is accurate and computational efficient. It is very suitable to be used in many satellite data applications.
5.3 DOAM Application in ARMS
ARMS is a fast radiative transfer model developed by Chinese Meteorological Administration. In ARMS, the gaseous optical depth is calculated following [52,53] and the look-up table of cloud/aerosol optical properties is established based on Lorenz-Mie theory, Invariant imbedding T-matrix method [54,55] and parameterization of [22,23]. The superiority of ARMS has already been proved in [8,17].
In this section, we introduce DOAM into ARMS and simulate the brightness temperature at MWHS as well as MWRI frequencies. The results are compared to ones simulated by DA [25]. We use a U.S. standard atmosphere profile [50] and position a water cloud from 814.8 hPa to 450.8 hPa. The liquid water content and effective radius of cloud are set to 0.5 kg/m$^{2}$ and 50 $\mu$m, respectively. The angle-dependent ground emissivity is provided by the FAST microwave Emissivity Model (FASTEM) [56] under sea surface temperature of 275 K, southerly wind of 5 m/s and water salinity of 33 ppmv. The viewing zenith angle and viewing azimuth angle are set to 30$^{\circ }$ and 26.37$^{\circ }$. The simulated brightness temperatures at MWHS and MWRI frequencies are shown in the Table 1 and Table 2, respectively. As DOAM and DA apply different expressions to represent the vertical variation of Planck function within a layer (e.g. Eq. (72)), the maximum difference between two schemes is 0.523 K at MWHS frequencies. The maximum differences can reduce to 0.0138 K at MWRI frequencies.

Table 1. The Brightness Temperature (Unit: K) Simulated by DOAM as well as Doubling Adding method at Each Channel of MWHS. The Brightness Temperature Differences between DOAM and Doubling Adding method are Listed in Parentheses.

Table 2. Same as Table 1 but for MWRI Channel.
6. Summary and conclusions
The Discrete Ordinate Adding Method (DOAM) is proposed in this study to solve the scalar radiative transfer equation. In DOAM, source terms related to solar energy and thermal emission are considered separately. The discrete ordinate method is introduced to solve the single-layer radiative transfer equation, while it is replaced by the single scattering approximation in the layer with an optical depth less than $10^{-8}$ or a single scattering albedo less than $10^{-10}$. The adding method which is based on four principles of invariance is applied to extend the single-layer solution to the multiple-layer inhomogeneous atmosphere.
Two idealized triple layer cases and two multiple layer cases are investigated to benchmark the accuracy of DOAM. The results of DISORT with 64 streams are used as a truth. In the triple layer case of visible spectrum, the relative bias of reflectance at TOA calculated by the 4-stream DOAM ranges from -5.03 % to 5.92 % while the maximum bias of the 8-stream DOAM is only about 1 % within a wide range of single scattering albedo, optical depth, viewing zenith angle and solar zenith angle. Above results can be significantly improved by the single scattering correction. Compared to the visible spectrum, the accuracy of the 4-stream DOAM is much higher in the infrared spectrum. The maximum bias of the 4-stream DOAM is only -1.69 % in the triple layer case of infrared spectrum. Similar results can also be found in two multiple layer cases. In the Rayleigh scattering case, the relative bias of the 2-stream DOAM ranges from -30.17 % to -27.28 % while the maximum bias of the 16-stream DOAM is only 0.21 %. In the aerosol case, the RMSE of DOAM is gradually reduced with increasing stream numbers. The RMSE of the 2-, 4-, 16- stream DOAM is 0.0118, 0.0067 and 0.00007, respectively.
Computational efficiency is another key issue for any radiative transfer scheme which is applied into satellite data assimilation and remote sensing. In our small Mac laptop, the 2-stream DOAM only takes 1.68 seconds when calculating azimuthally independent radiance of 3000 profiles while the 4-stream DOAM takes 4.06 seconds and the 16-stream DOAM takes 45.93 seconds. In the azimuthally averaged radiance calculation, the time of the 2-, 4- and 16- stream DOAM are 0.86 seconds, 1.09 seconds and 4.34 seconds, respectively. DISORT with 16 streams takes 1521.56 seconds and 127.64 seconds under the same condition.
We also introduce DOAM into ARMS and simulate brightness temperatures at MWHS as well as MWRI frequencies. Since DOAM and DA apply different expressions to represent the vertical variation of Planck function within a layer, the two schemes can result in the difference of brightness temperature on the order of 0.5 K at MWHS frequencies. The differences are much smaller than 0.0138 K at MWRI frequencies. All the results show that DOAM is accurate and computational efficient. It is highly suitable to be used in NWP data assimilation and satellite remote sensing. It should be mentioned that DOAM can also be applied into irradiance calculations in climate models through integrating the hemispheric radiances into the flux. With an increasing technology in space measuring the radiative properties through polarization instruments [57–59], we will extend the DOAM to handle the full polarimetric atmospheric radiative transfer in the near future. For users who is interested in using DOAM, its codes are available by contacting shiyn@cma.gov.cn after signing an user license agreement.
Funding
National Key Research and Development Program of China (2018YFC1506501); National Natural Science Foundation of China (41675030); Chinese Academy of Meteorological Sciences (2019Z001); .
Acknowledgments
Authors appreciate three anonymous reviewers for giving constructive comments. Linear Algebra PACKage (LAPACK) (available at http://www.netlib.org/lapack/) provide routines for solving systems of linear equations and eigenvalue problems in DOAM.
Disclosures
The authors declare that there are no conflicts of interest related to this article.
References
1. H. Hu, F. Weng, Y. Han, and Y. Duan, “Remote sensing of tropical cyclone thermal structure from satellite microwave sounding instruments: Impacts of background profiles on retrievals,” J. Meteorol. Res. 33(1), 89–103 (2019). [CrossRef]
2. K. L. Chan, M. Wiegner, J. van Geffen, I. De Smedt, C. Alberti, Z. Cheng, S. Ye, and M. Wenig, “MAX-DOAS measurements of tropospheric NO2 and HCHO in Munich and the comparison to OMI and TROPOMI satellite observations,” Atmos. Meas. Tech. 13(8), 4499–4520 (2020). [CrossRef]
3. K. L. Chan, P. Valks, S. Slijkhuis, C. Köhler, and D. Loyola, “Total column water vapor retrieval for Global Ozone Monitoring Experience-2 (GOME-2) visible blue observations,” Atmos. Meas. Tech. 13(8), 4169–4193 (2020). [CrossRef]
4. W. Li, F. Zhang, Y.-N. Shi, H. Iwabuchi, M. Zhu, J. Li, W. Han, H. Letu, and H. Ishimoto, “Efficient radiative transfer model for thermal infrared brightness temperature simulation in cloudy atmospheres,” Opt. Express 28(18), 25730–25749 (2020). [CrossRef]
5. W. Kan, Y. Han, F. Weng, L. Guan, and S. Gu, “Multisource assessments of the FengYun-3D MicroWave Humidity Sounder (MWHS) on-orbit performance,” IEEE Trans. Geosci. Remote Sensing 58(10), 7258–7268 (2020). [CrossRef]
6. A. Doicu, D. S. Efremenko, and T. Trautmann, “A spectral acceleration approach for the spherical harmonics discrete ordinate method,” Remote Sens. 12(22), 3703 (2020). [CrossRef]
7. P. Hedelt, D. S. Efremenko, D. G. Loyola, R. Spurr, and L. Clarisse, “Sulfur dioxide layer height retrieval from Sentinel-5 Precursor/TROPOMI using FP_ILM,” Atmos. Meas. Tech. 12(10), 5503–5517 (2019). [CrossRef]
8. F. Weng, X. Yu, Y. Duan, J. Yang, and J. Wang, “Advanced Radiative transfer Modeling System (ARMS): A new-generation satellite observation operator developed for numerical weather prediction and remote sensing applications,” Adv. Atmos. Sci. 37(2), 131–136 (2020). [CrossRef]
9. J. Yang, Z. Zhang, C. Wei, F. Lu, and Q. Guo, “Introducing the new generation of Chinese geostationary weather satellites, FengYun-4,” Bull. Am. Meteorol. Soc. 98(8), 1637–1658 (2017). [CrossRef]
10. J. Veefkind, I. Aben, K. McMullan, H. Förster, J. de Vries, G. Otter, J. Claas, H. Eskes, J. de Haan, Q. Kleipool, M. van Weele, O. Hasekamp, R. Hoogeveen, J. Landgraf, R. Snel, P. Tol, P. Ingmann, R. Voors, B. Kruizinga, R. Vink, H. Visser, and P. Levelt, “TROPOMI on the ESA Sentinel-5 Precursor: A GMES mission for global observations of the atmospheric composition for climate, air quality and ozone layer applications,” Remote. Sens. Environ. 120, 70–83 (2012). [CrossRef]
11. S. English, C. Prigent, B. Johnson, S. Yueh, E. Dinnat, J. Boutin, S. Newman, M. Anguelova, T. Meissner, M. Kazumori, F. Weng, A. Supply, L. Kilic, M. Bettenhausen, A. Stoffelen, and C. Accadia, “Reference-quality emission and backscatter modeling for the ocean,” Bull. Am. Meteorol. Soc. 101(10), E1593–E1601 (2020). [CrossRef]
12. R. Saunders, J. Hocking, E. Turner, P. Rayer, D. Rundle, P. Brunel, J. Vidot, P. Roquet, M. Matricardi, A. Geer, N. Bormann, and C. Lupu, “An update on the RTTOV fast radiative transfer model (currently at version 12),” Geosci. Model Dev. 11(7), 2717–2737 (2018). [CrossRef]
13. F. Weng and Q. Liu, “Satellite data assimilation in numerical weather prediction models. part I: Forward radiative transfer and Jacobian modeling in cloudy atmospheres,” J. Atmos. Sci. 60(21), 2633–2646 (2003). [CrossRef]
14. Y. Han, P. van Delst, Q. Liu, F. Weng, Y. Banghua, R. Treadon, and J. Derber, “JCSDA Community Radiative Transfer Model (CRTM) - Version 1,” Tech. rep., U.S. Department of Commerce, National Oceanic and Atmospheric Administration, National Environmental Satellite, Data, and Information Service (2006).
15. B. Mayer and A. Kylling, “Technical note: The libRadtran software package for radiative transfer calculations - description and examples of use,” Atmos. Chem. Phys. 5(7), 1855–1877 (2005). [CrossRef]
16. C. Emde, R. Buras-Schnell, A. Kylling, B. Mayer, J. Gasteiger, U. Hamann, J. Kylling, B. Richter, C. Pause, T. Dowling, and L. Bugliaro, “The libRadtran software package for radiative transfer calculations (version 2.0.1),” Geosci. Model Dev. 9(5), 1647–1672 (2016). [CrossRef]
17. J. Yang, S. Ding, P. Dong, L. Bi, and B. Yi, “Advanced Radiative transfer Modeling System developed for satellite data assimilation and remote sensing applications,” J. Quant. Spectrosc. Radiat. Transfer 251, 107043 (2020). [CrossRef]
18. S. Clough, M. Shephard, E. Mlawer, J. Delamere, M. Iacono, K. Cady-Pereira, S. Boukabara, and P. Brown, “Atmospheric radiative transfer modeling: a summary of the AER codes,” J. Quant. Spectrosc. Radiat. Transfer 91(2), 233–244 (2005). [CrossRef]
19. M. Matricardi, The generation of RTTOV regression coefficients for IASI and AIRS using a new profile training set and a new line-by-line database (ECMWF, 2008).
20. L. L. Strow, S. E. Hannon, S. D. Souza-Machado, H. E. Motteler, and D. Tobin, “An overview of the airs radiative transfer model,” IEEE Transactions on Geoscience and Remote Sensing 41(2), 303–313 (2003). [CrossRef]
21. L. M. McMillin, L. J. Crone, and T. J. Kleespies, “Atmospheric transmittance of an absorbing gas. 5. Improvements to the OPTRAN approach,” Appl. Opt. 34(36), 8396–8399 (1995). [CrossRef]
22. B. Yi, X. Huang, P. Yang, B. A. Baum, and G. W. Kattawar, “Considering polarization in MODIS-based cloud property retrievals by using a vector radiative transfer code,” J. Quant. Spectrosc. Radiat. Transfer 146, 540–548 (2014). [CrossRef]
23. B. Yi, P. Yang, Q. Liu, P. van Delst, S.-A. Boukabara, and F. Weng, “Improvements on the ice cloud modeling capabilities of the Community Radiative Transfer Model,” J. Geophys. Res.: Atmos. 121, 13577–13590 (2016). [CrossRef]
24. 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]
25. Q. Liu and F. Weng, “Advanced doubling adding method for radiative transfer in planetary atmospheres,” J. Atmos. Sci. 63(12), 3459–3465 (2006). [CrossRef]
26. K. Evans and G. Stephens, “A new polarized atmospheric radiative transfer model,” J. Quant. Spectrosc. Radiat. Transfer 46(5), 413–423 (1991). [CrossRef]
27. X. He, Y. Bai, Q. Zhu, and F. Gong, “A vector radiative transfer model of coupled ocean–atmosphere system using matrix-operator method for rough sea-surface,” J. Quant. Spectrosc. Radiat. Transfer 111(10), 1426–1448 (2010). [CrossRef]
28. F. Weng, “A multi-layer discrete-ordinate method for vector radiative transfer in a vertically-inhomogeneous, emitting and scattering atmosphere—I. theory,” J. Quant. Spectrosc. Radiat. Transfer 47(1), 19–33 (1992). [CrossRef]
29. F. Schulz, K. Stamnes, and F. Weng, “VDISORT: An improved and generalized discrete ordinate method for polarized (vector) radiative transfer,” J. Quant. Spectrosc. Radiat. Transfer 61(1), 105–122 (1999). [CrossRef]
30. R. Spurr, T. Kurosu, and K. Chance, “A linearized discrete ordinate radiative transfer model for atmospheric remote-sensing retrieval,” J. Quant. Spectrosc. Radiat. Transfer 68(6), 689–735 (2001). [CrossRef]
31. A. Doicu and T. Trautmann, “Discrete-ordinate method with matrix exponential for a pseudo-spherical atmosphere: Scalar case,” J. Quant. Spectrosc. Radiat. Transfer 110(1-2), 146–158 (2009). [CrossRef]
32. F. Zhang, Z. Shen, J. Li, X. Zhou, and L. Ma, “Analytical delta-four-stream doubling-adding method for radiative transfer parameterizations,” J. Atmos. Sci. 70(3), 794–808 (2013). [CrossRef]
33. 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]
34. H. Zhang, Z. Wang, F. Zhang, and X. Jing, “Impact of four-stream radiative transfer algorithm on aerosol direct radiative effect and forcing,” Int. J. Climatol. 35(14), 4318–4328 (2015). [CrossRef]
35. X. Jing, H. Zhang, J. Peng, J. Li, and H. W. Barker, “Cloud overlapping parameter obtained from CloudSat/CALIPSO dataset and its application in AGCM with McICA scheme,” Atmos. Res. 170, 52–65 (2016). [CrossRef]
36. J. Li, H. Barker, P. Yang, and B. Yi, “On the aerosol and cloud phase function expansion moments for radiative transfer simulations,” J. Geophys. Res.: Atmos. 120, 12128–12142 (2015). [CrossRef]
37. K. Liou, Radiation and Cloud Processes in the Atmosphere: Theory, Observation and Modeling (Oxford University, 1992).
38. S. Chandrasekhar, Radiative transfer (Dover Publications, 1960).
39. J. F. de Haan, P. B. Bosma, and J. W. Hovenier, “The adding method for multiple scattering calculations of polarized light,” Astron. Astrophys. 183, 371–391 (1987).
40. K.-N. Liou, Q. Fu, and T. P. Ackerman, “A simple formulation of the delta-four-stream approximation for radiative transfer parameterizations,” J. Atmos. Sci. 45(13), 1940–1948 (1988). [CrossRef]
41. F. Zhang, K. Wu, J. Li, H. Zhang, and S. Hu, “Radiative transfer in the region with solar and infrared spectra overlap,” J. Quant. Spectrosc. Radiat. Transfer 219, 366–378 (2018). [CrossRef]
42. F. Zhang and J. Li, “Doubling-adding method for delta-four-stream spherical harmonic expansion approximation in radiative transfer parameterization,” J. Atmos. Sci. 70(10), 3084–3101 (2013). [CrossRef]
43. F. Zhang, J. Li, H. Zhang, L. Ma, and X. Zhou, “On the relationship between direct and diffuse radiation,” J. Quant. Spectrosc. Radiat. Transfer 115, 60–65 (2013). [CrossRef]
44. T. Nakajima and M. Tanaka, “Algorithms for radiative intensity calculations in moderately thick atmospheres using a truncation approximation,” J. Quant. Spectrosc. Radiat. Transfer 40(1), 51–69 (1988). [CrossRef]
45. W. J. 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]
46. O. B. Toon, C. P. McKay, T. P. Ackerman, and K. Santhanam, “Rapid calculation of radiative heating rates and photodissociation rates in inhomogeneous multiple scattering atmospheres,” J. Geophys. Res.: Atmos. 94(D13), 16287–16301 (1989). [CrossRef]
47. 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]
48. R. D. M. Garcia and C. E. Siewert, “Benchmark results in radiative transfer,” Transp. Theory Stat. Phys. 14(4), 437–483 (1985). [CrossRef]
49. L. G. Henyey and J. L. Greenstein, “Diffuse radiation in the galaxy,” The Astrophysical J. 93, 70–83 (1941). [CrossRef]
50. R. McClatchey, R. Fenn, J. A. Selby, F. Volz, and J. Garing, Optical Properties of the Atmosphere (Third Edition) (Air Force Cambridge Research Laboratories, 1972).
51. M. Hess, P. Koepke, and I. Schult, “Optical Properties of Aerosols and Clouds: The software package OPAC,” Bull. Am. Meteorol. Soc. 79(5), 831–844 (1998). [CrossRef]
52. W. Kan, P. Dong, Z. Zhang, and S. Ding, “Development and application of ARMS fast transmittance model for GIIRS data,” J. Quant. Spectrosc. Radiat. Transfer 251, 107025 (2020). [CrossRef]
53. E. Turner, P. Rayer, and R. Saunders, “AMSUTRAN: A microwave transmittance code for satellite remote sensing,” J. Quant. Spectrosc. Radiat. Transfer 227, 117–129 (2019). [CrossRef]
54. L. Bi, P. Yang, C. Liu, B. Yi, B. A. Baum, B. van Diedenhoven, and H. Iwabuchi, “Assessment of the accuracy of the conventional ray-tracing technique: Implications in remote sensing and radiative transfer involving ice clouds,” J. Quant. Spectrosc. Radiat. Transfer 146, 158–174 (2014). [CrossRef]
55. L. Bi, W. Lin, D. Liu, and K. Zhang, “Assessing the depolarization capabilities of nonspherical particles in a super-ellipsoidal shape space,” Opt. Express 26(2), 1726–1742 (2018). [CrossRef]
56. Q. Liu, F. Weng, and S. J. English, “An improved fast microwave water emissivity model,” IEEE Trans. Geosci. Remote Sensing 49(4), 1238–1250 (2011). [CrossRef]
57. P. J. Werdell, M. J. Behrenfeld, P. S. Bontempi, E. Boss, B. Cairns, G. T. Davis, B. A. Franz, U. B. Gliese, E. T. Gorman, O. Hasekamp, K. D. Knobelspiesse, A. Mannino, J. V. Martins, C. R. McClain, G. Meister, and L. A. Remer, “The Plankton, Aerosol, Cloud, ocean Ecosystem mission: status, science, advances,” Bull. Am. Meteorol. Soc. 100(9), 1775–1794 (2019). [CrossRef]
58. Y. Liu and D. J. Diner, “Multi-angle imager for aerosols: A satellite investigation to benefit public health,” Public Health Rep. 132(1), 14–17 (2017). [CrossRef]
59. T. Marbach, J. Riedi, A. Lacan, and P. Schlüssel, “The 3MI mission: Multi-viewing-channel-polarisation imager of the EUMETSAT polar system: second generation (EPS-SG) dedicated to aerosol and cloud monitoring,” in Polarization Science and Remote Sensing VII, J. A. Shaw and D. A. LeMaster, eds. (SPIE, 2015).