1. Introduction

The radiative transfer community requires multifunctional, rapid, and accurate radiance and radiance gradient models for satellite data assimilation, sensor design and specification, calibration, and validation of remote sensing data, research, and education. To achieve these goals, the Joint Center for Satellite Data Assimilation (JCSDA) in the United States has developed a framework for a Community Radiative Transfer (RT) Model (CRTM). The framework is prepared based on the scientific and software requirements from the analysis system. The approach is known as modular program development (Weng et al. 2005); it breaks down the RT model into components, each of which is encapsulated in one or several program modules and can be developed independently of the others. The core module or engine of the CRTM is the solver of the radiative transfer equation.

Several solvers of the radiative transfer equation have been proposed for CRTM. Those solvers are independently developed either on site or off site of the Joint Center and can be plugged into the CRTM framework. A polarized Delta-4-stream model for thermal and microwave radiative transfer has been developed byLiou et al. (2005) at the University of California, Los Angeles. The successive order of interaction (SOI) radiative transfer model was developed byHeidinger et al. (2006) at the University of Wisconsin. Direct ordinate tangent linear radiative transfer (DOTLRT) has been developed byVoronovich et al. (2004) at the NOAA Environmental Technology Laboratory in Colorado. The vector discrete-ordinate radiative transfer method (VDISORT:Weng 1992) has been enhanced byWeng and Liu (2003) at the Office of Research and Applications of NOAA. This study improves the computation efficiency of the doubling–adding (DA) method. The principle of the doubling–adding method was stated byStokes (1862).Hansen (1971) applied the doubling–adding method to the multiple scattering of polarized light in the planetary atmosphere. A well-documented code of the doubling–adding method is available online and the description of the method is provided byEvans and Stephens (1991). The doubling–adding method has not been considered for operational retrievals or data assimilation owing to its huge demand on computational resources. The method has rather been used for accurate and detailed radiative transfer calculations in the field of research and education. The advanced doubling–adding (ADA) method is a recent improvement to the DA method that has been demonstrated to be a powerful tool for multiple-scattering calculations. ADA is about 60 times faster than DA on our PC and alleviates the problem of large demand on computational resources. Although we only introduce the RT models proposed for CRTM, it should be mentioned that there are many other very useful radiative transfer models such as the two-stream model (e.g.,Schmetz 1984), Eddington approximation (Kummerow 1993), and a linearized discrete ordinate radiative transfer model (Spurr et al. 2001). In the following, ADA is derived and the intercomparison with other models is carried out.

2. Radiative transfer solver

Using the discrete ordinate form fromWeng and Liu (2003), the nonpolarized radiative transfer equation in the infrared and microwave ranges can be written as

View Expanded

whereμ_i (noteμ_−i = −μ_i) andw_i are Gaussian quadrature points and weights, respectively. Hereμ_i andμ_−i represent the cosine of the viewing zenith angle in upward and downward directions, respectively. The repeated subscriptj involves a summation. The phase matrix elementsP(μ_i,μ_j) andP(μ_i,μ_−j) are the azimuth-averaged forward and backward parts, respectively;P(μ_i,μ_j) =P(μ_−i,μ_−j) andP(μ_−i,μ_j) =P(μ_i,μ_−j) due to the symmetry conditions of the phase function for spherical scatterers or for randomly oriented particles with a symmetric plane (Hovenier 1969). Then, Eq.(1) can be rewritten in a matrix–vector form as

View Expanded

whereα andβ areN ×N matrices and

View Expanded

View Expanded

whereδ_ij is the Kronecker delta. The subscripts u (d) indicate upward (downward) directions, respectively, and 𝘂 is anN ×N matrix that has nonzero elements in its diagonal such as

View Expanded

Ξ is a vector ofN elements as

View Expanded

For an infinitesimal optical depthδ₀, multiple scattering can be neglected and the reflection matrix can be expressed as (Plass et al. 1973)

View Expanded

and the transmission matrix can be written as

View Expanded

where 𝗘 is anN ×N unit matrix.

Using the doubling procedure fromVan de Hulst (1963), the reflection and transmission matrices for a finite optical depth (δ =δ_n = 2ⁿδ₀) can be computed by doubling the optical depth (i.e.,δ_i+1/δ_i = 2) recursively:

View Expanded

and

View Expanded

fori = 0,n − 1. We denote 𝗿(k) = 𝗿(δ_n) and 𝘁(k) = 𝘁(δ_n) for the reflection and transmission matrices ofkth layer.

There exist formulas for building the layer source functions (Heidinger et al. 2006) depending on the Planck function of the temperature at the top of the layer and the gradient of Planck function over the layer optical depth. However, the formulas are complicated and computationally expensive. In this study, we found a very simple and strict expression for the layer source function using the existing layer reflection and transmission matrices. For an atmospheric layer of an optical depthδ and having the top temperature ofT₁ and the bottom temperature ofT₂, the upward layer source function can be derived as (see theappendix)

View Expanded

and the downward source of the layer can be written as

View Expanded

where ϖ andg are the single scattering albedo and asymmetry factor of the layer, respectively. The new expressions for the layer source functions take very little extra computation time. The expression can also be applied for other radiative transfer models such as the Matrix Operator Method (Fischer and Grassl 1984).

Equations(6)–(8) give the layer reflection and transmission matrices as well as the source vectors at the upward and downward directions. For a planetary atmosphere, the atmosphere may be divided inton optically homogeneous layers. The optical properties (e.g., extinction coefficient, single scattering albedo, and phase matrix) are the same within each layer although the temperature may vary within the layer. The adding method is for integrating the surface and multiple atmospheric layers. The method was applied to flux calculation using a two-stream approximation (Selby and Clough 1988). The method was also used in radiance calculations with multiple scatterings using a two-stream approximation (Liu and Weng 2002). In the following, we briefly describe the methodology. We denote 𝗥_u(k) for reflection matrix and I_u(k) for radiance vector at the levelk in the upward direction, andk =n andk = 0 represent the surface level and the top of the atmosphere, respectively. The adding method starts from the surface without atmosphere. At the surface, 𝗥_u(n) is the surface reflection matrix andI_u(n) equals the surface emissivity vector multiplied by the Planck function at the surface temperature. The upward reflection matrix and radiance at the new level can be obtained by adding one layer from the present level:

View Expanded

View Expanded

The physical meaning of Eq.(9) is obvious. The first term on the right side of Eq.(9a) is the reflectance of the layer to be added. The second term on the right side of Eq.(9a) is the reflectance due to the radiation from the new level transmitted to and multiply reflected by the present level and then transmitted back to the new level. The three terms on the right side of Eq.(9b) represent the upward layer source, from the present level reflected layer downward source, and from the present level transmitted upward radiance, respectively. The upward radianceI_u at the top of the atmosphere can be obtained by looping the index fromk =n tok = 1 and adding the contribution from cosmic background radiance (Planck function at the temperature of 2.7 kelvin) vectorI_sky, that is,

View Expanded

Equations(6)–(10) give the necessary and sufficient formulas for the advanced doubling–adding method. It needs to be mentioned that the above procedure is for the upward radiance at the top of the atmosphere. It is sufficient for the satellite data assimilation. However, an additional loop from the top to the surface is necessary in order to obtain the vertical profiles of radiances at both upward and downward directions.

For the viewing angle departure from the angles at Gaussian quadrature points, an additional stream as an extra Gaussian quadrature point associated with an integration weight of zero may be inserted to haveN quadrature points in total in either upward or downward directions. For this case, the upward intensity vector will contain the upward solutions atN −1 quadrarture points and at a specified viewing angle. The result by inserting the additional stream is exactly the same as inserting the multiple scattering solutions atN − 1 quadrarture points back to the integration equation for the specified viewing angle [see Eqs. (24) and (25) ofStamnes et al. 1988]. However, the present procedure avoids using extra codes for the specified viewing angle so that it much simplified all forward, tangent linear, and adjoint coding.

3. Results

The advanced doubling–adding method is tested for forward model comparisons with VDISORT and the doubling–adding method under cloudy and clear-sky conditions and for intercomparisons with current operational model under clear-sky conditions. For gaseous absorption, this study uses the Optical Path Transmittance (OPTRAN) model to compute gaseous optical depth at the nadir. OPTRAN is a regression-based fast radiative transmittance model (Xiong and McMillin 2005). OPTRAN is used in the Global Data Assimilation System by NOAA. It may result in some errors for cloudy cases since OPTRAN is optimized under clear-sky conditions and the channel transmittance may not be monochromatic. Spherical scatterers are used for liquid and ice clouds in the infrared and microwave ranges. A lookup table contains the extinction coefficients, single scattering albedos, asymmetry factors, delta truncation factors for removing forward peaks, and expansion coefficients from Mie calculations (Simmer 1994). The lookup table contains also the extinction coefficients, single scattering albedos, and asymmetry factors for nonspherical particles of cirrus clouds in the infrared range (Yang et al. 1997). For surface emissivity and reflectivity, the land microwave emissivity model (Weng et al. 2001) is implemented into CRTM. Ocean microwave emissivity models (Hollinger 1971;English 1999) are a function of the sea surface temperature and wind speed and the center frequency of the microwave channel. An ocean microwave polarimetric model based on airborne microwave radiometer measurements (St. Germain et al. 2002) is applied (Liu and Weng 2003). For infrared land surface emissivity, the measured spectral emissivity for 24 surface types used in net heat flux studies (Carter et al. 2002) is adapted. Infrared water emissivity (van Delst and Wu 2000) is applied to account for the effect of surface wind speed.

The model intercomparisons are carried out between the doubling–adding model (Evans and Stephens 1991), VDISORT (Weng and Liu 2003), and the advanced doubling–adding method. We use the CRTM platform, which allows us to insert various solvers for radiative transfer calculations. Three solvers mentioned above share the same atmospheric optical data, the same surface emissivity and reflectivity, and the same Planck function for an atmosphere, surface, and cosmic background. The differences of results from the three solvers are purely from the differences in the solvers. For 24 000 simulations with various clear and cloudy cases, computation times on our personal computer are 1041, 29, and 17 s for the DA, VDISORT, and ADA models, respectively. ADA is about 1.7 times faster than VDISORT and 61 times faster than DA. The huge gain of ADA to DA is partly contributed by the efficiency of matrix and vector manipulation in FORTRAN 95 because the DA code is still in FORTRAN 77. The maximum difference of the simulated brightness temperatures between using the three solvers for Advanced Microwave Sounding Unit A (AMSU-A) channels and 281 selected Atmospheric Infrared Sounder (AIRS) channels is less than 0.01 K. The subset of AIRS data used in National Centers for Environmental Prediction (NCEP) data assimilation contains necessary information on atmospheric temperature and water vapor (Goldberg et al. 2003).Tables 1 and2 list the comparison of the brightness temperatures for AMSU-A water vapor channel at 23.8 GHz and the infrared window channel of AIRS at 10.88μm computed from the ADA, VDISORT, and DA methods, respectively. A profile containing temperature and water vapor as well as ozone from our test dataset for OPTRAN is selected. For the microwave calculation, a rain cloud having an effective particle size of 200μm and 0.5-mm rainwater content was put at 850 hPa. One layer ice cloud having the same effective particle size and 0.1-mm ice water path is located at 300 hPa. A wind speed of 5 m s⁻¹ over ocean is used. The maximum difference of the brightness temperature computed from the three models is less than 0.01 K (seeTable 1). For the infrared calculation, an ice cloud having an effective particle size of 20μm and 0.1-mm ice water path was located at 300 hPa, and a liquid water cloud at 850 hPa having an effective particle size of 10μm and 0.5 mm are chosen. The results computed from the three models agree very well (seeTable 2).

4. Summary

ADA can be used in satellite radiance assimilation and retrieval systems for infrared and microwave observations. The newly derived analytical formulas for the layer source function significantly improve the performance of the doubling–adding method. The formulas can also be used in other models, for example, the matrix operator method, to take account for the temperature variation within the layer. With the improvement, the doubling–adding method can be applied to calculate cloud radiance in operational satellite data assimilation. Since the coding of ADA for the forward model is relatively simple, tangent linear and adjoint parts are also simple; that is very beneficial for the code maintenance and consistency in the forward and adjoint satellite data assimilation.