Development of a Rapid Retrieval Method for Cloud Optical Thickness and Cloud-Top Height Using Himawari-8 Infrared Measurements

We have developed a rapid simplified algorithm to retrieve cloud optical thickness and cloud-top height from measurements of the infrared split-window bands of Himawari-8. The method is based on a rapid calculation model for clear-sky brightness temperatures and empirical equations for cloud, for which the coeffi- cients are determined by a fit to a more rigorous radiative-transfer model. This method can be applied regardless of regions excluding the polar regions and season by taking into account the temperature, humidity, sea surface temperature, and surface emissivity. In this study, we have demonstrated that this method captures well the diurnal cycle of cloud amounts of different cloud types in the warm-pool region around Indonesia. With an accelerated retrieval process by a factor of around 1,000 compared with the the phys- ics-based retrieval, our rapid cloud retrieval algorithm yielded cloud amounts that agree quantitatively with those from a more rigorous, physics-based cloud retrieval method.


Introduction
The radiative effect of clouds varies greatly depending on cloud properties such as cloud-top height (CTH) and cloud optical thickness (COT), and future changes in different cloud types may result in different warming/cooling impacts on the Earth's climate. For a better understanding on the cloud-climate feedback, it is important to understand the physical and dynamical mechanisms of variations of the cloud-precipitation systems on various scales. Geostationary meteorological satellites are useful for examining temporal variations including diurnal and life cycles of clouds. Himawari-8 began operating in July 2015 and is equipped with the Advanced Himawari Imager (AHI) visible-infrared radiometer (Bessho et al. 2016). The temporal resolution is 10 minutes for the full disk observation. The AHI has six and ten bands in the solar and infrared spectra, respectively, and the horizontal resolution for the infrared bands is roughly 2 km at nadir.
In previous studies of cloud variations, cold clouds were often identified by brightness temperature (equivalent blackbody temperature) thresholds (Chen and Houze 1997;Machado et al. 1998;Tobin et al. 2012;Inoue et al. 2008). The brightness temperature is close to the cloud-top temperature (CTT) if the cloud is optically thick. However, in cases of optically thin cloud, radiation from the atmosphere below the cloud and the background surface is easily transmitted, making the brightness temperature considerably higher than the actual CTT (Protopapadaki et al. 2017). For this reason, it is difficult to distinguish optically thin high clouds and lower clouds by using only the brightness temperature of a single band. Inoue (1985Inoue ( , 1987 pioneered the use of the 11-µm band and the brightness temperature difference (BTD) between the bands at 11 and 12 µm (split-window bands) for ice-cloud remote sensing. Those two bands have different absorptivities by water, making them suitable for detecting different types of cloud, particularly opaque and transparent ice clouds. Hamada and Nishi (2010) have derived a look-up table for the CTH of tropical clouds as a function of brightness temperature in the 11-µm band and BTD using the reference CTH determined by the CloudSat cloud-profiling radar.
Modern methods for cloud remote sensing can retrieve more cloud properties (e.g., COT, effective particle radius, and cloud thermodynamic phase) based on a rigorous radiative-transfer model and the use of three or more spectral bands and instantaneous atmospheric and surface data. The Integrated Cloud Analysis System (ICAS) performs cloud analysis using eight AHI bands (Iwabuchi et al. 2016(Iwabuchi et al. , 2018. Cloud retrieval in the ICAS uses the optimal estimation method (Rodgers 2000), which fits the physics model to the measurements using optimization. Although this physics-based cloud retrieval (PCR) is highly accurate and reliable, the computational cost tends to be extremely high when analyzing a large amount of data for a wide area. To benefit from geostationary meteorological satellites, rapid cloud retrieval is desirable.
In this study, a rapid simplified cloud retrieval (SCR) algorithm for estimating COT and CTH is developed using the splitwindow bands. Section 2 describes (i) a rapid radiative-transfer calculation for clear-sky conditions and (ii) the SCR algorithm. Section 3 presents an example result of the diurnal cycle of clouds as retrieved by SCR and its comparison with the result from PCR. Finally, conclusions are given in Section 4.

Rapid clear-sky calculation
The present algorithm is based on a rapid radiative-transfer calculation for clear sky. The vertical profile of the radiative source function is represented as a third-order polynomial. For an absorbing (non-scattering) atmosphere, the radiative-transfer equation can be written as (1) where µ is the cosine of the emission angle, τ is the optical depth from the top of the atmosphere, I is the band-mean spectral radiance, and a n (n = 0, …, 3) are the coefficients of the third-order polynomial. A least-squares fit to the Planck-function profile determines the polynomial. Solving (1) analytically gives the downward radiation at the bottom of the atmosphere. The sum of the thermal emission from the surface and the reflection of downward radiation at the surface is then calculated. The upward radiance at the top of the atmosphere is obtained analytically from (1) and transmission of the upward radiance at the bottom of the atmosphere. The solution of (1) is presented in Supplementary Materials. Finally, the radiance is converted to the brightness

Development of a Rapid Retrieval Method for Cloud Optical Thickness and Cloud-Top Height
Using Himawari-8 Infrared Measurements T T  T  y T T   T  T   obs  obs  eq   obs  obs  eq   clr  e where T clr and T obs denote clear-sky BT13 and observed BT13, respectively, ΔT obs is the observed BTD, T eq is BT13 assuming optically thick cloud with COT of 10, and θ is the satellite zenith angle. The empirically determined function y is defined as where ΔT clr is BTD for clear sky and a, b, c, d, and e are constants. The role of y is to make the relationship between u and v close each other for different conditions as shown in Fig. 1b. Here, y originates from the point (u, v) = (0, 0), and other terms were empirically designed after testing sensitivities on the view zenith angle and water vapor amount. The relationship between u and v is mostly unique (hereinafter referred to as the u-v curve) as shown in Fig. 1b, irrespective of CTP, µ, and water vapor. A single approximation is made with a fourth-order polynomial as v Au Au where A i (i = 1, …, 4) are determined by independent fits for each cloud phase (i.e., ice or liquid). T clr and ΔT clr are calculated using the method in Section 2.1. u is derived from (2) after solving T eq from simultaneous equations of u and v as formulated by (2)−(5), and T eq is calculated by substituting u into (2). Note that T eq is higher than the CTT (T top ) because of transmission from the lower (warmer) part within the cloud. To estimate T top , the following empirical equation is determined by comparing SCR-derived T eq and PCR-derived T top using granules in different seasons and for daytime/nighttime data: The root-mean-square-error of T top for the fit is 14.5 K.
Since radiance generally decreases exponentially with increasing COT (τ c ), a relationship can be written as (8) using u¢ defined as temperature.
The data used are basically the same as those used in the ICAS PCR. Land and sea surface temperatures are obtained from the Moderate Resolution Imaging Spectroradiometer eight-day mean land and ocean products (Wan et al. 2015). Land surface emissivity is obtained from the baseline-fit emissivity database monthly average product (Seemann et al. 2008). Sea surface emissivity is calculated using the Fresnel equations for a flat water surface. Water vapor, ozone, and carbon dioxide are considered as the absorbing gases. The atmospheric temperature, specific humidity, and ozone mixing ratio are obtained from the assimilated meteorological fields of the Modern-Era Retrospective Analysis for Research and Application, Version 2 (MERRA-2) data (Gelaro et al. 2017). Carbon dioxide concentrations are obtained from the global monthly molar ratio data provided by the World Data Centre for Greenhouse Gases of the World Meteorological Organization Global Atmosphere Watch program. The absorption coefficients of each gas and band were predetermined by a fit to the ICAS forward model (Iwabuchi et al. 2016). The accuracy of brightness temperature in the present method is summarized in Supplementary Materials. The root-mean square errors in brightness temperature of bands 13 (10.4 µm) and 15 (12.4 µm) are 0.65 and 0.91 K, respectively.

Characteristics of split-window bands
In the SCR algorithm, AHI infrared bands 13 (10.4 µm) and 15 (12.4 µm) are mainly used. Figure 1a shows the relationship between the brightness temperature of AHI band 13 (BT13) and the BTD of bands 13 and 15 for various COT and cloud-top pressure (CTP). This relationship was obtained with the ICAS forward model, which calculates the infrared brightness temperatures with a radiative-transfer model assuming the cloud properties, the atmospheric temperature and humidity, the surface temperature, and the satellite view zenith angle. The effective particle radius varies depending on CTT and COT, and the details are described in Supplementary Materials 3. For a constant CTP, the relationship between BT13 and BTD is an arch-like curve (hereinafter referred to as the BT13-BTD curve). The maximum BT13 is shown for clear sky (i.e., zero COT), and BT13 decreases with increasing COT because of decreasing radiation from the lower part of atmosphere and surface. The arches of the BT13-BTD curves become higher with larger BTD values for lower CTP (colder cloud top) because of increasing difference between surface temperature and cloud-top temperature.

Simplified cloud retrieval
In SCR, brightness temperatures are represented as functions of COT, CTT, surface temperature, satellite zenith angle, and amount of water vapor. For the formulation, the two equations for normalized BT13 and BTD, namely, u and v, respectively, are introduced first: Fig. 1. (a) Brightness temperature of Advanced Himawari Imager (AHI) band 13 (BT13) and brightness temperature difference (BTD) between bands 13 and 15 assuming various values of cloud optical thickness (COT). Circles correspond to COT = 0, 0.1, 0.3, 1, 3, 10. The colors blue, green, and red correspond to cloud-top pressures (CTPs) of 200, 320, and 440 hPa, respectively. The BT13 and BTD were calculated by ICAS forward model given the COT and CTPs. A mid-latitude atmospheric model in summer is assumed with an ice cloud over the surface with a temperature of 290 K. The effective particle radius varies depending on the cloud-top temperature (CTT) and COT. The satelite zenith angle is assumed to be 0°. (b) Relationship between u and v defined in (2) and (3), calculated from the results corresponding to (a).
where B is the Planck function at AHI band 13. F (u¢) is defined as 2nd-order polynomial: where f, g, and h are constants. The constants in the above equations are determined for each cloud phase by fitting to the ICAS forward model and are summarized in Tables 1 and 2. The constants a, b, c, d, e, and A i (i = 1, …, 4) were determined to minimize the error of the u-v curve represented by (5) under various conditions. In this formulation, we considered various atmospheric profiles in different regions (tropics and mid-latitude) and seasons (summer and winter), COTs, CTPs, surface temperatures, satellite zenith angles, and amounts of water vapor. The cloud-particle effective radius was parameterized by using statistical values obtained from ICAS PCR and is represented as a function of τ c and T top , as shown in Supplementary Materials (Fig. S1). The root-mean-square error of v derived from (5) is about 0.003 for ice clouds and 0.008 for liquid clouds when u = 0.5.
Because cloud properties are retrieved for both ice and liquid clouds, the optimal cloud phase is selected depending on T top and the uncertainty of the modeled brightness temperature of band 13. This way is based on the viewpoint that cloud phase can be reasonably determined if the model better explains the observed brightness temperature. Look-up tables for the root-mean-square error of the brightness temperature uncertainty are shown in Fig. 2 and are calculated from T obs retrieved by substituting the retrieved u into (2) and T obs of the ICAS forward model. The modeling uncertainties for the two cloud phases are compared to determine the optimal cloud phase if T top is between −38 and 0°C. Multilayer and mixed-phased clouds are not supported in this work.
Finally, T top is converted to CTP and CTH by interpolation using the MERRA-2 atmospheric profile. If there is a temperatureinversion layer, the lowest level that corresponds to CTT is chosen as the cloud top.

Comparison with physics-based cloud retrieval
The calculation speed of SCR is accelerated by a factor of around 1,000 compared to PCR. As an example of applying SCR, we analyzed the cloud diurnal cycle in the warm-pool region around Indonesia. The region corresponds to 90°E−160°E and 15°S to 15°N with a spatial resolution of around 2 km. The analysis period was one month of January 2016, and hourly data were used. To investigate the diurnal cycle of different cloud types, we defined the cloud types in terms of COT and CTP obtained by SCR. Low-level clouds have CTP above 680 hPa, mid-level clouds have CTP between 440 and 680 hPa, and high-level clouds have CTP below 440 hPa. High-level clouds are further classified into thick clouds with COT above 6, cirrostratus (Cs) with COT between 1 and 6, and cirrus (Ci) with COT between 0.2 and 1. Our SCR is not applicable to very optically thin clouds with COT below 0.2. Figure 3 shows the time series of pseudo-color composite images constructed from AHI infrared bands 11, 12 and 15 and cloud types categorized based on the aforementioned criteria around New Guinea Island. In the pseudo-color images, high and thick clouds are shown in white and cirrus clouds in purple. In RGB composite images, high and thick clouds are shown white and cirrus clouds looks like blue shade around white part in (a). This figure shows the typical diurnal cycle of convection, particularly over land. Thick clouds developed from 1700 LT over land accompanied by small amounts of Cs and Ci, grew in the evening, and disappeared in the early morning.
The same analysis is performed by PCR. A comparison between SCR and PCR is shown in Fig. 4. The results are shown separately for land and open ocean 100 km from the coast. In Fig.  4, the SCR and PCR results show similar diurnal variations. The intensity and phase of the diurnal variations of the cloud fraction for each cloud type closely match those from PCR. Quantitatively, SCR underestimates low clouds and cirrus compared with PCR. SCR did not be used for data with T clr − T obs less than 5 K because it is difficult to estimate cloud property in SCR when CTH is very low or optically-thin, high cloud is present. On the other hands, SCR overestimates thick clouds compared with PCR as shown in Figs. 4c and 4d. This could be ascribed that PCR tends to underestimate COT greater than 6 (Iwabuchi et al. 2018). According to Figs. 4a and 4c, high-level clouds exhibit maximum coverage around 1800 LT over land, with thick cloud contributing considerably. This agrees with the fact that convection becomes active in the afternoon because of the sea breeze induced by solar heating on the ground over land. Furthermore, there is a large extent of Cs clouds in later hours. This delay suggests that thick cloud spreads widely and become thinner over time. The amplitude of diurnal variation of the cloud fraction is weaker over ocean (Figs. 4b and 4d) than over land (Figs.4a and 4c). In Fig. 4d, thick clouds increase slightly between 2100 LT and 0600 LT over the ocean. This agrees with the finding of a previous study that deep convective clouds become active by longwave cooling from the cloud top during the night over the ocean (Randall and Dazlich 1991). Table  3 shows the relative error of COT and the error of CTH. These errors are calculated from data in which the cloud phase in SCR  matches with that in PCR in this analysis. The error is reasonable.
The main objective of this study is development of a rapid calculation method which well match with the PCR by using the same data (infrared measurement, MERRA-2 atmospheric and MODIS surface products). Iwabuchi et al. (2018) showed comparisons between the PCR and radar and lidar cloud product.

Summary
We have developed a rapid SCR calculation algorithm to estimate COT and CTH. Our SCR uses the split-window bands (AHI bands 13 and 15) of Himawari-8. The SCR is performed by constructing an empirical model of brightness temperatures as functions of COT, CTT, satellite zenith angle, surface temperature, and amount of water vapor. The SCR can be applied regardless of regions excluding the polar regions and season since we have considered temperature, humidity, sea surface temperature, and surface emissivity for various conditions in our formulation. The method is based on a rapid calculation of brightness temperatures for clear sky.
Here we have analyzed the cloud diurnal cycle in the warmpool region around Indonesia. SCR and PCR gave almost the same results for the diurnal variation of the cloud fraction of each cloud type. Because SCR allows a rapid calculation of cloud retrieval, in future work we will perform large-scale analysis of cloud variations, taking advantage of the high spatial and temporal resolutions of Himawari-8.