Long-term variability in black carbon emissions constrained by gap-filled absorption aerosol optical depth and associated premature mortality in China

Black carbon (BC) plays an important role in air quality, public health, and climate, but its long-term variations in emissions and health effects have been insufficiently understood for China. Here, we present the spatiotemporal evolution of BC emissions and the associated premature mortality in China during 2000–2020 based on an integrated framework combining satellite observations from the Ozone Monitoring Instrument (OMI), an extreme gradient boosting (XGBoost) algorithm, a “top-down” inversion approach, and an exposure–response model. We found that the “bottom-up” approach most likely underestimated BC emissions, particularly in less developed western and remote areas. Pollution controls were estimated to reduce the annual BC emissions by 26 % during 2010–2020, reversing the 8 % growth during 2000–2010. BC emissions in the main coal-producing provinces declined by 2010 but rebounded afterwards. By contrast, provinces with higher economic and urbanization levels experienced emission growth (0.05–0.10 Mg km⁻² yr⁻¹) by 2010 and declined greatly (0.07–0.23 Mg km⁻² yr⁻¹) during 2010–2020. The national annual BC-associated premature mortality ranged between 733 910 (95 % confidence interval: 676 790–800 250) and 937 980 cases (864 510–1 023 400) for different years. The changing BC emissions contributed to an increase of 78 590 cases (72 520–85 600) during 2000–2005 and a decrease of 133 360 cases (123 150–145 180) during 2010–2015. Strategies differentiated by region are needed for further reducing BC emissions and their health and climate impacts.

1 Introduction

Black carbon (BC), commonly emitted during incomplete combustion of fossil fuels (Bond et al., 2013; Y. Liu et al., 2022; Shindell et al., 2012), is an important species in airborne fine particulate matter (PM_2.5). BC poses greater health risks than total PM_2.5 due to its absorption of harmful matters (polycyclic aromatic hydrocarbons and volatile organic compounds) and its penetration abilities (Y. Li et al., 2016; Pani et al., 2020; Wang et al., 2014; Xue et al., 2021), and it is a crucial short-lived climate forcer (Harmsen et al., 2020; Samset et al., 2020). With a large population and high energy consumption, China has become a major contributor of global BC emissions (Lu et al., 2019; Wang et al., 2012) and has suffered from BC-associated climate and health effects since the 2000s (Gu et al., 2020; Y. Liu et al., 2022). Compared with widely measured total PM_2.5 across the country (Liang et al., 2020; Zhang et al., 2019), fewer BC data are available from ground observations, and the spatiotemporal coverage of BC concentrations is far less sufficient (Cui et al., 2015; Tao et al., 2017). As a result, the long-term evolution of BC pollution and its associated health burden remain unclear.

Alternatively, satellite observations provide broader spatiotemporal coverage of aerosol-related variables (Schutgens et al., 2021), e.g., aerosol absorption optical depth (AAOD) that reflects light extinction due to light-absorbing aerosols including BC. However, most sensors can only monitor total aerosol information rather than individual components (S. S. Li et al., 2016), and cloud cover and surface reflectance lead to considerable missing values and uncertainty (Liang et al., 2020; Zhang et al., 2015). For example, the Ozone Monitoring Instrument (OMI; Stein-Zweers and Veefkind, 2012; https://disc.gsfc.nasa.gov/datasets/OMAEROe_003/summary, last access: 10 March 2022) and the POLarization and Directionality of the Earth's Reflectance instrument (POLDER; https://www.grasp-open.com, last access: 4 May 2022) provided long-term national average AAOD coverage of 9 %–22 % (2005–2020) and 8 %–12 % (2006–2013) in China, respectively. Satellite-derived AAOD needs to be comprehensively processed to fill gaps in the data and to improve its representativeness of BC before it can be effectively applied.

Complete and reliable emission estimates are essential for diagnosing pollution sources and evaluating the benefits of pollution controls. Compared with species generated largely from industrial and energy infrastructures (e.g., SO₂ and NO_x), BC emissions are more challenging to estimate, as they are commonly from residential and commercial sources that are more difficult to track (Bond et al., 2013; Li et al., 2017; Zhu et al., 2019). Existing “bottom-up” estimates varied between 0.9 and 2.5 Tg yr⁻¹ during 2000–2020, with inconsistent interannual changing patterns (the Multiresolution Emission Inventory for China (MEIC; Tsinghua University, 2023), the Emissions Database for Global Atmospheric Research (EDGAR; European Commission, 2022), Community Emissions Data System (CEDS; McDuffie et al., 2020), the Peking University Fuel Inventory (PKU-Fuel; Wang et al., 2014), Regional Emission inventory in ASia (REAS; Kurokawa and Ohara, 2020), and others (Lu et al., 2011; Lei et al., 2011; Klimont et al., 2009; Qin and Xie, 2012)). The uncertainty of those estimates reached ±360 % due to diverse and quickly changing manufacturing technologies and emission controls (Streets et al., 2003; Wang et al., 2016). Consequently, chemical transport models (CTMs) often underestimate BC concentrations and AAOD, particularly in Asia, ranging from a factor of 2 to a factor of 10 (L. Chen et al., 2019; Hu et al., 2016; Wang, 2015). To overcome this limitation, “top-down” approaches constraining BC emissions with available observations have been developed and applied to correct BC emissions in China (Cohen and Wang, 2014; Evangeliou et al., 2018; Fu et al., 2012; Guerrette and Henze, 2017; Wang et al., 2013; Zhao et al., 2019). However, restricted by insufficient spatiotemporal coverage, studies were usually conducted for individual years/months and showed considerable discrepancies (Wang et al., 2018; Zhang et al., 2015). Incomplete and inconsistent information could hardly be combined to provide full knowledge of long-term BC emissions. Based on the bottom-up emission estimates with great uncertainty and on CTMs, previous studies have evaluated the BC-associated premature mortality in China for limited years (2000, 2013, and 2016; Cui et al., 2022; Qin et al., 2019; Saikawa et al., 2009; Wang et al., 2021). Moreover, the results showed magnitude discrepancy (50 100–1 436 960 cases), and few analyses have evaluated the long-term spatiotemporal variations and driving forces of BC-associated health effects. The influence of human activities on quickly changing BC emissions and their associated health impact is inadequately or inaccurately understood, weakening science-based decision-making for air pollution control.

Herein, we developed an integrated framework combining available satellite observations, an improved machine learning technique, a top-down inversion approach, and an exposure–response model to obtain a panoramic perspective of China's BC emissions and the associated mortality for the past 2 decades. We first predicted full-coverage monthly AAOD for mainland China during 2000–2020 using an extreme gradient boosting (XGBoost) model. Combining this new dataset with air quality and BC light absorption empirical models, we then improved the top-down inversion technique to estimate the interannual changes in BC emissions. We further calculated the BC-associated premature mortality and attributed its interannual changes to individual driving factors. The outcomes highlight an improved BC emission estimation and the influence of human activities on the long-term evolution of BC emissions and the associated health effects, thereby supporting policies coordinating air quality, health, and climate issues.

2 Materials and methods

2.1 Filling gaps in AAOD data using a machine learning algorithm

We applied the XGBoost model to fill gaps in satellite-derived AAOD data at the monthly level during 2000–2020. XGBoost has been widely used in predicting air pollution and shown to outperform various statistical and machine learning models (Liang et al., 2020; S. Liu et al., 2022; Wang et al., 2023; Xiao et al., 2018). The XGBoost algorithm is an additive model based on hundreds of decision tree models. It first builds multiple classification and regression trees and then integrates these trees as a new tree model using an additive function (Liu et al., 2021). The model continues to iteratively improve, and the new tree model generated in each iteration will fit the residual of the previous tree. The complexity of the ensemble model will gradually increase until the training achieves the best results. Different from the boosting approach of XGBoost, the random forest model fits a set of decision trees, and then a majority vote method is taken for final prediction (Lyu et al., 2019). Generally, the XGBoost model requires less training and prediction time and has better performance than the random forest model. The target domain included mainland China at a horizontal resolution of 0.25° × 0.25° (Fig. S1 in the Supplement).

OMI takes advantage of the near-UV algorithm based on the sensitivity of radiances measured in the atmosphere to the varying aerosol species in order to derive AAOD (Zhang et al., 2017). The algorithm excluded the very low AOD values to reduce the uncertainty of the AAOD retrieval at low AOD values. Previous studies have proven the good agreement between OMI AAOD and Aerosol Robotic Network (AERONET) ground observations, as well as other satellite observations (Ahn et al., 2008; Zhang et al., 2017). Here we used OMI-derived AAOD at 483 nm, obtained from the OMAEROe L3 global aerosol product at a horizontal resolution of 0.25° × 0.25° (Stein-Zweers and Veefkind, 2012; https://disc.gsfc.nasa.gov, last access: 10 March 2022), as the dependent variable for model training and validation. For each grid cell, daily AAOD values for no less than 7 d in a given month were averaged as the monthly AAOD value. Owing to its long service time and damage to the satellite sensor, the original spatial coverage of monthly OMI-derived AAOD ranged from 1 % to 53 % and was commonly lower for later years than earlier years. The multiyear average coverage was relatively low in southern China (< 25 %), attributed to cloud cover and high surface reflectance (Fig. S2).

A total of 24 interpretation variables were selected for model training, including aerosol optical, meteorological, geographic, and temporal parameters (Table S1 in the Supplement). Aerosol optical and meteorological parameters were extracted from the Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2) dataset at a horizontal resolution of 0.625° × 0.5° (https://disc.gsfc.nasa.gov) to reflect the optical properties, transport, and diffusion of pollutants. Here we use the daily MERRA-2 data because the daily average MERRA-2 data were proved to be more reliable than the hourly data when compared with the observation (Xu et al., 2020). As ancillary variables associated with BC emission sources and transport conditions, land-use data were obtained from the Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences, at a horizontal resolution of 1 × 1 km (https://www.resdc.cn/DOI/DOI.aspx?DOIID=54​​​​​​​​​​​​​​, last access: 5 June 2024). The elevation data were obtained from the Shuttle Radar Topography Mission at a horizontal resolution of 1 × 1 km (https://www.resdc.cn/data.aspx?DATAID=123, last access: 5 June 2024). These parameters were resampled to the 0.25° × 0.25° grid system by averaging the 1 km resolution data. Details of specific feature selection can be found in Sect. S1 in the Supplement.

To evaluate model performance, we applied a 10-fold cross-validation (CV) to evaluate out-of-sample accuracy. The CV process randomly split training data records into 10 subsets, in which 9 subsets were used to train models and the remaining 1 was used to examine the performance. Through 10 CV repetitions, all of the data records were tested once. Besides, we collected monthly AAOD data from four aerosol monitoring networks to further verify the model reliability, including the Aerosol Robotic Network (AERONET, https://aeronet.gsfc.nasa.gov, last access: 10 March 2022), China Aerosol Remote Sensing Network (CARSNET; Che et al., 2015), Campaign on Atmospheric Aerosol Research Network of China (CARE-China; Xin et al., 2015), and Sun Sky Radiometer Observation Network (SONET; Li et al., 2018), as shown in Fig. S1a. Detailed site descriptions can be found in corresponding studies. Given the complicated technologies and high costs required for measurement operation, instrument maintenance, and calibration, current aerosol monitoring sites are rare and unevenly located in the country. There are clear missing values in the time series, and most measurements we collected focused on the period 2015–2019. All the ground-level AAOD data were interpolated to 483 nm using the Angstrom exponent to independently evaluate the performance of machine learning predictions. Model performance was evaluated with selected statistical indicators including the correlation coefficient (R), normalized mean error (NME; Eq. 1), normalized mean bias (NMB; Eq. 2), and root mean square prediction error (RMSE; Eq. 3):

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathrm{NME}={\displaystyle \frac{\sum _{i=\mathrm{1}}^N|X_i-O_i|}{\sum _{i=\mathrm{1}}^N{O}_i}},\end{array}$$

$$\begin{array}{}\text{(2)}& {\displaystyle }\mathrm{NMB}={\displaystyle \frac{\sum _{i=\mathrm{1}}^N\left(X_i-O_i\right)}{\sum _{i=\mathrm{1}}^N{O}_i}}\times \mathrm{100}\mathit{\%},\text{(3)}& {\displaystyle }\mathrm{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N{\left(X_i-O_i\right)}^{\mathrm{2}}},\end{array}$$

where X and O indicate the results from XGBoost prediction and OMI observation, respectively, and N is the number of data points.

2.2 Constraining BC emissions with gap-filled AAOD and CTM

We developed a top-down inversion approach to estimate the monthly BC emissions in China during 2000–2020 (see conceptual diagram in Fig. 1). To avoid abundant calculations, 5-year intervals were adopted in the simulation, and January, April, July, and October were selected as representative months of different seasons (widely applied in previous inversion studies; Zhang et al., 2015; Zhao et al., 2019).

Figure 1 The top-down inversion approach to estimate monthly BC emissions.

2.2.1 Integrated model for AAOD simulation

We simulated AAOD using the Community Multi-scale Air Quality (CMAQ; US EPA Office of Research and Development, 2015) model version 5.1 and an empirical BC light absorption model. AAOD is defined as the integrated absorption coefficient (m⁻¹) over the atmospheric column, and the absorption coefficient is the product of the mass concentration (g m⁻³) and mass absorption efficiency of BC (MAE, m² g⁻¹) (Bond et al., 2013). Prior anthropogenic BC emissions during 2000–2020 were obtained from MEIC (http://www.meicmodel.org, last access: 25 May 2022), and prior BC emissions from open biomass burning (OBB) were obtained from the Global Fire Emissions Database version 4.1s (GFED V4.1s, https://www.geo.vu.nl/~gwerf/GFED/GFED4/, last access: 5 June 2024). First, BC concentrations at different vertical layers were simulated using the CMAQ model at a horizontal resolution of 27 × 27 km (see Sect. S2 for model settings). Based on our previous measurements, we then adopted an empirical model to quantify the enhanced light absorption of the coating on BC particles (D. Chen et al., 2019) and obtained the simulated BC AAOD:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{\mathrm{MAE}}_{i,m,n}& =\mathrm{6.83}-\mathrm{0.0007}\times {\left({\displaystyle \frac{[\text{NA-PM}{]}_{i,m,n}}{[\mathrm{BC}{]}_{i,m,n}}}\right)}^{\mathrm{2}}\\ & +\mathrm{0.08}{\displaystyle \frac{[\text{NA-PM}{]}_{i,m,n}}{[\mathrm{BC}{]}_{i,m,n}}}\end{array}\end{array}$$

$$\begin{array}{}\text{(5)}& \mathrm{AAOD}\mathrm{\_}\mathrm{BC}\mathrm{\_}{\mathrm{sim}}_{i,m,n}=\underset{z_{\mathrm{0}}}{\overset{z}{\int }}{\mathrm{MAE}}_{i,m,n}\times [\mathrm{BC}{]}_{i,j,m,n}\times \mathrm{d}z,\end{array}$$

where MAE, AAOD_BC_sim, and [BC] represent the simulated BC MAE, BC AAOD, and BC concentration, respectively; z and z₀​​​​​​​ represent the simulated top and bottom of the atmosphere (0), respectively; [NA-PM] represents the simulated concentration of total non-absorbing matter (i.e., SO${}_{\mathrm{4}}^{\mathrm{2}-}$, NO${}_{\mathrm{3}}^{-}$, and organic carbon, OC); i and j represent the number of grids and vertical layers, respectively; dz represents the height of the vertical layer; and m and n represent the year and month, respectively. All parameters and variables shown in the equations are summarized in Table S2.

2.2.2 Inversion system for BC emissions

We developed an inversion system based on the spatiotemporal-dependent relationship between BC emissions and BC AAOD.

As AAOD is attributed to all light-absorbing aerosols, including BC, dust, and brown carbon (BrC), we first separated the contribution of BC to the XGBoost-predicted AAOD obtained in Sect. 2.1 using the fraction of BC in AAOD obtained from MERRA-2:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}\mathrm{AAOD}\mathrm{\_}\mathrm{BC}\mathrm{\_}{\mathrm{xgb}}_{i,m,n}& =\mathrm{AAOD}\mathrm{\_}{\mathrm{xgb}}_{i,m,n}\\ & \times {\displaystyle \frac{\mathrm{AAOD}\mathrm{\_}\mathrm{BC}\mathrm{\_}\mathrm{merra}{\mathrm{2}}_{i,m,n}}{\mathrm{AAOD}\mathrm{\_}\mathrm{merra}{\mathrm{2}}_{i,m,n}}},\end{array}\end{array}$$

where AAOD_BC_xgb represents the separated XGBoost BC AAOD, AAOD_xgb represents XGBoost AAOD, AAOD_BC_merra2 represents MERRA-2 BC AAOD, and AAOD_merra2 represents MERRA-2 AAOD. The hourly MERRA-2 data at 0.625° × 0.5° were reallocated to the horizontal resolution of the CMAQ model (27 km) and averaged to a monthly level.

With XGBoost BC AAOD, we inferred monthly BC emissions using Eq. (7):

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{E}_{\mathrm{posterior}i,m,n}& =E_{\mathrm{prior}i,m,n}\\ & \times \left(\mathrm{1}+{\displaystyle \frac{\begin{array}{c}\mathrm{AAOD}\mathrm{\_}\mathrm{BC}\mathrm{\_}{\mathrm{xgb}}_{i,m,n}-\\ \mathrm{AAOD}\mathrm{\_}\mathrm{BC}\mathrm{\_}{\mathrm{sim}}_{i,m,n}\end{array}}{\mathrm{AAOD}\mathrm{\_}\mathrm{BC}\mathrm{\_}{\mathrm{xgb}}_{i,m,n}}}\times {\mathit{\alpha }}_{i,m,n}\right),\end{array}\end{array}$$

where E_posterior and E_prior represent posterior and prior BC emissions, respectively; AAOD_BC_sim represents the simulated BC AAOD based on prior BC emissions; and α is a unitless factor representing the sensitivity of changes in BC AAOD to those in BC emissions in each model grid. We carried out a perturbation simulation to obtain α:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}{\mathit{\alpha }}_{i,m,n}& ={\displaystyle \frac{\mathrm{\Delta }{E}_{\mathrm{perturbed}i,m,n}}{E_{\mathrm{prior}i,m,n}}}\\ & ÷{\displaystyle \frac{\begin{array}{c}\mathrm{AAOD}\mathrm{\_}\mathrm{BC}\mathrm{\_}{\mathrm{sim}}_{\mathrm{perturbed}i,m,n}-\\ \mathrm{AAOD}\mathrm{\_}\mathrm{BC}\mathrm{\_}{\mathrm{sim}}_{\mathrm{prior}i,m,n}\end{array}}{\mathrm{AAOD}\mathrm{\_}\mathrm{BC}\mathrm{\_}{\mathrm{sim}}_{\mathrm{prior}i,m,n}}},\end{array}\end{array}$$

where prior and perturbed represent prior and perturbation simulations, respectively; $\mathrm{\Delta }{E}_{\mathrm{perturbed}}/E_{\mathrm{prior}}$ represents a 10 % reduction in prior BC emissions; and AAOD_BC_sim_perturbed and AAOD_BC_sim_prior represent the simulated BC AAOD with the perturbation and prior simulation, respectively.

We adopted the posterior BC emissions as the new BC emission input, repeating the simulation until the NME of BC AAOD from CTM and XGBoost was reduced to < 30 %. We evaluated the CMAQ model performance based on our gap-filled AAOD dataset and available observations of surface BC concentrations collected from 64 published studies (Table S3). The collected BC concentration observations cover various sampling regions in China and the study period from 2000 to 2020. Most studies analyzed BC using well-acknowledged reliable and widely used analyzers, for example, a DRI carbon analyzer or Sunset carbon analyzer. Besides this base case as mentioned above, we conducted four sensitivity tests to recalculate posterior BC emissions in order to explore the uncertainty in the inversion (Sect. S3 and Table S4). We applied the estimated AAOD at a longer wavelength of 865 nm (Test 1), the different dust AAOD fractions (Test 2), the adjusted MAE based on observations (Table S5; Test 3), and the adjusted simulated BC lifetime (Test 4).

2.3 Estimating the associated mortality burden and determining its drivers

A log-linear model was applied to estimate the attributable fraction (AF) of premature mortality to BC exposure (Wang et al., 2021):

$$\begin{array}{}\text{(9)}& \mathrm{AF}\left(C_{i,m}\right)=\mathrm{1}-e^{-{\mathit{\beta }}_{\mathrm{BC}}\times \mathrm{\Delta }{C}_{i,m}},\end{array}$$

where C_i,m represents the posterior simulated average BC concentration of 4 months in year m of grid i, and ΔC_i,m represents the difference between C_i,m and the health impact threshold. Due to the lack of reported BC concentration thresholds in current epidemiological studies (Cui et al., 2022), we applied the 1.25 percentile of BC concentrations (suggested by Pani et al., 2020; Wang et al., 2021) as the threshold (0.02 µg m⁻³). β_BC represents the concentration–response coefficient. Here, we used a β_BC value of 0.0204 (95 % confidence interval (CI): 0.0187–0.0224) based on a unique cohort study conducted in eastern China (Chen et al., 2021).

The all-cause premature mortality (M) attributable to BC exposure was calculated using Eq. (10) (Wang et al., 2021):

$$\begin{array}{}\text{(10)}& M_{i,s,m}=P_{i,m}\times {\mathrm{PS}}_{s,m}\times B_{s,m}\times \mathrm{AF}\left(C_{i,m}\right),\end{array}$$

where s represents the population subgroup, P represents the population, PS represents the proportion of the population subgroup to the national population, and B represents the national baseline mortality rate of all-cause diseases. The gridded population data were aggregated from the 1 km population density dataset in WorldPop (WorldPop, 2018; https://doi.org/10.5258/SOTON/WP00674). We corrected the annual total population using Chinese census data obtained from the State Statistics Bureau (https://data.stats.gov.cn/english/​​​​​​​, last access: 15 February 2023). The national average population age structure and baseline mortality rate of all-cause diseases were collected from the Global Burden of Disease study (Cohen et al., 2017; https://vizhub.healthdata.org/gbd-results/, last access: 5 June 2024).

We evaluated the impact of each of the four factors in Eq. (10) (three vulnerability factors and BC exposure) on the changing mortality through a series of sensitivity analyses (Cohen et al., 2017; Geng et al., 2021). We established 24 decomposition sequences with the four factors and calculated the mean changing mortality of each factor through all sequences. Moreover, mortality changes attributed to BC exposure were further disaggregated into contributions from BC emissions and meteorological factors using a direct proportion approach with the CMAQ model (Sect. S4).

3 Results and discussion

3.1 Gap-filled AAOD during 2000–2020

3.1.1 Evaluation of XGBoost model performance

By filling the missing values with the XGBoost model, we obtain full coverage of monthly AAOD for China from 2000 to 2020. When evaluated by 10-fold CV, the predicted AAOD shows good agreement with OMI observations, with R = 0.92, RMSE = 0.013, and NMB = −4 % (Fig. 2). A comparison with individual ground measurements further verifies the reliability and robustness of the model across regions and periods, with an RMSE value of 0.017 and NMB of 5 % for all available observations (Table S8). The comparisons with observations at typical individual sites are shown in Fig. S5. Overall, better performance is shown by our predictions compared with other AAOD datasets. As a reference, the RMSE and NMB between MERRA-2 and the same observations are 0.021 and −19 % (Table S8), respectively. In addition to the evaluation for the whole period, Table S9 shows the model performance for each year. The performance of 10-fold CV becomes moderately poorer for more recent years, accompanied by the reduced sampling size of OMI observations. Regarding the ground observations, a better performance appeared in recent years, indicated by the smaller RMSE, and the underestimation in earlier years could probably be attributed to fewer sites being available and more difficulty in data quality control.

Figure 2 Comparison between XGBoost-predicted and OMI-observed AAOD for the period 2005–2020. Different colors of dots represent values for different years. The dashed red line indicates the 1 : 1 line. The dashed blue line indicates the regression line. The interval of bins of the marginal histograms is 0.02.

Besides, the spatial coverage of OMI AAOD influences the training data size for the XGBoost model and thereby the spatiotemporal pattern of gap-filled AAOD. Table S10 shows the XGBoost performance against OMI and ground measurements by OMI coverage. In general, a poorer performance was found for areas with more missing values. When evaluated against ground measurements, for example, the RMSE and NMB for areas with OMI coverage of less than 20 % are 0.019 and 9 %, respectively, and they decline to 0.015 and 3 % for areas with OMI coverage of more than 60 %. All the biases are kept within a limited range. The analysis indicates the satisfactory quality of our gap-filled AAOD dataset, with full spatiotemporal coverage for the research domain and period.

Table S1 summarizes the importance levels of model predictors in the XGBoost model, expressed with three indicators; “Gain” is the most decisive predictor. MERRA-2 AAOD, longitude, latitude, and temporal parameters are identified as the most dominant variables for filling the gap of OMI AAOD. Furthermore, crop coverage is another dominant variable, which reflects the contribution of OBB to light-absorbing aerosol emissions and thereby to AAOD. Specific meteorological parameters – e.g., surface pressure (PS), shortwave radiation flux (SWGDN), evaporation from turbulence (EVAP), and planetary boundary layer height (PBLH) – reflect the surface energy budget, transport, and diffusion of air pollutants and thus play an important role in AAOD prediction.

Figure 3 (a) Spatial distribution of multiyear average AAOD during 2000–2020 and (b–h) interannual variations of AAOD for China and six key regions during 2000–2020. The gray and white areas represent phase 1 (2000–2012) and phase 2 (2013–2020), respectively. The red dots and dashed line represent time series of monthly AAOD (left vertical axis). The solid black lines represent the interannual variability after removing the seasonal change through time-series decomposition (right vertical axis). The straight red and blue lines present the linear trends of AAOD for different phases (right vertical axis). The annual variation rates (1 yr⁻¹) during different phases with significance levels (^∗ p < 0.05, ${}^{\ast \ast }$ p < 0.01, ${}^{\ast \ast \ast }$ p < 0.001) are presented. The map data provided by the Resource and Environment Data Cloud Platform are freely available for academic use (https://www.resdc.cn/DOI/DOI.aspx?DOIID=122, last access: 30 May 2024). © Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences.

3.1.2 Spatiotemporal patterns of XGBoost-predicted AAOD during 2000–2020

Figure 3a illustrates the spatial distribution of the averaged XGBoost-predicted AAOD during 2000–2020. Hot spots mainly existed in eastern China, with a regional average of 0.05, which was higher than the national average of 0.03. The AAOD values of Beijing–Tianjin–Hebei (BTH), Fenwei Plain (FWP), Yangtze River Delta (YRD), Sichuan Basin (SCB), and Northeast China (NE) were 1.4–1.8 times those of the national average, while the value for the Pearl River Delta (PRD) was much closer to the national average (see Fig. S1a for the locations of regions). The relatively small proportion of the rural population to the total population (34 %) and the highly developed economy in PRD might have resulted in limited light-absorbing aerosol emissions.

Figure 3b–h illustrate the interannual variability in the predicted monthly AAOD for China and the key regions during 2000–2020, which were divided into two temporal phases. In phase 1 (2000–2012), the AAOD of China experienced a slight decline (−6.94 × 10⁻⁵ yr⁻¹), while that of key regions showed a moderate increase, except for FWP. The annual AAOD growth of the most economically developed regions (BTH, YRD, and PRD) ranged 1.13 × 10⁻⁴–5.20 × 10⁻⁴ yr⁻¹, larger than that in other regions (SCB and NE, 2.84 × 10⁻⁷–7.54 × 10⁻⁵ yr⁻¹), reflecting the influence of regional differences based on human activities and meteorological conditions. Increasing industrial production and residential combustion led to elevated anthropogenic emissions of light-absorbing aerosols, thereby increasing AAOD values in key regions. However, nationwide increasing precipitation during phase 1 (0.3 g m⁻² s⁻¹ yr⁻¹ indicated by MERRA-2) may have enhanced the removal of pollutants and offset the effects of increasing emissions. In addition, strengthened afforestation policies in northern China led to a decline in dust aerosol emissions and thus also a decline in the AAOD values during this period (Middleton, 2019).

In phase 2 (2013–2020), the AAOD of China demonstrated a clear decline (1.99 × 10⁻⁴ yr⁻¹), with faster rates in key regions (3.61 × 10⁻⁴–1.42 × 10⁻³ yr⁻¹). The AAOD decline in key regions in phase 2 was much faster than its growth in phase 1, indicating the benefits of China's air pollution controls, such as implementation of the strictest-ever emission standards in the industrial sector and promotion of clean energy use for household heating and cooking. AAOD growth clearly occurred in northern regions from 2018 onwards, which differed from its continuous decrease in southern regions. The increasing surface wind speed in northern regions elevated the frequency of sandstorms, resulting in greater regional emissions of light-absorbing dust aerosols that partly contributed to AAOD growth (Yang et al., 2021).

3.2 Long-term evolution of constrained BC emissions during 2000–2020

3.2.1 Verification of constrained BC emissions

We find that application of posterior BC emissions constrained by XGBoost-predicted AAOD (described in detail in Sect. 3.2.2) largely improved the model performance of BC concentration and AAOD simulation compared with the use of prior BC emissions.

Showing an acceptable performance in meteorological simulation (Table S11), the CMAQ model presented a clear underestimation of surface BC concentrations based on the prior BC emissions, with the NMB and NME calculated to be −46 % and 53 %, respectively (Fig. 4). Besides, a larger underestimation appeared for the very early and most recent year, with NMB calculated at −59 % and −60 % for 2000 and 2020, respectively (Fig. S6), which may be caused by the larger underestimation of BC emissions in these years (described in detail in Sect. 3.2.2). Application of posterior BC emissions greatly reduced the NMB and NME during the research period to −14 % and 36 %, respectively (Fig. 4). As a reference, the performance meets the benchmark of BC simulation (NMB < ±20 % and NME < 45 %) proposed by Huang et al. (2021). Moreover, improved model performance was achieved for all the years, with largely reduced NMB and NME compared with simulations with the prior emissions (Fig. S6). The larger uncertainty for 2000 may be caused by limited observations and less controlled data quality for earlier years.

Figure 4 Correlation between simulated and observed monthly surface BC concentrations in China. Simulations were conducted based on prior and posterior BC emissions, while observations were collected from publications.

We also compare the simulated and observed BC concentrations by land-use type (Fig. S7). Application of the prior emissions resulted in a greater underestimation of BC concentrations for the forest and grassland regions (NMB: from −51 % to −77 %; NME: from 51 % to 77 %) compared with urban and rural regions (NMB: from −35 % to −51 %; NME: from 47 % to 51 %). The model performance was clearly improved for all the land-use types when the posterior emissions were applied. In particular, the NMB and NME were calculated to range, respectively, from −2 % to −40 % and from 31 % to 52 % for the forest and grassland regions. The evaluation supported our estimates of posterior emissions, not only for areas with insensitive human activities (e.g., urban regions) but also for remote regions.

Similarly, simulation of AAOD based on the prior emissions presented a clear underestimation compared with the OMI-derived observations, with the monthly NMB and NME ranging from −85 % to −29 % and from 34 % to 85 %, respectively (Table S12). As pointed out by Bond et al. (2013), the incorrect assumption of the mixing state of BC in the CTM could result in the general underestimation of MAE and thereby of AAOD. A clear improvement in the performance of AAOD simulation was found when the posterior emissions were applied. The NMB, NME, and RMSE were calculated to range from −11 % to 14 %, from 15 % to 28 %, and from 0.01 to 0.03, respectively, i.e., much smaller values than with the prior emissions, while R was largely elevated from 0.15–0.86 to 0.73–0.95.

Figure 5 Comparisons between posterior and various bottom-up BC emission estimates in China during 2000–2020. (a) Multiyear average spatial distribution of prior BC emissions (MEIC + GFED), (b) posterior BC emissions, and (c) their relative differences. (d) Long-term variability in the relative differences between posterior and various bottom-up BC emission estimates with 5-year intervals. Note that OBB emissions from GFED are added to each anthropogenic emission estimate as the total bottom-up estimate, except for PKU-Fuel, Lu et al. (2011), and Qin and Xie (2012) that include their own OBB emission estimate. (e) Long-term variability in normalized posterior and various bottom-up BC emission estimates (relative to 2000). The gray area indicates the period with declining national BC emissions. (f) Long-term variability in the relative differences between prior (MEIC + GFED) and posterior BC emissions by region and land-use type.​​​​​​​

3.2.2 Spatiotemporal patterns of posterior BC emissions and differences between prior and posterior BC emissions

Figure 5a–c present the spatial distribution of multiyear averages of the prior and posterior BC emissions and their relative differences (see details for individual years in Fig. 6). The annual total posterior BC emissions were estimated to increase from 6.48 Tg in 2000 to 7.03 Tg in 2010 and to decrease to 5.21 Tg in 2020. Compared with the prior BC emissions, the posterior BC emissions presented a clear enhancement, with a multiyear average factor (the ratio of the posterior emission difference to prior BC emissions, i.e., (posterior − prior) $/$ prior) of 3.3 for the entire country (Fig. 5c). The value declined from 3.7 to 2.7 during 2000–2015 but rose again to 4.1 in 2020 (Fig. 5d). The posterior BC emissions also presented an enhancement compared with other bottom-up estimates of China's BC emissions (sum of anthropogenic and OBB emissions), with the lowest factor of 1.7 for PKU-Fuel (Wang et al., 2014) and the highest factor of 4.1 for EDGAR + GFED (https://edgar.jrc.ec.europa.eu/dataset_ap61, last access: 10 May 2023) (Fig. 5d and Table 1). The comparisons between bottom-up and top-down estimates of BC emissions suggested a possible underestimation of the former, resulting partly from the under-reporting of activity levels and lack of local measurements for specific BC emission factors (EFs; emissions per unit of activity level) (Fu et al., 2012; Guan et al., 2012). In addition, the omission of small fires from satellite observations and the application of global EFs led to an underestimation of biomass burning emissions (Yang and Zhao, 2019). Along with improved energy and economic statistics and the increased number of EF data obtained through field observations, the discrepancy between prior and posterior BC emissions was gradually reduced until 2020. The increased uncertainty in prior BC emissions in 2020 may have resulted partly from an underestimation of increased fuel use owing to residential heating and cooking during the COVID-19 lockdown and quarantine period (Zheng et al., 2020).

Figure 6 The annual (a) prior BC emissions, (b) posterior BC emissions (3 times total emissions of January, April, July, and October), and (c) their relative differences during 2000–2020 (with a 5-year interval).

Table 1 The multiyear average relative differences between the posterior and various bottom-up estimates of BC emissions by region (unitless). Note that OBB emissions from GFED are added to each anthropogenic emission estimate as the total bottom-up estimate, except for PKU-Fuel which includes its own OBB emission estimate.

The posterior emissions presented a smaller interannual variability compared with the prior and other bottom-up estimates, with a net growth of 8 % during 2000–2010 (the analogous numbers are 12 %–55 % for various bottom-up estimates including 24 % for the prior used in this work, MEIC + GFED) and a decline of 26 % during 2010–2020 (41 % for MEIC + GFED; Fig. 5e). Besides residential sources, prior emission estimate bias may also have occurred in the transportation sector, such as extra emissions derived from inadequately eliminated vehicles with relatively old standards and the use of specific after-treatment technologies (e.g., diesel particulate filters) causing the release of ultrafine particles (Louis et al., 2016).

Relatively smaller differences between posterior and prior BC emissions were found in eastern China, with a multiyear average ratio of posterior to prior BC emissions estimated at 2.4, while that for the rest of China reached 8.0 (Fig. 5f). The relative differences between posterior and other bottom-up BC emission estimates were smaller (1.1–2.1) in more economically developed regions (BTH, FWP, YRD, and PRD) but larger (3.5–5.6) in SCB, NE, and other regions (Table 1). To further explore the impact of human activities on BC emissions, we divided the country into different land-use types (Fig. S1b). The multiyear average BC emission intensity in urban areas was estimated at 1.86 Mg km⁻² yr⁻¹, i.e., higher than 1.47 Mg km⁻² yr⁻¹ in rural areas (Table 2). Industrial production, transportation, and commercial activities generate abundant emissions in urban areas, while straw burning, residential cooking, and heating are important emission sources in rural areas. Smaller BC emission intensity was estimated in regions less influenced by human activities, i.e., 0.84, 0.33, and 0.11 Mg km⁻² yr⁻¹ for forest, grassland, and unused regions, respectively. As shown in Fig. 5f, the relative differences between posterior and prior BC emissions were the smallest in urban areas, with a multiyear average enhancement factor of 1.6, followed by rural areas (a factor of 2.6), forest (a factor of 4.3), and grassland (a factor of 5.5). In general, the bottom-up approach could capture information about energy consumption and pollution controls more easily and accurately in regions with more intensive human activities. Such advantages helped reduce the uncertainty in emission estimates for developed urban areas compared with the uncertainty in remote areas. However, current official statistics do not sufficiently report biofuel consumption and are believed to greatly underestimate raw coal consumption in rural areas (Zhi et al., 2017; Zhu et al., 2019). Limited small-fire detection ability via satellite also led to an underestimation of OBB in forest and grassland areas (Schroeder et al., 2008; Yang and Zhao, 2019). Such limitations resulted in a greater underestimation of prior BC emissions in rural and remote regions compared with urban regions (Fig. S7).

Table 2 The annual posterior BC emission intensity of different land-use types during 2000–2020 (with a 5-year interval). The annual emission intensity was estimated as 3 times the sum of BC emission intensity of January, April, July, and October. “Urban” includes city and building categories, and “Rural” includes cropland and countryside categories (unit: Mg km⁻² yr⁻¹).

Previous limited studies employing top-down approaches estimated China's BC emissions at 5.7, 3.1, and 2.5 Tg yr⁻¹ for 2000, 2006, and 2008 (Cohen and Wang, 2014; Fu et al., 2012; Wang et al., 2016), respectively. These studies presented a 0.6–5 times enhancement of posterior BC emissions compared with prior BC emissions (C. Chen et al., 2019; Cohen and Wang, 2014; Fu et al., 2012; Wang et al., 2016; Zhang et al., 2015), showing a wider range than the results presented in the current study (2.7–4.1). This discrepancy may have resulted from the application of different inversion approaches, observational data, and prior BC emissions. Moreover, the posterior BC emissions applied in previous studies were lower than those used in the current study for corresponding years by a factor of 12 %–64 %. In general, incomplete spatial coverage in previous studies limited the emission inversion capability in regions far from the observational sites.

3.2.3 Influence of social and economic development on BC emissions

To explore the influence of social and economic development on BC emissions, we analyzed the diverse changing patterns in posterior BC emission intensity (Mg km⁻²) by province, based on the rural population fraction, provincial proportion of coal production to the national total, and industrial gross domestic product (GDP).

Figure 7 Changes in provincial BC emission intensity (annual BC emissions per km²) in posterior BC estimates. (a, b) Spatial distribution of the interannual change rate of BC emission intensity by province during 2000–2010 and 2010–2020. (c, d) Relationships between interannual BC emission intensity change rate and rural population fraction, proportion of coal production, and industrial GDP for each province. The x and y axes present the rural population fraction of each province and the provincial proportion of coal production to the national total for the middle year of the period concerned (i.e., 2005 for (c) and 2015 for (d)), respectively. Circle size represents provincial industrial GDP level. The colors in (c) and (d) are the same as those in (a) and (b). Population, coal production, and GDP statistics were obtained from the National Bureau of Statistics (https://data.stats.gov.cn/english/).

Shanxi, Inner Mongolia, Henan, and Shaanxi were identified as the main coal-producing provinces, collectively contributing 50 % and 68 % to the national total coal production in China during 2000–2010 and 2010–2020, respectively. Shanxi, Henan, and Shaanxi experienced a decline in posterior BC emission intensity during 2000–2010, with annual average decreasing rates estimated at 0.13, 0.03, and 0.02 Mg km⁻² yr⁻¹, respectively (Fig. 7a). These provinces have long suffered air pollution from coal burning, and the coal industry's structure and technology were improved earlier than in other provinces. For example, Shanxi eliminated over 7000 small coal mines and an outdated production capacity of 3.85×10⁸ t of coal during 2000–2010 (The Central People's Government of the People's Republic of China, 2011; Han and Wang, 2015). These measures significantly improved coal consumption efficiency, resulting in a sharp decline in BC emission intensity.

In comparison, a slower decline or even growth in BC emission intensity was demonstrated in Shanxi and Henan during 2010–2020, with changing rates of −0.007 and 0.004 Mg km⁻² yr⁻¹, respectively (Fig. 7b). Merging and reorganization of the coal industry in these two provinces have enhanced coal production and consumption in recent years, thereby reducing the benefits of BC emission controls. Comparatively faster declines were demonstrated in Inner Mongolia and Shaanxi, with decreasing rates of 0.01 and 0.03 Mg km⁻² yr⁻¹, respectively (Fig. 7b). The proportion of coal production in Inner Mongolia and Shaanxi increased from 8 % and 2 % in 2000 to 28 % and 17 % in 2020, respectively, demonstrating the increasingly important role of national coal production in these two provinces (Fig. 7c and d). Compared with Shanxi and Henan, relatively later but greater efforts were made to improve the coal industry's structure in Inner Mongolia and Shaanxi, leading to considerable BC emission reductions after 2010.

Increasing emission intensities were found for most other provinces during 2000–2010, particularly those with high industrial GDP (larger circles in Fig. 7c). Prominent BC emission intensity growth was demonstrated in Jiangsu, Shandong, and Beijing, with the annual average growth rates ranging from 0.05 to 0.10 Mg km⁻² yr⁻¹. Intensive industrial production activities in these provinces/cities resulted in quickly increasing emissions. Second to industry-developed regions, Anhui, Guangxi, and Yunnan experienced BC emission growth rates of 0.03–0.05 Mg km⁻² yr⁻¹. In these provinces with higher rural population fractions, the enhanced consumption of household solid fuel and limited progress in air pollution controls led to fast BC emission growth, along with an increased demand for living standard improvements.

During 2010–2020, BC emission intensities were estimated to decline for most provinces (Fig. 7b), and faster declines commonly occurred in provinces with higher urbanization and industrial GDP levels (darker blue and larger circles from right to left on the x axis in Fig. 7d). The greatest reductions were demonstrated in Shanghai, Liaoning, Chongqing, Jiangsu, and Fujian, with the decreasing rate ranging from 0.07 to 0.23 Mg km⁻² yr⁻¹. During this period, stringent pollution controls in the industrial and transportation sectors took effect, particularly in economically developed and highly urbanized regions, resulting in a faster decline in BC emissions compared with less developed regions. Moreover, the reduced rural population proportion and increased clean energy use jointly restrained household BC emissions.

3.2.4 Comparison of emission and concentration trends for multiple species

We compare the interannual changes in posterior BC emissions with those in national PM_2.5, OC, and BC emissions (i.e., prior BC emissions) derived from MEIC, PM_2.5 concentrations derived from Tracking Air Pollution in China (http://tapdata.org.cn/, last access: 31 January 2023), and CMAQ-simulated BC concentrations based on posterior BC emissions (Fig. 8). During 2005–2020, the annual BC emissions were estimated to decline by 26 % (posterior BC emissions) or 43 % (prior BC emissions), which was slower than PM_2.5 (56 %), and the relative reduction in BC concentrations (14 %) was less than that in PM_2.5 (35 %). Compared with total PM_2.5, for which the health effects are widely recognized, more attention should be paid to the health effects and control of BC emissions, given their relatively slower decline in ambient concentrations and well-acknowledged higher health risks (Wang et al., 2021; see our estimate on the mortality attributable to BC exposure in Sect. 3.3). Moreover, the comparison between emission trends in warming (BC) and cooling (e.g., OC) species reveals a climate challenge. A faster decline was estimated in OC emissions (47 % during 2005–2020) than in BC emissions, resulting mainly from greatly reduced biofuel use and biomass burning. By contrast, development of the transportation and industrial sectors makes further reductions in BC emissions challenging, and more effective strategies for BC emission controls are urgently needed to restrain climate warming in the future.

Figure 8 Relative changes in the prior and posterior BC, OC, and PM_2.5 emissions (left vertical axis) and the posterior BC and PM_2.5 concentrations (right vertical axis) during 2000–2020 compared with those of 2000 (with a 5-year interval). OC and PM_2.5 emissions are obtained from the Multi-resolution Emission Inventory for China (MEIC, http://www.meicmodel.org). PM_2.5 concentrations are obtained from Tracking Air Pollution in China (TAP, http://tapdata.org.cn/, last access: 31 January 2023).

3.3 Mortality attributable to BC exposure and its drivers

The all-cause premature deaths attributed to BC increased from 733 910 cases (95 % CI: 676 790–800 250) in 2000 to 903 030 cases (832 830–984 530) in 2005, decreased to 857 510 cases (790 500–935 370) in 2015, and finally reached the highest level at 937 980 cases (864 510–1 023 400) in 2020 (Fig. 9). All-cause premature deaths in China attributed to BC exposure were estimated by previous studies to be within the wide range of 50 100–1 436 960 cases, with different BC exposure levels and β_BC values used (Cui et al., 2022; Qin et al., 2019; Saikawa et al., 2009; Wang et al., 2021). More premature deaths in eastern China (Fig. S8) were attributed mainly to the relatively high population density and BC exposure from developed industrial and commercial activities (Fig. S9). The highest multiyear average of premature mortality was 1482 cases per 1000 km² (the all-cause premature deaths attributed to BC per area of 1000 km²) in Shanghai, followed by 793, 761, 520, 450, and 442 cases per 1000 km² in Beijing, Tianjin, Jiangsu, Henan, and Shandong, respectively (Table 3). These values were much higher than the national average of 86 cases per 1000 km².

Figure 9 Total all-cause premature deaths attributed to BC exposure in China (gray bars) and drivers of changing premature mortality (colored bars) during 2000–2020 (with a 5-year interval). Error bars show the 95 % confidence interval of estimates in this study. The numbers of total all-cause premature deaths are at the bottom of the figure. The contributions of major factors to the national changing mortality are above the colored bars.

Table 3 The annual all-cause premature mortality associated with BC exposure by province in mainland China during 2000–2020 with a 5-year interval (unit: cases per 1000 km²). Locations of provinces are shown in Fig. S1a.

Also shown in Fig. 9 are the contributions of major factors to the national changing mortality during 2000–2020 (provincial-level results are presented in Table S13). The changing emission levels played an important role in premature mortality. Along with swift growth in the economy and energy consumption from 2000 to 2010, increasing BC emissions enhanced health risks in China overall and in most provinces. The largest increase in national annual mortality was demonstrated during 2000–2005 (78 590 cases, 95 % CI: 72 520–85 600). However, BC emission reductions in the main coal-producing provinces, i.e., Shanxi, Henan, and Sichuan, led to a decrease in regional health risks during 2005–2010, with declines in annual mortality of 14 320, 13 580, and 8410 cases, respectively. Benefiting from improved air pollution controls from 2010 to 2020, declining health risks associated with BC emission reductions were demonstrated in China overall and in most provinces. The largest decline in national annual mortality associated with BC emission reductions was demonstrated during 2010–2015 (133 360 cases, 123 150–145 180). BC emission rebound in Shanxi and Henan led to elevated regional health risks during 2015–2020, with increases in annual mortality of 7170 and 6190 cases, respectively.

Varying meteorological conditions also affected the health burden. With the exception of 2005–2010, meteorological conditions were estimated to increase BC exposure and the associated mortality for most of the research period, particularly in eastern China (Fig. S10). Among vulnerability factors, population aging contributed most to the increased BC-related health burden, elevating the annual national mortality by 99 800–142 310 cases for the various 5-year intervals. Population growth contributed modestly to the increased health burden in most provinces, with annual national total changes of 10 940–21 230 cases. Exceptions included Sichuan and Hubei, where reduced populations resulted in declining mature mortality. Improved healthcare partly offset the adverse effect, with the annual avoided deaths increasing from 70 100 cases during 2000–2005 to 120 440 cases during 2010–2015. This positive effect shrank to 56 690 cases during 2015–2020. With increased population aging and reduced potential of medical care improvement, greater BC emission abatement will be needed to further prevent health damage.

3.4 Uncertainty analysis

The uncertainties in our results are mainly attributed to (a) separating the contribution of BC to AAOD, (b) simulating BC AAOD based on WRF-CMAQ and the empirical light absorption model, (c) lacking diurnal information of observed BC AAOD, and (d) estimating the premature mortality attributable to BC. Due to the limited number of studies reporting spatiotemporal variability in BrC and dust light absorption in China, we separated BC AAOD based on the MERRA-2 dataset, which may underestimate BrC light absorption at 483 nm (Buchard et al., 2017). Herein, the multiyear average BrC share in AAOD was estimated at 16 % (MERRA-2) in eastern China (Fig. S11), lower than that based on observations (18 %–44 %; S. Chen et al., 2019; Li et al., 2019; Zhu et al., 2021). Notably, the multiyear average of posterior BC emissions using AAOD at a longer wavelength (865 nm) with little BrC effect was estimated to be 11 % lower than that using 483 nm (Test 1 in Text S3 and Fig. 10). The posterior BC emissions based on the dust light absorption fraction in AAOD from Copernicus Atmosphere Monitoring Service (CAMS), which was only half that of MERRA-2 (Fig. S3), were estimated to be very close to those of Test 1 (Test 2 in Text S3 and Fig. 10). In addition, we adopted an empirical model based on observations from one city to simulate MAE. Even with improved performance compared with the average level of multiple CTMs at the global scale (Gliß et al., 2021), our model underestimated MAE and presented smaller regional heterogeneity than existing observations, most likely due to the limited spatial extrapolation ability. When MAE was modified according to available observations across the country, the multiyear average of posterior BC emissions was estimated to be 10 % lower than that without MAE modification (Test 3 in Text S3 and Fig. 10). Besides, the removal processes in CMAQ also affect the BC lifetime simulation and thereby its atmospheric column concentration and emission inversion. The simulated lifetime of 4.7 d in the base case is within the range of 3.8–11.4 d reported by previous studies (Fig. S4; Bond et al., 2013; Vignati et al., 2010). By adjusting the simulated BC lifetime to the multi-model average level (5.5 d; Gliß et al., 2021), the posterior emissions were estimated to be 4 % smaller than those of the base case (Text S3 and Fig. 10). Although modest uncertainties were revealed by the above sensitivity tests, they did not change the main findings of this study, with similar interannual variabilities between the base and sensitivity test cases (Fig. 10). These uncertainties should be reduced with improved spatiotemporal coverage of BC, BrC, and dust light absorption observations. Another limitation of the inversion system was the missing diurnal information of observed AAOD. Since OMI could only provide AAOD at 13:30 LT (local time)​​​​​​​, we could not capture the diurnal distribution of BC emissions. The Geostationary Environment Monitoring Spectrometer (GEMS; Kim et al., 2020) was launched in 2020 and provided hourly daytime observations of aerosols. This can potentially be helpful for improving the temporal accuracy of BC emission inversion in the future.

Figure 10 (a) Long-term variability in the normalized prior and posterior BC emissions in China during 2000–2020 (compared with 2000, solid lines with left vertical axis, with a 5-year interval) and relative difference between the posterior and prior emissions (dashed lines with right vertical axis) for the base case, Test 1, Test 2, Test 3, and Test 4. (b) Relative difference in eastern China and the rest of China for the base case (red lines and marks), Test 1 (green lines and marks), Test 2 (purple lines and marks), Test 3 (blue lines and marks), and Test 4 (yellow lines and marks).

Furthermore, the health impact estimation could be biased by rare domestic β_BC values in China. Previous studies commonly adopted the same functions as PM_2.5 (Saikawa et al., 2009) or β_BC values obtained from American or European studies (Wang et al., 2021), resulting in large uncertainty. Herein, we relied on a unique cohort study in China and calculated the all-cause premature deaths attributed to BC at 733 910–937 980 per year. The β_BC values obtained from national-scale studies in the US and Europe indicate a 10-fold difference (220 980–2 386 060 per year, Table S14), similar to the estimation conducted in the US (Y. Li et al., 2016). More domestic epidemiological studies focusing on BC emissions are expected to further reduce the uncertainty.

4 Concluding remarks

Compared with previous studies with nonconsecutive or incomparable estimates, this study provides a panoramic view of the spatiotemporal patterns of AAOD, BC emissions, and the associated mortality in China for the past 2 decades. We found that the bottom-up approach most likely underestimated BC emissions, particularly in less developed western and remote areas. Our findings also reveal the influence of human activities on the evolution of BC emissions and the remarkable emission abatement resulting from the implementation of national pollution controls, particularly in developed regions. Pollution controls were estimated to reduce the annual BC emissions by 26 % during 2010–2020, reversing the 8 % growth during 2000–2010. However, the benefits were smaller than those previously estimated employing the bottom-up approach, which most likely overestimated progress in pollution controls for certain sources, such as the transportation sector and residential solid fuel burning. The long-term BC emission trends in this study address both health and climate risks combined with the effects of other short-lived aerosol species.

The energy transition path to achieve China's goal of peak emissions and carbon neutrality provides an opportunity to further reduce BC emissions. Compared with developed regions, the energy transition and emission abatement is more challenging in coal-producing and less urbanized regions; thus, region-specific emission controls should be formulated. For the main coal-producing provinces, BC emissions have declined much slower than those in economically developed provinces or even rebounded along with increased industrial production capacity and energy demand in recent years. As China's traditional energy base, these provinces need to accelerate energy infrastructure adjustments and reduce their dependence on coal through, for example, the development of photovoltaic and wind power. In addition, aggressively promoting advanced manufacturing technology is recommended. For example, expanding the coal chemical industry chain could shift the role of coal consumption from traditional fuel to raw material, thus achieving its clean utilization (e.g., coal liquefaction and gasification technique). For less urbanized regions, solid fuel, including coal, firewood, and crop residues, remains the major energy source, and actual BC emissions in rural areas could be greatly underestimated compared with those in urban areas. Expansion of natural gas and electricity use for cooking and heating could effectively limit BC emissions in these regions. These efforts can be supported through better infrastructure development and subsidy policy design in the future.