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.

1 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⁻². The report also highlighted that the range of uncertainties in these radiative forcing estimations remains large (−1.9 W m⁻² to −0.1 W m⁻²). 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–$\mathrm{7.16}\times {\mathrm{10}}^{-\mathrm{8}}i$ for small-mode particles and 1.445–$\mathrm{1.00}\times {\mathrm{10}}^{-\mathrm{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., 2016, 2018; Dai 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.

2 Methodology

The aerosol retrieval algorithm in this study is based on Yoshida et al. (2018). 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 Yoshida et al. (2018), 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 Flowchart of data processing for aerosol retrieval at time T₁.

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₀ as per Yoshida et al. (2018). The Level-3 (L3) τ and α at T₀ are then estimated using L2 products in 1 h by an hourly-combined algorithm (Kikuchi et al., 2018). 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 Kikuchi et al., 2018, for more details).

The L3 τ at T₀ is then assimilated into a global aerosol transport model by the 2D-Var assimilation system (Yumimoto et al., 2018). For the aerosol transport model, we use MASINGAR (Model of Aerosol Species IN the Global AtmospheRe; 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₂), 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 sigma–pressure coordinates from the ground to 0.4 hPa. The integration time step is set to 600 s. Anthropogenic emissions of SO₂ and black and organic carbon are taken from the MACCity 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 volume concentration (then τ) of each aerosol component at the next time (T₁) is then forecasted using the assimilated aerosol transport model.

Figure 2 Forecast of aerosol transport model used for retrieval at time T₁. Solid and dashed lines show the assimilation period (1 d) and forecast run, respectively.

Figure 3 Mean (upper) and standard deviation (lower) of τ for free-run model from 2011 to 2015 in (a) March, April, and May, (b) June, July, and August, (c) September, October, and November, and (d) December, January, and February.

In the new retrieval process, we retrieve the L2 aerosol properties (τ, α, and ω) from AHI-observed radiance at T₁ using these L4 forecasts for a priori estimates of the retrieval. In this way, the information from previous observations at T₀ is used for the next aerosol retrievals at T₁ through the aerosol transport model. The retrieval obtained at T₁ is further used in the same way to derive the retrieval at the following time step (T=T₂, not shown) by using L4 forecasts for an a priori estimate. Figure 6 compares the improved retrieval with the original retrieval at T₁ 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 ${\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}$) 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 ${\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}$ 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 ${\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}$ (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; Nakajima and Tanaka, 1986, 1988; Stamnes et al., 1988). The detailed aerosol setting is explained in Yoshida et al. (2018) 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 ${\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}$). We retrieve the aerosol parameters (τ,η_f, and ${\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}$) using an optimal estimation method (Rodgers, 2000). The state vector of a set of aerosol parameters $\mathit{x}=\left\{\mathit{\tau },{\mathit{\eta }}_f,{\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}\right\}$ is derived by minimizing object function J (Eq. 1). It uses the measurement vector of a gas-corrected observed reflectance set $\mathit{R}=\mathit{\{}{\mathit{\rho }}_{\mathrm{i}}^{{\mathrm{obs}}^{\prime }}$, $i=\mathrm{1},\mathrm{\dots },N\mathit{\}}$ and simulated TOA reflectance $F\left(\mathit{x}\right)=\mathit{\{}{\mathit{\rho }}_{\mathrm{i}}^{\mathrm{sim}},i=\mathrm{1},\mathrm{\dots },N\mathit{\}}$, where N is the channel number.

$$\begin{array}{}\text{(1)}& J={\left[\mathit{R}-F\left(\mathit{x}\right)\right]}^T{\mathbf{S}}_{\mathrm{e}}^{-\mathrm{1}}\left[\mathit{R}-F\left(\mathit{x}\right)\right]+{\left[\mathit{x}-{\mathit{x}}_{\mathrm{a}}\right]}^T{\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}\left[\mathit{x}-{\mathit{x}}_{\mathrm{a}}\right],\end{array}$$

where ${\mathit{x}}_{\mathrm{a}}=\left\{{\mathit{\tau }}_{\mathrm{a}},{\mathit{\eta }}_{{\mathrm{f}}_{\mathrm{a}}},{{\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}}_a\right\}$ 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 component of covariance matrices was set to 0 (Yoshida et al., 2018).

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. ${\mathit{\eta }}_{{\mathrm{f}}_{\mathrm{a}}}$ is derived from the ratio of total aerosol optical thickness between 500 and 870 nm. As the selection of ${{\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}}_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 (Yumimoto et al., 2018). 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 non-diagonal 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 itself. For that reason, in order to estimate total S_a, we add a model absolute error (${\mathbf{S}}_{\mathrm{a}}^{\mathbf{A}}$) to the error estimated from the ensemble (${\mathbf{S}}_{\mathrm{a}}^{\mathrm{E}}$) as follows:

$$\begin{array}{}\text{(2)}& {\displaystyle }{\mathbf{S}}_{\mathrm{a}}={\mathbf{S}}_{\mathrm{a}}^{\mathrm{E}}+{\mathbf{S}}_{\mathrm{a}}^{\mathrm{A}}=\left[\begin{array}{ccc}{\mathit{\sigma }}_{{\mathit{\tau }}_{\mathrm{a}}}^{\mathrm{2}}& {\mathit{\sigma }}_{{\mathit{\tau }}_{\mathrm{a}}{\mathit{\eta }}_{{\mathrm{f}}_{\mathrm{a}}}}& {\mathit{\sigma }}_{{\mathit{\tau }}_{\mathrm{a}}{{\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}}_a}\\ {\mathit{\sigma }}_{{\mathit{\tau }}_{\mathrm{a}}{\mathit{\eta }}_{{\mathrm{f}}_{\mathrm{a}}}}& {\mathit{\sigma }}_{{\mathit{\eta }}_{{\mathrm{f}}_{\mathrm{a}}}}^{\mathrm{2}}& {\mathit{\sigma }}_{{\mathit{\eta }}_{{\mathrm{f}}_{\mathrm{a}}}{{\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}}_a}\\ {\mathit{\sigma }}_{{\mathit{\tau }}_{\mathrm{a}}{{\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}}_a}& {\mathit{\sigma }}_{{\mathit{\eta }}_{{\mathrm{f}}_{\mathrm{a}}}{{\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}}_a}& {\mathit{\sigma }}_{{{\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}}_a}^{\mathrm{2}}\end{array}\right],\text{(3)}& {\displaystyle }{\mathit{\sigma }}_{{\mathit{\tau }}_{\mathrm{a}}}={\mathit{\sigma }}_{{\mathit{\tau }}_{\mathrm{a}}}^{\mathrm{E}}+{\mathit{\sigma }}_{{\mathit{\tau }}_{\mathrm{a}}}^{\mathrm{A}},\text{(4)}& {\displaystyle }{\mathit{\sigma }}_{{\mathit{\eta }}_{{\mathrm{f}}_{\mathrm{a}}}}={\mathit{\sigma }}_{{\mathit{\eta }}_{{\mathrm{f}}_{\mathrm{a}}}}^{\mathrm{E}}+{\mathit{\sigma }}_{{\mathit{\eta }}_{{\mathrm{f}}_{\mathrm{a}}}}^{\mathrm{A}},\text{(5)}& {\displaystyle }{\mathit{\sigma }}_{{{\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}}_a}={\mathit{\sigma }}_{{{\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}}_a}^{\mathrm{E}}+{\mathit{\sigma }}_{{{\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}}_a}^{\mathrm{A}},\end{array}$$

where ${\mathit{\sigma }}_{{\mathit{\tau }}_{\mathrm{a}}}^{\mathrm{E}}$, ${\mathit{\sigma }}_{{\mathit{\eta }}_{{\mathrm{f}}_{\mathrm{a}}}}^{\mathrm{E}}$, and ${\mathit{\sigma }}_{{{\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}}_a}^{\mathrm{E}}$ are the standard deviations of τ_a, ${\mathit{\eta }}_{{\mathrm{f}}_{\mathrm{a}}}$, and ${{\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}}_a$, respectively, estimated from the ensemble. ${\mathit{\sigma }}_{{\mathit{\tau }}_{\mathrm{a}}}^{\mathrm{A}}$, ${\mathit{\sigma }}_{{\mathit{\eta }}_{{\mathrm{f}}_{\mathrm{a}}}}^{\mathrm{A}}$, and ${\mathit{\sigma }}_{{{\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}}_a}^{\mathrm{A}}$ are the same as those of the model absolute error.

Figure 4 Aerosol optical thickness at 500 nm τ (upper), external mixing ratio of dry-volume concentration of fine particles η_f (middle), and external mixing ratio of the dry-volume concentration of dust particles for the coarse model ${\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}$ (lower) that are (a) retrieved from the original algorithm (i.e., using constant an a priori estimate), (b) retrieved from this algorithm, (c) retrieved from this algorithm but without model absolute error (${\mathbf{S}}_{\mathrm{a}}^{\mathrm{A}}$), and (d) of the model forecast at 02:00 UTC on 19 May 2016. (e) Standard deviations of model forecast (${\mathit{\sigma }}_{{\mathit{\tau }}_{\mathrm{a}}}$, ${\mathit{\sigma }}_{{\mathit{\eta }}_{{\mathrm{f}}_{\mathrm{a}}}}$, and ${\mathit{\sigma }}_{{{\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}}_a}$) used for retrieval in (b).

Figure 5 Same as Fig. 4, except for the case at 05:00 UTC on 7 May 2017.

${\mathit{\sigma }}_{{\mathit{\tau }}_{\mathrm{a}}}^{\mathrm{A}}$ is set to ${\mathit{\sigma }}_{{\mathit{\tau }}_{\mathrm{a}}}^{\mathrm{G}}$ or ${\mathit{\sigma }}_{{\mathit{\tau }}_{\mathrm{a}}}^{\mathrm{M}}$ (whichever is larger) as follows:

$$\begin{array}{}\text{(6)}& {\mathit{\sigma }}_{{\mathit{\tau }}_{\mathrm{a}}}^{\mathrm{A}}=\left\{\begin{array}{ll}{\mathit{\sigma }}_{{\mathit{\tau }}_{\mathrm{a}}}^{\mathrm{G}}& \text{ if }{\mathit{\sigma }}_{{\mathit{\tau }}_{\mathrm{a}}}^G\ge {\mathit{\sigma }}_{{\mathit{\tau }}_{\mathrm{a}}}^{\mathrm{M}}\\ {\mathit{\sigma }}_{{\mathit{\tau }}_{\mathrm{a}}}^{\mathrm{M}}& \text{ else}\end{array}\right\},\end{array}$$

where ${\mathit{\sigma }}_{{\mathit{\tau }}_{\mathrm{a}}}^{\mathrm{G}}$ is the root mean square error (RMSE) of the model's τ from ground observation (0.399 in Fig. 6c) and ${\mathit{\sigma }}_{{\mathit{\tau }}_{\mathrm{a}}}^{\mathrm{M}}$ 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). ${\mathit{\sigma }}_{{\mathit{\eta }}_{{\mathrm{f}}_{\mathrm{a}}}}^{\mathrm{A}}$ (0.093) is calculated from RMSE of the model's α (0.223 in Fig. 6f) at α of 1.2 and ${\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}$ of 0.5. ${\mathit{\sigma }}_{{{\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}}_a}^{\mathrm{A}}$ is set to 0.5 because ${\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}$ takes a value from 0. to 1, and ω (which is uniquely determined by ${\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}$) has little correlation with the ground observation data (Fig. 6i). As the non-diagonal component of ${\mathbf{S}}_{\mathrm{a}}^{\mathrm{A}}$ cannot currently be calculated from our limited database, we use the non-diagonal components of ${\mathbf{S}}_{\mathrm{a}}^{\mathrm{E}}$ as those of S_a.

3 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 ${\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}$ 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 Yoshida et al. (2018), 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, ${\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}$ 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 ${{\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}}_a$ as ${\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}$ over the ocean leads to a worse estimation of τ (not shown). However, using non-absorbing 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 ${{\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}}_a$ as ${\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}$ over the ocean after obtaining a better model's ${{\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}}_a$ in the future. Note that ${\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}$ over land is properly retrieved from satellite data (i.e., not set to 0) using the model's ${{\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}}_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).

Figure 6 Frequency distribution of τ (a, b, c), α (d, e, f), and ω (g, h, i) retrieved from AHI and those from AERONET. Panels (a), (d), and (g) show the results from the original algorithm (i.e., using constant a priori), (b), (e), and (h) show the results from this algorithm, and (c), (f), and (i) are an a priori estimate used for (b), (e), and (h), respectively. E, B, R, and N above the figures show the root mean square error, mean bias, correlation, and total number, respectively. Red asterisks in (a), (b), and (c) indicate the results at the red circles in Fig. 7a and b.

Figures 4 and 5 compare the retrieval results from the new algorithm using ${\mathbf{S}}_{\mathrm{a}}^{\mathrm{E}}+{\mathbf{S}}_{\mathrm{a}}^{\mathrm{A}}$ for S_a with the original algorithm (Figs. 4, 5a and b). The retrieval results from the new algorithms using only ${\mathbf{S}}_{\mathrm{a}}^{\mathrm{E}}$ for S_a are also shown to evaluate the effect of ${\mathbf{S}}_{\mathrm{a}}^{\mathrm{A}}$ on the retrieval result (Figs. 4, 5c). Figure 4 depicts the retrieved τ, η_f, and ${\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}$ 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 ${\mathit{\sigma }}_{{\mathit{\tau }}_{\mathrm{a}}}$, ${\mathit{\sigma }}_{{\mathit{\eta }}_{{\mathrm{f}}_{\mathrm{a}}}}$, and ${\mathit{\sigma }}_{{{\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}}_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 (Yumimoto et al., 2018). 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 ${\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}$ 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 ${\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}$ 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.

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 and King, 2000; Dubovik et al., 2000; Dubovik et al., 2002a, b, 2006; Sinyuk 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 7 Aerosol optical thickness at 500 nm τ (upper), external mixing ratio of dry-volume concentration of fine particles η_f (middle), and external mixing ratio of the dry-volume concentration of dust particles for the coarse model ${\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}$ (lower) that are (a) retrieved from the original algorithm (i.e., using constant an a priori estimate), (b) retrieved from this algorithm, and (c) from the model forecast at 06:40 UTC on 29 June 2018. (d) Standard deviations of model forecast (${\mathit{\sigma }}_{{\mathit{\tau }}_{\mathrm{a}}}$, ${\mathit{\sigma }}_{{\mathit{\eta }}_{{\mathrm{f}}_{\mathrm{a}}}}$, and ${\mathit{\sigma }}_{{{\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}}_a}$) used for retrieval in (b). Red circles in (a) τ and (b) τ indicate the results for the red asterisks in Fig. 6a, b, and c.

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 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 ${\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}$ even where the original algorithm failed to retrieve them. Thus, integrating the model and satellite information resulted in an improvement of the τ estimations.

Figure 8 Frequency distribution of τ retrieved from AHI and those from AERONET. The results retrieved from this algorithm in the case of χ² 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.

Figure 9 Frequency distribution of α retrieved from AHI and those from AERONET. The results retrieved from this algorithm in the case of χ² less than 20, 0.5, and 0.2 and uncertainties of the retrieved η_f (${\mathbf{S}}_{{\mathit{\eta }}_{\mathrm{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.

Figure 10 Frequency distribution of ω retrieved from AHI and those from AERONET. The results retrieved from this algorithm in the case of χ² less than 20, 0.5, and 0.2 and uncertainties of the retrieved ${\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}$ (${\mathbf{S}}_{{\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}}$) 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.

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.

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 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 validation 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 (χ²) and the uncertainties of the retrieved three parameters (τ, η_f, and ${\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}$) are smaller than a threshold. The chi-square value (χ²) is calculated as follows:

$$\begin{array}{}\text{(7)}& {\mathit{\chi }}^{\mathrm{2}}\left[\mathit{R}-F\left(\mathit{x}\right)\right]={\left[\mathit{R}-F\left(\mathit{x}\right)\right]}^T{\mathbf{S}}_{\mathrm{e}}^{-\mathrm{1}}\left[\mathit{R}-F\left(\mathit{x}\right)\right]/N.\end{array}$$

It shows the closeness of the retrieved value to the observed value. The covariance matrix of the uncertainties of the retrieved parameters ${\mathbf{S}}_{\widehat{x}}$ is calculated using the law of error propagation, as follows:

$$\begin{array}{}\text{(8)}& {\mathbf{S}}_{\widehat{x}}={\left({\mathbf{A}}^T{\mathbf{S}}_{\mathrm{e}}^{-\mathrm{1}}\mathbf{A}\right)}^{-\mathrm{1}},\end{array}$$

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 (Yoshida et al., 2018). 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 ${\mathbf{S}}_{\widehat{x}}$ is mostly caused by the surface reflectance uncertainty. Figure 8 shows that RMSE for τ decreases as the threshold of χ² or S_τ becomes strict (i.e., decreases). On the other hand, RMSE for α (in Fig. 9) is not dependent on the threshold of ${\mathbf{S}}_{{\mathit{\eta }}_{\mathrm{f}}}$ but decreases as the χ² threshold decreases. RMSE for ω (in Fig. 10) is almost independent on the threshold of ${\mathbf{S}}_{{\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}}$ and χ². 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 and the assimilation method and increasing the assimilation frequency may further improve the retrieval accuracy in the future.

Figure 12 Same as Fig. 4, except for using the forecast on 27 April 2018 as an a priori estimate.

3.3 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 ${\mathbf{S}}_{\mathrm{a}}^{\mathrm{E}}$ is used as S_a (Fig. 12c), all parameters (especially η_f and ${\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}$) 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 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).

4 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 assimilation 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 (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). ${\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}$ 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:

$$\begin{array}{}\text{(A1)}& {\displaystyle \frac{\mathrm{d}V\left(r_{\mathrm{d}}\right)}{\mathrm{d}\ln r_{\mathrm{d}}}}={\displaystyle \frac{C_{\mathrm{v}}}{\sqrt{\mathrm{2}\mathit{\pi }}\ln \mathit{\sigma }}}\exp \left[-{\displaystyle \frac{(\ln r_{\mathrm{d}}-\ln r_{\mathrm{v}}{)}^{\mathrm{2}}}{\mathrm{2}{\ln }^{\mathrm{2}}\mathit{\sigma }}}\right],\end{array}$$

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 non-spherical 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 $\mathrm{3.0}\times {\mathrm{10}}^{-\mathrm{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 ${\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}$ such that the fine and coarse models exhibit the same ω at 500 nm.

Appendix B: Relationship of α and ω with η_f and ${\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}$

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 (${\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}$). The ω at 500 nm can be uniquely determined by the ${\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}$ (Fig. B1a), since ${\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}$ 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 ${\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}$ but also η_f. The α is mainly determined by η_f but also depends slightly on ${\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}$ (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 (${\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}$). Each color represents a different ${\mathit{\eta }}_{\mathrm{c}}^{\mathrm{dst}}$.