Assimilation and Forecasting Experiment for Heavy Siberian Wildfire Smoke in May 2016 with Himawari-8 Aerosol Optical Thickness

The Japan Meteorological Agency (JMA) launched a next-generation geostationary meteorological satellite (GMS), Himawari-8, on October 7, 2014, which began its operation on July 7, 2015. The Advanced Himawari Imager (AHI) onboard Himawari-8 has 16 observational bands that enable the retrieval of full-disk maps of aerosol optical properties (AOPs), including aerosol optical thickness (AOT) and the Ångström exponent (AE), with unprecedented spatial and temporal resolutions. In this study, we combined an aerosol transport model with the Himawari-8 AOT using the data assimilation method and performed aerosol assimilation and forecasting experiments on smoke from an intensive wildfire that occurred over Siberia between May 15 and 18, 2016. To effectively utilize the high observational frequency of Himawari-8, we assimilated 1-h merged AOTs generated through the combination of six AOT snapshots taken over 10-min intervals, three times per day. The heavy smoke originating from the wildfire was transported eastward behind a low-pressure trough and covered northern Japan from May 19 to 20. The southern part of the smoke plume then traveled westward, in a clockwise flow associated with high pressure. The forecast without assimilation reproduced the transport of the smoke to northern Japan; however, it underestimated AOT and the extinction coefficient compared with observed values mainly because of errors in the emission inventory. Data assimilation with the Himawari-8 AOT compensated for the underestimation and successfully forecasted the unique C-shaped distribution of the smoke. In particular, the assimilation of the Himawari-8 AOT in May 18 greatly improved the forecast of the southern part of the smoke flow. Our results indicate that the inheritance of assimilation cycles and the assimilation of more recent observations led to better forecasting in this case of a continental smoke outflow. ©The Author(s) 2018. This is an open access article published by the Meteorological Society of Japan under a Creative Commons Attribution 4.0 International (CC BY 4.0) license (http://creativecommons.org/license/by/4.0). Corresponding author: Keiya Yumimoto, Research Institute for Applied Mechanics, Kyushu University, 6-1 Kasugakoen, Kasuga-city, Fukuoka 816-8580, Japan E-mail: yumimoto@riam.kyushu-u.ac.jp J-stage Advance Published Date: 8 April 2018 Journal of the Meteorological Society of Japan, Vol. 96B, pp. 133−149, DOI:10.2151/jmsj.2018-035, 2018 Special Issue on Meteorology and Climate Change Studies by Using the Geostationary Meteorological Satellite Himawari-8 Journal of the Meteorological Society of Japan Vol. 96B 134


Introduction
Airborne aerosols play an important role in air quality and human health. The World Health Organization (WHO) has suggested that 6.5 million deaths per year may be associated with air pollution exposure and its health impact (WHO 2016). It is well known that aerosols affect Earth's radiation budget by scattering and absorbing radiation (direct effect) and modifying the cloud properties (indirect effect). This aerosol effect on the radiation balance is crucial not only for estimating climate change (i.e., long-range effects; IPCC 2013) but also for weather prediction (i.e., short-range effects). Various studies have reported the influence of accurate simulation of aerosols and their radiative effects on the accuracy of numerical weather prediction (NWP; e.g., Rodwell and Jung 2008;Grell and Baklanov 2011;Morcrette et al. 2011;Mulcahy et al. 2014;Toll et al. 2016). For example, Pérez et al. (2006) took the radiative effects of African mineral dust into account and showed an improvement in temperature and mean sea-level pressure forecasts.  applied the chemistry version of the Weather Research and Forecasting (WRF-Chem) model to intense wildfires during the 2004 Alaska wildfire season and found that the interaction of aerosols from wildfires with atmospheric radiation led to significant modifications in vertical temperature and moisture profiles.
To monitor air quality, provide early warning of heavy pollution events, and improve the accuracy of NWP, operational forecasting centers have been developing and operating aerosol forecasting systems that forecast the three-dimensional (3D) distribution of aerosols. To share best practices and discuss pressing issues facing the operational aerosol community, multi-model inter-comparison projects, including the International Cooperative for Aerosol Prediction (ICAP; Sessions et al. 2015), and the World Meteorological Organization (WMO) Sand and Dust Storm Warning Advisory and Assessment System (SDS-WAS; Terradellas et al. 2015) were established by operational forecasting centers. A variation in multimodel forecasting indicated that large uncertainties remained in aerosol forecasting, requiring further research for improved accuracy (Sessions et al. 2015). To this end, data assimilation techniques have been incorporated into aerosol transport models to obtain better initial conditions for forecasting. Both variation-based methods (e.g., Benedetti et al. 2009;Zhang et al. 2008;Wang and Niu 2013) and ensemble-based methods (e.g., Di Tomaso et al. 2016;Rubin et al. 2016;Yumimoto et al. 2016b) have been used by forecasting centers. More recently, Rubin et al. (2017) have performed assimilation and forecasting experiments to determine aerosol optical thickness (AOT) from ground-and satellite-based observations with the Navy Aerosol Analysis Prediction System (NAAPS) using two-dimensional variation (2D-Var) methods and ensemble-based assimilation systems; their study results showed that ensemble-based assimilation is a more suitable method for incorporating spatially sparse observations (i.e., ground-based networks) than 2D-Var assimilation.
The Meteorological Research Institute (MRI) of the Japan Meteorological Agency (JMA) has developed an online aerosol transport model coupled with an atmospheric general circulation model (AGCM), the Model of Aerosol Species IN the Global AtmospheRe (MASINGAR; Tanaka and Chiba 2005). JMA has been using MASINGAR for operational sand and dust storm forecasting since January 2004. JMA launched Himawari-8, the eighth Japanese geostationary meteorological satellite (GMS), on October 7, 2014, and started its operation on July 7, 2015. The Advanced Himawari Imager (AHI) onboard Himawari-8 has 16 observational bands over the visible and infrared wavelength ranges, with higher spatial (0.5 -2 km) and temporal (full-disk image every 10 min) resolutions compared with those of previous GMSs (Bessho et al. 2016). In particular, the three visible bands (0.47, 0.51, and 0.64 μm) and one near-infrared (0.86 μm) band are sensitive to aerosol scattering and absorption and enable the retrieval of aerosol optical properties (AOPs), including AOT and the Ångström exponent (AE), over wide areas of Asia and Oceania with unprecedented frequency from geostationary orbit.
Our research team previously developed an aerosol data assimilation and forecasting system based on MASINGAR mk-2 and an ensemble-based assimilation technique called the Local Ensemble Transform Kalman Filter (LETKF; Hunt et al. 2007) and reported successful results with AOPs from the Moderate Resolution Imaging Spectroradiometer (MODIS) and Himawari-8 (Yumimoto et al. 2016a, b). However, the ensemble-based assimilation system has a relatively Keywords aerosol; aerosol transport model; data assimilation; Himawari-8; wildfire smoke high computational cost, and the development of operational and quasi-real-time aerosol forecasting remains unrealistic. Therefore, in addition to the ensemble-based system, we developed an aerosol data assimilation and forecasting system based on the 2D-Var method called MASINGAR/2D-Var (Yumimoto et al. 2017) for operational assimilation and forecasting.
In this study, we present the integration of MASINGAR mk-2 and the Himawari-8 AOPs through the 2D-Var assimilation to improve aerosol forecasting. Assimilation and forecasting experiments were performed, targeting smoke from a severe wildfire that occurred over Siberia from May 15 to 18, 2016 (the east side of Lake Baikal). The aerosol model, the data assimilation system, the observational data used in assimilation and validation, and the experimental conditions are outlined in Section 2. In Section 3, we briefly describe the emission and transport of smoke. Section 4 presents the assimilation and forecasting results and the results of their validation. Finally, conclusions are summarized in Section 5.

Aerosol transport model
We used MASINGAR mk-2, which is the first major update of the global aerosol transport model MASINGAR developed at MRI/JMA. MASINGAR mk-2 covers the major tropospheric aerosol components (i.e., black and organic carbon, mineral dust, sea salt, and sulfate aerosols) and their precursors (e.g., sulfur dioxide and dimethyl sulfide) and is coupled online with an atmospheric general circulation model (MRI-AGCM3; Yukimoto et al. 2012). MASINGAR mk-2 took over from the original MASINGAR for operational sand and dust forecasts in November 2014, and it serves as a member of the ICAP multi-model ensemble (Sessions et al. 2015). Emissions of mineral dust and sea salt aerosols are calculated online using schemes developed by Shao et al. (1996) and Gong (2003), respectively. We used the MACCity (MACCC/ City ZEN EU projects) emission inventory (http:// ether.ipsl.jussieu.fr/eccad) and the Global Fire Assimilation System (GFAS) data set (http://www.gmesatmosphere.eu/about/project_structure/input_data/ d_fire) for anthropogenic and biomass burning emissions, respectively. More detailed information about MASINGAR mk-2 can be obtained in Yumimoto et al. (2017) and Yukimoto et al. (2012).

Data assimilation system
The cost function in the 2D-Var method is defined where τ denotes the modeled AOT vector (a twodimensional diagnostic variable); the suffix f represents the first guess (a priori); τ o denotes a vector that contains the observed AOT used for the assimilation; P τ and R are the background error covariance matrix for AOT and the observation error covariance matrix, respectively; and H I is the interpolation into observation space and a linear operator. The analysis increment of AOT (δτ a ) is derived from Eq. (1) as where τ a is the analysis (a posteriori) ΑΟΤ, which provides the minimum of the cost function; τ τ τ τ where δ x a l, k is the analysis increment of the mixing ratio of the l-th aerosol component and the k-th vertical layer, and x f l, k and α f l, k represent the first-guess mixing ratio and the extinction coefficient of the l-th aerosol component and k-th vertical layer, respectively. Finally, the analysis (a posteriori) mixing ratio (x a l, k ) is given by To include observed information, forecasting is initialized by the analysis (a posteriori) aerosol field. The timeline for the operational assimilation and forecasting is shown in Fig. 1. One cycle consists of a one-day analysis run and a five-day forecasting run. During the analysis run, the model experiences data assimilation with observations. We performed assimilations three times during the daytime of East Asia (at 03:00, 06:00, and 00:00 Universal Time (UTC)) with the Himawari-8 AOT in the analysis run. The five-day forecasting was then initialized with the results of the analysis run. In the subsequent cycle, the assimilation run used the results from the previous analysis run as the initial condition.
Background error covariance is an important parameter that weights the first term of the cost function (Eq. 1) and propagates observed information to regions far from the observed point through its non-diagonal components (i.e., correlation). In this study, we developed a method similar to the National Meteorological Center (NMC) method (Parrish and Derber 1992) and defined the background error covariance between model grids m and n (i.e., the (m, n) element of the covariance matrix) as where τ m and τ n are the AOTs at model grids m and n, f refers to the first guess, i refers to the i th ensemble member, N is the number of the ensemble, and C m, n is the smooth weighting function. To prevent a spurious large covariance between distant model grids, we introduced the smooth weighting function C m, n , whose value diminishes as the distance between model grids increases, and a localization technique (Yumimoto et al. 2017). We employed the second-order autoregressive (SOAR) approximation (Daley and Barker 2001) to calculate the value of C m, n : where R m, n denotes the distance (km) between model grids m and n, and L is the horizontal error correlation length. The localization technique divides the model space into local regions using a prescribed localization scale to reduce spurious error covariance with distance and enables parallel processing to be used to reduce computational cost. The AOT ensemble (τ m, i and τ n, I , i = 1, …, N ) was collected from AOT values within ± 5 h of the targeted hour of the five previous forecasts. For example, for the assimilation at 06:00 UTC on May 18, AOT values at 01:00, 02:00, 03:00, 04:00, 05:00, 06:00, 07:00, 08:00, 09:00, 10:00, and 11:00 UTC on May 18 from forecasts started on May 14, 15, 16, 17, and 18 were assigned as the AOT ensemble. In this case, the ensemble number is 55 (11 × 5). Figure 2 shows examples of variance and covariance at 00:00 UTC on May 18 and 19. Distributions of covariance (Figs. 2c, f) were similar to those of the wildfire smoke represented by the first guess AOT (Figs. 2a,d). In the case of May 18 (Fig. 2, upper panels), the aerosol plume around 45°N, 136°E has a relatively large correlation with another aerosol plume around 47°N, 126°E (approximately 800 km apart). This result implies that both plumes originated from the same source.

Observation data a. Himawari-8 data
The AOT retrieved from AHI onboard Himawari-8 was used for the assimilation. The Earth Observing Research Center (EORC) of the Japan Aerospace Exploration Agency (JAXA) processes Himawari-8 data and provides Himawari-8 AOPs via the JAXA Himawari Monitor (http://www.eorc.jaxa.jp/ptree/ index.html). Through this system, two levels of AOP products (L2 and L3 products) are available. The L2 product includes full-disk AOT and AE at a wavelength of 500 nm every 10 min. The details of the AOP retrieval algorithm for the L2 product are provided in Yoshida et al. (2017). The L3 product consists of hourly AOPs derived from a combination of the 10-min AOPs (L2 product). Taking advantage of Himawari-8's high-frequency observations, the combination enables minimization of the number of pixels missing because of cloud and sunglint and removes degradation due to cloud contamination. Figure 3 shows a map of AOT from the L2 product and merged AOT from the L3 product at 03:00 UTC on May 19, 2016. The details of the combined algorithm for the L3 product are provided in Kikuchi et al. (2018). In this study, to efficiently utilize high-frequency Himawari-8 observations, we applied merged AOT in the L3 product (Version 1.0) as assimilation data. Because only AOT at 550 nm is not available in the current version of the system, we extrapolated the merged AOTs at 500 nm into AOTs at 550 nm using the AE of the L3 product. The extrapolated AOTs were averaged in the model grid and then used in the assimilation. We also used wildfire detection data provided by the JAXA Himawari Monitor to identify the source region of the wildfire smoke.

b. Validation data
We used observation data from satellite-based and ground-based measurements that were not used in the data assimilation to analyze the wildfire smoke and validate the forecasting results. MODIS onboard Terra and Aqua low Earth orbit (LEO) satellites provide AOPs once per day (at approximately 10:30 local time (LT) for Terra and approximately 13:30 LT for Aqua).  We used the "optical depth land and ocean mean" contained in the Collection 6.0 L2 aerosol data (Levy et al. 2013(Levy et al. , 2015. Lidar observations provide vertical profiles of suspended layers in the wildfire smoke.

Experimental conditions
To estimate the impact of data assimilation with the Himawari-8 AOPs on the forecast of wildfire smoke, we performed one free run without assimilation (hereafter referred to as noDA) and two cycles of assimilation and forecast runs (Fig. 1). The first cycle (DA1) started the analysis run at 00:00 UTC on May 17, 2016, and then executed the forecast run with the analyzed aerosol fields from 00:00 UTC on May 18. We performed a subsequent cycle (DA2) to investigate the impact of the number and timing of assimilations on forecasting. DA2 was initialized by the analyzed DA1 results, from the analysis run on May 18, and we then ran a five-day forecast.
We performed assimilations with the merged AOTs three times (at 03:00, 06:00, and 00:00 UTC) in analysis runs (Fig. 1). The merged AOTs used for assimilation were limited to ocean extents (merged AOTs over the land were not used for assimilation) because Himawari-8 AOPs over land have larger retrieval uncertainties than those over the ocean (see Fig. 3c). We assigned the uncertainty estimate provided by the L3 product to the observation error. The uncertainty estimate includes the errors from the random noise of the sensor and surface reflectance that varies in time and space. The model horizontal resolution was set to TL479 (approximately 50 km) with 40 vertical layers from the ground to 0.4 hPa using the hybrid sigma pressure coordinate system. The operational global analysis provided by JMA (GANAL/JMA; Japan Meteorological Agency 2002) at 6-h intervals was used for AGCM nudging. The nudging scheme is applied to the horizontal wind components and air temperature. The nudging term is applied at each time step of the integration by temporary interpolating the variables by linear interpolation. For the assimilation procedure, we set the localization scale and horizontal error correlation length to 750 and 200 km, respectively, based on the results from Zhang et al. (2008).

Wildfire smoke from Siberia
During May 15 -18, 2016, a heavy wildfire occurred to the eastern side of Lake Baikal. Figure 4 shows the GFAS wildfire organic carbon flux and Himawari-8 wildfire detection on May 16. The Himawari-8 wildfire detection was consistent with the GFAS wildfire flux, which was based on MODIS observations. Smoke emitted from the wildfires traveled eastward to the Sea of Japan behind of a low-pressure trough. The CALIOP observations show that the smoke was transported over the east coast of Siberia at an altitude of 2 -8 km ( Fig. 5; CALIOP observation path shown in Fig. 4). Satellite observations captured the spread of the smoke behind a cold front located along the Japanese archipelago on May 18, as shown in the upper panels of Fig. 6. The merged AOT from Himawari-8 is consistent with the AOT from MODIS. The observed AOT exceeds 1.8 around the core of the smoke plume. The AERONET AOT at Hokkaido University (Fig. 7a) observed a rapid increase of AOT corresponding to the motion of the smoke on May 18; the observed AOT reached 1.4, which is consistent with satellite observations. On the following day (May 19), the smoke covered the island of Hokkaido and the northern part of Honshu (the main island of Japan; Fig. 8, upper panels). The Himawari-8 AOT (Fig. 8a) shows a C-shaped smoke distribution. Ground-based observations also captured the smoke. The AD-net Lidar at Sendai observed the dense aerosol layer at 1 -5 km altitude on May 19 (Fig. 9a). At Niigata, the AD-net Lidar (Fig. 9e) detected the main body of the smoke; however, the extinction coefficient was lower than that of Niigata because of dispersion. On May 19 and through May 20 at 0.5 -4.5-km altitude (thickness: ~4 km), consistent with the results observed by AER-ONET, the Niigata AERONET showed an increase in AOT because of the smoke on May 19 (Fig. 7b).
On May 20 -21, smoke transport proceeded in two directions. The northern part of the smoke traveled eastward to the northern Pacific Ocean, whereas the southern part was transported from east to west across Honshu by clockwise flow associated with a highpressure system centered around Hokkaido (approximately 43°N, 142°E; Fig. 10d). This is a plausible explanation for AD-net Lidar at Matsue station (in western Japan) capturing a suspended smoke layer at 1 -3 km on May 21 (Fig. 9i). The blue line in Fig. 4 shows the backward trajectory initiated at 18:00 UTC on May 21 at Matsue at an altitude of 2 km (red circle in Fig. 9i). The air mass that traveled over the source region on May 16 was transported eastward over Siberia, reached Hokkaido on May 18, and was then transported westward over Honshu by a clockwise flow. The arriving times of the smoke predicted by the model simulation and the backward trajectory analysis are consistent with the observation results. Figure 6 compares the horizontal AOT distributions on May 18 between satellite observations and noDA model results. The root-mean-square error (RMSE), linear Pearson correlation coefficient (R), mean fractional error (MFE), mean fractional bias (MFB), and index of agreement (IOA) calculated over the land are also shown in Fig. 6 (the formulations of these statistical measures are given in Appendix A). The transport of the smoke was captured by the model (R and IOA are 0.66 and 0.60, respectively); however, noDA (Fig. 6d) underestimated the observed AOT level and the extent of smoke in the Sea of Japan and the Sea of Okhotsk (Fig. 6f) Fig. 4). White/blue colors in (a) indicate relatively high/low total attenuated backscatters.

Hokkaido Island
The main island of Japan  forecast for several days. On May 19 (Fig. 8), noDA underestimated the extent of the smoke over Hokkaido and the northern part of Honshu and did not reproduce the C-shaped distribution observed by satellites (see Figs. 8d, g); the modeled AOTs were 0.2 -0.6, whereas the observed AOTs exceeded 0.8, and the MFB of noDA shows a large negative value (−73.2 %). On May 20 (Fig. 10), the AOT forecast by noDA was much lower over Hokkaido and the Sea of Japan than that observed by satellite (see Fig. 10g); the MFB is −66.4 %. Comparisons between the AERONET AOT (Fig. 7) and AD-Net observations (Fig. 9) also showed that noDA could forecast the advent of the smoke, but with very low AOT and extinction coefficient values; the maxima of the AERONET (noDA) AOT at Hokkaido University and Niigata are 1.47 (0.51) and 0.86 (0.22), respectively. We also performed a free run with lower spatial resolution (TL319; approximately 60 km) and found that the forecast smoke distribution was quite similar between the fine and coarse models and that the coarse model also underestimated the observed AOT and extinction coefficient (data not shown). This fact implies that underestimation was mainly attributed to uncertainty associated with emis-sion rather than uncertainty associated with transport (i.e., meteorological fields). Uncertainty in the spatial and temporal distributions and intensity of emissions for aerosols and their precursors is a major source of error in aerosol forecasting. This is particularly true for aerosols from natural origins (e.g., mineral dust from desert, smoke from wildfire, and ash and sulfate aerosols from volcanic eruptions); their episodic and explosive occurrences make the accurate estimation of emission difficult and result in large errors in forecasting (e.g., Uno et al. 2006). In this case, the introduction of observed information in the model simulation via assimilation can be effective in improving forecasting performance. The analyzed AOT from DA1 at 00:00 UTC on May 18 is shown in Fig. 6e. It is clear that data assimilation with Himawari-8 AOT improved the smoke distribution model. Figure 8e shows a 24-h DA1 forecast (FT = 24) of AOT at 00:00 UTC on May 19 (the 24-h forecast started from the analysis aerosol field at 00:00 UTC on May 18; see Fig. 1); the assimilation improved the smoke forecast (all statistical measures improved compared with noDA). Specifically, the assimilation compensated for the underestimation and reproduced the C-shaped distribution (the MFB value was reduced from −73. 2 % to −38.1 %); however, the lower bias remained in the southern part of the smoke plume (around 38 -40°N, 140 -150°E) and persisted beyond the 24-h forecast (see Fig. 8h). Although the 48-h DA1 forecast (FT = 48; 00:00 UTC on May 20; Fig. 10e) shows better agreement and statistics with the satellite observations than noDA, the forecast AOT underestimated observed AOT values (the MFB is −45.4 %). This underestimate was the result of the lower bias in the southern part of the smoke plume on May 19. Although the AERONET AOT (Fig. 7) and AD-Net observation (Fig. 9) confirmed the improvement in the 24 -72-h forecasts due to assimilation, the DA1 forecast AOT and extinction coefficient were still lower than the observed values.
We performed a subsequent assimilation and forecasting cycle from the DA1 analysis to introduce the Himawari-8 AOT May 18 observations into the forecast (Fig. 1). The results of the subsequent cycle are labeled "DA2". Figure 8f shows the analyzed AOT from DA2 at 00:00 UTC on May 19. Subsequent assimilation improved the modeled AOT; in particular, the assimilation compensated for the underestimation in the southern part of the C-shaped smoke distribution (the MFB value was further reduced from −38.1 % to −13.7 %; compare Figs. 8h, i). Because the Himawari-8 AOTs over the land were excluded from    Fig. 9. Time-height plots of extinction coefficients (1 km −1 ) at Sendai, Niigata, and Matsue: (a, e, and i) AD-Net Lidar observations, (b, f, and j) modeled results from the noDA run, (c, g, and k) modeled results from the DA1 run, and (d, h, and l) modeled results from the DA2 run. Red vertical lines denote the data assimilation times. Red/blue colors indicate relatively high/low extinction coefficients.   (approximately 45°N, 132°E) in the Himawari-8 observation (Fig. 8a) were not reflected in the DA2 analysis and forecast. On the following day (May 20; Fig. 10), the 24-h DA2 forecast (FT = 24;Figs. 10f,i) improved the underestimation observed in the noDA forecast (Figs. 10d, g) and the 48-h DA1 forecast (Figs. 10e,h) and resulted in much better agreement with satellite-observed AOTs. DA2 achieved better scores in all statistics than both noDA and DA1. This improvement was verified in a comparison with AERONET observations (Fig. 7); the DA2 forecast quantitatively agreed with the AERONET AOT at both Hokkaido University and Niigata. The high AOT values early on May 20 at Hokkaido University ( Fig. 7a) may be attributed to another smoke plume transported after the smoke examined in this study and not captured by the Himawari-8 AOT used in the assimilation. This is a plausible reason for the DA1 and DA2 assimilations providing a slight improvement in the high-AOT forecast. Compared with AD-Net observations (Fig. 9), the DA2 forecast exhibited improvement with respect to the underestimation and much better agreement with the observations than noDA and DA1 forecasts. In particular, a better agreement in the 48 -96-h DA2 forecast at Matsue on May 20 -22 (Fig. 9l) indicates that the assimilation improved the forecast of the southern part of the smoke plume transported across Honshu by a clockwise flow and improved the performance in relatively long-term (48 -96-h) forecasting.
The improved performance of DA2 over DA1 suggests that the inheritance of assimilation cycles and the assimilation of more recent observations provide a more accurate forecast. Additionally, our results imply that the benefits of the assimilation were limited to short-term (i.e., 1 -24-h) forecasts; this is particularly true for DA1. The DA1 assimilation run on May 17 improved the smoke forecast on May 18 (i.e., a 1 -24-h forecast). However, in the forecasts performed on the following days, the benefit of assimilation was insufficient, and the concentration of the modeled smoke underestimated that of the observed smoke. We excluded the Himawari-8 observations over the land from the assimilation. Therefore, the assimilation data were obtained only around the Japanese archipelago (i.e., downwind of the smoke). This is a plausible explanation for the limited benefit provided by the assimilation to long-term forecasts over Japan. Including the Himawari-8 AOT observed over the land (i.e., closer to the source region) in the assimilation will further improve the forecasts, especially long-term forecasts.
We performed additional assimilation and forecasting experiments (DA1′ and DA2′) in which the Himawari-8 AOTs were assimilated only at 00:00 UTC (i.e., one assimilation in one assimilation run) and found a little difference between DA1 and DA1′ and between DA2 and DA2′ (data not shown). This result implies that the most recent assimilations (in this case, 00:00 UTC) had a significant impact on the forecast and that multiple assimilations (in this case, three times per day) had a little effect on the forecast performance for the smoke modeled in this study. One reason for this result may be the schedule of the assimilation run. In DA1 and DA2, we ran the assimilations at 03:00 UTC and 06:00 UTC, took an 18-h interval, ran the assimilation at 00:00 UTC, and then began the forecasts. The relatively long interval between the assimilations at 06:00 UTC and 00:00 UTC increased the importance of the assimilations at 00:00 UTC in these forecasts. The limitation of the ocean assimilation data (i.e., the downwind region of the smoke) may provide an alternative explanation for this. Utilizing the Himawari-8 AOTs over both land and ocean in the analysis will enable multiple assimilations of the smoke in both the downwind and upwind regions (i.e., near the source region) to further improve forecasting accuracy.

Conclusion
The AHI onboard Himawari-8, the eighth Japanese GMS put into operation on July 7, 2015, can produce AOPs, including AOT and AE, with wide vision and unprecedented temporal resolution (every 10 min). In the present study, we developed an aerosol assimilation and forecasting system, with which we assimilated the Himawari-8 AOT multiple times per day, and performed assimilation and forecasting experiments for a heavy wildfire smoke event that occurred on May 15 -18, 2016, over Siberia. To take advantage of the high observational frequency of the Himawari-8, we used merged AOT, which was derived from a combination of 10-min AOTs, minimizing the number of pixels missing because of cloud and sunglint and degradation by cloud contamination. Our results are summarized below. 1. The dense smoke emitted by a wildfire during May 15 -18, 2016, over the eastern side of Lake Baikal traveled eastward behind a low-pressure trough, covered northern Japan (Hokkaido and the northern part of the main island of Japan) at an altitude of 1 -5 km on May 19, and then branched into two parts. The northern part of the smoke was transported further eastward to the northern Pacific Ocean. The southern part traveled across the main island of Japan in a clockwise flow associated with a high-pressure system and was observed in the western part of Japan. 2. The free run without assimilation (noDA) forecast the observed transport of the smoke, but underestimated its concentration, mainly because of errors in the emission inventory. The data assimilation during May 17 (DA1) compensated for underestimation in the forecast and reproduced the unique C-shaped distribution that noDA failed to predict. However, the AOT and extinction coefficient forecast by DA1 remained lower than those observed. The subsequent assimilation and forecast cycle (DA2), which assimilated the Himawari-8 AOT during May 18, was greatly improved and successfully predicted the transport of the smoke to the western part of Japan in its 48 -96-h forecast. 3. That DA2 performed better than DA1 in forecasting the smoke movement indicates that the succeeding assimilation cycles and the assimilation of the most recent observations resulted in better forecasting.
Our results also implied that the assimilation had an impact only on relatively short-range (i.e., 1 -24hour) forecasts of the smoke in this study. The limitations of the assimilation data over the ocean (i.e., data only downwind of the smoke) may be a plausible explanation for this result. Further studies targeting other aerosol events than wildfire smoke (e.g., dust storms and trans-boundary pollutions) and including Himawari-8 data over land are required.
The assimilation system has been developed for the sand and dust storm forecasting operated in JMA. Our results clearly show that data assimilation with AOTs from Himawari-8 successfully improved the initial condition and then brought much better forecast. Although this study targeted an intensive wildfire smoke case, our results indicate that AOTs from Himawari-8 are promising not only for monitoring but also for improved forecasting of other aerosol events, including dust storms, through data assimilation. In future studies, we will perform assimilation and forecast experiments for Asian dust and storms and continental anthropogenic pollutions. AOT over land (closer to continental outflow source regions) and AE that are sensitive to aerosol size distribution will be utilized in assimilation during the development process of these forecast models.

Appendix A: Statistical metrics
We used the root-mean-square error (RMSE), the Pearson correlation (R), the mean fractional error (MFE), the mean fractional bias (MFB; also known as the fractional gross error), and the index of agreement (IOA) for validation. The statistical metrics are defined as follows: , respectively. RMSE represents the standard deviation of the discrepancies between modeled and observed values. MFE can range from 0 % to 200 % and is a measure of overall modeling error without emphasizing outliers. MFE = 0 is a perfect score. MFB is a measure of the estimation bias error that allows symmetric analysis of over-or underestimation by the model relative to observed values. The best value of MFB is 0 with ±200 % of the minimum and maximum values. Boylan and Russell (2006) proposed that MFE should be ≤ +50 % and MFB should be ≤ ±30 % to meet model performance goals for particulate matter (PM) and light extinction. IOA was developed by Willmott (1981) as a standard measure of the degree of model prediction accuracy, and it ranges from 0 to 1. IOA = 1 indicates perfect agreement.
investigators and those who contributed to the establishment and maintenance of the Himawari-8 AOPs, the AERONET sites, CALIOP/CALIPSO, the NOAA HYSPLIT trajectory model, and the NASA MODIS L2 aerosol products.