Satellite retrieval of aerosol combined with assimilated forecast

We developed a new aerosol satellite retrieval algorithm combining a numerical aerosol forecast. In the retrieval algorithm, the short-term forecast from an aerosol data assimilation system was used as an a priori estimate instead of spatially and temporally constant values. This method was demonstrated using observation of the Advanced Himawari Imager onboard the Japan Meteorological Agency’s geostationary satellite Himawari-8. Overall, the retrieval results incorporated strengths of the observation and the model and complemented their respective weaknesses, showing spatially finer distributions than the model forecast and less noisy distributions than the original algorithm. We validated the new algorithm using ground observation data and found that the aerosol parameters detectable by satellite sensors were retrieved more accurately than an a priori model forecast by adding satellite information. Further, the satellite retrieval accuracy was improved by introducing the model forecast instead of the constant a priori estimates. By using the assimilated forecast for an a priori estimate, information from previous observations can be propagated to future retrievals, leading to better retrieval accuracy. Observational information from the satellite and aerosol transport by the model are incorporated cyclically to effectively estimate the optimum field of aerosol.

Abstract. We developed a new aerosol satellite retrieval algorithm combining a numerical aerosol forecast. In the retrieval algorithm, the short-term forecast from an aerosol data assimilation system was used as an a priori estimate instead of spatially and temporally constant values. This method was demonstrated using observation of the Advanced Himawari Imager onboard the Japan Meteorological Agency's geostationary satellite Himawari-8. Overall, the retrieval results incorporated strengths of the observation and the model and complemented their respective weaknesses, showing spatially finer distributions than the model forecast and less noisy distributions than the original algorithm. We validated the new algorithm using ground observation data and found that the aerosol parameters detectable by satellite sensors were retrieved more accurately than an a priori model forecast by adding satellite information. Further, the satellite retrieval accuracy was improved by introducing the model forecast instead of the constant a priori estimates. By using the assimilated forecast for an a priori estimate, information from previous observations can be propagated to future retrievals, leading to better retrieval accuracy. Observational information from the satellite and aerosol transport by the model are incorporated cyclically to effectively estimate the optimum field of aerosol.

Introduction
Aerosols have a fundamental influence on the energy budget of the earth's climate system through the scattering and absorption of solar radiation. The fifth assessment report of the Intergovernmental Panel on Climate Change (IPCC, 2014) stated that radiative forcing of the total aerosol effect in the atmosphere, including cloud adjustments due to aerosols, is −0.9 W m −2 . The report also highlighted that the range of uncertainties in these radiative forcing estimations remains large (−1.9 W m −2 to −0.1 W m −2 ). Identifying the frequency and properties of aerosols over the globe by satellite measurements is essential in estimating the radiation budget and the impacts of aerosols on climate systems.
In satellite aerosol remote sensing, not all aerosol properties can be accurately detected by satellite sensors, as there are more unknown aerosol parameters (e.g., particle size distributions, vertical density distribution, shape, refractive index) than the actual information obtained by the sensors. Most studies use assumptions or information about aerosol parameters and limit the number of parameters retrieved. For example, Higurashi and Nakajima (1999) and Fukuda et al. (2013) assumed fixed complex refractive indices (1.5-0.005i in Higurashi and Nakajima (1999), and 1.503-7.16 × 10 −8 i for small-mode particles and 1.445-1.00 × 10 −8 i for coarse-mode particles based on sulfate and sea spray models, respectively, in Fukuda et al., 2013) and retrieved the aerosol optical thickness and Ångström exponent. Some studies assumed aerosol particle models according to the location and season. For example, Kaufman et al. (1997), Remer et al. (2005), and Levy et al. (2007) estimated aerosol optical thickness and fine-mode fraction over a dark target using the Moderate Resolution Imaging Spectroradiometer (MODIS), by selecting the fine-dominated particle models as a function of geography and season. Hsu et al. (2004) retrieved the optical thickness and type of aerosols over desert regions using blue channels (< 500 nm), where the surface reflectance was relatively low, assuming dust or a mixture of dust and smoke depending on the region and season. Jeong et al. (2016) used a priori information according to the location and season. However, these studies did not take temporal changes into account. Because it is impossible to completely fix aerosol type as a function of geographical location and season, the unrealistic assumptions hence implemented lead to one of the major causes of retrieval error.
Aerosol data assimilation methods using satellite data have also been developed to obtain better initial conditions for the aerosol transport model. The aerosol data assimilation study was first developed with low earth orbit (LEO) satellites (Benedetti et al., 2009;Saide et al., 2013;Dai et al., 2014;Rubin et al., 2016;Yumimoto et al., 2015). In recent years, assimilation studies have been extended to using geostationary satellites with large spatial coverage and fine observation frequencies (Saide et al., 2014;Lee et al., 2016;Yumimoto et al., 2016Yumimoto et al., , 2018Dai et al., 2019;Jin et al., 2019).
Due to the development of such assimilation studies, the satellite data have contributed to improving aerosol forecast simulations. However, no studies have utilized assimilated model forecast as an a priori estimate of the retrieval. Since satellite sensors cannot accurately detect all parameters and unrealistic assumptions of aerosol parameters are a major cause of retrieval errors as mentioned above, adding the model information is expected to improve the retrieval accuracy. Therefore, in this study, we utilize the forecast of an aerosol transport model for a priori estimates of the retrieval. This allows the aerosol information in the aerosol transport model to be used for retrieval. By using the assimilated forecast, information from previous satellite observations can be propagated to future satellite retrievals through the aerosol transport model.
The sections in this study are organized as follows: Sect. 2 explains the retrieval methodology in detail. Section 3.1 presents the results of application to the Advanced Himawari Imager (AHI) onboard Himawari-8. Section 3.2 describes the validation of the estimations using ground observations, and Sect. 3.3 tests the worst-case scenario. Finally, Sect. 4 summarizes our findings.

Methodology
The aerosol retrieval algorithm in this study is based on . As an a priori estimate of the retrieval, the algorithm introduces aerosol forecast from a transport model that has assimilated previous satellite observations. Given the general applicability of the retrieval algorithm by , the methodology explained in this section can also be applied to various sensors. Here, we demonstrate the algorithm using the Himawari-8/AHI whose assimilation system is operationally available. The AHI has six observation bands from visible to near-infrared wavelength ranges and observes the top-of-atmosphere (TOA) radiance at a resolution of 0.5-2.0 km over Asia and Oceania at 10 min intervals (Bessho et al., 2016). Figure 1 depicts an overview of the algorithm, showing the process of using forecast data for a priori estimates of the retrieval. In the original retrieval process, the Level-2 (L2) aerosol optical thickness at 500 nm (τ ), the Ångström exponent between 400 and 600 nm (α), and the single-scattering albedo at 500 nm (ω) are retrieved using Level-1 (L1) AHI-observed radiance every 10 min around time T 0 as per . The Level-3 (L3) τ and α at T 0 are then estimated using L2 products in 1 h by an hourly-combined algorithm . The hourly-combined algorithm is a method that (1) minimizes cloud contamination using the difference between aerosol and cloud spatiotemporal variability characteristics and (2) interpolates the aerosol retrievals using 1 h of data and the movement of clouds within the hour (see , for more details).
The L3 τ at T 0 is then assimilated into a global aerosol transport model by the 2D-Var assimilation system . For the aerosol transport model, we use MASINGAR (Model of Aerosol Species IN the Global At-mospheRe; Tanaka et al., 2003;Tanaka and Chiba, 2005) developed at the Meteorological Research Institute (MRI) of the Japan Meteorological Agency (JMA). MASINGAR covers the major tropospheric aerosol components (i.e., black and organic carbon, mineral dust (10 size bins), sea salt (10 size bins), sulfate aerosols) and their precursors (e.g., sulfur dioxide (SO 2 ), dimethyl sulfide, terpenes), and is coupled online with an atmospheric general circulation model (MRI-AGCM3; Yukimoto et al., 2012). The model's grid resolution is set to horizontal Gaussian TL479 (960 × 480 grids, about 0.375 • ) and 40 vertical layers in hybrid sigmapressure coordinates from the ground to 0.4 hPa. The integration time step is set to 600 s. Anthropogenic emissions of SO 2 and black and organic carbon are taken from the MAC-City emission inventory (Granier et al., 2011). Daily biomass burning emission flux is taken from the Global Fire Assimilation System (GFAS, Kaiser et al., 2012) version 1.2 provided by the European Centre of Medium Range Forecast (ECMWF). The horizontal wind components and sea surface temperature are nudged toward the global analyses and forecasts of JMA (GANAL). The forecast from the assimilation system serves as the operational sand and dust forecasting by JMA, the aerosol property model product in the JAXA Himawari Monitor (https://www.eorc.jaxa.jp/ptree/index.html, last access: 20 January 2021) and a member of the ICAP multi-model ensemble (MME) (Xian et al., 2019). The vol- ume concentration (then τ ) of each aerosol component at the next time (T 1 ) is then forecasted using the assimilated aerosol transport model.
In the new retrieval process, we retrieve the L2 aerosol properties (τ , α, and ω) from AHI-observed radiance at T 1 using these L4 forecasts for a priori estimates of the retrieval. In this way, the information from previous observations at T 0 is used for the next aerosol retrievals at T 1 through the aerosol transport model. The retrieval obtained at T 1 is further used in the same way to derive the retrieval at the following time step (T = T 2 , not shown) by using L4 forecasts for an a priori estimate. Figure 6 compares the improved retrieval with the original retrieval at T 1 as later described in Sect. 3.2. The methodology for using the forecast as a priori estimates of the retrieval is detailed as follows: in the retrieval process, the final retrieval parameters (τ , α, and ω) are calculated from the set of aerosol parameters (τ , external mixing ratio of dry-volume concentration of fine particles η f , and external mixing ratio of the dry-volume concentration of dust particles for the coarse model η dst c ) defined by Yoshida et al. (2018). Here, the imaginary part of the refractive index (m i ) for the fine-aerosol model varies with change in η dst c such that the fine and coarse models exhibit the same ω at 500 nm (see Yoshida et al., 2018, for more details). The α and ω are calculated from the retrieved η f and η dst c (i.e., m i for fine-aerosol model) using the tables previously calculated by radiative transfer code called the System for the Transfer of Atmospheric Radiation, whose development was initially led by the University of Tokyo (STAR; Tanaka, 1986, 1988;Stamnes et al., 1988). The detailed aerosol setting is explained in  and is outlined in Appendix A. Appendix B shows the relationship of the final retrieval parameters (α and ω) with the set of aerosol parameters (η f and η dst c ). We retrieve the aerosol parameters (τ, η f , and η dst c ) using an optimal estimation method (Rodgers, 2000). The state vector of a set of aerosol parameters x = τ, η f , η dst c is derived by minimizing object function J (Eq. 1). It uses the measurement vector of a gascorrected observed reflectance set R = {ρ obs where x a = τ a , η f a , η dst c a is the vector of a prior estimate of x and S e and S a are the covariance matrices of R and x a , respectively. The calculations of R, F (x), and S e are the same as those of Yoshida et al. (2018), but we apply canonical correlation analysis to find the optimal coordinate system and converted R, F (x), and S e to the coordinate system whose dimension is reduced to the number of retrieved parameters (i.e., three). In the original retrieval process, we used spatially and temporally constant values of x a and S a that were derived from climatology analysis and assumed that the non-diagonal  To introduce a more realistic a priori estimate and covariances into the retrieval process, we employ the forecast from the aerosol assimilation system instead of the constants. The model forecast includes the total aerosol optical thickness at 500 and 870 nm and the absorption aerosol optical thickness at 500 nm derived from the modeled volume concentration and extinction cross section of each aerosol component (Yumimoto et al., 2017). We assign an a priori estimate x a as follows: the model's total aerosol optical thickness at 500 nm is used for τ a . η f a is derived from the ratio of total aerosol optical thickness between 500 and 870 nm. As the selection of η dst c a , we use the model's ω as calculated from the total and absorption aerosol optical thickness at 500 nm.
The assimilation system uses an ensemble method to calculate the background error covariance matrix . In the method, the ensemble was collected from forecast values within ±2 h of the targeted hour of the five previous forecasts (Fig. 2). We employ this method to define S a . The model ensemble enables S a to include the nondiagonal component and express the error of aerosol transport. However, S a from model ensemble may become too small when the model does not predict the aerosol event it-self. For that reason, in order to estimate total S a , we add a model absolute error (S A a ) to the error estimated from the ensemble (S E a ) as follows: where are the standard deviations of τ a , η f a , and η dst c a , respectively, estimated from the ensemble. σ A τ a , σ A η fa , and σ A η dst c a are the same as those of the model absolute error.
σ A τ a is set to σ G τ a or σ M τ a (whichever is larger) as follows: where σ G τ a is the root mean square error (RMSE) of the model's τ from ground observation (0.399 in Fig. 6c) and σ M τ a is the standard deviation of τ for 5 years as calculated by the free-run model without assimilation. The free-run model's spatial resolution is around 1.2 • , and the standard deviation is calculated for MAM (March, April, and May), JJA (June, July, and August), SON (September, October, and November), and DJF (December, January, and February) (Fig. 3).
is set to 0.5 because η dst c takes a value from 0. to 1, and ω (which is uniquely determined by η dst c ) has little correlation with the ground observation data (Fig. 6i). As the non-diagonal component of S A a cannot currently be calculated from our limited database, we use the non-diagonal components of S E a as those of S a .

Results and discussion
3.1 Results of application to Himawari-8 We applied the methodology described in Sect. 2 to the Himawari-8/AHI. We retrieved τ , η f , and η dst c and then derived ω and α at 10 min intervals from the calibrated L1 data subsampled at 0.05 • using the method described in Sect. 2. The channels used for the retrieval are the same as those used by , which are channels 1 (0.46 µm), 2 (0.51 µm), 3 (0.64 µm), 4 (0.86 µm), and 5 (1.6 µm) over land and channels 4 and 5 over the ocean. As the number of satellite channels (two) used over the ocean is less than the number of retrieval parameters (three), not all parameters are stably retrieved by satellite data. Therefore, η dst c over the ocean, which is the least sensitive to satellite observation, is set to 0 (i.e., non-absorbing aerosol) at this time because the aerosol over the ocean is generally less absorbing than that over land and using the model's η dst c a as η dst c over the ocean leads to a worse estimation of τ (not shown). However, using nonabsorbing aerosol over the ocean causes a big problem in the case of dust or smoke transported over the ocean, so we will use the model's η dst c a as η dst c over the ocean after obtaining a better model's η dst c a in the future. Note that η dst c over land is properly retrieved from satellite data (i.e., not set to 0) using the model's η dst c a as an a priori estimate, since the number of satellite channels (five) used over land is greater than the number of retrieval parameters (three).  (Figs. 4, 5a and b). The retrieval results from the new algorithms using only S E a for S a are also shown to evaluate the effect of S A a on the retrieval result (Figs. 4, 5c). Figure 4 depicts the retrieved τ , η f , and η dst c at 02:00 UTC on 19 May 2016 when aerosols originating from wildfires near Lake Baikal in Russia reached Japan. The model's x a used for retrieval in Fig. 4b and c is indicated in Fig. 4d. The σ τ a , σ η fa , and σ η dst c a used for retrieval in Fig. 4b are shown in Fig. 4e. The white regions indicate the area where retrieval is not executed due to the presence of clouds, etc. The 2 h forecasts starting from 00:00 UTC on 19 May were assimilated with L3 merged τ at 03:00, 06:00, and 09:00 UTC on 18 May and at 00:00 UTC on 19 May and then used for an a priori estimate (Fig. 4d). Figure 5 is the same as Fig. 4 except for another case at 05:00 UTC on 7 May 2017, when Asian dust was observed in Japan. The 5 h forecast starting from 00:00 UTC on 7 May (and assimilated at 03:00, 06:00, and 09:00 UTC on 6 May and at 00:00 UTC on 7 May) is used for an a priori estimate (Fig. 5d). These short-term forecasts with data assimilation are considered relatively realistic compared to long-term forecasts or a free run without assimilation ). If only model ensemble error is used for S a (Figs. 4c, 5c), that is, the absolute error is not included in S a , all retrieved parameters (especially η f and η dst c over land) are highly dependent on an a priori estimate (Figs. 4d, 5d). However, when using an appropriate S a containing absolute error, the retrieved τ , η f , and η dst c are updated by satellite data or remain close to an a priori estimate depending on the location (Figs. 4b, 5b). Specifically, spatially finer τ distributions than the model forecast are retrieved for cases of both wildfire aerosol (Fig. 4b) and Asian dust (Fig. 5b) due to the relatively coarser model horizontal resolution. Similar η f is retrieved over open ocean in Figs. 4b and Fig 5b, and the large η f (i.e., small particle) and small η f (i.e., large particle) are successfully retrieved in areas corresponding to wildfire aerosol (Fig. 4b) and Asian dust (Fig. 5b), respectively. This distribution is also expressed in the forecast model in Fig. 5d but cannot be expressed sufficiently in the forecast model in Fig. 4d because information about the aerosol particle size (e.g., α, η f ) is not assimilated into the model. That is, by using an appropriate S a , both the model and satellite data are used for estimating the aerosol parameters. In addition, the local noise in τ and η f is apparently reduced for this algorithm (Figs. 4b, 5b) as compared with the original algorithm (Figs. 4a, 5a). This will be discussed in Sect. 3.2.

Validation
We conducted a preliminary validation of our method by comparing the retrieved τ , α, and ω from the Himawari-8/AHI with those from ground observation known as the Aerosol Robotic Network (AERONET). AERONET's τ and α were derived from Level 2.0 quality-assured Version 3 direct sun algorithm data (Giles et al., 2019;O'Neill et al., 2003), and ω was derived from Version 3 AERONET inversion products Dubovik et al., 2002aDubovik et al., , b, 2006Sinyuk et al., 2007). AERONET's ω at 500 nm was calculated from linear interpolation of ω at 440 and 675 nm. In this study, the 60 AERONET sites on the full disk of Himawari-8 were used for the validation. We used the AERONET data averaged over 10 min of AHI observation time. For our retrieval data, we used τ , α, and ω estimated from AHI L1 radiance data subsampled at 0.05 • nearest to the AERONET sites. Initial validation was conducted for 6 months (March, April, May, June, and July 2018 and February 2019). Long-term validation will be required in future studies. Figure 6 compares the τ , α, and ω retrieved from the AHI with those from AERONET. The validation of α and ω is limited to cases where the simultaneously retrieved τ are greater than 0.3 because there is little information of α and ω from satellite observation for thin aerosol layers. The total number of validation points (14 711, 14 031, and 521) from this algorithm (Fig. 6b, e, and h) is about 6 %-7 % higher than those (13 714, 13 137, and 493) from the original algorithm (Fig. 6a, d, and g), which means that the new algorithm successfully retrieved the aerosol in more cases than the original algorithm. Here, the total number of validation points for ω is less than those for τ and α because the number of ω data  Fig. 6a, b, and c. from AERONET inversion products is less than those of τ and α from the direct sun measurements.
For the τ estimations ( Fig. 6a, b, c), the root mean square error (RMSE), mean bias (MB), and correlation (0.290, −0.099, and 0.758) from this algorithm (Fig. 6b) are all better than those (0.399, −0.224, and 0.572) from the model forecast (i.e., an a priori estimate) in Fig. 6c, which means that satellite information is very effective for the retrieval of τ . In addition, the RMSE (0.290) in Fig. 6b is better than that (0.307) without the forecast model (Fig. 6a), which means that the model information is also effective and the improved algorithm shows better performance than the original algorithm. The MB (−0.099) in Fig. 6b is worse than that (−0.023) in Fig. 6a, probably because the large outlier in Fig. 6a is improved in Fig. 6b. Figure 7 shows an example of the retrieval results of the outlier (red asterisks in Fig. 6a, b, and c) and shows that the outlier of the original algorithm (Fig. 7a) is improved in the new algorithm by constraining τ to the model's τ a . In addition, the retrieval results around the red circles show that the new algorithm successfully retrieved the τ , η f , and η dst c even where the original algorithm failed to retrieve them. Thus, integrating the model and satellite information resulted in an improvement of the τ estimations.
For the α estimations (Fig. 6d, e, f), large variance in the original method is considerably reduced by this method. The RMSE and correlation (0.271 and 0.581) from this algorithm (Fig. 6e) are much better than those (0.429 and 0.353) from the original algorithm without the forecast model (Fig. 6d), which indicates that the new algorithm could improve the precision of α estimations by adding more accurate α (RMSE of 0.223) information from the model. In addition, the MB (−0.052) from this algorithm (Fig. 6e) is better than that (−0.057) from the model forecast (Fig. 6f), due to the improvement of negative bias in the large α in the model forecast. Thus, the new α can be retrieved with good accuracy by Figure 8. Frequency distribution of τ retrieved from AHI and those from AERONET. The results retrieved from this algorithm in the case of χ 2 less than 20, 0.5, and 0.2 and uncertainties of the retrieved τ (S τ ) less than 20, 1.0, and 0.5 are plotted in each panel. E, B, and R above the figures show the root mean square error, mean bias, and correlation, respectively.
utilizing the relatively accurate model's α and correcting the bias by using the satellite data.
For the ω estimations (Fig. 6g, h, i), the RMSE, MB, and correlation (0.035, 0.000, and 0.550) from this algorithm (Fig. 6h) are better than those (0.048, −0.002, and 0.176) from the model forecast (Fig. 6i), which indicates the effectiveness of satellite information for ω retrieval. In addition, this algorithm improved RMSE, MB, and correlation by introducing the model forecast. Note that the current system assimilates only total τ (total amount of aerosols), and information about the fraction of fine and absorbing aerosols from the retrieval is not fed back to the model forecast. The improved retrieval accuracy of ω can be expected if the model's ω becomes more realistic in the future, such as by assimilating the satellite's ω to the model. Considering the vali-dation results of τ , α, and ω, this new algorithm effectively improved the retrieval accuracy using information from both the model and the satellite by setting appropriate S a and S e .
We also investigated the cause of the possible large deviation between the retrieved parameters from the new algorithm and the ground observation. Figures 8, 9, and 10 show the validation results of τ , α, and ω, respectively, when the chi-square value (χ 2 ) and the uncertainties of the retrieved three parameters (τ , η f , and η dst c ) are smaller than a threshold. The chi-square value (χ 2 ) is calculated as follows: It shows the closeness of the retrieved value to the observed value. The covariance matrix of the uncertainties of the retrieved parameters Sx is calculated using the law of error Figure 9. Frequency distribution of α retrieved from AHI and those from AERONET. The results retrieved from this algorithm in the case of χ 2 less than 20, 0.5, and 0.2 and uncertainties of the retrieved η f (S η f ) less than 20, 0.5, and 0.2 are plotted in each panel. E, B, and R above the figures are the same as in Fig. 8. propagation, as follows: where A is the Jacobian matrix. S e is the covariance matrix of R and calculated from sum of sensor noise and the uncertainty in TOA reflectance that results from surface reflectance uncertainty . In reality, the S e is almost determined by the uncertainty in TOA reflectance that results from surface reflectance uncertainty because sensor noise is much smaller. Therefore, the Sx is mostly caused by the surface reflectance uncertainty. Figure 8 shows that RMSE for τ decreases as the threshold of χ 2 or S τ becomes strict (i.e., decreases). On the other hand, RMSE for α (in Fig. 9) is not dependent on the threshold of S η f but decreases as the χ 2 threshold decreases. RMSE for ω (in Fig. 10) is almost independent on the threshold of S η dst c and χ 2 . Next, in Fig. 11 we investigated how the retrieved accuracy (difference between aerosol parameters retrieved from AHI and those of AERONET) depends on the model's (i.e., a priori) accuracy. The retrieved accuracy of α and ω has strong linear relationships (a correlation of 0.801 and 0.739, respectively) to the model's accuracy, while that of τ has a moderate linear relationship (a correlation of 0.622). Summarizing these results, the retrieved accuracy of τ depends on all of the closeness to the observed value, accuracy of the surface reflectance estimation, and accuracy of an a priori estimate, while the accuracy of an a priori estimate is critical for the retrieved accuracy of α and ω. Thus, introducing more realistic a priori estimates into the new retrieval algorithm instead of the constant values in the original algorithm led to the improvement of RMSE. It is also shown that the improvement of a numerical aerosol forecast by improving the aerosol transport model Figure 10. Frequency distribution of ω retrieved from AHI and those from AERONET. The results retrieved from this algorithm in the case of χ 2 less than 20, 0.5, and 0.2 and uncertainties of the retrieved η dst c (S η dst c ) less than 20, 0.5, and 0.2 are plotted in each panel. E, B, and R above the figures are the same as in Fig. 8. and the assimilation method and increasing the assimilation frequency may further improve the retrieval accuracy in the future.

Worst-case scenario
We have shown that the new retrieval algorithm using the forecast of an aerosol transport model improves the retrieval accuracy. However, in order to use this algorithm constantly (such as in an operational system), the effects of the model forecast (a priori estimate) that deviate from reality must be examined because the model forecast may miss an aerosol event. Therefore, we conducted a sensitivity test to investigate the impact on the retrieval results of using unrealistic forecast as an a priori estimate. Figure 12 shows the retrieval results on the same day as in Fig. 4, except for using the forecast on another day (27 April 2018) as an a priori estimate of the retrieval (Fig. 12d). If only S E a is used as S a (Fig. 12c), all parameters (especially η f and η dst c ) are retrieved unrealistically by being dependent on the unrealistic a priori estimate. However, when using an appropriate S a (Eq. 2), the retrieved parameters are updated well by satellite data with less dependence on an unrealistic a priori estimate (Fig. 12b). Even in Figure 11. Frequency distribution of the difference between τ (a), α (b), and ω (c) retrieved from AHI and those from AERONET, as a function of the difference between τ (a), α (b), and ω (c) of an a priori estimate and AERONET. R shows the correlation. such an extremely worst-case scenario, this new algorithm is apparently not significantly worse than the current algorithm, especially where the model forecast is missing an aerosol event, which may occur in the model forecast for natural aerosols (e.g., mineral dust and smoke from biomass burning).

Summary
We developed a new satellite aerosol retrieval algorithm combining a numerical aerosol forecast. In the retrieval algorithm, the short-term forecast from an aerosol data assimi-lation system was used for an a priori estimate instead of spatially and temporally constant values. This is the first study that utilizes the assimilated model forecast of aerosol as an a priori estimate of the satellite retrieval. We applied this new algorithm to the Himawari-8/AHI and confirmed that the aerosol parameters detectable by satellite sensors were retrieved more accurately (RMSE of 0.290 for τ and 0.035 for ω) than an a priori model forecast (RMSE of 0.399 for τ and 0.048 for ω) by adding satellite information. Moreover, the satellite retrieval accuracy was improved (RMSE of 0.290 for τ , 0.271 for α, and 0.035 for ω) by using the model forecast as compared with those using constant a priori estimates M. Yoshida et al.: Satellite retrieval of aerosol combined with assimilated forecast 1809 (RMSE of 0.307 for τ and 0.429 for α and 0.039 for ω). As a result, aerosol retrievals were improved by effectively incorporating both model and satellite information, depending on each covariance. By using the assimilated forecast as an a priori estimate, information from previous observations can be propagated to future retrievals, thereby leading to better retrieval accuracy. In this way, satellite observation and model simulation are used synergistically to continuously estimate the optimum field of aerosol. Future work would include applying the methodology proposed in this study to polar-orbiting satellites and combining them with geostationary satellite measurements, in order to offer consistent geostationary and polar-orbiting estimates and thereby improve aerosol properties over the globe.

Appendix A: Aerosol setting
We assume that the aerosol model is an external mixture of fine and coarse particles (η f is the external mixing ratio of the dry-volume concentration of fine particles). We set the fine-aerosol model based on the average properties of the fine mode for categories 1-6 by Omar et al. (2005). For the coarse-aerosol model, we set the external mixture of the pure marine aerosol on the basis of the model illustrated by Sayer et al. (2012) and the dust model based on the coarse model of category 1 (dust) as illustrated by Omar et al. (2005). η dst c is the external mixing ratio of the dry-volume concentration of dust particles for the coarse model.
Regarding each aerosol size, we use a monomodal lognormal volume size (r d ) distribution, which is defined as follows: where C v is the particle volume concentration, r v is the volume median radius, and σ is the standard deviation. r v is set to 0.143, 2.59, and 2.834 (σ is 1.537, 2.054, and 1.908) for fine, coarse marine, and coarse dust, respectively, based on the observations by Omar et al. (2005) and Sayer et al. (2012). Regarding the aerosol shape, we assume a spherical model for the fine and coarse marine models and a nonspherical model for the coarse dust model (Nakajima et al., 1989). The aerosol vertical distribution is set to the same distribution that was used for rural (dominant at 0-2 km), sea spray (below 2 km), and yellow sand (4-8 km) for fine, coarse marine, and coarse dust in the STAR code, respectively. The real part of the refractive index is set to 1.439, 1.362, and 1.452 for fine, coarse marine, and coarse dust, respectively, and the imaginary part of the refractive index (m i ) is set to 3.0 × 10 −9 and 0.0036 at all wavelengths for coarse marine, and coarse dust, respectively, based on Sayer et al. (2012) and Omar et al. (2005). The m i for the fine-aerosol model is perturbed to represent non-absorbing and absorbing aerosols. To decrease the number of derived parameters, the m i for the fine-aerosol model varies with change in η dst c such that the fine and coarse models exhibit the same ω at 500 nm.
Appendix B: Relationship of α and ω with η f and η dst c Figure B1 shows the relations of the final retrieval parameters α and ω with the external mixing ratio of dry-volume concentration of fine particles (η f ) and external mixing ratio of the dry-volume concentration of dust particles for the coarse model (η dst c ). The ω at 500 nm can be uniquely determined by the η dst c (Fig. B1a), since η dst c for the coarse aerosol changes in conjunction with m i for the fine aerosol so that the ω at 500 nm has the same value without depending on the η f . Note that the ω at wavelengths other than 500 nm are dependent on not only η dst c but also η f . The α is mainly determined by η f but also depends slightly on η dst c (Fig. B1b). Figure B1. The relations of (a) single-scattering albedo at 500 nm (ω) and (b) Ångström exponent between 400 and 600 nm (α) with the external mixing ratio of dry-volume concentration of fine particles (η f ) and the external mixing ratio of the dry-volume concentration of dust particles for the coarse model (η dst c ). Each color represents a different η dst c .
Author contributions. MY developed the retrieval code and analyzed the data with significant conceptual input from KY and critical technical support from TMN, MK, and HM. KY and TYT prepared the assimilated forecast data for test cases and long time validations, respectively. MY prepared the paper with contributions from all co-authors.
Competing interests. The authors declare that they have no conflict of interest.
Acknowledgements. The authors are grateful to the Open CLASTER project for allowing us to use the RSTAR package for this research. We would like to thank the AERONET project and its staff for establishing and maintaining the AERONET sites considered in this investigation. We also thank Haruma Ishida for providing the cloud detection algorithm (CLAUDIA). Finally, we appreciate the valuable discussions and support provided by Takashi Maki, Tsuyoshi Sekiyama, Makiko Hashimoto, and Teruyuki Nakajima.
Financial support. This research has been supported by the JSPS Grants-in-Aid for Scientific Research (grant no. JP16H02946) and the JAXA Earth Observation Priority Research (grant no. PI.ER2GCF212).
Review statement. This paper was edited by Stelios Kazadzis and reviewed by Alexander Kokhanovsky and two anonymous referees.