Development of Gas Absorption Tables and an Atmospheric Radiative Transfer Package for Applications Using the Advanced Himawari Imager and Change Studies by Using the

We developed an atmospheric gas absorption table for the Advanced Himawari Imager (AHI) based on the correlated k -distribution (CKD) method with an optimization method, which was used to determine quadrature weights and abscissas. We incorporated the table and band information of the AHI into a multi-purpose atmospheric radiative transfer package, RSTAR. We updated the package so that users could easily specify the satellite and band number. Use of this update made it possible to carry out calculations rapidly and accurately using the optimized CKD method. RSTAR is easy for beginners to use and facilitates comparison of results. Cloud retrieval tests using different numbers of quadrature points showed that cloud retrievals could be significantly affected by the CKD model’s accuracy.


Introduction
Radiative transfer code is fundamental to analyzing satellite observations. Several modern satellite remote sensing algorithms (Watts et al. 2011;Poulsen et al. 2012;Sourdeval et al. 2015Sourdeval et al. , 2016 are continually being developed using an optimal estimation method (Rodgers 2000). Iwabuchi et al. (2016Iwabuchi et al. ( , 2018 developed the Integrated Cloud Analysis System (ICAS), which is an optimal estimation-based algorithm used to investigate global distributions of cloud properties.
In addition, Hashimoto and Nakajima (2017) developed the Multi-Wavelength and multi-Pixel Method (MWPM), which is also based on an optimal estimation approach, to retrieve aerosol optical properties over heterogeneous surfaces. Although computational efficiency is not very important for creating look-up tables (often used to analyze satellite imagery data), these forward models needed radiative transfer code to rapidly and accurately calculate radiative transfer many times and in many cases.
Several atmospheric radiative transfer codes have been developed for satellite analysis. For example, Radiative Transfer for TOVS (RTTOV; Saunders et al. 1999) and Second Simulation of the Satellite Signal in the Solar Spectrum (6S; Vermote et al. 1997) are well known open source codes. In Japan, the authors manage the OpenCLASTR (Open Clustered Libraries for Atmospheric Science and Transfer of Radiation) project, from which packages and libraries for atmospheric radiation are developed and distributed. The System for Transfer of Atmospheric Radiation (STAR) series, which plays the main role in this project, contains the packages RSTAR Tanaka 1986, 1988) for radiance calculations, PSTAR (Ota et al. 2010) for polarized radiance calculations, FSTAR for radiative flux calculations, and MCSTAR (Okata et al. 2017) for three-dimensional Monte Carlo calculations. RSTAR is a well-known radiative transfer package that has been introduced in several satellite retrieval algorithms. RSTAR was created in 1988 and has been under continuous development, with its latest version being version 7. However, some algorithms are incompatible with recent satellite sensors. RSTAR's developmental policy is to provide a general and versatile package, but such a package requires complicated and detailed specifications to be used for a particular purpose. For a satellite analysis, users must set the wavelength at the band center, the bandwidth, and the spectral response function of the band they want to analyze. It may be difficult for beginners to appropriately specify these settings to facilitate comparison of results using different settings, and it is time-consuming to perform calculations many times with a single setting (i.e., the same set of values for many parameters). Most radiative transfer codes for satellite analyses already include information about sensors; users only indicate the index numbers of the band and sensor as input data. For example, 6S, a radiative transfer code for clear sky and solar wavelengths, contains spectral response functions at a 0.025 µm resolution and calculates radiative quantities using an approximation method involving successive orders of scattering at each wavelength. RSTAR uses the discrete-ordinate method to calculate radiative transfer and accurately treats particle scattering, but it takes a relatively long time to achieve the same resolution as 6S.
The Himawari-8 satellite carries the Advanced Himawari Imager (AHI), which has greatly improved spatial and temporal resolutions compared to Japan's previous meteorological satellites, the Geostationary Meteorological Satellite (GMS) and Multi-Functional Transport Satellite (MTSAT) series (Bessho et al. 2016). Figure 1 shows the spectral distribution of outgoing radiance at the top of the atmosphere (TOA) and at the Earth's surface multiplied by the response function within the AHI band 16 wavelength range. We assumed a satellite zenith angle of 0° and that atmospheric conditions were typical of an Air Force Geophysics Laboratory (AFGL) standard atmosphere at mid-latitudes during the summer (Anderson et al. 1986); however we assumed a CO 2 concentration of 360 ppm. This spectral response function is a transmittance, which is provided by the Meteorological Satellite Center of the Japan Meteorological Agency (JMA). The spectrally integrated value in Fig. 1 is equal to satellite-observed radiance in the above described atmospheric conditions. Integrating such a spiky spectral distribution requires many quadrature abscissas.
The correlated k-distribution (CKD) method (Lacis and Oinas 1991;Fu and Liou 1992) is a rapid method for evaluating atmospheric gas absorption, and has been incorporated into many broadband models.  In this study, we used this method for each of the sensors' wavelength bands to increase the efficiency of the gas absorption process. Moreover, by adopting an optimization method to determine quadrature abscissas and weights, we reduced the number of required radiative transfer calculations. Furthermore, we introduced an optimization method into RSTAR and updated its package to make it suitable for satellite retrieval analysis. Section 2 describes the models and datasets we used. In Section 3, we explain the optimized CKD method and its evaluations. Section 4 shows the results the AHI retrieved using this method. In Section 5 we summarize this paper.

Models and datasets
Rstar7 is a narrow-band model containing two standard gas absorption tables that covers the wavenumber spectrum from 10 to 54,000 cm −1 , divided log-linearly into either 3732 or 7464 bands. The base-10 logarithm bandwidths are 0.001 and 0.0005, thus, the bandwidths tend to be wider at shorter wavelengths. These bandwidths may be incompatible with the AHI's resolution. When the CKD method is applied, the number of quadrature points is fixed at two per band. The abscissas and weights for integration are determined by a squared Gaussian quadrature with abscissas that are doubled abscissas of a Gaussian quadrature and weights that are products of abscissas and weights. We assumed perfectly correlated overlapping in bands involving multiple gases. We tabulated the absorption coefficients of seven major gases (H 2 O, CO 2 , O 3 , N 2 O, CO, CH 4 , and O 2 ) in each band for 26 log-linear pressure levels and three temperature levels. We used the atmospheric gas absorption database HITRAN 2004 (Rothman et al. 2005) for line absorptions, and MT_CKD_1 code (Mlawer et al. 2012) for continuous absorption. As discussed above, users set quadrature points and weights for spectral integration by RSTAR.
To make calculations rapid and accurate, we created a new gas absorption table corresponding to each AHI band. We applied the CKD method, and reduced the number of quadrature points using an optimization method to determine the abscissas and weights for integration. This method combining the CKD method and optimization was originally developed for Model Simulation radiation TRaNsfer code (MSTRNX; Sekiguchi and Nakajima 2008, hereafter called SN08), which is a broadband model adapted to the Model for Interdisciplinary Research on Climate (MIROC; Watanabe et al. 2010) and the Nonhydrostatic Icosahedral Atmospheric Model (NICAM; Satoh et al. 2008). MSTRNX is known to be a fast and stable radiative transfer code. We used HITRAN 2012 (Rothman et al. 2013), the latest version of HITRAN data and MT_ CKD_2.1 code, to obtain absorption coefficients. The MT_CKD_2.1 code we used to calculate continuous absorption spectra was developed by AER Inc., and has been updated to version 3, which corresponds to HITRAN 2012. We plan to continue updating with the latest version.

Methods
When using the CKD method, the optical characteristics of particle scattering, solar irradiance, and the Planck function for blackbody radiation should be independent of wavelength within each band. However, optical properties vary significantly with wavelength in two of AHI's bands. In AHI band 3, absorption by oxygen is present in part of the spectrum but weak in other parts. AHI band 8 is very wide compared to the other bands, and the Planck function within it varies greatly with wavelength. We divided these two bands, AHI band 3 and 8, into two sub-bands each. We initially specified the bandwidths and gases in each AHI band as indicated in Table 1. In this study, we defined the bandwidths as the full width at half maximum (FWHM) of the spectral response function (SRF) and applied the CKD method to each bandwidth. To determine the important absorption gases in each band, we compared calculated radiative fluxes with and without a target gas using reference calculations with resolutions the same as that of the SRF distributed by JMA. The resolutions were 0.1 cm −1 (bands 7 -16) or 1.0 cm −1 (bands 1 -6), respectively. We took a target gas into consideration if the error caused by not considering it was greater than 1 %. We sorted the spectral absorption coefficients within each bandwidth to generate a k-distribution, and stored target gas species k-distributions at 26 pressure levels from 0.01 to 1013.25 hPa and three temperatures (200, 260, and 320 K) for each band.
The k-distribution method often uses a trapezoidal or Gaussian quadrature i to carry out numerical integrations. However, in this study, we determined the quadrature abscissas and weights using sequential quadratic programming, which is an interactive nonlinear optimization method. This method is almost identical to SN08, but with two differences. One is the setting of the objective function and the other is the way of selecting appropriate quadrature abscissas and weights.
In this optimization process, quadrature abscissas and weights are set to minimize an objective function, which is defined as the square root of the sum of the squared differences between the calculations using this method and the reference results of the radiative fluxes at the TOA and surface, and profiles of heating rates. This setting of the objective function is the same as in SN08, but SN08 instead used a Line-By-Line method as a reference calculation. The reference results are the calculated radiative fluxes and heating rate, with high-spectral resolutions, multiplied by the spectral response function and integrated over a bandwidth. The spectral resolutions of the reference calculation were set to the same value as the SRF. We used six AFGL standard atmospheric conditions (with CO 2 concentrations modified to 360 ppm) in these calculations. We assumed clear sky conditions because scattering due to clouds and aerosols impacted the reference results and made optimization difficult. To account for optical path length variations due to multiple scattering and differences in sun and satellite position, we assumed optical path lengths with satellite zenith angles equivalent to either 0° or 60°. This is also different from SN08.
When the objective function decreased, the optimization process sometimes identified a local minimum rather than a global one. To avoid this problem, we used two initial conditions to start the optimization process if the band included more than two gas species. The first condition was a completely correlated overlapping and the second was a completely uncorrelated overlapping. We calculated the initial abscissas and weights using a Gaussian quadrature. Subsequently, we used two processes to change the number of quadrature points, N c . One process was to optimize for each N c separately and the other process was used to sequentially decrease N c . In the latter case, the optimization process started from N c = 8 (the initial condition corresponding with completely correlated overlapping for all gas species except 1 or 3, in which case it corresponded with completely uncorrelated overlapping) or N c = 9 (the initial condition corresponding with completely uncorrelated overlapping for gas species 2). When the optimization process with N c quadrature points was completed, the initial condition corresponding to N c -1 quadrature points was defined. We then removed the quadrature point that contributed the least from the optimized quadrature with N c quadrature points. Using this method, we were able to obtain optimized results in each case of N c quadrature points, and we selected the set of abscissas and weights that gave the best results. Finally, we determined the CKD parameters for N c values from 1 to 6. In SN08, N c was not selectable and was already set to perform properly with all bands in a general circulation model. On the other hand, in this study, users could select the N c best suited for their purposes. Figure 2 shows the differences of radiative fluxes calculated using the optimized CKD method for N c = 1, 2, 4, and 6 from the reference calculation. Panels corresponding to the visible (VIS) and near infrared (NIR) bands (1 -7) show net flux differences, and panels corresponding to the thermal infrared (TIR) bands (8 -16) show upward flux differences. We assumed that atmospheric conditions corresponded to those of the mid-latitude summer model of the AFGL standard atmosphere. We also assumed the solar zenith angle to be 60°, and the surface to be Lambertian with an albedo 0.1 in the VIS and NIR bands and surface emissivity of 1.0 in the TIR bands. In bands 3 and 8, which were divided into two sub-parts, we used the same N c in each sub-part and summed the results. The satellite received radiances were observed at TOA; however, the radiance profile is important for cloudy sky cases. In general, the larger the value of N c , the better the results. In the case of bands treated as a single gas species, the difference from the reference was smaller than in the case of multiple gases (see Table 1) because the optimization process converged easily. The altitude of maximum difference would indicate a maximum dependency of a main target gas in each band. We also checked the angular dependencies of satellite-received radiance. The differences between results from the combined CKD method and those from the reference method did not change with satellite zenith angle in the TIR bands, but in the VIS and NIR bands, the differences became large with larger satellite zenith angles, and their variance was about same as the method difference. We plan to increase the satellite zenith angles of the objective functions in the next update.  Fig. 2. Differences in radiative fluxes calculated from the reference calculations using the optimized correlated k-distribution (CKD) method. Solid, broken, dotted, and bold lines indicate results with a number of quadrature points (N c ) = 1, 2, 4, and 6, respectively. Panels corresponding to bands 1 -7 show net flux differences and panels corresponding to bands 8 -16 show upward flux differences. Atmospheric conditions were assumed to be equal to those of the mid-latitude summer, Air Force Geophysics Laboratory (AFGL) standard atmosphere. The solar zenith angle was assumed to be 60°. In bands 3 and 8, which were divided into two sub-parts, the results for each N c case were added.

Calculation check
We used the CKD table corresponding to N c = 6 in a forward model of the ICAS (Iwabuchi et al. 2018) to simulate brightness temperatures in eight TIR bands (9 -16) for clear sky pixels over the ocean observed by the AHI. In the calculations, we obtained land and ocean surface temperatures and emissivity from the Moderate Resolution Imaging Spectroradiometer (MODIS) 8-day mean land and ocean products, and spatially and temporally interpolated atmospheric temperature profiles and water vapor and ozone concentrations from the Modern-Era Retrospective analysis for Research and Applications (MERRA) meteorological field product (Rienecker et al. 2011). Concentrations of CO 2 , N 2 O, and CH 4 were assumed to be equal to monthly mean values of the global mean. The forward model is based on a two-stream radiative transfer solution in a plane-parallel multilayered atmosphere. Details of error evaluation for the AHI longwave bands are presented in Iwabuchi et al. (2018). Figure 3 shows scatter plots between the ob- served and model-simulated brightness temperatures for each AHI band. We used full disk data for 19 -28 August 2015. We discriminated clear sky pixels based on confi dently clear sky pixels identifi ed by the collocated MODIS cloud mask product. Table 2 shows the means and standard deviations of differences between model simulations. We attribute these differences to errors in the assumed temperatures and water vapor profi les, sea surface temperatures (SST), and the gas absorption table calculated by the optimized CKD method. The standard deviation of the SST error was estimated to be about 0.4 -0.5 K and the radiometric calibration accuracy ranged from 0.20 -0.29 %, which we converted to 0.11 -0.18 K when assuming brightness temperature was 300 K. Given this estimated standard deviation, the trends of the model results and observations were in good agreement, except for band 11, in which the model results were overestimates. One of the reasons for this may be that SO 2 was not considered in band 11 in this version, whereas the other three gas species (H 2 O, N 2 O, and CH 4 ) were taken into consideration (Table 1). The transmittance, which is considered SO 2 in this band, is estimated at about 90.2 % with the AFGL concentration profi le of SO 2 . On the other hand, transmittance when it is not considered is estimated to be about 90.7 %. SO 2 is considered a main gas species in this band, and as it is the target of this band in the AHI design, it should be included in the gas absorption model and will be introduced in next update. In the water vapor bands (9 and 10), the standard deviation of the error was larger than 1 K. These fi gures are used in clear sky pixels over the  Table 2. ocean, and their error sources are limited to errors in atmospheric gas profi les and SST measurements. Because these bands are sensitive to atmospheric profi les in the middle and upper troposphere, SST error is not estimated as the main reason in these band. Compared to temperatures, water vapor amounts in the middle and upper troposphere in the atmospheric reanalysis product are more uncertain. A main source of uncertainty in these bands is error in the assumed amounts of water vapor in the middle and upper troposphere. For the VIS and NIR bands, we considered a similar analysis that uses the TOA refl ectance (not shown), but for that analysis, we would have needed more information such as aerosol optical properties and sea surface emissivity. Therefore, that analysis is left for a future study. Infrared measurements for cloudy pixels are sensitive mainly to cloud top temperature, cloud optical thickness, and particle effective radius, and secondly are sensitive to cloud geometrical thickness and vertical inhomogeneity along with surface properties and atmospheric profi les. Not all cloud parameters are reliably available from observation data; however, to evaluate the CKD models, it would be interesting to Fig. 4. Histograms of brightness temperature (BT) differences between measurements and model calculations for the cloud retrieval results for full disk data at 12 UTC on August 20, 2015. Cloud retrieval was performed with correlated k-distribution (CKD) sets with different numbers of quadrature points (N c ). The means and standard deviations of the BT differences are shown in the legend in each panel.
test model-measurement consistency and the effects of the CKD model on cloud retrieval. Cloud properties were retrieved using the ICAS, in which different values of N c (2, 3, and 6) were used for forward brightness temperature calculations. We obtained a set of retrieved cloud properties and simulated brightness temperatures for optimal solutions. The optimal estimation framework used in ICAS attempts to estimate cloud properties that best fi t the measurements. Thus, the model calculations should fi t well with measurements if the forward model has smaller errors, and vice versa. Figure 4 shows histograms of brightness temperature differences between model calculations and measurements for the cloud retrieval results. The mean difference between the measurements and the model with N c = 2 was different from zero, whereas the mean bias with N c = 6 was almost zero except in band 12. In band 12, the BT difference means were larger with larger N c values; however, the histogram was more symmetric with larger N c values. The mean biases were generally smaller than in the clear sky cases shown in Fig. 3 because model calculations were fi tted to the measurements in the optimal estimation-based cloud property inversion. In the water vapor bands (9 and 10), the standard deviations of the differences were larger than in other bands. This was because of the same reason as in the clear sky cases. The standard deviations of the differences for N c = 2 were signifi cantly larger than those for N c = 3 and 6, indicating that the model using N c = 2 did not fi t the measurements well. Figure 5 shows global distributions of cloud top height (CTH) and cloud optical thickness (COT) retrieved using the optimized CKD table with different numbers of quadrature points (N c = 2, 3, and 6). These results are shown only for pixels with solutions optimized by the ICAS. The number of pixels was signifi cantly smaller when N c = 2 than when N c = 3 and 6, and the spatial distributions of cloud properties retrieved with N c = 3 and 6 seemed reasonable. Figure 6 quantitatively shows differences in COT and CTH as a function of N c . Although a difference between N c = 3 and 6 is not clearly apparent in Fig. 5, a difference in COT for low clouds is apparent in Fig.  6b. By comparing only pixels with optimal solutions, we better estimated high-cloud COT than low-cloud COT, irrespective of N c . Because the ICAS uses eight thermal infrared AHI bands, the COT estimate has greater certainty for high clouds than low clouds. This is because of the larger difference between cloud and underlying-surface temperatures. The difference between N c = 2 and 6 was widely distributed (Figs. 6a, c), which is primarily due to misinterpretation of low clouds as high clouds in the retrieval using N c = 2. This misinterpretation is found over the Indian Ocean near the western coast of Australia, as shown in Fig. 5. Because brightness temperatures are generally comparable for optically thick low clouds and optically thin high clouds, an inter-band consistency (i.e., brightness temperature differences among bands) is required to discriminate the two cloud types. It is suggested that the inter-band consistency in the N c = 2 case is not very good. Calculations using the modifi ed RSTAR package were rapid. Figure 7 shows an example of cloud retrieval fi elds over the Sea of Japan at 01 UTC and 04 UTC on 7 April 2017 that used the satellite analysis version of RSTAR we developed in this study. Cirrus clouds were present in the upper left portion of the images at 01 UTC, and low clouds were in the center of the same images. Cirrus and low clouds overlapped at 04 UTC. The top, second, third, and bottom panels show true-color composite images, COT retrieval results, cloud effective radii (CER), and CTHs, respectively. The panels on the left were retrieved using VIS, NIR, and TIR bands (4, 6, and 13, respectively). The panels on the right were retrieved using TIR bands (11, 13, and 15). Compared with two results of CERs and CTHs at overlapping pixels, it is apparent that the VIS-NIR-TIR method retrieves mainly low clouds, whereas the TIR method retrieves mainly cirrus clouds. The various bands of the AHI provide much information about many targets, thus the modifi ed RSTAR package would work for high-resolution satellite analysis.

Summary
We developed gas absorption tables using an optimized CKD method to produce rapid and accurate satellite measurement simulations. The number of quadrature points, N c , directly affected computational efficiency. We made CKD parameter tables in which N c varied from 1 to 6 and could be selected by the user. We have checked cloud retrieval results with different values of N c . In cases where N c = 2, cloud retrieval results were significantly different from cases where N c > 2. For high accuracy, we recommend N c > 3.
In this study, we used radiative flux and heating rate for the objective function as was the case in SN08; however, for satellite analysis it is suitable to use radiance in various solar and satellite angles and surface conditions. In addition, we considered it difficult to adopt cloud and aerosol for the objective function. This was because it is required to divide error sources from gas absorption and particle scattering, and the objective function is not easy to converge if the number of parameters increases. In future work, we VIS-NIR-TIR method (bands #4, #6, and #13) 01UTC 04UTC TIR method (bands #11, #13, and #15) 01UTC 04UTC Fig. 7. Sample of low clouds and overlapping cirrus clouds over the Sea of Japan at 01 and 04 UTC on 7 April 2017. The top, second, third, and bottom panels are true-color composite images and retrieved results of cloud optical thickness (COT), cloud effective radii (CER), and cloud top height (CTH), respectively. The panels on the left were retrieved using visible (VIS), near infrared (NIR), and thermal infrared (TIR) bands (4, 6, and 13, respectively). The panels on the right were retrieved using TIR bands (11, 13 and 15) with a version of RSTAR that was updated for this study.
plan to study the effects of several parameters on the objective function, update the latest version of continuum program MT_CKD_3, and introduce SO 2 absorptions in bands 10 and 11. The spectral responsivity of AHI-9, mounted on Himawari-9, is slightly different from that of AHI-8, thus an extension to AHI-9 is also planned. We incorporated the CKD tables and band information from the AHI into a multi-purpose atmospheric radiative transfer package, RSTAR. We updated the package for satellite analysis so users could easily specify the satellite and band number. We also developed tables for Aqua/MODIS, CALIPSO/IIR, and Landsat-7/ETM+. This package makes high-speed and high-precision cloud and aerosol retrievals that are suitable for high-frequency and high-resolution observations made by satellites such as Himawari-8 possible.