Diagnosing drivers of PM_2.5 simulation biases in China from meteorology, chemical composition, and emission sources using an efficient machine learning method

Chemical transport models (CTMs) are widely used for air pollution modeling, which suffer from significant biases due to uncertainties in simplified parameterization, meteorological fields, and emission inventories. Accurate diagnosis of simulation biases is critical for the improvement of models, interpretation of results, and management of air quality, especially for the simulation of fine particulate matter (PM_2.5). In this study, an efficient method with high speed and a low computational resource requirement based on the tree-based machine learning (ML) method, the light gradient boosting machine (LightGBM), was designed to diagnose CTM simulation biases. The drivers of the Community Multiscale Air Quality (CMAQ) model biases are compared to observations obtained by simulating PM_2.5 concentrations from the perspectives of meteorology, chemical composition, and emission sources. The source-oriented CMAQ was used to diagnose the influences of different emission sources on PM_2.5 biases. The model can capture the complex relationship between input variables and simulation bias well; meteorology, PM_2.5 components, and source sectors can partially explain the simulation bias. The CMAQ model underestimates PM_2.5 by −19.25 to −2.66 µg m⁻³ in 2019, especially in winter and spring and during high-PM_2.5 events. Secondary organic components showed the largest contribution to the PM_2.5 simulation bias for different regions and seasons (13.8 %–22.6 %) of all components. Relative humidity, cloud cover, and soil surface moisture were the main meteorological factors contributing to PM_2.5 bias in the North China Plain, Pearl River Delta, and northwestern China, respectively. Primary and secondary inorganic components from residential sources showed the two largest contributions to this bias (12.05 % and 12.78 %), implying large uncertainties in this sector. The ML-based methods provide valuable complements to traditional-mechanism-based methods for model improvement, with high efficiency and low reliance on prior information.

1 Introduction

Fine particulate matter (PM_2.5) is a complex mixture of primary particulate matter (PPM) and secondary inorganic and organic components (SIAs/SOAs), with adverse effects on public health and ecosystems. Ambient levels of PM_2.5 are influenced by meteorological conditions, emissions from different sources, and atmospheric chemical processes (World Health Organization, 2021; Xiao et al., 2022; Yang et al., 2016; J. Liu et al., 2021b; Zhai et al., 2019). China has experienced severe PM_2.5 pollution over the past 2 decades (Bai et al., 2022; F. Liang et al., 2020). For effective air quality management, accurate PM_2.5 modeling is essential. Chemical transport models (CTMs), like the Community Multiscale Air Quality (CMAQ) model, have been widely developed and applied to PM_2.5 simulations through the atmospheric processes of dispersion and deposition, as well as chemical reactions (Qiao et al., 2018; Wang et al., 2021; Hu et al., 2017a). Application of CTM simulations is often limited by the biases due to uncertainties in simplified parameterization, meteorological prediction, emission inventories, and initial and boundary conditions (Binkowski and Roselle, 2003; Hu et al., 2014, 2016; Wang et al., 2023a, 2021). Thus, it is essential to diagnose the specific sources of simulation biases according to specific model applications, including grid resolution, parameterization, mechanisms, meteorological inputs, and emission inventories.

Traditional bias diagnosis approaches for CTMs usually rely on empirical and prior assumptions with extensive sensitivity testing and high demands on computational resources such as Monte Carlo methods or Latin hypercube sampling (Beekmann and Derognat, 2003; Hanna et al., 2005; Aleksankina et al., 2019). Recently, machine learning (ML) methods, such as random forest and eXtreme Gradient Boosting (XGBoost), have been widely used in environmental science research due to their simple structure, high speed, and ability to deal with non-linear relationships (Liu et al., 2022). Many studies used ML to predict air pollutant concentrations like those of PM_2.5 and ozone (Wei et al., 2021a; Sun et al., 2021; Zhu et al., 2022; Bai et al., 2022), improve the accuracy of CTM simulations (Wang et al., 2023a; Wei et al., 2020), and explain the prediction results using interpretable ML techniques (Hou et al., 2022; Li et al., 2023; Stirnberg et al., 2021). To date, few studies have used ML to diagnose the simulation bias of CTMs. One study showed the potential of machine learning in explaining the simulation bias of ozone (Ye et al., 2022). However, as it is a complex multi-phase mixture, it is still challenging to diagnose biases in PM_2.5 simulations using ML methods (Liu and Xing, 2022). Moreover, given the significant impact of emissions, it is instructive to diagnose CTM biases of PM_2.5 based on a source apportionment perspective.

In this study, we use the light gradient boosting machine (LightGBM) model, an efficient ensemble ML approach, to diagnose the drivers of CMAQ biases in simulating PM_2.5 concentrations. A source-oriented version of CMAQ is used to track sectoral source contributions to PM_2.5. Model biases are diagnosed by observations from multiple perspectives, including meteorology, chemical components, and emission sources.

2 Materials and methods

2.1 Surface PM_2.5 observations

This study specifically targets the year of 2019 due to the extensive availability of observational data, the reliability of emission inventories, and the absence of COVID-19-related disruptions. Hourly PM_2.5 observations for 2019 are taken from the China National Environmental Monitoring Centre (CNEMC; http://www.cnemc.cn/, last access: 23 April 2024). The daily observation data < 0.1 % quantile and > 99.9 % quantile, data showing PM_2.5 exceeding PM₁₀, and days with fewer than 20 valid hourly records are excluded. For observation sites located on the same CMAQ simulation grid (36 km×36 km), their average PM_2.5 concentrations were calculated to be compared with CMAQ simulation. Approximately 350 000 observations, which met the quality control criteria, were selected from the entire time series data points collected from various monitoring stations. The distribution of (about 1200) observation sites is shown in Fig. S1 in the Supplement. The stations are unevenly distributed, with dense stations in eastern populated areas and sparse stations in western regions of Xinjiang and Tibet. Analysis has been carried out in several haze-polluted regions and the whole country (Fig. S1), including the Beijing–Tianjin–Hebei region (BTH), the Yangtze River Delta (YRD), the Pearl River Delta (PRD), the Sichuan Basin (SCB), and the region of northwestern China (NWCHN).

2.2 CMAQ simulation

The CMAQ simulation (36 km×36 km) was carried out to simulate PM_2.5 components in mainland China and surrounding regions in 2019. The list of PM_2.5 components simulated by CMAQ is shown in Table S1 in the Supplement. The Weather Research and Forecasting Model (WRF v4.2) was used to generate meteorological fields, driven by the National Centers for Environmental Prediction (NCEP) Final (FNL) Operational Model Global Tropospheric Analysis dataset (http://rda.ucar.edu/datasets/ds083.2/, last access: 23 April 2024) (NCEP, 2000). Several meteorological factors (Table S1) that are highly relevant to aerosol concentrations are selected for ML model building (Xiao et al., 2021; Z. Y. Chen et al., 2020; Meng et al., 2019). CMAQ v5.0.2 with a modified SAPRC-11 photochemical mechanism and an AERO6 aerosol module was applied to aerosol simulations (Carter and Heo, 2013; Ying et al., 2015; Binkowski and Roselle, 2003). The Multi-resolution Emission Inventory for China (MEIC) was used as a source of anthropogenic emissions (http://meicmodel.org/, last access: 23 April 2024), and the Model of Emissions of Gases and Aerosols from Nature (MEGAN) v2.1 was used to generate biogenic emissions (Guenther et al., 2012, 2006). The Fire INventory from the National Center for Atmospheric Research (FINN) based on satellite data was used to generate open burning emissions (Wiedinmyer et al., 2011).

The source apportionment method was used to quantify the contributions of the industry, energy, residential, transportation, agriculture, open burning, and biogenic sources to PPM and SIA concentrations using a modified version of CMAQ (Zhang et al., 2012; Ma et al., 2021; Qiao et al., 2018). PPM from different source sectors is tracked by non-reactive tracers (which account for 10⁻⁵ of the PPM emission rates), and source-specific PPM concentrations are then calculated by multiplying the tracer by 10⁵. The contributions of source sectors to SIAs are quantified using specific reactive tagged tracers. Specifically, NO_x, SO₂, and NH₃ from different sources were tracked separately through a series of chemical and physical processes involved in SIA formation. The source of SOAs was not traced due to the complex and currently imperfect mechanism of SOA formation and the high uncertainties in the precursor VOC emissions (J. Liu et al., 2021b; Hu et al., 2017b). Details on source apportionment can be found in previous studies (Zhang et al., 2012; Ma et al., 2021; Qiao et al., 2018; Ying et al., 2014). The contributions of source sectors to SOAs were not tracked due to insufficient knowledge of its precursors and incomplete formation mechanisms (Yang et al., 2019; Carlton et al., 2007; Zhang et al., 2011).

2.3 Machine learning method

Tree-based ML models typically outperform deep learning approaches in tabular data (e.g., air pollutant observation datasets) and thus have been widely developed (Grinsztajn et al., 2022). Wei et al. (2021a) compared several models when reconstructing PM_2.5 data records in China and found that the tree model showed superior performance. The LightGBM model is an optimized gradient boosting decision tree (GBDT) algorithm (Ke et al., 2017) and has shown accurate performance in many fields (Wei et al., 2021b; Yan et al., 2021; Sun et al., 2020; W. Liang et al., 2020). Compared to XGBoost, a widely used GBDT, LightGBM uses the histogram decision tree algorithm along with gradient-based one-side sampling (GOSS), which saves memory and computation time (Ke et al., 2017). Three tree-based models, random forest, XGBoost, and LightGBM, were compared in our previous study (Wang et al., 2023b). Using the same input data and hyperparameters, LightGBM is as accurate as XGBoost but faster and less susceptible to overfitting (the difference in accuracy between training and testing). Nevertheless, multicollinearity between features such as pollutant concentrations and meteorological factors can greatly affect the performance of traditional linear models. When two independent variables are correlated, changes in one variable are associated with changes in the other, making it difficult for the model to independently estimate the relationship between each independent and dependent variable. However, these collinearities do not affect the performance of tree-based models like random forest and LightGBM because they do not require the assumption of feature independence (Belgiu and Drăguţ, 2016; Chen and Guestrin, 2016; Ke et al., 2017). Therefore, the LightGBM model was used to diagnose PM_2.5 simulation biases in this study. Two metrics were calculated to evaluate model performance, including the coefficient of determination (R²) and the root mean square error (RMSE) (Wei et al., 2020).

$$\begin{array}{}\text{(1)}& {\displaystyle }{R}^{\mathrm{2}}=\mathrm{1}-{\displaystyle \frac{{\sum }_i(y_i-f_i{)}^{\mathrm{2}}}{{\sum }_i(y_i-\widehat{y}{)}^{\mathrm{2}}}}\text{(2)}& {\displaystyle }\text{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}{\sum }_{i=\mathrm{1}}^n(y_i-f_i{)}^{\mathrm{2}}}\end{array}$$

Cross-validation (five-fold) combined with the RMSE was used to select hyperparameters. The dataset was randomly divided into five parts; one was taken in turn as the test set, and the rest were used for training, which was repeated five times, and the average test RMSE was calculated by looping to increase model complexity, ending the loop, and returning to the hyperparameters when the model test RMSE does not decrease significantly (< 0.01) or the gap between the training and test RMSE increases significantly (< 0.05). The separate test sets (not involved in the training and CV hyperparameter selection process) were chosen by randomly sampling 20 % of the data from all stations in the region of interest.

The target variable was defined as the difference between observed and simulated daily PM_2.5 concentrations, and the key contributors to the simulation bias were identified through the relative importance (calculated by gain) of the input features (Ye et al., 2022; Loyola-González et al., 2023). Three categories of input variables were designed to separately fit the simulation biases: meteorological factors, chemical components, and emission sources. Meteorological factors, including wind speed, wind direction, temperature, humidity, surface pressure, cloud fraction, and boundary layer height, are used to investigate the impact of meteorology on the CMAQ simulation biases. The components of PM_2.5 are divided into SIAs (sulfate, nitrate, and ammonium), primary and secondary organic aerosols (POAs/SOAs), elemental carbon (EC), and other components. The contributions to the simulation bias were quantified using seven sectoral sources: industry, energy, residential sources, transportation, agriculture, open burning, and biogenic emissions.

3 Results and discussion

3.1 Observation and simulation of PM_2.5

Figure 1a shows the time series of observed and simulated daily surface PM_2.5 concentrations in China as a whole and in five specific regions (BTH, YRD, PRD, SCB, and NWCHN) during 2019. Observed PM_2.5 concentrations were the highest in the BTH region (51.172 µg m⁻³) and lowest in the PRD region (28.273 µg m⁻³). The CMAQ model underestimates PM_2.5 concentrations of −8.59, −2.66, −6.21, and −19.25 µg m⁻³ in the BTH, YRD, PRD, and NWCHN regions, respectively (Fig. 1b). Moreover, the underestimation occurred mainly in winter and spring (Fig. 1c) and during high-PM_2.5 events (Fig. 1d) (Hu et al., 2016; Huang et al., 2017).

Figure 1 (a) The time series of observed (black) and CMAQ-simulated (red) daily surface PM_2.5 concentrations in China and five regions. Mean concentrations of the observed and simulated PM_2.5 and MNB are also shown in the inset. (b) Box plots of CMAQ-simulated biases (simulated minus observed) for different regions. Crosses indicate average values and outliers are determined to be > 1.5 times the upper quartile and < 1.5 times the lower quartile. (c) Same as in panel (b) but for four seasons. Spring, summer, autumn, and winter are defined as March to May; June to August; September to November; and December, January, and February, respectively. (d) Same as in panel (b) but for different PM_2.5 concentration levels (L1: [0,35], L2: [35,75], L3: [75,115], L4: [115,150], and L5: [150,1000]).

Table S2 in the Supplement shows the validation of CMAQ simulations against observations in different regions. Four indicators (mean normalized bias, MNB; mean normalized error, MNE; mean fractional bias, MFB; mean fractional error, MFE) were used to systematically evaluate the performance of the CMAQ simulations. The PM_2.5 simulations in the BTH, YRD, and PRD regions were in better agreement with observations, with average MNB of −0.08, −0.07, and −0.08, respectively (within the standard of 0.66). The PM_2.5 simulations in SCB and NWCHN regions show large biases with MNB of 0.46 and −0.42, respectively. The differences in CMAQ performance between regions can be attributed to multiple factors, including emission inventories, dominant mechanisms of PM_2.5 generation, and topographic and meteorological conditions (Ma et al., 2021; Xue et al., 2019; Hu et al., 2014).

Annual and monthly mean PM_2.5 components (SIAs, POAs, SOAs, EC, and other components) were calculated for China as a whole and for five key regions (Fig. 2). PM_2.5 and its components show a similar spatial distribution, with high concentrations occurring in the eastern regions (SCB, BTH, and central YRD). SOAs showed high concentrations in summer over China (6.80 µg m⁻³), which could be related to enhanced solar radiation and atmospheric oxidation capacity in summer (precursors of SOAs such as isoprene are highly dependent on temperature and light) (Yang et al., 2019; S. Chen et al., 2020; J. Liu et al., 2021b). Nitrate and POAs were the dominant components in winter (10.14 and 9.11 µg m⁻³, respectively). In the BTH and SCB regions, POAs account for a higher proportion of total particles than nitrate in winter, whereas nitrate has a higher proportion in the YRD region. Nitrate showed higher concentrations than sulfate in most regions and seasons due to the implementation of coal combustion control policies (Shang et al., 2021; J. Liu et al., 2021b; Xu et al., 2019).

Figure 2 Annual mean concentration map (a–g) and monthly mean concentrations (h–m) of PM_2.5 and its components (SIAs, POAs, SOAs, EC, and other components) for China as a whole and five key regions in 2019. Dotted lines in (h–m) indicate PM_2.5 observations.

The results of the PM_2.5 sectoral source apportionment (Figs. 3 and S2 in the Supplement) show that industries and residential sources were the main contributors to daily PM_2.5 concentrations for all regions and seasons, with seasonal fractional contributions of 25.31 %–31.92 % and 11.13 %–30.64 %, respectively. The seasonal average fractional contributions from energy, transportation, and agricultural NH₃ in China as a whole were 3.26 %–5.67 %, 6.82 %–11.26 %, and 7.50 %–8.67 %, respectively. The contributions from biogenic sources were negligible in all regions and seasons (< 1 %). In contrast to the contributions from energy, transportation, industrial, and agricultural sources, significant seasonal variations occurred in emissions from residential sources in all five regions, with high contributions in winter (17.60 %–30.90 %) and low contributions in summer (5.53 %–16.46 %).

Figure 3 Seasonal average fractional contributions of different sources to PM_2.5 concentrations (black circle in relation to the right-hand axis) in China as a whole and five key regions.

As the secondary component makes up a large proportion of the total PM_2.5, the source sectors of SIAs were analyzed for five regions (Fig. S2). High concentrations of SIAs were found in winter (12.36–34.08 µg m⁻³), with large contributions from industrial and agricultural sources and transportation (31.49 %–36.41 %, 20.40 %–22.40 %, and 19.77 %–22.46 %). The low contribution of the residential sector to secondary PM_2.5 but the high contribution to total PM_2.5 indicates that most residential emission sources emit PPM directly, with a small fraction of PPM coming from secondary generation. The contributions from biogenic and open burning sectors to SIAs were relatively low in all regions and seasons (< 10 %).

3.2 Drivers of PM_2.5 simulation bias

The ML models were trained separately using information on meteorology, PM_2.5 components, and source sectors for different regions and seasons, and separate test sets were used to evaluate the model performance (Fig. 4). All three feature combinations can partially explain the simulation bias. The mean test R² values for meteorology, PM_2.5 components, and source sectors were 0.64, 0.52, and 0.50, respectively, and the RMSEs were 17.41, 19.82, and 19.56 µg m⁻³, respectively. The model performed better in summer than in winter. This may be attributed to the high simulation biases in winter due to severe PM_2.5 pollution and complex causes, while PM_2.5 pollution in summer is lighter with a lower CMAQ simulation bias.

Figure 4 Test results of CMAQ bias model training according to meteorology (a), PM_2.5 components (b), and source sectors (c). RMSE unit: µg m⁻³.

Using PM_2.5 components as input features to fit the total simulation bias enables the identification of components with a large simulation bias. Among the PM_2.5 components (Fig. S4 in the Supplement), SOAs showed the largest contribution to the PM_2.5 simulation bias for different regions and seasons (13.8 %–22.6 %), which is consistent with previous studies (J. Liu et al., 2021b; Yang et al., 2019; Fry et al., 2014). The inorganic aerosols (e.g., sulfates) are produced mainly by chemical pathways, while the SOAs are produced by the condensation of numerous partially oxidized gases and are therefore influenced by complex precursor concentrations and multi-stage oxidation processes. The incomplete description of SOA formation pathways in CTMs (simplified SOA parameterization) leads to significant differences between simulations and observations (Carlton et al., 2007; Zhang et al., 2018; Yang et al., 2019). In addition, biogenic emissions play an important role in SOA formation, with biogenic SOAs accounting for more than 70 % of the total SOAs in China during summer (Hu et al., 2017b; Wu et al., 2020), so the uncertainties in biogenic emissions can further contribute to the uncertainties in SOA emissions. Nitrate showed a large contribution to the PM_2.5 simulation bias in winter in the BTH region, which is consistent with a previous study (Liu and Xing, 2022). Nitrate contribution to the simulation bias further implies the inaccuracy of nitrate simulations, which can be related to the imperfect pathways of nitrate production (e.g., non-homogeneous oxidation) in SAPRC-11 (that we used) and the uncertainties in nitrate precursor emission inventories in winter (Xu et al., 2019; Zhang et al., 2018; Carter and Heo, 2013).

The contribution of meteorological factors to the simulation bias varies across regions and seasons (Fig. 5). In the BTH region, surface pressure and relative humidity contribute the most to the simulation bias. In the PRD region, relative humidity, cloud cover, and wind direction were the main contributors during all four seasons.

Figure 5 Contribution (%) of each meteorological factor to CMAQ simulation biases by region and season.

Humidity positively or negatively influences PM_2.5 concentrations through several mechanisms. By enhancing the PM_2.5 hygroscopic increase, promoting the secondary formation, and facilitating the gas-to-particle partitioning, high humidity positively influences PM_2.5 concentrations and contributes significantly to haze pollution (Z. Y. Chen et al., 2020; Cheng et al., 2015; Zhang et al., 2011). The contribution of humidity to CMAQ simulation biases can be partly attributed to the uncertainties in the WRF simulation. The mean RMSE of relative humidity from WRF simulations compared to observations was 20.38 % in this study (Table S3 in the Supplement). In addition, imperfections in the mechanism of humidity-promoted secondary particle formation (e.g., non-homogeneous catalysis of SOA) can also lead to simulation biases (Zhang et al., 2011; J. Liu et al., 2021b). Atmospheric pressure indirectly influences PM_2.5 concentrations through other meteorological factors (e.g., humidity and wind). High-pressure systems are connected to stationary weather, which is unfavorable for PM_2.5 dispersion. On the other hand, low pressure is usually accompanied by high humidity, influencing PM_2.5 nucleation, condensation, and coagulation and leading to increased PM_2.5 concentrations (Z. Y. Chen et al., 2020). Therefore, the influence of atmospheric pressure on the CMAQ simulation biases in the BTH region may be attributed to the uncertainties in the meteorological field (Bei et al., 2017; Zhang et al., 2015). The contribution of wind direction in the YRD region may also be related to the uncertainties in the WRF simulation (mean RMSE of 90.39°). Aerosols have feedback on meteorology (Qu et al., 2021). In addition to directly changing the radiation received by the Earth through scattering and absorbing (direct radiation effect), PM_2.5 can also influence radiation through aerosol–cloud interactions (indirect radiation effect) (Zhao et al., 2017; Yang et al., 2016). Moreover, PM_2.5 can act as cloud condensation and nucleation sites, contributing to clouds' microphysical development and precipitation formation process (Wang et al., 2020). However, the aerosol-to-meteorology feedback mechanism is missing in the CMAQ used in this study. A previous study showed the dominant role of cloud chemistry in aerosol–cloud interactions in southern China (Zhao et al., 2017). Therefore, the influence of cloud cover on simulation biases in YRD can be attributed to the lack of an aerosol feedback mechanism.

In the NWCHN region, soil surface moisture and stomatal conductance contributed significantly to the simulation bias. These factors can be associated with ground-level sand rise and dust emission (S. Liu et al., 2021). Underestimation of dust aerosol in NWCHN can be attributed to emission (both natural and anthropogenic sources), and an accurate emission inventory (empirical- or physical-based numerical models) should be established in northwest China by detailed activity data and emission factors (Hu et al., 2016; J. Liu et al., 2021a). In addition, the parameterization and mechanism for dust aerosol simulation in CMAQ should be further examined and improved.

Dry and wet days were separated to analyze the influence of humidity on the simulation biases (Table S4 in the Supplement). In most areas of China, CMAQ underestimates PM_2.5 more severely on dry days than on wet days, with larger absolute biases (−14.56, −7.09, −7.11, and −27.87 µg m⁻³ in spring, summer, autumn, and winter, respectively). On dry days, the BTH region showed severe underestimation in winter (−22.86 µg m⁻³), while the PRD region showed large simulation bias in spring (−21.55 µg m⁻³). Severe underestimation of PM_2.5 was observed on both wet and dry days in NWCHN. These underestimates of PM_2.5 on dry days can be related to the dry deposition process. Dry deposition is a critical but highly uncertain sink for aerosols which depends on the chemical and physical properties of aerosols and is influenced by land surface properties and meteorological conditions (Shu et al., 2022). Different land-use types (e.g., vegetation, deserts, and snow) possess markedly different capacities to capture particulate matter. The CMAQ model in this study used the dry deposition scheme PR11 from Pleim and Ran (2011). This study shows that the PR11 scheme underestimates PM_2.5 concentrations in China. Recent studies in the United States also showed an underestimation of PM₁₀ concentrations (Shu et al., 2022). Therefore, it is necessary to further develop and optimize the dry deposition scheme, especially for PM_2.5. PM_2.5 underestimation on wet days may be attributed to the biases in wet deposition and secondary organic aerosol formation under high humidity conditions (Wu et al., 2018; Ryu and Min, 2022; J. Liu et al., 2021b; Zhang et al., 2011).

Source sector contributions of PPM and SIAs (obtained from the source-oriented CMAQ) were used to build the model and diagnose the influences of different emission sources on PM_2.5 simulation biases (Fig. 6). The PPM and SIAs in residential emissions showed the largest contribution (12.05 % and 12.78 %) to the PM_2.5 simulation bias. The same conclusion was obtained when building a model with total PM_2.5 concentrations from different source sectors (Table S5 in the Supplement). PM_2.5 from residential emissions is the main contributor to the CMAQ simulation bias, accounting for 20.2 % of the total bias.

Figure 6 Contribution (%) of each source sector to CMAQ biases by region and season (res: residential, ene: energy, tra: transportation, agr: agriculture, ind: industry, AEC: elemental carbon, OTHER: other components).

In China, the residential sector consumed fossil fuels (coal, oil, and natural gas) and biofuels (wood and crop straw) with low combustion efficiency for cooking and heating and emitted large amounts of air pollutants (Li et al., 2017). However, due to the lack of reliable data (a locally accurate emission factor and fuel consumption data), the residential sector has been recognized as a major uncertainty source in anthropogenic emission inventories (X. Liu et al., 2021; Shen et al., 2021), which is consistent with the results identified by a machine learning model in this study. Therefore, developing an accurate residential sector emission inventory is essential for accurate PM_2.5 modeling, which requires reliable data regarding fuel consumption and emission factors based on fuel type, fuel characteristics, and combustion conditions (X. Liu et al., 2021).

3.3 Comparisons and uncertainties

Huang et al. (2019) used a new reduced-form model coupled with a higher-order decoupled direct method and a stochastic response surface model to identify the sources of uncertainty in CMAQ simulations. An analysis of the PRD region in China in the spring of 2013 revealed a systematic underestimation of SOA and identified wind speed and primary PM_2.5 emissions as the key sources of uncertainties in PM_2.5 simulations, which is consistent with the results identified using LightGBM in this study. Aleksankina et al. (2019) identified the PM_2.5 simulation bias in Europe using optimized Latin hypercube sampling and also demonstrated the important impact of primary emissions on PM_2.5 simulation uncertainties. Liu and Xing (2022) used a fully connected neural network to identify PM_2.5 simulation biases and found that NO₂ is the main contributor in BTH during heavy-pollution events in the winter, which is consistent with the main contribution of nitrate that we found in the BTH region (Fig. S4).

Although we filtered the features according to their relative importance and prior knowledge, collinearity still exists among the input features. Multicollinearity among features does not affect the performance of tree-based models like random forest and LightGBM (Belgiu and Drăguţ, 2016; Chen and Guestrin, 2016; Ke et al., 2017), but the contribution of a single feature might be slightly influenced by the other features. Previous studies (Hou et al., 2022; Ye et al., 2022) used ML to explain the causes of air pollution and model bias, and although there was multicollinearity between the input features they used, they got reliable conclusions, showing the minimal impact of multicollinearity and the reliability of tree-based machine learning methods.

The main objective of this study is to diagnose the contributors to CMAQ simulation biases using machine learning rather than for prediction. Since meteorology or emissions can only partially explain the simulation bias, a low R² is justified when fitting the model with only meteorology or emission variables (Fig. 4). We designed a complementary experiment to measure the impact of the model itself by comparing popular regression models (including multiple linear regression, polynomial regression (second-degree), random forest, XGBoost, and LightGBM) with the same features (PM_2.5 components). All models show similar performance (Table S6 in the Supplement); for example, all models show lower R² in winter in the BTH region (0.16–0.4) and higher R² in the SCB region (0.7–0.8). This is also evidence that the low R² is more affected by the features than the model itself, as the commonly used regression models can hardly obtain high R² with insufficient explanatory features (e.g., chemical component features in winter in BTH). Nevertheless, LightGBM shows comparable accuracy to XGBoost but is faster and shows smaller accuracy gaps between training and testing data with less overfitting.

Previous pollution prediction studies based on tree models usually added time-related features to describe the temporal pattern of pollutant changes to further improve the prediction ability; for example, Wei et al. (2021a) improved the model performance by adding temporal features of the day of year and Unix timestamps. However, the inclusion of temporal features cannot provide any useful information about contributors to simulation biases; instead, it is difficult to attribute them to meteorological or emissions contributions. Therefore, temporal features were not included in our model. Nevertheless, the ML bias diagnosis model constructed in this study is based entirely on local data. Some temporal and regional processes influencing PM_2.5 concentrations are not considered in this study, such as vertical transport and long-distance transport, which should be better diagnosed in future work, and the main bias contributors identified by variable importance are in good agreement with the current findings.

4 Conclusion

Based on artificial intelligence technology, this study systematically diagnoses the possible drivers of biases in PM_2.5 simulations from the perspectives of meteorology, chemical components, and emission sources. The relative importance of multiple factors helps to understand the sources of simulation bias and the deficiencies of the CMAQ mechanisms. SOA is the main contributor to simulation biases among chemical components. PM_2.5 is more underestimated in dry weather. Among source sectors, residential sectors contributed the most simulation bias for both PPM and SIAs. These results provide valuable information for CMAQ model improvement in terms of SOA and dust aerosol underestimation, meteorological field uncertainties, imperfection of the dry deposition scheme, and inaccurate residential emission inventories. As efficient bias diagnosis methods, machine-learning-based methods provide valuable complements to traditional-mechanism-based methods. This approach also greatly reduces the amount of prior information needed for diagnosing simulation bias and efficiently identifies the important contributors, so it can be easily extended to other CTMs and other pollutants.