Assessment of AOD Simulation by Assimilating Himawari and MODIS Products

1.   Introduction

In recent years, with increases in economic growth and urbanization, air pollution events have become more frequent. Aerosols, a key atmospheric pollutants, consist of fine solid particles or liquid droplets suspended in the air, with aerodynamic diameters typically ranging from 0.001 to 100 micrometers (Bao et al., 2019). East Asia and its surrounding regions are some of the world’s primary sources of aerosol emissions, with aerosol concentrations generally higher than those in most other regions (Cheng et al., 2022). Due to their complex distribution and chemical composition (Wu et al., 2013), aerosols have significantly affected weather, climate, the environment, and public health (Qian et al., 2009; Zhang, 2014; Li et al., 2019; Marina-Montes et al., 2020; Zhao et al., 2020, 2023, 2024). China's air quality has significantly improved over the past decade because of the implementation of stringent air pollution control policies (Zhang et al., 2019; Fan et al., 2020). However, the annual average concentration of fine particulate matter in regions such as Beijing Tianjin Hebei and Yangtze River Delta still exceeded the second level concentration limit of the ambient air quality standards (35 μg m⁻³) (Zheng et al., 2015; Wu et al., 2018).

Aerosol optical depth (AOD) is a fundamental index used to describe the optical properties of atmospheric aerosols. It plays a crucial role in assessing atmospheric pollution levels and conducting climate-related studies (Quaas et al., 2008). AOD observations are primarily obtained through ground-based measurements and satellite remote sensing retrievals (Holben et al., 1998). While ground-based AOD observations offer continuous measurements at high temporal resolutions and high accuracy, their spatial coverage is limited. In contrast, satellite remote sensing methods can provide extensive and continuous AOD data. With advancements in remote sensing software and hardware, the accuracy of satellite-derived AOD products has steadily improved (Yang et al., 2019). As a result, satellite remote sensing AOD products are now widely used in air pollution monitoring and forecasting. In addition to observational methods, numerical models can generate high-resolution meteorological and atmospheric composition data. Researchers have conducted extensive studies on air quality, cloud‒aerosol interactions, aerosol optical properties, and their climate effects using numerical models (Li et al., 2009; Luo et al., 2014). However, the spatial and temporal distributions of aerosols simulated by numerical models are prone to errors due to uncertainties in emission inventories, chemical mechanisms, parameterization schemes, etc. (Cheng et al., 2022). Improving the accuracy of aerosol and AOD simulations in numerical models remains a critical and challenging research focus in atmospheric science.

Data assimilation (DA) is a technique used to integrate observed data with model-simulated background fields, aiming to produce the best possible estimate of the atmospheric state (Xia et al., 2019). It is one of the key methods for improving the accuracy of AOD forecasts and reanalysis products (Wang et al., 2004). Common DA methods include optimal interpolation (OI), three-dimensional variational (3DVAR), four-dimensional variational (4DVAR), and the ensemble Kalman filter (EnKF) (Collins et al., 2001; Benedetti et al., 2009; Liu et al., 2011; Dai et al., 2014). Collins et al. (2001) developed a system for assimilating AOD data from the Advanced Very High Resolution Radiometer (AVHRR) satellite using the OI method, successfully simulating the aerosol distribution during the Indian Ocean Experiment (INDOE). The Gridpoint Statistical Interpolation (GSI) system, developed by the National Centers for Environmental Prediction (NCEP), integrates the 3DVAR assimilation technique for Moderate Resolution Imaging Spectroradiometer (MODIS) AOD products, significantly increasing the accuracy of AOD simulations (Liu et al., 2011). Dai et al. (2014) further advanced DA by developing a global aerosol assimilation system using a complex icosahedral grid configuration based on the local ensemble transform Kalman filter (LETKF), optimizing the accuracy of aerosol optical properties. Xia et al. (2023) and Cheng et al. (2022) employed 3DVAR to assimilate data from two satellites simultaneously, achieving better results than those obtained using a single satellite. While more complex assimilation methods, such as 4DVAR and LETKF, have demonstrated excellent performance in various applications, 3DVAR is still a widely adopted DA method in the aerosol assimilation community. Its popularly is attributed to its relatively low computational cost, simpler auxiliary models, and easier integration with complex aerosol processes (Wang et al., 2022).

The limited spatial and temporal coverage of a single satellite can restrict the effectiveness of assimilation efforts. Further evaluation is needed to assess the effectiveness of assimilating multiple satellite data simultaneously. In this study, the Weather Research and Forecasting model coupled with Chemistry (WRF-Chem), combined with the 3DVAR method implemented within the GSI assimilation system, was applied to assimilate AOD products from the Himawari 8 and MODIS satellites over East Asia during the summer of 2017 and the winter 2017/2018. The assimilated results were then compared with reanalysis data from Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA-2), and the accuracy of the WRF-Chem simulations was evaluated using ground-based AOD observations from the Aerosol Robotic Network (AERONET) stations.

2.   Data and Methods

2.1   Data Sources

Himawari 8 is the first next-generation geostationary satellite launched by the Japan Meteorological Agency (JMA) in October 2014. It is equipped with the Advanced Himawari Imager (AHI), which provides cloud and aerosol observations across the Asia–Pacific region (80°E–160°W/60°S–60°N) (Bessho et al., 2016). In this study, we utilized its AOD product at 550 nm (level 3; Version 031) for the summer of 2017 and the winter of 2017/2018 (Kikuchi et al., 2018), which was released by the Japan Aerospace Exploration Agency (JAXA). These data feature a temporal resolution of 1 hour and a spatial resolution of 0.05°, with a detection time range from 20:00 UTC each day to 11:00 UTC the following day. The aerosol optical characterization applied in this study uses a newly developed algorithm with more stringent cloud contamination rejection than the Level 2 data (Yoshida et al., 2018; Tan et al., 2022). Additionally, we used the Level 2 AOD product at 550 nm, provided by MODIS sensors aboard the Terra and Aqua polar-orbiting satellites, for the same months. MODIS aerosol products leverage the combined results of the Dark Target and Deep Blue retrieval algorithms, with a spatial resolution of 10 km and a temporal resolution of 5 minutes (Hsu et al., 2004, 2006; Remer et al., 2005). However, ground-based and satellite AOD should differ with some bias due to their different detection methods (depending on the aerosol scattering and absorption properties) (Li et al., 2014). Previous studies have demonstrated good agreement between the MODIS AOD and ground-based observations (Bilal et al., 2017; Gupta et al., 2018). Compared with other newer satellites, MODIS offers distinct advantages in aerosol remote sensing owing to its dual-satellite constellation configuration and ultrawide scanning swath of 2330 kilometers, enabling superior spatial coverage. Moreover, Himawari-8’s AOD products stand out among geostationary satellites in the Eastern Hemisphere because of their relatively high spatial and temporal resolutions, whereas its multispectral collaborative retrieval technology provides enhanced aerosol retrieval accuracy. Furthermore, both MODIS and Himawari-8 have undergone extensive validation, with their AOD products repeatedly evaluated through ground-based station networks, ensuring higher precision and reliability. In contrast, the uncertainties associated with newer satellite datasets still require more rigorous assessments to establish comparable credibility.

Although both Himawari-8 and MODIS perform well in AOD retrieval, they still exhibit some notable differences. These differences are influenced not only by their distinct orbital configurations but also by variations in their retrieval algorithms. For example, Yumimoto et al. (2016) noted that MODIS has a systematic error of 5% over land and 15% over the ocean, whereas Remer et al. (2005) noted that Himawari 8 has a systematic error of approximately 6%.

In addition, Himawari 8 and MODIS differ significantly in terms of spatial coverage. There are significant instances of missing satellite-derived AOD data due to cloud cover and other factors. Fig. 1 shows the ratio of assimilated AOD information from Himawari 8 and MODIS. Himawari 8, as a geostationary satellite, covers most of the simulation domain; however, the proportion of valid assimilated data averages only 8.3% of the total area. In contrast, while MODIS observations offer nearly full coverage of the domain, their polar-orbiting nature limits each observation to a smaller range, resulting in an average coverage of only 6.5% during the assimilation process. Furthermore, these observations are primarily concentrated in the southern part of the simulation domain. The decrease in effective observations from both Himawari 8 and MODIS in the northern part of the modeled study area during winter may be attributed to the effects of the morning terminator line. At this time, the sun remains over the Southern Hemisphere, leading to insolation deficiency in regions near the Arctic Circle. As a result, these satellites discard a significant number of low-quality observations because the received radiant energy falls below the threshold required for effective data collection and accurate retrieval calculations (Levy et al., 2013).

Download

Larger image

PowerPoint slide

Fig  1.  Ratios of assimilated AOD data from Himawari 8 (a-c) and MODIS (d-f). The areas simulated by WRF-Chem are indicated with purple borders, while regions with no observed AOD data are shown with a gray background. AERONET site locations are marked with a "+" symbol in panel a. The selected AERONET sites for this study included Beijing (BJ), Xuzhou (XZ), Yonsei University, Korea (YS), Baotou (BT), and Kanpur (KP). The Tibetan Plateau, Sichuan Basin, and central and eastern China are highlighted in red, yellow and blue boxes, respectively, in each panel.

In this study, AOD at 550 nm from the MERRA-2 reanalysis data was compared against that simulated by the WRF-Chem model. Additionally, AOD data observed at five AERONET ground stations were employed to evaluate the model’s performance after AOD assimilation. MERRA-2 is the latest generation of reanalysis products developed by NASA's Global Modeling and Assimilation Office (GMAO) (Jiang et al., 2015). It provides AOD reanalysis information from 1980 to the present, with a grid resolution of 0.625° × 0.5° (Gelaro et al., 2017). AERONET is a globally distributed network of sun-sky photometers (Holben et al., 1998). In this study, the latest version 3 data were used, which ensures high precision through fully automated cloud screening and instrumental outlier quality control (Giles et al., 2019). AERONET data exhibit AOD errors of only 0.01–0.02. Thus, they are widely used for assessing aerosol simulation and assimilation (Bao et al., 2019; Giles et al., 2019). The AERONET AOD data at 440 and 675 nm were interpolated to derive the AOD at 550 nm using the Ångström relation (Ångström, 1929). To maintain the diversity between the eastern and western regions, as well as between coastal and inland areas, we selected AERONET sites for this study, including Beijing (BJ), Xuzhou (XZ), Yonsei University, Korea (YS), Baotou (BT), and Kanpur (KP), as shown in Fig. 1a.

Meteorological observations from regional automatic weather stations in China were utilized to assess the accuracy of the meteorological fields simulated by the WRF-Chem model. The data were obtained from the National Meteorological Information Center and included hourly measurements of 2 m temperature (T₂), 2 m relative humidity (RH₂), and 10 m wind speed (WS₁₀).

2.2   3DVAR Data Assimilation Method

In this study, the 3DVAR assimilation system of the GSI was employed to assimilate AOD observations from MODIS and Himawari 8. The 3DVAR assimilation technique involves fusing data by optimizing an objective function (Lorenc, 1986), which can be expressed as follows:

where $\mathit{x}$ represents the analytical field to be solved, which is the state matrix of the optimally estimated variables. ${\mathit{x}}^{\mathit{b}}$ denotes the background field, reflecting the state of the variables on the basis of the forecast prior to assimilation. ${\mathit{B}}_{\mathit{c}\mathit{o}\mathit{v}}$ is the background error covariance matrix, which captures the uncertainties between the variables within the background field, whereas ${\mathit{R}}_{\mathit{c}\mathit{o}\mathit{v}}$ is the observation error covariance matrix, which describes the uncertainty associated with the observational data. Additionally, $H$is the observation operator responsible for mapping the aerosol mass mixing ratio to the satellite-observed AOD in this study.

To construct the background error covariance matrix ${\mathit{B}}_{\mathit{c}\mathit{o}\mathit{v}}$, this study employed the GENerate Background Errors (GENBE) model, which utilizes the National Meteorology Center (NMC) method (Parrish and Derber, 1992). The background error covariance is described as follows:

where ${\boldsymbol{x}}^{\boldsymbol{b}}$ represents the forecast field, which is the state of the variable predicted in real time, while ${x}^{t}$ denotes the true field, reflecting the actual state of the aerosol. ${\boldsymbol{x}}^{\boldsymbol{T}1}$ and ${\boldsymbol{x}}^{\boldsymbol{T}2}$ represent two prediction fields at the same point in time but with different initialization times. In this study, T1 is set to 12 hours, and T2 is set to 24 hours, indicating that these fields correspond to the prediction results for 12 and 24 hours after the onset of the forecast, respectively.

Given that most observations are obtained independently, it is assumed that the observation error covariance matrix ${\mathit{R}}_{\mathit{c}\mathit{o}\mathit{v}}$ is a diagonal matrix. This assumption implies that the errors between different observations are independent of one another (Wang et al., 2022). On the basis of the satellite error parameters mentioned in Section 2.1, in this study, we adopted the observation errors for MODIS AOD data proposed by Remer et al. (2005), with an AOD uncertainty of 5% over the ocean and 15% over land. For the Himawari 8 AOD, we set the observation error at 6% (Yumimoto et al., 2016).

The GSI system uses the Community Radiative Transfer Model (CRTM) as the observation operator, which consists of a forward operator and a tangential/critical operator. The forward operator calculates the corresponding AOD value on the basis of the input aerosol mass concentration, whereas the tangential/critical operator is essential for the optimization process; it minimizes the objective function (Eq. 1) to achieve the best estimate of the atmospheric state. The official version of the GSI 3DVAR assimilation system can only assimilate AOD data retrieved from MODIS. We further extended the system to include 550 nm AOD retrieved from Himawari 8. In regions covered by both satellites, the respective contributions of the two satellites to the AOD analysis field are determined by the observation errors in those areas. All missing observations and those flagged as low quality by the satellite retrieval algorithms were excluded from assimilation.

2.3   Experiment Settings

In this study, the spatial and temporal distributions of aerosols were simulated using the WRF-Chem model version 4.3. The simulation domain, with a horizontal resolution of 25 km and 45 vertical layers, covering China and its neighboring regions (Fig. 1). Two time periods were modeled: Jun. 1^st, 2017 to Aug. 29^th, 2017, and Dec. 1^st, 2017, to Feb. 28^th, 2018, representing summer and winter, respectively. The lateral boundary conditions and initial meteorological fields for the WRF-Chem model were acquired from NCEP FNL reanalysis with a resolution of 2.5° × 2.5°. The flowchart of the assimilation experiment is presented in Fig. 2, where panel (a) outlines the overall process of multiple assimilation cycles of the chemical field, and panel (b) shows the workflow of a single assimilation cycle. A segmented integration method was applied, with the model starting at 18:00 UTC each day and running for 36 hours. The first 12 hours (from 18:00 UTC to 6:00 UTC the next day) were treated as the spin-up time and were excluded from the analysis, and the forecasts of the next 24 hours (from 6:00 UTC to 6:00 UTC the next day) were used for analysis. For first day simulation, the initial chemical field used the default chemical profile in WRF-Chem (Zhao et al., 2012). AOD was assimilated after 12 hours simulation. For the following simulations, only meteorological model was run during first 12 hours. After spin-up time, meteorological and chemical models were run simultaneously. The chemical initial field (06:00 UTC) was derived from the previous outputs: for the unassimilated experiments, we used the background field from the previous day’s 24 h chemical field forecast, whereas for the assimilated experiments, we used the assimilation analysis field from the previous day’s 24 h chemical field forecast. Anthropogenic emissions in China were sourced from the 2017 Multiresolution Emission Inventory for China (MEIC) inventory (Li et al., 2017a; Zheng et al., 2018), whereas emissions for regions outside China were obtained from the 2017 Model Intercomparison Study for Asia (MIX) inventory (Li et al., 2017b).

Download

Larger image

PowerPoint slide

Fig  2.  Flowchart of WRF-Chem model simulation and GSI 3DVAR assimilation process, including multiple assimilation cycles for the chemical field (a) and a single assimilation cycle (b).

The setting of spin-up time has a certain impact on simulation results. A few hours was enough for meteorological spin-up (Chu et al., 2018; Liu et al., 2023). For chemical field simulation, about a week was used as spin-up time with default chemical initial conditions, and grid nudging was used to decrease the accumulated error of meteorological simulation (Ritter et al., 2012; Kim et al., 2024). A simple evaluation also revealed that relatively long spin-up can improve AOD simulation (Fig. S1). Different from simulation, forecast cannot adopt grid nudging. As the forecast time extends, the uncertainty of the forecast will gradually increase. Prolonged spin-up time increases computational cost and may also lead to a decrease in forecast accuracy. In addition, different assimilation roadmaps also have a certain impact on the results. We evaluated three assimilation roadmaps, i.e., assimilating AOD (06:00 UTC) with previous forecast, today’s forecast using default chemical initial conditions, and today’s forecast using previous forecast as chemical initial conditions, and found that assimilating AOD with previous forecast had best performance (Fig. S2). These preliminary evaluations confirm that our model settings are effective.

The parameterization schemes utilized in the model are provided in Table 1. The model configuration includes the Thompson microphysics scheme for simulating cloud microphysical processes such as the formation and evolution of hydrometeors and the Grell‒Freitas cumulus parameterization to represent subgrid-scale convective processes, which are critical for capturing deep convection and precipitation. For boundary layer processes, the Mellor‒Yamada‒Nakanishi‒Niino (MYNN) scheme was used to model turbulence and vertical mixing, whereas land‒atmosphere interactions were simulated using the Noah land surface model, which accounts for soil moisture, surface energy fluxes, and vegetation effects. The revised MM5 surface layer scheme was employed to calculate surface fluxes of heat, moisture, and momentum. For aerosol processes, the Global Ozone Chemistry Aerosol Radiation and Transport (GOCART) scheme was adopted to simulate aerosol emissions, transport, and chemical transformations. Additionally, gas-phase chemical reactions were represented using the regional acid deposition model (RADM2) mechanism.

Table  1.  WRF-Chem parameterization scheme

Download as CSV

	Parameterization scheme	Refs
Microphysics	Thompson	(Thompson et al., 2008)
Cumulus Parameterization	Grell–Freitas	(Grell and Freitas, 2014)
Boundary Layer	MYNN	(Nakanishi and Niino, 2004, 2006)
Land Surface	Noah land surface model	(Chen et al., 1996)
Surface Layer	Revised MM5 surface layer scheme	(Jiménez et al., 2012)
Aerosol Scheme	GOCART	(Chin et al., 2000)
Gas Chemical	RADM2	(Stockwell et al., 1990)

 | Show Table

DownLoad: CSV

To compare the performance in terms of assimilating the AOD data derived from different satellites, four simulation experiments were conducted (Table 2). Although our study area spans six time zones (UTC+4 to UTC+9), MODIS data availability is limited. The Terra satellite (morning orbit) observes the region from approximately 01:30 UTC to 07:30 UTC, whereas the Aqua satellite (afternoon orbit) covers the area from approximately 04:30 UTC to 11:30 UTC. For most of the day, MODIS cannot provide simultaneous data from both satellites for assimilation. To ensure sufficient AOD observations for our study domain, we chose to assimilate at 06:00 UTC. At this time, the Terra satellite is nearing the end of its pass over the domain, and the Aqua satellite has just begun its coverage. Given the constrained coverage of individual observations from the polar-orbiting MODIS satellite, the AOD assimilation window is designated between 05:30 and 06:30 UTC to ensure the consistent and reliable generation of valid observation data. This time window allows for the inclusion of data from both MODIS satellites, maximizing observation availability for assimilation; however, since the data within the time window are not strictly aligned with the assimilation time points, this discrepancy could introduce temporal mismatch errors.

Table  2.  WRF-Chem Simulation setup

Download as CSV

Simulation	Assimilated data	Assimilation time window
CT	None	None
DAM	MODIS AOD	5:30-6:30 UTC
DAH	Himawari 8 AOD	6:00 UTC
DAC	MODIS AOD + Himawari 8 AOD	5:30-6:30 UTC + 6:00 UTC

 | Show Table

DownLoad: CSV

2.4   Statistical parameters

In this study, statistical metrics such as the Pearson correlation coefficient (R), root mean square error (RMSE), and mean bias (MB) were employed to evaluate the performance of the WRF-Chem model:

where $n$ is the total number of samples, ${\boldsymbol{S}}_{i}$ and ${\boldsymbol{O}}_{i}$ represent the simulated and observed AODs for the $i$^th group, respectively, and $\overline{\boldsymbol{S}}$ and $\overline{\boldsymbol{O}}$ are the mean values of the simulations and observations, respectively. R reflects the strength of the linear relationship between the simulation and observation, and a large R indicates a strong linear correlation. The RMSE value measures the overall difference between the simulated and observed values, with lower values indicating closer agreement and better simulation performance. The MB value quantifies the average error between the simulation and observation data and is used to assess the overall accuracy of the AOD simulations.

For spatial correlation coefficients, the method proposed by He et al. (2017) was adopted to calculate the effective degrees of freedom, ensuring the statistical significance of the results. The significance of the differences between the statistical indices derived from simulations and observations was tested by comparing their 95% confidence intervals. For R, the confidence intervals were calculated using the asymptotic normal distribution of $0.5\times \mathrm{log}\left(\dfrac{1+R}{1-R}\right)$. For other statistical indices, bootstrap tests were employed to determine the confidence intervals.

3.   Results and Discussions

3.1   Assessment of the Meteorological Field

The accuracy of meteorological field simulations is crucial to accurate chemical field simulations. In this study, the WRF-Chem simulated meteorological fields were assessed by automatic weather station observations (Table 3). The R values between the model-simulated T₂ and observed values were greater than 0.8, whereas for RH₂, they exceeded 0.6, and for WS₁₀, the coefficients were approximately 0.4. These results suggest that the WRF-Chem model effectively captured the temporal and spatial characteristics of the near-surface temperature and humidity fields in China. However, the model exhibited larger errors in simulating WS₁₀ than in simulating T₂ and RH₂. Previous studies have highlighted simulation uncertainties in wind speed, especially in cases of coarse grid resolution, where an inaccurate calculation of surface drag due to insufficient terrain depiction is a key source of error (Zhang et al., 2013; Liu et al., 2016; Mar et al., 2016; Song et al., 2017; Zhang et al., 2021; Feng et al., 2022). Research has shown that the RMSE values for T₂ simulations at similar or finer resolutions typically fall between 3.1 and 3.5 °C (Song et al., 2017), those for RH₂ are between 13.95% and 18.77% (Zhang et al., 2021), and those for WS₁₀ are between 2.1 and 2.3 m s⁻¹ (Feng et al., 2022). The simulation performance in this study was comparable to that of previous studies.

Table  3.  Performance for WRF-Chem simulations of meteorological fields^*

Download as CSV

	T2			RH2			WS10
	Summer	Winter		Summer	Winter		Summer	Winter
R	0.89	0.94		0.82	0.66		0.36	0.48
RMSE	2.87 °C	3.48 °C		14.24%	18.44%		2.49 m s−1	2.60 m s−1
MB	2.14 °C	2.64 °C		−1.70%	1.33%		1.92 m s−1	2.03 m s−1
*The correlation coefficients all passed the significance test with 95% confidence intervals.

 | Show Table

DownLoad: CSV

Overall, the WRF-Chem model performed better at simulating T₂ and RH₂, whereas the simulation error for WS₁₀ was relatively large. A previous study highlighted that although the model’s near-surface wind speed and direction simulations tend to produce relatively low correlations, they could still effectively drive the atmospheric chemistry model (Song et al., 2017). Despite the poor wind field simulation, the aerosol simulation generated by the model were considered reliable with a high level of confidence.

3.2   Comparison of MERRA-2 AOD reanalysis data

The distributions of the simulated monthly mean AOD data in the CT experiment and the three assimilation experiments, along with a comparison with the MERRA-2 monthly mean AOD, are shown in Fig. 3. In terms of spatial distribution, the AOD in central and eastern China is significantly greater in summer than in winter. This seasonal contrast may be attributed to several factors. First, a higher boundary layer height (BLH) during summer enables more efficient transport and mixing of aerosol particles to higher altitudes, reducing the near-surface concentration while increasing the columnar AOD. Second, the inshore monsoon circulation during this period leads to elevated relative humidity, which enhances aerosol extinction and contributes to a higher AOD. Finally, regional open stalk burning during the summer harvest season, as observed in MODIS fire products, introduces significant aerosol emissions, further amplifying AOD levels (Qu et al., 2016). In contrast, the AOD in the Sichuan Basin is greater in winter. The decreasing AOD trends in summer in the Sichuan Basin are associated with frequent rainfall because rainwater can remove the majority of ambient respirable suspended particulates and act as an effective way to shorten their lifetime. Another reason for the lower AOD values in summer is the weather patterns of the summer monsoon, under which these regions are seldom affected by dust from large-scale transport originating in the northern part of China. The winter increase in the AOD in the Sichuan Basin is likely associated with poor dispersion conditions and heavy local industrial emissions (Li et al., 2003). In the control experiment (CT), higher simulated AOD values are concentrated over the Tibetan Plateau and central and eastern China, a distribution that aligns well with the MERRA-2 outputs but is significantly lower overall. Several factors likely contributed to this underestimation of AOD in the CT experiment. Biases in the simulated meteorological fields, such as overestimation of atmospheric dispersion due to higher-than-expected wind speeds (as shown in Table 3), could have contributed to the lower AOD results. The exclusion of natural source emissions during the simulation might have further contributed to the simulated AOD values. The spatial correlation coefficients between the CT experiment and MERRA-2 monthly mean AOD were 0.44 in summer and 0.65 in winter (Table 4). After assimilating the AOD from Himawari 8 or MODIS, the spatial distribution of the WRF-Chem simulations more consistently aligned with the spatial distribution of the MERRA-2-output AOD. The Himawari 8-derived AOD was less abundant in the western part of the simulation domain (Fig. 1c, 1f), limiting the improvement in that region in the DA_H experiment. In contrast, in the eastern part of the simulation domain, where Himawari 8 AOD data were more abundant, a significant improvement in the AOD simulation was observed. In regions covered by both satellites, such as central and eastern China and the southern Tibetan Plateau, the DA_H simulations produced lower AOD values than the DA_M simulations did, likely due to differences in the observational instruments and retrieval algorithms used by the satellites. By combining data from both satellites, most of the study area was effectively covered. Relative to MERRA-2, the DA_C experiment achieved spatial correlation coefficients of 0.57 in summer and 0.70 in winter (Table 4). In summer, the DA_C experiment significantly improved the spatial distribution of AOD, whereas in winter, despite an improvement in correlation, the results did not pass the significance test. Overall, WRF-Chem simulated the spatial distribution of AOD more accurately in winter, whereas the improvement in assimilation was more pronounced in summer.

Download

Larger image

PowerPoint slide

Fig  3.  Monthly mean AOD distributions simulated from four different numerical experiments for the summer (a-e) of 2017 and winter (f-j) of 2017/2018 and comparison with the MERRA-2 monthly mean AOD.

Table  4.  Spatial correlation coefficients between the WRF-Chem-simulated monthly mean AOD and the MERRA-2 monthly mean AOD^*

Download as CSV

	CT	DAH	DAM	DAC
Summer 2017	0.44	0.56	0.60	0.57
Winter 2017/2018	0.65	0.69	0.72	0.70
*The correlation coefficients all passed the significance test with 95% confidence intervals

 | Show Table

DownLoad: CSV

Figs. 4 and 5 present probability density scatter plots comparing the hourly WRF-Chem simulated AOD with the MERRA-2 reanalysis AOD for the summer of 2017 and winter of 2017/2018. In the plots, most of the data are concentrated near the origin (coordinate 0), indicating a high density, whereas the scatter points in the peripheral areas are much sparser. Moving outward from the origin, the probability density decreases sharply, suggesting that in most cases, the reanalysis and simulation results of the AOD within the study domain remain relatively low. However, in reality, high AOD values are often associated with severe pollution, making accurate forecasting of high-AOD events particularly critical (Zheng et al., 2017). The CT experiment showed a significant deviation between the simulated AOD and MERRA-2 reanalysis data, especially in summer (Fig. 4a), where most simulated AOD values were noticeably lower. The MB values in the CT experiment were −0.17 for summer and −0.12 for winter, indicating underestimation in both seasons. In contrast, the DA_H and DA_M experiments demonstrated significant improvements in the AOD simulation. These experiments yielded more reasonable probability density distributions than did the CT experiment. In summer, the MB value for the DA_M experiment was −0.082, with an RMSE value of 0.223 (Fig. 4c), whereas in winter, the MB value was -0.055, with an RMSE value of 0.184 (Fig. 5c). The DA_H experiment yielded slightly larger errors, with an MB value of −0.112 and an RMSE value of 0.240 in summer (Fig. 4b) and an MB value of −0.092 and an RMSE value of 0.194 in winter (Fig. 5b). On the basis of these statistics, the DA_M experiment's simulated AOD aligned more closely with the MERRA-2 reanalysis data than did that of the DA_H experiment. One possible explanation for the better performance of the DA_M experiment is that MERRA-2 assimilated both MODIS and AVHRR AOD observations, creating some level of homology with the DA_M experiment, which also incorporated MODIS data. This similarity could explain why the DA_M experiment more accurately captured the MERRA-2-derived AOD than did the DA_H experiment, which used Himawari 8 data. When AOD data from both satellites were assimilated, the WRF-Chem model showed pronounced improvements over single-satellite assimilation (Figs. 4d, 5d and Table 4). As shown in Table 4, regardless of whether assimilation was applied, WRF-Chem tended to simulate AOD more accurately in winter than in summer. However, the improvements induced by assimilation are more pronounced during the summer months. This phenomenon may result from larger simulation errors in summer meteorological fields (e.g., T₂ and WS₁₀) and precipitation biases. Specifically, frequent summer precipitation events coupled with the model's pronounced deficiencies in simulating rainfall patterns likely contribute to amplified errors in aerosol optical depth (AOD) simulations during this season (Gao et al., 2014). Moreover, the greater initial uncertainties in summer meteorological fields reduce the effectiveness of data assimilation techniques compared with those in winter conditions, when more stable atmospheric states allow assimilation methods to achieve better error correction.

Download

Larger image

PowerPoint slide

Fig  4.  Scatter plot of probability densities for CT (a), DA_H (b), DA_M (c), and DA_C (d) versus MERRA-2 reanalyzed AOD, based on hourly forecast data for summer 2017.

Download

Larger image

PowerPoint slide

Fig  5.  Scatter plot of probability densities for CT (a), DA_H (b), DA_M (c), and DA_C (d) versus MERRA-2 reanalyzed AOD, based on hourly forecast data for winter 2017–2018.

Central and eastern China and the Sichuan Basin, two areas known for heavy pollution (as shown in Fig. 1 and their locations marked in yellow and blue boxes), and the Tibetan Plateau, a representative clean region (as shown in Fig. 1 and the location marked in the red box), were selected for further evaluation of the performance of the WRF-Chem simulated AOD (Fig. 6). For central and eastern China and the Sichuan Basin, although the CT experiment simulated AOD values were significantly lower than those from MERRA-2, it was still able to capture the day-to-day variations in AOD, with R values of 0.53 for summer and 0.82 for winter. This finding showed that even without assimilation, the model could represent the temporal variation in AOD reasonably well. The assimilation experiments significantly improved the AOD simulations in these heavily polluted areas. However, as shown in Fig. 6, the daily average AOD fluctuations in the MERRA-2 experiment are significant, whereas the CT experiment, although it partially reflects the trend of AOD variation, has a relatively smoother rate of change. A comparison of the results across the experiments revealed that the CT experiment yielded notably lower AOD than MERRA-2 did, whereas the AOD levels and variations in the assimilated experiments were comparable to those in the latter. In summer (winter), the RMSE value of the DA_C experiment's AOD simulation was reduced by 42% (55%) in the central and eastern China and 37% (42%) in the Sichuan Basin compared to the CT experiment. This finding indicates that assimilating satellite-derived AOD greatly enhanced the model’s ability to simulate aerosol concentrations in these regions. The sample size of the assimilated data was also a key factor affecting the simulation performance. In summer, the assimilation rate of the MODIS-derived AOD was higher than that of Himawari 8, which contributed to the DA_M experiment simulating AOD values closer to those of MERRA-2 than those of the DA_H experiment. In central and eastern China, where both Himawari 8 and MODIS AOD data were well covered, the DA_C experiment-produced AOD data were more closely aligned with the MERRA-2 data. However, in the Sichuan Basin, the situation was different. The region's complex topography led to low wind speeds, weak horizontal diffusion, and high humidity, which caused aerosols to accumulate within the basin. These factors resulted in persistent cloud cover throughout the year, limiting the availability of AOD observations. As a result, the DA_C experiment's AOD simulations for the Sichuan Basin showed larger deviations from the MERRA-2 reanalysis than those for central and eastern China. This finding demonstrates the challenge of accurately simulating AOD in areas with complex meteorological and geographical conditions.

Download

Larger image

PowerPoint slide

Fig  6.  Comparison of MERRA-2 and WRF-Chem 24-hour forecasts of daily mean AOD in Central and Eastern China (a-b), the Sichuan Basin (c-d) and the Tibetan Plateau (e-f).

For the Tibetan Plateau, as shown in Fig. 3e and 3f, the AOD values were significantly lower than those in other regions. In this area, the correlation coefficients between the CT experiment and MERRA-2 were 0.57 and 0.60 for summer and winter, respectively, indicating that the CT experiment reasonably captured the spatial distribution and temporal variation in AOD even without assimilation. However, the CT experiment tended to underestimate AOD levels compared with those of MERRA-2. After assimilation, the RMSEs of AOD in the DA_C experiment significantly improved by 53% in summer and 58% in winter compared with those in the CT experiment, demonstrating that data assimilation substantially enhanced AOD simulations over the Tibetan Plateau. Notably, this improvement in winter was more pronounced than that in the two heavily polluted regions studied. This finding may be attributed to the complex topography of the Tibetan Plateau, which introduces greater simulation uncertainties (Zhao et al., 2017). Although satellite observations in this region are relatively sparse compared with those in other areas, they are sufficient to reflect the AOD distribution characteristics.

3.3   Comparison of AOD observations with AERONET sites

Fig. 7 shows the scatter plot of the WRF-Chem model simulated AOD versus the station observed AOD. The CT experiment effectively captured the variation characteristics of the AOD, with an R value of 0.57 (Fig. 7a); however, it exhibited significant simulation errors and a general underestimation of the AOD values. In comparison, the DA_M experiment demonstrated smaller MB and RMSE values along with a higher R than the DA_H experiment did, indicating that the overall assimilation of the MODIS AOD led to better model performance than the assimilation of the Himawari 8 AOD. Furthermore, the DA_C experiment showed improved AOD simulations with smaller errors and a higher R than single satellite assimilation did, suggesting that assimilating both the MODIS and Himawari 8 AODs can enhance model simulations.

Download

Larger image

PowerPoint slide

Fig  7.  Scatter plot of the 24-hour forecast AOD versus the observed AOD for CT (a), DA_H (b), DA_M (c), and DA_C (d).

Table 5 presents a comparison of the WRF-Chem simulated AOD with the AERONET site-observed AOD. Consistent with the results from the MERRA-2 reanalysis, the CT experiment simulated lower AOD values than the site observations did, with MB values ranging from -0.64 to -0.05. Assimilating satellite-derived AOD significantly mitigated the underestimation in simulating AOD. After both Himawari 8 and MODIS AODs were incorporated, the MB for the DA_C simulated AOD ranged from -0.19 to -0.01. In the eastern part of the study area, where there is a higher proportion of Himawari 8 data coverage, the assimilation of Himawari 8 AOD (DA_H experiment) resulted in a smaller RMSE and a higher R value. For example, the DA_H experiment performed better at the BJ site in both summer and winter and at the YS site in summer. Conversely, in the western region, where only MODIS provided available AOD data, the DA_M experiment yielded better AOD simulations. At the KP site, the assimilation of Himawari 8 data (DA_H) had less improvement for AOD simulations in both summer and winter, whereas assimilating MODIS AOD (DA_M) significantly enhanced AOD simulation in this region. The DA_C experiment, which integrated data from both satellites, yielded simulation results comparable to those obtained from single satellite assimilation when only one satellite provided observational data, as observed at the KP site. Nonetheless, at specific sites, such as the YS site in summer and the BJ site in winter, the DA_C simulation outperformed the individual satellite assimilations. Significance tests indicated that the DA_C significantly improved the AOD simulation across all sites (except for the BT site in summer and the YS site in winter) compared with the CT simulation of the AOD (Table 6).

Table  5.  Statistical performance of the simulated AOD for the summer of 2017 and the winter of 2017/2018^*

Download as CSV

	CT				DAM				DAH				DAC
	MB	RMSE	R		MB	RMSE	R		MB	RMSE	R		MB	RMSE	R
Summer
BJ	−0.40	0.61	0.62		−0.25	0.53	0.41		−0.19	0.42	0.69		−0.19	0.43	0.67
BT	−0.08	0.10	0.71		0.04	0.10	0.59		−0.05	0.09	0.65		−0.02	0.09	0.59
XZ	−0.54	0.72	0.44		−0.19	0.54	0.53		−0.19	0.41	0.69		−0.09	0.42	0.66
KP	−0.19	0.4	0.15		0.15	0.33	0.78		−0.45	0.54	0.13		0.06	0.29	0.68
YS	−0.39	0.55	0.54		−0.20	0.43	0.43		−0.09	0.30	0.73		−0.10	0.30	0.74
Winter
BJ	−0.11	0.21	0.70		−0.10	0.20	0.68		−0.06	0.15	0.82		−0.05	0.15	0.83
BT	−0.05	0.08	0.69		−0.05	0.08	0.67		−0.03	0.06	0.77		−0.02	0.06	0.76
XZ	−0.25	0.39	0.47		−0.18	0.34	0.53		−0.08	0.28	0.64		−0.08	0.28	0.63
KP	−0.64	0.73	0.21		0.05	0.39	0.52		−0.5	0.61	0.29		−0.08	0.41	0.59
YS	−0.08	0.15	0.72		−0.05	0.14	0.68		−0.01	0.14	0.67		−0.01	0.14	0.67
*The correlation coefficients all passed the significance test with 95% confidence intervals

 | Show Table

DownLoad: CSV

Table  6.  Significance tests for differences in the statistical indices of the simulated AODs in the summer of 2017 and winter of 2017/2018^*

Download as CSV

DAC vs CT
Time	Indicator	BJ	BT	XZ	KP	YS
Summer	RMSE	Y	N	Y	Y	Y
	MB	Y	Y	Y	Y	Y
	R	N	n	Y	Y	Y
Sinter	RMSE	Y	Y	Y	Y	N
	MB	Y	Y	Y	Y	Y
	R	Y	Y	Y	Y	n
*Y (y) and N (n) indicate significant and nonsignificant differences in the statistical indices, respectively. Upper (lower) case letters denote better (worse) simulation statistics.

 | Show Table

DownLoad: CSV

To evaluate the impact of satellite data availability on assimilation performance, we analyzed the spatial coverage of Himawari-8 observations within a 2° radius centered at the BJ site under the DA_H framework. Seasonal thresholds were established to categorize coverage rates: during summer, data coverage was classified as low (<5%), moderate (5–15%), and high (>15%), whereas winter thresholds were adjusted to <10%, 10–20%, and >20%, respectively, to account for reduced satellite visibility under frequent cloud cover. As quantified in Table 7, the 24-hour correlation coefficients between the DA_H and independent AERONET measurements exhibited a progressive increase with increasing data coverage. The summer correlations improved from 0.62 under low-coverage conditions to 0.73 at high coverage levels, whereas the winter correlations exhibited a more pronounced increase from 0.79 to 0.88.

Table  7.  Statistics of the correlation coefficients between the DA_H-simulated AOD and observations in BJ under different assimilation data ratios^*

Download as CSV

	Low	Moderate	High
Summer 2017	0.62	0.66	0.73
Winter 2017/2018	0.79	0.81	0.88
*The correlation coefficients all passed the significance test with 95% confidence intervals. Low, moderate, and high assimilation data ratios represent less than 5% (10%), 5%-15% (10%-20%), and more than 15% (20%) of the data, respectively, in summer (winter).

 | Show Table

DownLoad: CSV

4.   Conclusion and Prospect

In this study, we employed the WRF-Chem model to simulate the AOD in China and its surrounding areas during the summer of 2017 and winter of 2017/2018. We assimilated Himawari 8 and MODIS AOD data using the 3DVAR method within the GSI assimilation system and compared the simulation results with the AOD data acquired from both the MERRA-2 reanalysis and the AERONET sites. We analyzed the impacts of single satellite assimilation and simultaneous assimilation of satellite AOD data on the WRF-Chem model's simulated AOD.

The results indicated that the WRF-Chem model effectively simulated meteorological fields for aerosol modeling studies but significantly underestimated the AOD. This underestimation was attributed to the biases in the meteorological field, and the exclusion of natural source emissions. After assimilating the AOD data from Himawari 8 and MODIS, the model improved the spatial and temporal variability characteristics of the AOD. However, significant differences in the accuracy of the assimilated simulated AOD were observed between the eastern and western regions of the study area owing to disparities in the spatial and temporal scales of the valid observed data from the two satellites. Significance tests revealed that the simultaneous assimilation of both the MODIS and Himawari 8 AOD data resulted in better WRF-Chem simulated AOD data than did single satellite assimilation, with improvements being particularly pronounced in summer. An analysis of key polluted areas revealed that dual satellite assimilation in the middle and eastern regions of China yielded better improvements than the improvements observed in the Sichuan Basin. The corresponding improvement in the lightly polluted area of the Tibetan Plateau was greater than that in the above two heavily polluted areas in winter. An analysis of data availability and assimilation quality near the BJ site revealed that greater data accessibility significantly enhanced assimilation performance. Overall, the simultaneous assimilation of MODIS and Himawari 8 AOD data can effectively address the limitations of single-source information and enhance model simulations. However, substantial improvements in assimilation effects are context dependent.

This study focused primarily on the differences between assimilating MODIS or Himawari 8 AOD data and their simultaneous assimilations in WRF-Chem AOD simulations. Nevertheless, there remains significant potential to enhance the assimilation effectiveness of WRF-Chem, given the considerable uncertainty in meteorological conditions that affects aerosol simulations, as well as the unclear impact of aerosol chemical composition on these simulations. In future work, we aim to improve the accuracy of meteorological field simulations through nudging and other methods, focus on the effects of aerosol and chemical composition simulations to better understand the impact of assimilation on aerosols, and further enhance the simulation of AOD. Furthermore, by leveraging the enhanced data assimilation potential from improved satellite observation coverage, we propose employing deep learning methodologies for optimizing the spatial interpolation of satellite-derived AOD products.

Download

Larger image

PowerPoint slide

  S1.  Scatter plot of WRF-Chem model simulated AOD with default chemical initial conditions for 12 hours’ spin up (a) and 24 hours’ spin up (b) versus observed AOD at the AERONET site.

Download

Larger image

PowerPoint slide

  S2.  Scatter plots of DAH vs. AERONET-observed AOD: assimilating AOD (06:00 UTC) with previous forecast (a), today’s forecast using default chemical initial conditions (b), and today’s forecast using previous forecast as chemical initial conditions (c).