The meteorological fields were simulated using the WRF (version 4.2.1) model with the FNL reanalysis dataset. The FNL data were obtained from the U.S. National Centre for Atmospheric Research (NCAR), with a spatial resolution of 1.0° × 1.0° (http://rda.ucar.edu/datasets/ds083.2/, last accessed on 15 November 2021). The physical parameterizations used in this study include the Thompson microphysical process, RRTMG longwave/shortwave radiation scheme; Noah land-surface scheme; MYJ boundary layer scheme; and modified Tiedtke cumulus parameterization scheme. The detailed configuration settings could be found in the work of Hu et al. (2016) and Wang et al. (2021).

The CMAQ version 5.2 (CMAQv5.2) model (Fahey et al., 2017), configured with the gas-phase mechanism of SAPRC07tic and the aerosol module of AERO6i, was employed in this study to simulate the air quality over China during 2013–2019. Air quality simulations were performed for a period of 7 years (2013–2019) using a horizontal resolution of 36 km. The corresponding domain covered China and the surrounding countries and regions with 197 × 127 grids (Figure 1A). The vertical resolution had 18 layers. The initial and boundary conditions were provided by the default profiles of the CMAQ model. The simulated results of the first three days were not included in the model analysis, thus serving as a spin-up and reducing the effects of the initial conditions on the simulated results.

The Multi-resolution Emission Inventory for China version 1.3 (MEICv1.3) (http://www.meicmodel.org) and Regional Emission inventory in ASia (REASv3.2) (https://www.nies.go.jp/REAS/) were used to provide the anthropogenic emissions. MEIC served as the anthropogenic emissions from China, while REAS served as the anthropogenic emissions from neighboring countries and regions. The MEICv1.3 emissions of the year 2017 were used for the simulations of the years 2018 and 2019, as no reliable sources for emission changes in China were currently available. For REAS, since no emission inventory was released for the years after 2015, we used the emission inventory in the year 2015 for 2016–2019. Although emission inventories are usually released 2–3 years behind, we acknowledge that this may cause additional uncertainties in the simulation for 2018 and 2019. Biogenic emissions were generated using the Model for Emissions of Gases and Aerosols from Nature (MEGANv2.1) (Guenther et al., 2012) for the whole simulation period. The open biomass burning emissions were processed using the Fire Inventory for NCAR (FINN) during the entire study period (Wiedinmyer et al., 2011). The spatiotemporal variations of the total emission of PM_2.5, SO₂, NOx, NH₃, and VOC across the five regions are shown in Supplementary Table S3.

The daily observation data of meteorological variables (wind speed, wind direction, relative humidity, and temperature) for the selected regions were downloaded from the Chinese Meteorological Agency (http://data.cma.cn/en, last accessed on 30 November 2021). There were 18, 11, 14, 14, and 11 meteorological stations considered in the NCP, YRD, PRD, CY, and FW regions, respectively. In addition, the hourly observation data of air pollutants (PM_2.5, O₃, NO₂, and SO₂) were obtained from the Chinese Ministry of Ecology and Environment (https://www.mee.gov.cn/, last accessed on 20 December 2021). In this study, five regions in China were selected as target areas (Figure 1B), and they include the North China Plain (NCP, with 70 air quality monitoring stations), Yangtze River Delta (YRD, with 107 air quality monitoring stations), Pearl River Delta (PRD, with 54 air quality monitoring stations), Chengyu (CY, with 40 air quality monitoring stations), and Fenwei Plain (FW, with 60 air quality monitoring stations) regions. Besides the NCP and YRD, most of the cities in other regions had a month observation data in 2013. The cities with data of more than 3 months in each region were selected for citywide analysis in 2013, while all of the data in each region were selected for regional analysis. In the subsequent years (2014–2019), the observation data of the monitoring stations in each of the cities were used citywide, while the observation data of all the cities located in each region were used to estimate the average observation value of the region (Hu et al., 2016; Sulaymon et al., 2021c). The 21 cities selected in 2013 in the five regions are as follows: Beijing, Tianjin, Shijiazhuang, Qinhuangdao, Chengde, and Zhangjiakou in the NCP; Shanghai, Wuxi, Nanjing, Suzhou, Xuzhou, and Hangzhou in the YRD; Guangzhou, Dongguan, and Shenzhen in the PRD; Chengdu, Mianyang, and Chongqing in the CY; and Xi’an, Xianyang, and Baoji in the FW. The details about the selected cities in each region are shown in Supplementary Table S1. Furthermore, the predicted major chemical components of PM_2.5 (SO₄ ^2-, NO₃ ⁻, and NH₄ ⁺) were evaluated using the daily observation data in nine cities (Beijing, Shijiazhuang, Nanjing, Suzhou, Xuzhou, Hangzhou, Guangzhou, Shenzhen, and Chengdu) during the study period.

Previous studies have investigated and documented the impacts of meteorological conditions on the formation, transportation, and dissipation of air pollutants (Hu et al., 2016; Hua et al., 2021; Sulaymon et al., 2021c, 2021d). In addition, the influences of some meteorological parameters (such as wind speed, wind direction, temperature, and relative humidity) on air quality modeling have been elucidated (Hu et al., 2016; Sulaymon et al., 2021b; Wang et al., 2021). Therefore, the evaluation of the WRF model performance was carried out prior to the usage of its meteorological fields in the air quality simulations. The evaluation of the WRF model was achieved by comparing the predicted wind speed (WS) and wind direction (WD) at 10 m above the surface, as well as the simulated relative humidity (RH) and temperature (T2) at 2 m above the ground level to their corresponding observed values in each region during the entire study period. The statistical indices used in evaluating the WRF model were the mean bias (MB), mean error (ME), and root mean square error (RMSE) (Table 1). The benchmarks of statistical indices employed in this study were suggested by Emery et al. (2017). T2 was generally over-predicted in all the regions except the NCP, whose MB value fell below the suggested benchmark (≤±0.5), while the ME values in all the regions, except the CY (the southeast basin of the Tibetan Plateau, with poor terrain and complicated weather conditions), were found below the suggested benchmark (≤2.0). With low ME indices (≤2.0) in four out of the five regions, it is shown that T2 was well simulated in the four regions. Previous studies have reported over-prediction of T2 in the YRD (Ma et al., 2021; Sulaymon et al., 2021a; Wang et al., 2021) and PRD (Wang N. et al., 2016). The results of this study are consistent with previously reported ones in the studied regions. However, no benchmarks were suggested for the MB and ME values of RH, and RH was underestimated in the NCP and FW, while it was overestimated in the other three regions. The MB values of WS in all the regions except the CY greatly exceeded the recommended criterion (≤±0.5), while the ME values in all the five regions were below the benchmark (≤2.0). In addition, the RMSE values of WS met the benchmark (≤2.0) in all the regions except the PRD. Considering the ME and RMSE values, the simulated WS reasonably captured the observations in all the regions. Over-prediction of WS has been previously found in the PRD (Wang N. et al., 2016; Qing Chen et al., 2019), NCP (Hanyu Zhang et al., 2019; Sulaymon et al., 2021b; Mengmeng Li et al., 2021), and YRD (Sulaymon et al., 2021a; Ma et al., 2021; Wang et al., 2021; Yu et al., 2021). Except in FW, the MB indices of WD in other regions were greater than the suggested benchmark (≤±10), while the ME values in all the regions greatly exceeded the recommended criterion (≤±30), especially in the PRD (68.98), YRD (48.83), CY (47.63), and FW (45.84). The model performance of WD in this study was consistent with previous studies in the YRD (Sulaymon et al., 2021a; Wang et al., 2021; Yu et al., 2021) and NCP (Sulaymon et al., 2021b) regions. Generally, the WRF model in this study performed better when compared to previous studies across China (Hu et al., 2016, 2017; Wang H. L. et al., 2016; Hanyu Zhang et al., 2019; Qing Chen et al., 2019; Sulaymon et al., 2021a, 2021b; Ma et al., 2021; Mengmeng Li et al., 2021; Yu et al., 2021), and the simulated meteorological fields were further utilized in driving the CMAQ model.

Supplementary Figures S1-S5 show the comparison of the simulated daily mean concentrations of O₃, NO₂, SO₂, and PM_2.5 with the observations in 21 cities during 2013–2019. The time series of pollutants’ concentrations in six major cities in the NCP region are illustrated in Supplementary Figure S1. In the NCP, the observed O₃ was about 200 μg/m³ in summer and 50 μg/m³ in winter; NO₂ was 100 μg/m³ in winter and 20 μg/m³ in summer; SO₂ was 100 μg/m³ in winter and 10 μg/m³ in summer before 2016 and 10 μg/m³ without seasonal change after 2016; and PM_2.5 was 200 μg/m³ in winter and 10 μg/m³ in summer. The simulated O₃ and PM_2.5 were captured well with the observation data; SO₂ was well simulated on monthly trends at Beijing, Tianjin, and Shijiazhuang sites than in other cities such as Qinhuangdao, Chengde, and Zhangjiakou, where it was high in winter and low in summer before 2016. NO₂ was underestimated by over 20 μg/m³, without obvious seasonal change after 2016. Supplementary Figure S2 shows the six megacities in the YRD region. The monthly trend and simulated results of O₃, NO_2, and PM_2.5 were the same as in the NCP; the predicted SO₂ was better in Shanghai, Wuxi, and Hangzhou compared to Nanjing, Suzhou, and Xuzhou, which was near the observed data within 10 μg/m³. Supplementary Figure S3 shows the three cities in the PRD region. The monthly trends of O₃ and PM_2.5 were the same as those of the NCP and YRD. SO₂ was well simulated in Guangzhou and Dongguan sites relative to Shenzhen, which was overestimated by more than 20 μg/m³; NO₂ was well simulated in Shenzhen with a little seasonal change. Supplementary Figure S4 shows the three urban stations in the CY region. The simulated O₃ was close to the observation data; NO₂ was well simulated in Chengdu without seasonal changes; PM_2.5 and SO₂ were generally overestimated. FW (Supplementary Figure S5) was similar to CY on monthly trend, and both O₃ and PM_2.5 were well predicted, while NO₂ and SO₂ were underestimated.

Figure 2 shows the comparison of simulated pollutants with observed data in the five regions. The observed data of all stations and the simulated results in the five regions were grouped into three (2013–2014, 2015–2016, and 2017–2019) for model evaluation. The statistical metrics of normalized mean bias (NMB), normalized mean error (NME), and the correlation coefficient (R) were calculated to evaluate the pollutant predictions in each region (Table 2). The model performance criteria for O₃ and PM_2.5 were suggested by Emery et al. (2017). In the NCP, O₃ was slightly overpredicted during 2013–2014 (NMB: 0.19), while it was well predicted during the 2015–2016 and 2017–2019 periods (NMB less than the benchmark), and the model performance improved with an increase in years as observed in R-values. The model performance in predicting PM_2.5 and SO₂ also improved significantly with an increase in years with higher R-values during 2015–2016 and 2017–2019 compared to the 2013–2014 period and with the statistical metrics of PM_2.5 meeting the suggested benchmarks (Table 2). NO₂ was underestimated, but the model performance also improved with increased R-value as the years increased. Similar to the NCP, O₃ was slightly overpredicted in the YRD during 2013–2014 (NMB: 0.16), while it was well predicted during the 2015–2016 (NMB: 0.02) and 2017–2019 (NMB: 0.02) periods. The model performance improved with higher R-values with an increase in years. In terms of NMB and R-values, the model performance of SO₂ in the YRD decreased with an increase in years. PM_2.5 was well estimated with the NMB and R-values of 0.22–0.30 and 0.83–0.87, respectively, while NO₂ was underestimated with fluctuating R-values, a similar scenario to what was found in the NCP. The trio of O₃, SO₂, and PM_2.5 was well predicted in the PRD region with their NMB values found below the suggested benchmarks during the grouped study periods, except for O₃ with NMB value slightly higher than the criterion during 2015–2016. The model performance of O₃, SO₂, and PM_2.5 significantly improved with higher R-values as the years increased. NO₂ was underestimated during the three periods, and the bias increased with an increase in years. It should be noted that the four pollutants (O₃, NO₂, SO₂, and PM_2.5) had their highest R-values during the 2015–2016 period. In the CY region, O₃ was slightly overpredicted during 2013–2014 (NMB: 0.21), while it was well predicted during the 2015–2016 and 2017–2019 (NMB: 0.09) periods. The model performance of O₃ improved with higher R-values with an increase in years. NO₂ was underestimated while SO₂ and PM_2.5 were highly overestimated, with the NMB and NME values of PM_2.5 greatly exceeding the suggested criteria during the three periods. The R-values of PM_2.5 also increased with the increase in years, while fluctuation was noted in the R-values of SO₂. PM_2.5 in the FW region was well predicted with very low NMB values. The R-value decreased as the years increased, an indication that the best model performance occurred during the 2013–2014 period (Table 2). NO₂ was underestimated throughout the three periods, while O₃ was overestimated during the 2013–2014 and 2015–2016 periods. The overall model performance improved as the years increased. Above all, the model exhibited better performances in reproducing O₃, SO₂, and PM_2.5 in the NCP, YRD, and PRD regions.

Figure 3 shows the results of the model performance of SO₄ ^2-, NO₃ ⁻, and NH₄ ⁺ in the nine selected cities across the five regions. In Beijing and Shijiazhuang, SO₄ ^2- had an NMB value of ∼ -0.2, an indication of underprediction, while NO₃ ⁻ and NH₄ ⁺ were well simulated. The NMB value of NO₃ ⁻ (0.15) in Shijiazhuang was relatively the same as that of PM_2.5 in the NCP (Figure 2), indicating that NO₃ ⁻ was more dominant in the NCP. In Nanjing, Suzhou, Xuzhou, and Hangzhou (YRD), SO₄ ^2- was well simulated, while NO₃ ⁻ and NH₄ ⁺ were overestimated. The NMB of PM_2.5 in the YRD (Figure 2) was similar to the NMB of SO₄ ^2- in the selected four cities in the YRD (Figure 3), suggesting that SO₄ ^2- significantly dominated PM_2.5 in the YRD compared to NO₃ ⁻ and NH₄ ⁺. NO₃ ⁻ and NH₄ ⁺ in Guangzhou and Shenzhen were also well simulated, while SO₄ ^2- was underestimated. All of the three major components dominated PM_2.5 in the PRD. In Chengdu, NO₃ ⁻ and NH₄ ⁺ were overestimated with NMB values of 2.55 and 0.87, respectively. This feature was also observed in PM_2.5 in the CY (overprediction) (Figure 2), and this shows that NO₃ ⁻ and NH₄ ⁺ were the important components of PM_2.5 in the CY.

The analysis of the spatial and temporal variations shows that effective emissions control strategies are needed in eastern China. The control of PM_2.5 is always complicated as it is related to its major components. For instance, it was found that NO₃ ⁻ was significant in the NCP and SO₄ ^2- was dominant in the YRD, while all of the three components equally influenced PM_2.5 in the PRD region. In addition, NO₃ ⁻ and NH₄ ⁺ were more dominant than SO₄ ^2- in the CY. From the abovementioned analysis, it can be concluded that NO₃ ⁻ and NH₄ ⁺ were the main components of PM_2.5 in heavily polluted regions, where the maximum PM_2.5 concentration was above 250 μg/m³. However, in clean regions such as the PRD (with high temperature and precipitation), where the maximum PM_2.5 concentration was below 150 μg/m³, each component equally influenced PM_2.5.

The model performance during different pollution levels was also evaluated. The pollution levels used in this study were based on ambient air quality standards of the Chinese Ministry of Ecology and Environment. Table 3, Supplementary Table S2, and Figure 4 show the model performance during different pollution levels in different regions. O₃ was overestimated during the “Good” level and underestimated during the “Lightly Polluted” and “Moderately Polluted” levels in all regions (Table 3). PM_2.5 was overestimated in the NCP, YRD, and CY. NCP was overpredicted (<10 μg/m³) during low-pollution and underestimated (<10 μg/m³) during high-pollution events. In the YRD and CY, PM_2.5 was overestimated (>20 μg/m³). In the CY, PM_2.5 was highly overpredicted as all of the statistical metrics breached the suggested standards for all the pollution levels. This could be attributed to the poor terrain and complicated weather conditions in the CY region. In the PRD and FW, PM_2.5 was generally underpredicted (>20 μg/m³), especially during the four pollution stages. As shown in Supplementary Table S2, NO₂ was underestimated in all the five regions (>10 μg/m³) with different pollution levels except in the CY. SO₂ in the NCP, YRD, and PRD regions was underestimated (>5 μg/m³), while it was overestimated in the CY and FW (>10 μg/m³).

Supplementary Figure S6 illustrates the monthly trends and the highest and lowest concentrations of the pollutants (O₃, NO₂, SO₂, and PM_2.5), while Figure 5 shows the annual trends of the pollutants in the five regions during 2013–2019. In addition, Figure 6 shows the total trends of the pollutants, which were obtained as the difference between 2019 and 2013 annual average concentrations in each region. Considering the model performance of the pollutants, NO₂ was generally underpredicted on a monthly basis in all the regions during the whole study period except in the CY, while other pollutants exhibited better monthly predictions as illustrated in Supplementary Figure S6. It should be noted that SO₂ was also overpredicted in the CY region.

In the NCP, the changes in the four pollutants during the study period followed the same pattern (Figure 5). Overall, the observed and simulated SO₂ and PM_2.5 decreased by about 40 μg/m³ (Figure 6), while they both decreased by ∼8 μg/m³ per year. NO₂ and O₃ also decreased by 3 μg/m³ per year. It can be observed from the monthly average variation (Supplementary Figure S6) that O₃ had high concentrations in summer (∼200 μg/m³) and low concentrations in winter (∼20 μg/m³). Contrary to O₃, other pollutants exhibited high and low concentrations in winter and summer, respectively. The difference between the simulated and observed O₃, SO_2, and PM_2.5 was less than 10 μg/m³.

In the YRD, the trio of O₃, SO_2, and PM_2.5 followed the same trends as observed in the NCP (Figure 6). The observed value of NO₂ decreased to 31 μg/m³ in 2017 and later increased in the subsequent year (Figure 5). Similar to the NCP, O₃ had high concentrations in summer (∼140 μg/m³) and low concentrations in winter (∼50 μg/m³), while the reverse was the case for other pollutants. Regarding the model performance, NO₂ was underestimated (>10 μg/m³) while other pollutants were well simulated with better model performances. The annual trends of observed pollutants in the PRD were not different from those of NCP and YRD (Figure 6). It can be observed from the monthly average variation (Supplementary Figure S6) that O₃ had high concentrations in summer (∼140 μg/m³) and low concentrations in winter (∼80 μg/m³). Contrarily, other pollutants displayed high and low concentrations in winter and summer, respectively. The simulated O₃ was overestimated (∼10 μg/m³) (Figure 5) and increased during the study period (Figure 6). NO₂ was underestimated, while PM_2.5 and SO₂ were well simulated with minimum bias (Figure 5).

The trends of the observed SO₂ and PM_2.5 in the CY were similar to the previously discussed regions (Figure 6). SO₂ decreased by 55 μg/m³ between 2013 and 2014 (Figure 5). During the entire study period (Figure 5), SO₂ and PM_2.5 were overpredicted and NO₂ was underpredicted, while O₃ was well simulated with minimum bias. The simulated SO₂ during the first 4 years (2013–2016) exhibited apparent seasonal variations, while the levels of observed SO₂ were generally the same with no significant changes during the entire study period. In the FW region, the observed O₃ increased steadily (Figures 5, 6), while SO₂ and PM_2.5 showed a decreasing trend during 2013–2019. Similar to other regions, high and low concentrations of O₃ were observed in summer and winter, respectively, while the other three pollutants had their high and low concentrations during winter and summer, respectively. PM_2.5 and SO₂ were well simulated (Figure 5), NO₂ was underpredicted (Supplementary Figure S5, S6), and O₃ was highly overestimated during 2013–2015 (Figure 5), while the model performance for O₃ greatly improved during 2016–2019.

Generally, O₃ in the PRD and FW increased significantly during the study period and was about 110 μg/m³. No significant change was found in NO₂ during 2013–2019. SO₂ and PM_2.5 decreased on yearly basis. All of the pollutants were well predicted in 2019 except SO₂ and PM_2.5 in the CY region.

Figure 7 shows the correlation coefficients (R) between PM_2.5 and O₃, while the correlation coefficients between PM_2.5, NO_2, and SO₂ are illustrated in Supplementary Figure S8. Considering the simulated concentrations, there was an apparent seasonal change between PM_2.5 and O₃ across the five regions (Supplementary Figure S7). The seasonality and correlation gradually weakened from north to south, with the difference of the NCP being more obvious than that of the YRD in different seasons. In the NCP, there was no correlation between PM_2.5 and O₃ in spring (Figure 7). PM_2.5 was positively correlated to O₃ in summer (0.2) and negatively correlated in autumn (-0.3) and winter (-0.8). The observed concentrations also exhibited similar correlations, but the correlations were closer to zero compared to the simulated concentrations (Figure 7). In the YRD, the simulated PM_2.5 was positively correlated with simulated O₃ during spring and summer, while negative correlations were found between them in autumn and winter. The correlation coefficients were all below 0.5. For the observed concentrations, positive correlations (<0.5) were found between the two pollutants during the four seasons. The difference between the correlation coefficients of the simulated and observed concentrations might be attributed to the overprediction of PM_2.5 as O₃ was well simulated (Figure 4). In the PRD, positive correlations were found between PM_2.5 and O₃ during the four seasons for both simulated and observed concentrations. In the CY and FW regions, the relationships between the simulated PM_2.5 and O₃ were all negative during the four seasons except in summer in the FW, which was positive. The poor correlations found with the simulated concentrations might be attributed to the location of the CY and FW in an inland area. Considering the observed concentrations, positive correlations were noted in spring and summer, while negative correlations were found during autumn and winter in the CY. In the FW, however, the observed PM_2.5 and O₃ displayed positive relationships during summer and autumn and negative correlations in both spring and winter.

In addition to the relationship between PM_2.5 and O₃, the correlations between PM_2.5, NO_2, and SO₂ were also assessed (Supplementary Figures S7, S8). Low-to-medium correlation coefficients were noted between the simulated PM_2.5 and NO₂, while high coefficients were found between the observed PM_2.5 and NO₂ during the four seasons across the five regions, an indication of the strong relationship between PM_2.5 and NO₂ across China during the study period. The low correlation between the simulated PM_2.5 and NO₂ could be associated with the underestimation of NO₂ across the five regions during the four seasons. In addition, in both simulation and observation scenarios, there were strong and positive correlations between PM_2.5 and SO₂ during the four seasons in all five regions, and the correlation coefficients for the two scenarios were relatively similar.