Black carbon aerosols in China: spatial-temporal variations and lessons from long-term atmospheric observations

Black carbon (BC) significantly influences climate, air quality, and public health, and long-term observations are essential for understanding its adverse effects. While previous studies have primarily focused on spatiotemporal variations, deeper insights from such datasets remain uncovered. Using 13 years (2008–2020) of continuous measurements of equivalent black carbon (eBC) in China, this study reported the spatial-temporal variations of eBC and its sources, including solid fuel (eBC_sf) and liquid fuel combustion (eBC_lf). The results showed that eBC and its sources exhibited higher concentrations in eastern and northern China compared to western and southern China. Seasonal variations of eBC and eBC_sf generally showed lower values during summer and higher values during winter at most stations. Long-term trends indicated that eBC and eBC_lf decreased most rapidly at urban stations, while eBC_sf declined faster at rural stations. Comparisons of eBC concentrations and trends between this study and global observations revealed higher eBC levels but lower reduction rates in China. These long-term observations showed that the model simulations performed well in simulating spatial distribution but poorly in capturing inter-annual variations. The weather-normalized eBC concentrations showed potential for adjusting emission estimates. The normalized results also suggested that emission control was the dominant driver of the BC reduction. This decrease was primarily driven by reductions from solid fuel combustion at rural and background stations. This study provides insights for reducing uncertainties in black carbon emission inventories and improving model performance in simulating surface concentrations.

1 Introduction

Despite being a small fraction of the aerosol chemical composition (Bond et al., 2013; Huang et al., 2014; Tao et al., 2017), black carbon (BC) is important for climate, air pollution, and human health. Globally, BC can cause up to +0.41 W m⁻² effective radiative forcing (Szopa et al., 2021). At the regional scale, BC aerosol can suppress the deepening of the planetary boundary layer, thereby enhancing haze pollution (Ding et al., 2016; Wang et al., 2018) and radiative forcing (Peng et al., 2016). Personal exposure to BC can cause adverse health effects, such as carcinogenic risks and elevated blood pressure, and BC can reach the fetal side of the human placenta (Mordukhovich et al., 2009; Bové et al., 2019; Lin et al., 2019). As a short-lived climate forcer, the atmospheric lifetime of BC is about 5.5 d ± 35 % (median ± 1 standard deviation, hereafter), according to an ensemble of 14 models (Gliß et al., 2021). Therefore, the study of BC aerosol is important to reduce its adverse impacts.

To better understand the climate, environmental, and health effects of BC, it is important to model its atmospheric abundance. Although Chemistry Transport Models (CTMs) can capture the spatial-temporal variations of BC, their simulation accuracy is often poor, and there are discrepancies among different simulations. For instance, 17 global models overestimated the abundance of BC, with the ratio of simulated to observed values ranging from 0.5 to 10 in Europe (Koch et al., 2009). Similarly, Ikeda et al. (2022) evaluated the performance of six BC emission inventories and found an overestimation of modeled BC concentrations by a factor of 1.24 to 2.16 in China. In contrast, some simulations underestimated the annual mean BC levels in China, with normalized mean bias ranging from −83 % to −22 % (Mao et al., 2016). The high uncertainties in the simulation of BC arise from the emission inventory and physical mechanisms. One method to constrain the model results is to use long-term observation datasets to evaluate the emission inventory, transport, and removal processes (Bond et al., 2013; Wang et al., 2014; Alvarado et al., 2016; He et al., 2016; Evangeliou et al., 2018). For instance, Evangeliou et al. (2018) coupled BC observations with a top-down method to constrain emissions, and the posterior simulated BC concentrations showed better agreement with observations. Choi et al. (2020b) utilized long-term observations of BC to CO ratios to validate the regional emission inventory in Asia. They found that emissions in East and South Korea were overestimated, while emissions in North Korea were underestimated. Using long-term observations of BC and CO from 2009 to 2015 at Fukue Island, Kanaya et al. (2016) constrained the wet removal rate of BC and proposed a stretched exponential decay equation to describe the wet removal rate. Based on BC observations from 2010 to 2016 at three representative background sites in East Asia, Choi et al. (2020a) estimated the average transport efficiency of BC to be 0.73, which was lower than the mean rate of 0.91 from the FLEXible PARTicle Lagrangian transport model (FLEXPART, version 10.4).

Long-term observation of BC aerosol also contributes to a better understanding of the effectiveness of air quality management (Chen et al., 2016; Font and Fuller, 2016; Fuller and Font, 2019). For instance, a decreasing trend in BC from 2010 to 2014 in London was attributed to the application of diesel particle filters (Font and Fuller, 2016). Long-term observation of BC at 16 sites from 2009 to 2018 in Germany showed the greatest decreases at roadside and urban background sites, implying strong evidence of reduced traffic emissions in urban areas (Sun et al., 2020). The decreasing trend of BC at traffic sites in the Helsinki metropolitan area suggested the success of vehicle exhaust particle mitigation (Luoma et al., 2021). In China, similar conclusions were drawn. For instance, a 38 % reduction in BC concentration was observed in Beijing between 2005 and 2013, primarily due to the relocation of industries and a decrease in coal and coke consumption (Chen et al., 2016). Continuous observations of BC in Beijing from 2012 to 2020 revealed a 71 % reduction, with the largest decrease occurring at night. This suggests that pollution control measures implemented since 2013 have effectively reduced primary emissions (Sun et al., 2022). The long-term trend of BC in Wuhan suggested a positive role for pollution control measures across China in reducing BC at the observational site (Zheng et al., 2020). Continuous observation of BC at Mt. Waliguan showed a decreasing rate of 2 % yr⁻¹ from 2008 to 2017, reflecting the reduction of BC emissions in China (Dai et al., 2021).

As summarized above, long-term observation of BC is important for improving CTM simulations and evaluating the effectiveness of air pollution control measures on air pollutants. Long-term surface observations of BC aerosol have been widely reported (Hirdman et al., 2010; Boreddy et al., 2018; Kutzner et al., 2018). For instance, the earliest continuous in situ observation of BC was launched at 150 urban and rural stations in 1962 in the United Kingdom (Novakov and Hansen, 2004). Compared to long-term observations at a single or a small number of stations, a monitoring network consisting of many stations can provide more information about the spatial-temporal variations of air pollutants. However, long-term in situ observation networks are few. At present, the worldwide monitoring networks include the World Meteorological Organization Global Atmosphere Watch (WMO-GAW) (Bond et al., 2013), European Monitoring and Evaluation Programme (EMEP) (Yttri et al., 2007), German Ultrafine Aerosol Network (GUAN) (Birmili et al., 2016), UK Black Carbon Network (Singh et al., 2018), National Air Pollution Monitoring Network (NABEL) in Switzerland (Grange et al., 2020), Interagency Monitoring of Protected Visual Environments network (IMPROVE) in the US (Malm et al., 1994), Canadian Aerosol Baseline Measurement (CABM) (Chan et al., 2019), Aerosol Radiative Forcing over India Network (ARFINET) in India (Manoj et al., 2019), and China Black Carbon Observational Network (CBNET) (Zhang et al., 2019b). Research from these observational networks is mainly focused on the spatial-temporal variations of BC (Birmili et al., 2016; Manoj et al., 2019; Zhang et al., 2019b; Savadkoohi et al., 2023). More valuable information from long-term observations, however, remains untapped.

Therefore, this study presents a comprehensive BC dataset collected from CBNET covering the period from 2008 to 2020. Spatial-temporal characteristics, sources, and long-term trends of BC are reported. The information behind the results is explored with a focus on three questions: (1) What can we learn from the comparison between BC ground observations and CTM simulations? (2) Can the inter-annual variations of BC be used as an indicator of BC emissions? (3) Which factors dominated the variations of BC in China during the past 13 years? By analyzing the data from continuous measurements of BC in China, this study is expected to provide insights into reducing the uncertainties in BC emissions and simulations from CTMs.

2 Methodology

2.1 Observation sites

The CBNET, established in 2006, is a part of the China Atmosphere Watch Network (CAWNET). There are 68 stations across China, and the number of stations has increased recently. In this study, 48 sites with observations during 2008–2020 were selected (Fig. S1). According to the land type and surrounding environment, these 48 sites can be divided into baseline (N=7), urban (N=23), and rural (N=18) stations (Guo et al., 2020; Zhang et al., 2008, 2019b). The rural stations are situated around 50–100 km distant from local sources or nearby major cities, and they are positioned at a moderate height above the local elevation of the area. For the urban stations, the sampling heights are deliberately set to be 50–100 m higher than the average elevation of the city. Baseline stations are established in areas remote from strong emission sources and human activities to monitor the long-range transport of BC aerosols and their natural emissions. For instance, the Lin'an (LA) baseline station is situated at an altitude of 139 m a.s.l. in Zhejiang province, positioned approximately 150 km northeast of Shanghai and 50 km west of Hangzhou city. This strategic placement aims to collect samples that represent the entire region rather than just the immediate local surroundings (Zhang et al., 2008). More details about the selected stations can be found in Table S1.

2.2 Instruments and data prepossessing

The instruments used at each monitoring site are provided in Table S1 in the Supplement. The BC mass concentration was measured by filter-based absorption Aethalometers. Most stations utilized the AE31 to measure BC throughout the entire period (i.e., 2008–2020), and some stations used both the AE31 and AE33 (Table S1). Details about the principles, operations, and algorithms of AE31 and AE33 can be found elsewhere (Drinovec et al., 2015; Rajesh and Ramachandran, 2018). Briefly, ambient aerosol is drawn through a cyclone with a cut-off size of 2.5 µm at a flow rate of 2–5 L min⁻¹ and deposited on a quartz fiber tape (AE31) or tetrafluoroethylene-coated glass filter (AE33). Changes in attenuation before and after sampling are measured and converted into aerosol absorption coefficients (b_abs), then into mass concentrations using the pre-defined mass attenuation cross-section (for AE31) or mass absorption cross-sections (for AE33) at 7 channels (370, 470, 520, 590, 660, 880, 950 nm). It should be noted that we adopted the acronym MAC throughout this paper to collectively denote both the mass attenuation cross-section and mass absorption cross-section. The BC mass concentration reported by the Aethalometer is referred to as equivalent black carbon (eBC) (Petzold et al., 2013; Savadkoohi et al., 2024). However, due to multiple scattering and loading effects, the reported eBC mass concentration differs from that in ambient air. To correct the loading and multiple scattering effects, the parameterization scheme developed by Weingartner et al. (2003) was used for b_abs reported by AE31. For AE33, the filter loading effect is compensated by a new real-time loading effect compensation algorithm, and the multiple-scattering effect is corrected by a multiple-scattering enhancement correction (C₀) (Drinovec et al., 2015). The C₀ values differ between AE31 (C₀=2.14) and AE33 (C₀=1.39) due to differences in tape materials (quartz fiber for AE31 and tetrafluoroethylene-coated glass filter for AE33). Another factor that influences eBC mass concentration is the MAC, which is determined by particle size, the mixing state of aerosol components, and the morphology of BC particles(Bond and Bergstrom, 2006; Bond et al., 2013; Yuan et al., 2021), and it shows spatial-temporal variations (Pandolfi et al., 2014; Zanatta et al., 2016; Grange et al., 2020; Guo et al., 2020; Savadkoohi et al., 2024). Therefore, using the nominal MAC from default settings (e.g., 16.6 m² g⁻¹ for AE31 and 7.77 m² g⁻¹ for AE33 at 880 nm) without considering MAC spatial-temporal variability can result in misleading eBC estimations (Savadkoohi et al., 2024). To minimize the impact of MAC on eBC mass concentration, station-specific MAC values (Table S1) from previous studies (Zhang et al., 2008; Guo et al., 2020; Zheng et al., 2021), which utilized synchronous measurements of b_abs and elemental carbon (EC) to derive MAC values, were used.

Before data analysis, the raw data (5 min for AE31 and 1 min for AE33) acquired from the Aethalometers was checked and the outliers were removed according to the following procedures: (1) negative values of attenuation at time t and the next observation time (e.g., t+5 min) were removed at all 7 wavelengths; (2) the observations which failed to comply with the principle that the light absorption decreases with increasing wavelength were removed and hourly mean values of remaining observations were calculated; (3) the outliers of the hourly time series during the entire observation period were removed using the seasonal decomposition algorithm (Dancho and Vaughan, 2023). It should be noted that our method to process the negative value differed from the manufacturer's protocol, which recommends retaining negative values and applying smoothing techniques to exclude anomalous data points. In our method, the negative values were directly removed. This is because the Aethalometer employs the time-differential method to measure the optical attenuation (ATN) at the sampling point and subsequently converts ATN into BC mass concentration. The time-differential method implies that if the difference between the ATN in the current observation period (ATN_t) and that of the previous observation period (e.g., ATN_t−5 min) is negative, the calculated BC mass concentration for period t will be negative. Furthermore, the difference between the ATN in the next observation period (e.g., ATN_t+5 min) and that at time t (ATN_t) will be overestimated, consequently resulting in an overestimation of the BC concentration at time t+5 min. The absorption Ångström exponent of aerosol was calculated using the power law fitting between b_abs and wavelengths at seven wavelengths and reported as AAE_370–950. Observations with a p-value for fitting higher than 0.001 were considered as outliers and excluded from further analysis.

After data quality control, the data availability in each station of each year was calculated (Fig. S2). To get robust spatial-temporal variations of BC, the stations with annual data records exceeding 50 % were used. The 50 % threshold was selected to maximize the inclusion of monitoring stations while ensuring the collected data effectively captures the temporal variations of BC. To ensure data robustness, systematic screening criteria were applied: (1) long-term trend analysis required stations with more than 8 years of eBC observations, each year maintaining >50 % data availability; (2) spatial-temporal pattern analysis included stations with more than 2 years of valid observations during 2015 and 2017, where annual data completeness exceeded 50 %. According to these rules, 34 and 25 stations were selected for spatial-temporal and long-term trend analysis, respectively (Fig. 1). It should be noted that data from other stations with data availability higher than 50 % were also used (e.g., Sect. 4.1 for validation of the BC model simulations).

Figure 1 Spatial distribution of ground-level black carbon monitoring stations. The dots with white circular frame and white crosses represent the stations used for spatial-temporal variation analysis between 2015 and 2017 (N=34) and long-term analysis from 2008 to 2020 (N=25), respectively. The colors of the dots indicate different types of stations. The black diagonal line represents the “Hu-Line”, which separates China into eastern and western parts. The red polyline represents the “Qinling-Huaihe Line”, which separates China into northern and southern parts. The full names of acronyms can be found in Table S1.

2.3 Aethalometer model and eBC source apportionment

BC is formed via the combustion processes of carbonaceous matter and is mainly from residential, power, industrial, and transport sectors (Bond et al., 2013; Li et al., 2017; McDuffie et al., 2020). To apportion the sources of BC, an Aethalometer model using aerosol absorption measurements at two wavelengths was developed by Sandradewi et al. (2008). The Aethalometer model operates on the assumption that light-absorbing aerosols originate predominantly from two source types (Sandradewi et al., 2008), e.g., liquid fuel combustion and solid fuel combustion, which are characterized by distinct Absorption Ångström Exponent (AAE) values. Aerosols derived from solid fuel combustion (AAE_sf) generally exhibit higher and more variable AAE values compared to those from liquid fuel combustion (AAE_lf). The elevated AAE_sf values are attributed to the presence of significant amounts of organic carbon species, which absorb strongly in the ultraviolet and lower visible wavelength ranges (Sandradewi et al., 2008). In contrast, the lower AAE_lf values result from the dominant contribution of BC, which absorbs light broadly across the visible spectrum. These source-specific AAE characteristics are consistent with laboratory combustion experiments, which also report higher AAE values for solid fuel emissions and lower values for liquid fuel emissions (Olson et al., 2015). Consequently, measured ambient AAE values can serve as an indicator of dominant BC sources, e.g., lower AAE values are typically associated with environments strongly influenced by liquid fuel combustion, such as road tunnels (Blanco-Alegre et al., 2020).

This method has been widely used in eBC source apportionment (Harrison et al., 2012; Zheng et al., 2019, 2020, 2021; Savadkoohi et al., 2023), and the results of this method are generally consistent with other methods, e.g., receptor models (Favez et al., 2010; Herich et al., 2011; Liu et al., 2014). In this study, eBC was apportioned into solid fuel combustion (eBC_sf) and liquid fossil fuel combustion (eBC_lf) rather than biomass burning (eBC_bb) and fossil fuel combustion (eBC_ff) as reported by previous studies conducted in China (Zheng et al., 2020; Wang et al., 2023; Xu et al., 2024; Fan et al., 2025). The rationale for differentiating eBC into solid fuel (e.g., coal and biomass) combustion and liquid fossil fuel combustion (e.g., vehicle emissions) was based on the following reasons. (1) Bottom-up BC emission inventories showed that BC from coal, solid bio-fuel, and liquid fossil combustion contributed 28.6 %, 19.6 %, and 28.0 %, respectively, to total emissions (Fig. S3a). The discrepancy between the fraction of solid fuel combustion emissions and the attribution of eBC to eBC_sf results (48.2 % vs. 48.2 ± 7.14 %) was smaller than that between biomass burning emissions and the attribution of eBC to eBC_bb results (19.6 % vs. 48.2 ± 7.14 %). (2) Based on experimental measurements of the AAE conducted on various aerosol types (Fig. S3b), the discrepancy between AAE values for aerosols from biomass burning and coal combustion (1.76 ± 0.40 vs. 1.47 ± 0.13) was smaller than that between liquid fossil fuel and solid fossil fuel combustion (1.00 ± 0.04 vs. 1.47 ± 0.13). Given these two reasons, eBC was apportioned into eBC_sf and eBC_lf in this study.

The two key issues of using the Aethalometer model are the selection of the wavelengths and the AAE for aerosol from liquid fossil fuel combustion and solid fuel combustion. According to Zotter et al. (2017), wavelength of 470 nm as a near-ultraviolet wavelength is recommended while the choice between 950 and 880 nm in the near-infrared is less critical. Following previous studies (Sandradewi et al., 2008; Zotter et al., 2017), wavelengths of 470 and 950 nm were chosen in this study. The selection of AAE_lf and AAE_sf, however, is more complicated. The AAE is used to describe the wavelength dependence of aerosol absorption. Although the aerosol AAE is impacted by particle size (Gyawali et al., 2009; Liu et al., 2018), chemical composition (Flowers et al., 2010; You et al., 2016), and morphology (Liu et al., 2008; Li et al., 2016), it can act as a proxy of aerosol composition (Yang et al., 2009; Russell et al., 2010). As summarized in previous studies (Helin et al., 2021; Zheng et al., 2021), aerosol from liquid fossil fuel combustion has an AAE value of ∼ 1.0, while aerosol from solid fuel combustion is characterized by a higher AAE value. Different methods were used to constrain the optimal AAE combinations, and it showed spatial heterogeneity (Fuller et al., 2014; Becerril-Valle et al., 2017; Zotter et al., 2017; Helin et al., 2018; Mousavi et al., 2019; Zheng et al., 2021). Due to the lack of auxiliary measurements (e.g., radiocarbon, levoglucosan, and potassium) to constrain the station-specific AAE combination, the measured aerosol AAE frequency distribution (Tobler et al., 2021; Rovira et al., 2025; Navarro-Barboza et al., 2025; Savadkoohi et al., 2025) was used to select the station-specific AAE_lf and AAE_sf values here. Specifically, the AAE_lf and AAE_sf were determined as the 1st and 99th percentiles of aerosol AAE, respectively. It should be noted that aerosol AAE was determined by power law fitting between light absorptions and wavelengths, and the fittings with the determinate coefficient (r²) less than 0.99 were excluded (Savadkoohi et al., 2025). By this method, the station-specific AAE values for the Aethalometer model were determined (Table S2). Compared to the default AAE values (AAE_lf= 1.0, AAE_sf= 2.0), applying station-specific AAE values resulted in a mean change of 22.3 ± 8.10 % in the fractions of eBC_sf and eBC_lf across most stations (Table S2). The source apportionment using station-specific AAE values better matched the proportions of BC from different fuel types in China (Fig. S3a). Validation against dual-carbon isotope data from previous studies (Fang et al., 2018; Ni et al., 2018) supports the use of station-specific AAE values. As shown in Fig. S4a, compared to results using default AAE values, the eBC_sf fraction calculated with station-specific AAE values showed a stronger Pearson correlation (r= 0.37 vs. 0.20) with the solid-fuel combustion fraction (coal + biomass) determined by dual-carbon isotopes. Seasonal patterns of eBC_sf derived from the Aethalometer model were consistent with those from carbon isotope analysis (Fig. S4b). Therefore, station-specific AAE values were adopted for BC source apportionment in this study.

While using station-specific AAE values derived from percentiles improved agreement with dual-carbon isotope results, the Aethalometer model still has several uncertainties. These include sensitivity to wavelength selection, the choice between fixed and variable AAE values, measurement noise, interference from non-BC absorbing components, source-specific variability, and limitations in temporal resolution (Zheng et al., 2021; Savadkoohi et al., 2025). The percentile-based approach is especially sensitive to sources with high brown carbon content (such as coal combustion) and complex source mixtures (Savadkoohi et al., 2025). Due to the unknown of the “real” BC source apportionment results, the accuracy of the Aethalometer model was assessed through error propagation (Martinsson et al., 2017; Zheng et al., 2020). The primary sources of uncertainty considered were the absorption coefficient and the AAE values. The uncertainty of b_abs was reported as 5 % (Hansen, 2005) and this value was adopted for this study. For the uncertainty of AAE values, their uncertainties were defined as the difference between the values derived from the percentile method and the isotope-constrained values from Zotter et al. (2017), e.g., 0.90 for AAE_lf and 1.68 for AAE_wb. This approach was justified as the mean optimal AAE_lf and AAE_sf values obtained here (0.90 ± 0.05 and 1.70 ± 0.23, respectively) were closely aligned with the literature values. The error propagation analysis revealed station-specific relative uncertainties in the range of 19.8 %–53.8 % for eBC_sf and 23.7 %–54.0 % for eBC_lf (Table S3). The reported uncertainties here were comparable to previous studies, e.g., 41 % for fossil fuel combustion and 42 % for wood burning estimated by Martinsson et al. (2017). Despite the relatively high estimated uncertainty, the spatial and temporal patterns of BC sources resolved from the Aethalometer model were believed reliable here.

2.4 Weather normalization of eBC and its sources

To normalize the impact of meteorological conditions on variations in eBC and its sources, a machine learning method (Grange et al., 2018) was used in this study. According to previous research (Zheng et al., 2023b), the random forest (RF) model with hyperparameters tuned during model training and temporal variables excluded in weather normalization can reduce the bias in long-term trends of air pollutants. We adopted these recommendations to normalize the impact of weather on BC and its sources at each station. The entire observation dataset was divided into two parts: 70 % was used to train the model, and the remaining 30 % was used to evaluate model performance. The variables used to train the model included meteorological conditions, temporal variables, and a transport indicator. Specifically, temperature at 2 m (t2m), surface pressure (sp), wind components at 10 m (u10, v10), relative humidity (rh), accumulated precipitation along the trajectory (apt), boundary layer height (blh), surface downward solar radiation (ssr), and total cloud cover (tcc) were used. The temporal variables included Unix time (number of seconds since 1 January 1970) as a linear trend term, Julian day (day of the year) as the seasonal term, and day of the week. The transport indicator was the cluster category of air masses reaching the observational site, calculated by the Hybrid Single-Particle Lagrangian Integrated Trajectory (HYSPLIT) (Stein et al., 2015). Three key hyperparameters, including the number of trees (ntree), the number of variables that may split at each node (mtry), and the minimum size of terminal nodes (min.node.size), were tuned by random search with 5-fold cross-validation after 100 evaluations. The search space was 100–1000, 1–13, and 1–13 for ntree, mtry, and min.node.size, respectively. After tuning, the optimal parameter combination (Table S4) with the lowest root mean square error (RMSE) was used to train the RF model. Although the differences in statistical metrics between the RF model with tuned and untuned parameters were not significant (Fig. S5), the tuned model increased the Pearson correlation coefficient (r) and index of agreement (IOA), and reduced the RMSE. For instance, the r between observations and modeled eBC increased from 0.64 ± 0.07 (untuned) to 0.66 ± 0.07 (tuned). Therefore, the tuned parameters in Table S4 were used to train the RF model. After training, weather normalization of daily eBC and its sources was achieved by repeatedly resampling predictors and predicting with the trained RF model. For each prediction, the explanatory variables (meteorological conditions) in model training were randomly sampled from the entire dataset without replacement. This predictive process was repeated 1000 times, and the arithmetic mean of these predictions was considered the emission-related concentration (Zheng et al., 2023a, b).

2.5 Auxiliary dataset, data analysis and visualization

The auxiliary datasets used in this study included meteorological reanalysis data, surface BC mass concentrations from CTM simulations, and BC emission inventories. Meteorological conditions were obtained from the fifth generation of the European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis for global climate and weather (ERA5) (Hersbach et al., 2023). Monthly mean surface BC mass concentrations from 2008 to 2020 were sourced from the Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2) (Global Modeling And Assimilation Office and Pawson, 2015) and Tracking Air Pollution in China (TAP) (Liu et al., 2022). Additionally, monthly mean BC concentrations from AeroChemMIP, endorsed by the Coupled Model Intercomparison Project 6 (CMIP6), were used in this study. We utilized the historical experiment from CMIP6, which incorporates time-varying emissions of greenhouse gases, aerosols, and ozone, as well as volcanic and solar forcing, spanning the period from 2008 to 2014. Model simulations from the Community Earth System Model 2 (CESM2), Centre National de Recherches Météorologiques (CNRM), Geophysical Fluid Dynamics Laboratory (GFDL), NASA Goddard Institute for Space Studies (GISS), Japan Meteorological Research Institute (MRI), UK Earth System Model (UKESM1) were used here (see Table S5 for details). This dataset provides a comprehensive representation of the key drivers of climate change, allowing assessment of the impacts of both anthropogenic and natural factors on global climate patterns during this period. Yearly BC emission inventories in China were obtained from several datasets, including the Copernicus Atmosphere Monitoring Service (CAMS Granier et al., 2019), Community Emissions Data System (CEDS Hoesly et al., 2018), Emissions Database for Global Atmospheric Research (EDGAR, versions V4.3.2 (Crippa et al., 2018), V5 (Crippa et al., 2019), and V6 Monforti Ferrario et al., 2021), Hemispheric Transport of Air Pollution (HTAP, version V3 Crippa et al., 2023), and the Multi-resolution Emission Inventory model for Climate and air pollution research (MEIC, version 1.4 Geng et al., 2024).

Data analysis was conducted using the R language (R Core Team, 2024). Outlier detection was performed using the “anomalize” package (Dancho and Vaughan, 2023). The “rmweather” (Grange et al., 2018) and “mlr3” (Lang et al., 2019) packages and their extensions were used for weather normalization analysis. Model performance statistics were calculated using “openair” (Carslaw and Ropkins, 2012). Figures in this study were generated using “ggplot2” (Wickham, 2016) and its extensions. Other packages (e.g., “lubridate”, “plyr”) used in this study are available from the Comprehensive R Archive Network (CRAN, https://cran.r-project.org, last access: 3 June 2025).

3 Results

3.1 Spatial-temporal variations

Table S6 summarizes the statistics of eBC mass concentrations, source apportionment results, and AAE_370–950 values at 34 monitoring stations during the study period (2015–2017). The average mass concentrations of eBC, eBC_lf, and eBC_sf, as well as AAE_370–950, were 2.05 ± 2.85, 1.08 ± 1.73, 0.97 ± 1.52 µg m⁻³, and 1.33 ± 0.29, respectively. Fig. 2 shows the spatial distributions of eBC mass concentrations and their source apportionment results between 2015 and 2017. As shown in Fig. 2a, BC from MERRA-2 showed higher levels in the east of the “Hu-Line” compared to those west of the line. As shown in Fig. 2b, ground observations also showed higher levels in eastern China (2.25 ± 2.82 µg m⁻³) than in western China (1.51 ± 2.86 µg m⁻³). Similarly, the spatial distribution of eBC observations showed higher levels in northern China (2.19 ± 3.45 µg m⁻³) than in southern China (1.91 ± 2.05 µg m⁻³). The eBC concentrations at different types of stations also showed spatial heterogeneity, with the highest values at urban stations (2.67 ± 3.41 µg m⁻³), which were significantly higher than those at rural (1.25 ± 1.50 µg m⁻³, p<0.05) and baseline (1.07 ± 1.15 µg m⁻³, p<0.01) stations. The percentages of solid fuel combustion in eBC ranged from 33.5 ± 24.9 % to 65.2 ± 24.8 %, with the highest mean value at baseline stations (49.9 ± 10.9 %). Regarding the spatial distribution (Fig. 2c), eBC_sf in eastern China (49.7 ± 7.43 %) and northern China (50.0 ± 7.25 %) were higher than those in western China (44.4 ± 4.55 %) and southern China (46.5 ± 6.77 %) without statistical differences. As shown in Fig. 2d, the station-specific AAE_370–950 values were in the range of 1.11 ± 0.13–1.66 ± 0.40 and the urban stations had the lowest mean values (1.29 ± 0.11) compared to those in baseline stations (1.40 ± 0.09) and rural stations (1.40 ± 0.16). Similar to the spatial distribution of eBC_sf, AAE_370–950 in northern China (1.43 ± 0.12) was significantly higher (p<0.001) than those in southern China (1.24 ± 0.07) and AAE_370–950 in western China (1.44 ± 0.13) was also significantly higher than those in eastern China (1.30 ± 0.12) at 0.01 levels.

Figure 2 Spatial distributions of mean values of BC concentrations, sources of eBC including solid fuel (eBC_sf) and liquid fuel (eBC_lf) combustion, and absorption Ångström exponent between 2015 and 2017 (a), and their box plots (b)–(d) showing station-specific values for different types of stations, including baseline (B), rural (RU), and urban (U), as well as different regions in China: eastern China (EC), western China (WC), southern China (SC), and northern China (NC). The sizes and colors of the dots in panel (a) correspond to the mean values of AAE and eBC, respectively. The donuts in panel (a) represent the source apportionment results of eBC in each station. The data in the filled grid in panel (a) is from MERRA-2. NS, ^*, ${}^{**}$, and ${}^{***}$ in panels (b)–(d) indicate differences between two paired groups that are not significant (p>0.05) or significant at the 0.05, 0.01, and 0.001 levels, respectively.

The monthly variations of eBC, eBC_sf, and AAE_370–950 are shown in Fig. 3. Generally, the lowest values of eBC (1.46 ± 1.72 µg m⁻³), eBC_sf (36.3 ± 24.7 %), and AAE_370–950 (1.21 ± 0.28) occurred during summer (JJA), while the highest values were observed during winter (DJF). The higher values in winter were related to unfavorable meteorological conditions that hindered the dispersion of surface BC. Increased demand for heating also contributed to higher BC levels in winter, as evidenced by higher percentages of eBC_sf (60.2 ± 22.7 %) and higher AAE_370–950 values (1.45 ± 0.26). Some stations (e.g., 58370 and 58448), however, showed the highest eBC mass concentrations during May and July. These anomalously high eBC levels were related to biomass burning during the summer harvest season, as indicated by high eBC_sf percentages. Previous studies also showed that these stations or regions were impacted by biomass burning during the summer harvest (Zha et al., 2014; Zhou et al., 2019). Another station with high levels of AAE_370–950 during summer was Tazhong (51747), which was largely impacted by dusty weather from March to July (Zhou et al., 2023).

Figure 3 Monthly variations of black carbon concentrations (eBC, a, b), the fraction of BC from solid fuel combustion (eBC_sf, c, d), and the absorption Ångström exponent calculated by power law fitting at seven wavelengths (AAE_370–950, e, f). The filled grids in the left panels represent the scaled monthly values, ranging from 0 to 1 to better show the station-specific monthly variations. Note that the darker the color represents the higher the value.

To investigate the drivers of spatial-temporal variations in BC, we performed a correlation analysis between BC concentrations, emissions, and meteorological parameters. BC emissions from different fuel types were obtained from CEDS-MAPS, which divides sources into solid bio-fuel, coal, liquid fuel plus natural gas, and process emissions (McDuffie et al., 2020). Meteorological data for 2015–2017 were sourced from the ERA5 reanalysis dataset (Hersbach et al., 2023). As shown in Fig. S6, mean BC concentrations during 2015–2017 were positively correlated with emissions at most stations (N=32), with Pearson coefficients of 0.65 for eBC, 0.56 for eBC_sf, and 0.58 for eBC_lf. This significant correlation indicates that the spatial distribution of BC closely follows its emission patterns, which are higher in northern and eastern China (Wang et al., 2012). On a monthly scale, BC concentrations were positively correlated with emissions and negatively correlated with ambient temperature, boundary layer height, and wind speed at most stations (Fig. S6). Elevated BC levels in colder months were attributed to increased heating-related emissions and less favorable dispersion conditions (lower boundary layer height and wind speed) (Zhang et al., 2019b; Xie et al., 2025). Overall, both emission patterns and meteorological factors jointly controlled the spatial-temporal variability of BC in China.

3.2 Long-term trends from 2008 to 2020

Figure 4 shows the time series of eBC, eBC_lf, and eBC_sf from 2008 to 2020. Despite eBC and its sources showing overall decreasing trends from 2008 to 2020, the year with peak values differed among the different types of stations. As shown in Fig. 4a, the annual mean concentrations of eBC and its sources at urban stations increased from 2008 to 2012 and then showed a downward trend to 2020. The inter-annual variations of eBC, eBC_lf, and eBC_sf at rural stations showed peak values in 2009 and then continuously decreased until 2020. At baseline stations, eBC and its sources also showed downward trends with fewer variations compared to those at urban and rural stations. Using the Theil-Sen method, the slopes of eBC and its sources at each station were calculated. As listed in Table S7, most stations showed decreasing trends, while only a few sites (e.g., Yushe in Shanxi province) exhibited an increasing trend. The stations with decreasing trends had a mean slope of −0.17 ± 0.20 µg m⁻³ yr⁻¹ for eBC, −0.12 ± 0.14 µg m⁻³ yr⁻¹ for eBC_lf, and −0.06 ± 0.09 µg m⁻³ yr⁻¹ for eBC_sf, respectively. Regarding the trends in the different types of stations (Fig. 4b), eBC in urban stations (N = 14) showed the fastest decreasing rate with a slope of −0.20 ± 0.23 µg m⁻³ yr⁻¹, which was quicker than that in rural stations (−0.17 ± 0.21 µg m⁻³ yr⁻¹, N=6) and baseline stations (−0.10 ± 0.07 µg m⁻³ yr⁻¹, N=5). Similarly, as shown in Fig. 4d, eBC_lf also showed the fastest decreasing rates in urban stations (−0.15 ± 0.16 µg m⁻³ yr⁻¹), followed by rural stations (−0.09 ± 0.12 µg m⁻³ yr⁻¹) and baseline stations (−0.05 ± 0.04 µg m⁻³ yr⁻¹). The decreasing rates of eBC_sf in baseline (−0.05 ± 0.04 µg m⁻³ yr⁻¹) and urban stations (−0.06 ± 0.10 µg m⁻³ yr⁻¹) were larger than those in rural stations (−0.08 ± 0.09 µg m⁻³ yr⁻¹). The differences in trends of eBC, eBC_lf, and eBC_sf between southern/western China and northern/eastern China showed no statistical significance, but higher decreasing rates were found in southern China and western China (Fig. S7).

Figure 4 Annual concentrations of black carbon (a), black carbon from liquid fossil fuel combustion (eBC_lf, c), and solid fuel combustion (eBC_sf, e) from 2008 to 2020 and box plots (b, d, f) of their trends in different types of stations including baseline (B), rural (RU), and urban (U) stations. The solid line in panels (a), (c), and (e) represents the annual mean values, and the filled ribbons represent the 99 % confidence intervals of mean values.

Long-term observational studies across China have widely reported decreases in BC concentrations over various periods. For example, BC concentrations in Beijing declined at an average rate of 0.19 µg m⁻³ per year between 2013 and 2022 (Xie et al., 2025). Similarly, an annual decrease of 0.10 µg m⁻³ was observed in Nanjing from 2014 to 2021 (Abulimiti et al., 2025), while a trend of −0.12 µg m⁻³ yr⁻¹ was recorded at Mt. Lushan from 2008 to 2022 (Liu et al., 2025). These single-site observations consistently highlight the effectiveness of air pollution control policies in reducing BC concentrations. Our study reinforces these findings from a broader spatial perspective, confirming that the observed reduction in BC levels can be attributed to the implementation of national regulatory measures, such as the China Clean Air Action Plan and the Three-Year Action Plan.

3.3 Comparison with global observations

To better understand the levels of eBC and its characteristics in China, the results from this study were compared to those from other regions during the same observational period (2015–2017) on a global scale. Specifically, the observations were sourced from published papers, observation networks, and publicly accessible websites. Data for European countries and the US were obtained from EBAS (https://ebas-data.nilu.no/Default.aspx, last access: 18 July 2025) and the US EPA (https://aqs.epa.gov/aqsweb/airdata/download_files.html, last access: 18 July 2025), respectively. The location, station type, and instruments used in these observations are listed in Table S8. The data processing procedures for these observational datasets are described in Section S1 in the Supplement. Briefly, by harmonizing wavelengths, MAC, and AAE combinations, biases were minimized to ensure reliable comparisons.

Figure 5 shows the global distributions of mass concentrations of eBC, eBC_sf percentages, and AAE values from previous studies or observation networks during 2015 and 2017. Despite differences among the sampling inlets and monitoring instruments (Table S8), the mean values of eBC observations in China (2.04 ± 1.53 µg m⁻³, N=35) were significantly lower than those in India (6.30 ± 4.12 µg m⁻³, N=23) at the 0.001 levels, but significantly higher than those in the UK (1.22 ± 1.10 µg m⁻³, N=14, p< 0.05), Spain (1.07 ± 0.98 µg m⁻³, N=9, p< 0.05), US (0.83 ± 0.63 µg m⁻³, N=65, p< 0.001), Germany (0.74 ± 0.64 µg m⁻³, N=13, p< 0.001), and Switzerland (0.61 ± 0.47 µg m⁻³, N=5, p< 0.001). Regarding BC from solid fuel combustion derived from the Aethalometer model, the results are influenced by the choice of wavelength pair and AAE combination. Other factors such as station type (e.g., background vs. traffic), instruments, and sampling inlets also impact the Aethalometer model results. The comparison here used an AAE combination of 1.0 for AAE_lf, 2.0 for AAE_sf, and the wavelength pair of 370 nm and 880 nm to minimize the influence of AAE combination and wavelength pair selection on source apportionment results. As shown in Fig. 5e, eBC_sf had the highest contribution in the US (30.2 ± 12.8 %, N=25), which was higher than that in this study (28.1 ± 10.9 %, N=34). The lowest contribution of eBC_sf was found in Spain (13.3 ± 5.51 %, N=3), which was lower than that in China (p< 0.05). Similar to the spatial distribution of eBC_sf, AAE also showed the highest value in the US (1.40 ± 0.13, N=26) and the lowest value in Spain (1.21 ± 0.06, N=3). The comparison of eBC trends between this study and other observations from 2008 to 2020 is shown in Fig. 6. Results from 88 stations (34 in this study and 54 from other observations) showed widespread decreasing trends of eBC in most stations (N=83), with mean trends of −4.36 ± 3.06 % yr⁻¹. The slopes for eBC from 2008 to 2020 showed the fastest decreasing trend in the UK (−5.94 ± 2.26 % yr⁻¹, N=12) and the slowest decreasing trend in the US (−2.61 ± 2.04 % yr⁻¹, N=18). The eBC decreasing rates in Germany (−4.93 ± 2.51 % yr⁻¹, N=16) and Switzerland (−5.66 ± 5.44 % yr⁻¹, N=7) were also faster than those in this study for China (−3.98 ± 3.41 % yr⁻¹, N=25).

Figure 5 Spatial distributions and box plots in black carbon concentrations (a, d), black carbon from solid fuel combustion (b, e), and absorption Ångström exponent (c, f) during 2015–2017. The filled base maps of panels (a) and (c) are from MERRA-2 and the base map of panel (b) is from CEDS-MAPS (McDuffie et al., 2020). The dots in panels (a)–(c) represent the observation values (see Table S8 for details) and the size of each dot is mapped to its value. Panels (d)–(f) are box plots of observations in different countries including India (IN), Switzerland (CH), China (CN), United Kingdom (UK), Spain (ES), Germany (DE), and the United States (US).

Figure 6 Spatial distribution of black carbon concentration trends (a) and box plots of trends in different countries (b) from 2008 to 2020. The size of dots in panel (a) corresponds to the absolute value of slope. The fill color in panel (b) corresponds to the mean slope value in each country. Detailed information about slopes at each station can be found in Table S9.

This international comparison provides a global context for evaluating China's BC concentrations and trends, allowing for direct comparison with both developing countries (e.g., India) and developed countries (e.g., the United States, United Kingdom, Germany, Switzerland). The results highlight that while China achieved reductions in BC concentration, its concentrations still remain higher, and its rates of decline were slower than those observed in developed nations. In these countries, rapid BC reductions have been achieved through the implementation of stringent regulations (Luoma et al., 2021) and the adoption of cleaner technologies (Font and Fuller, 2016). These findings underscore the need for continued policy efforts in China, particularly targeting industrial and residential sources, and support the refinement of emission inventories to further reduce BC emissions.

4 Discussion

4.1 Comparison with CTM simulations

Due to the lack of long-term global-scale observations, knowledge of atmospheric BC abundance continues to rely on CTM simulations. Therefore, it is important to evaluate the performance of CTMs in simulating the spatial-temporal variations of BC. Owing to the availability of long-term BC observations in China, the BC simulation results from different models (see Sect. 2.5 for details) were evaluated in this study. Regarding model performance in simulating the spatial distribution of BC, these models showed high values in East China, which were consistent with ground observations (Fig. 7). The r between mean values of simulations and observations from 2008 to 2014 ranged from 0.51 (MRI) to 0.70 (TAP). Compared to the observations, the CTMs underestimated the surface BC mass concentrations by 80.8 % (GISS) to 42.1 % (UKEMS1) in China from 2008 to 2014. The underestimation of surface BC concentrations by CTMs has been widely reported elsewhere (Koch et al., 2009; Fu et al., 2012; Mao et al., 2016; Qi and Wang, 2019). Among the different types of stations, the CTMs showed the highest performance in reproducing the surface BC concentrations at baseline stations (r= 0.83 ± 0.07, slope = 0.88 ± 0.42), followed by urban (r= 0.65 ± 0.10, slope = 0.44 ± 0.26), and rural stations (r= 0.20 ± 0.10, slope = 0.13 ± 0.08). The higher performance of CTMs in simulating BC in urban areas compared to rural areas is related to uncertainty in BC emissions. For instance, BC is mainly from vehicle exhaust in urban areas and is dominated by biomass burning and coal combustion in rural regions. The uncertainty in BC emission factors for gasoline and diesel vehicles is lower than that for biomass burning (Rönkkö et al., 2023). As a result, BC emissions in urban areas show less uncertainty (Zhao et al., 2011; Zhang et al., 2020), and simulations from CTMs produce less deviation from observations in urban areas.

Figure 7 Spatial distributions of surface BC observations and simulations from different chemical transport models during 2008–2014 and the statistical metrics of model performance. The NMB and IOA represent normalized mean bias and index of agreement, respectively.

Regarding the temporal variations of BC in China, the simulations from CMIP6, MERRA-2, and TAP showed similar monthly variations, with higher concentrations in winter and lower values in summer (Fig. 8a). The r between CTMs and observations ranged from 0.63 (GISS) to 0.87 (MERRA-2) on a monthly resolution. On the annual scale (Fig. 8b), except for TAP, annual concentrations of BC from CMIP6 historical experiments and MERRA-2 showed increasing trends from 2008 to 2014, while the ground observations and TAP showed decreasing trends in China. Except for TAP and MERRA-2, the correlations between observations and CMIP6 results were negative at the annual scale (Fig. 8d), suggesting these models failed to reproduce the yearly variations of surface BC. The upward trends of BC from CMIP6 historical experiments and MERRA-2 were related to their inter-annual BC emissions, which also showed increasing trends from 2008 to 2014 in China (Fig. S8). Similarly, the downward trends of surface BC concentrations from TAP were also in accordance with their emissions (MEIC), which showed a decreasing trend from 2008 to 2014 (Fig. S8). The opposite BC trends between ground observations and model simulations from CMIP6 in this study were similar to previous studies (Zhang et al., 2019b; Ramachandran et al., 2022). For instance, aerosol optical depth from observations showed a decrease, while CMIP6 simulations showed increasing trends during 2002–2018 in China (Ramachandran et al., 2022). Therefore, caution should be exercised when interpreting the temporal variations of aerosols using simulations from CMIP6.

Figure 8 Monthly (a) and inter-annual (c) variations of BC observations and surface BC simulations from different chemical transport models during 2008–2014, and their Taylor Diagrams showing model performance at monthly (b) and yearly (d) scales.

Due to the chemical inertness of BC particles, the deviations between observations and simulations arise from uncertainties in two aspects: emissions and deposition (Fan et al., 2018, 2022). The underestimation of BC emissions in China results in lower simulations of BC by CTMs. For instance, BC simulations were underestimated by 35 %–60 % in different seasons using the bottom-up BC emissions in China, while the CTM showed better performance with BC underestimated only by 22 %–9 % using the top-down constrained emissions (1.59 times higher than bottom-up emissions) (Fu et al., 2012). Similarly, the RMSE for annual BC simulation decreased by 31.6% if the BC emissions in China increased by 1.8 times (Wang et al., 2016). Using the MEIC instead of the emissions inventory for the Intergovernmental Panel on Climate Change Fifth Assessment Report (BC emissions increased by 13.4 %), the CTM improved the BC concentrations by 42.6 % (Fan et al., 2018). Despite the increase of BC emissions in CTMs, the models still underestimated the BC concentrations, suggesting an overestimation of BC deposition within CTMs. For instance, a better agreement between observations and simulations was achieved by reducing the BC dry deposition velocity by 50 % (Huang et al., 2010). Similarly, using the lower dry deposition, the CTM simulations resulted in a lower global mean BC dry deposition flux and a higher correlation between simulations and observations (Wu et al., 2018). The observations also suggested an overestimation of the BC aging degree in the CTMs, and sensitive simulations showed that slowing down the BC aging degree in model simulation can reduce the large model bias in the simulation of surface BC concentration over China (Shen et al., 2023).

To identify which factor (uncertainty in emission or deposition) makes a greater contribution to the uncertainty of BC simulation, we used the importance of each variable in RF model building (Sect. 2.4). As shown in Fig. S9, generally, date_unix was the variable with the greatest importance for the prediction of eBC, eBC_lf, and eBC_sf in the RF model. Additionally, Julian day (seasonal term) had the third largest importance for the prediction of eBC and eBC_lf, while temperature showed higher importance for the prediction of eBC_sf. These explanatory variables can serve as proxies for emissions of BC, e.g., more biomass burning during colder months for heating. The explanatory variables related to deposition (e.g., wind speed, precipitation), however, had less importance for prediction. Therefore, emissions were more important than deposition when explaining BC concentrations.

The comparison between the observed BC concentration values and simulations suggests that to improve the model's accuracy, more attention should be given to updating the BC emission inventory. Here, we only discuss how to use BC observation data for updating the BC emission inventory, rather than attempting to complete this task. In practice, the BC emission inventory, observed concentrations, FLEXPART simulations, and Bayesian inversion framework can be integrated to update BC emissions (Evangeliou et al., 2018, 2021; Jia et al., 2021). For instance, Evangeliou et al. (2021) used BC observations, outdated prior BC emissions, and an inversion framework to estimate BC emissions during the COVID-19 lockdown in Europe. Their results showed a 23 % decline in BC emissions during the lockdown compared to the same period in the previous five years (2015–2019). Similarly, Jia et al. (2021) estimated a maximum weekly emission reduction of 70 % in eastern China and 48 % in northern China from 1 January to 10 March 2020. These studies showed improved accuracy of model simulations using posterior BC emissions compared to the prior emission inventory, highlighting the potential to reduce BC simulation uncertainty by incorporating BC observations and prior emissions.

4.2 Weather normalized concentrations: an indicator of BC emissions

Using the weather normalization method, the concentrations of eBC and its two sources (eBC_lf and eBC_sf) were normalized. Compared to the eBC observations, a higher correlation between weather-normalized eBC concentrations and BC emissions was found (Fig. S10), suggesting that the normalized eBC concentration can better reflect emission variations. Therefore, the weather-normalized eBC concentration was used here to discuss its relationship with BC emissions. As shown in Fig. 9a, except for MEIC, BC from other emission inventories showed increasing trends since 2008, reaching peak values in certain years and then decreasing. In contrast, BC emissions from MEIC showed a continuous decrease since 2008, and their inter-annual variation was more consistent with weather-normalized eBC concentrations. In 2010, China began to fight against air pollution, and the Action Plan on the Prevention and Control of Air Pollution and the Three-Year Action Plan were implemented during 2013–2017 and 2018–2020, respectively. As a result of these actions, emission standards and activity levels changed over time. These factors are considered in the MEIC framework, which employs a technology-driven methodology to track China's energy statistics and technological advancements, reflecting shifts in emission characteristics across various sectors, fuels, products, combustion/process technologies, and emission control technologies in recent years (Li et al., 2017; Zheng et al., 2018; Geng et al., 2024). Therefore, the BC emissions from MEIC coincided more closely with weather-normalized concentrations. The BC emissions from other inventories, such as CEDS, are not calibrated; instead, they are derived from default emissions using the Speciated Pollutant Emission Wizard, due to the limited availability of national BC inventory estimates in other countries. The similar reduction rates of surface BC concentrations and emissions (53.6 % vs. 51.9 %) and their similar inter-annual variations from 2008 to 2020 in China also suggest that the normalized BC concentrations can be used as a proxy to scale BC emissions in recent years. Updating bottom-up BC emissions is time-consuming, as it requires collecting activity data (e.g., energy consumption). As summarized in Table S8 and mentioned above, there are several BC observation networks worldwide (e.g., UK, US, Germany, and India); the normalized BC observations can be used to calibrate the total BC emissions in these countries.

Figure 9 The inter-annual variations of black carbon emissions from different emission inventories and observations (eBC) during 2008–2020 (a) and the BC emissions from solid fuel combustion (EMI_sf) and liquid fossil fuel combustion (EMI_lf) derived from MEIC and the BC sources including solid fuel combustion (eBC_sf) and liquid fossil fuel combustion (eBC_lf) derived from the Aethalometer model (b).

To further explore whether the weather-normalized BC source contributions can be used as an indicator of BC emissions from different fuel types, the BC emissions from MEIC were used here. The sectors in MEIC (industrial, power, residential, and transportation) were further broken down into fossil fuel (liquid fuel + coal) and biomass burning according to (Wang et al., 2012). As shown in Fig. 9b, the variations of eBC_lf concentrations from the Aethalometer model aligned with their emissions, showing a continuous reduction from 2008 to 2020. In contrast, the inter-annual variation of eBC_sf was not in line with its emissions, with a small peak in emissions in 2013, while its concentration continuously decreased since 2008. The higher correlation between eBC_lf and liquid fuel combustion emissions compared to that between eBC_sf and solid fuel emissions (Fig. S11) suggests higher uncertainty in BC emissions from solid fuel combustion, as discussed in Sect. 4.1 and previous studies (Zhao et al., 2013; Zhang et al., 2020). Therefore, to further improve the accuracy in the BC emission inventory, more attention should be paid to the solid fuel combustion sub-sector, e.g., focusing on surveys of activity rates and measurements of local emission factors (Li et al., 2017; Jiang et al., 2024).

Accuracy in sectoral emissions is important to improve the performance of CTMs. Previous studies have demonstrated that uncertainties in BC emission inventory, originating from incomplete activity data, unrepresentative emission factors , and sectoral inconsistencies (Li et al., 2017), will propagate significant biases into chemical transport models. For instance, sectoral uncertainties are pronounced for residential and industrial sources in developing countries due to inconsistent emission factors and activity data. From a modeling perspective, these uncertainties dominate over other processes (e.g., aging, deposition) in near-source regions, while transport and removal uncertainties amplify errors in remote areas (Vignati et al., 2010). Another study showed that BC simulation was overestimated by a factor of 2.22, which can be attributed to the highest domestic emissions during winter in East Asia (Ikeda et al., 2022). Results from a top-down study showed the largest bias in residential and transportation sectors during January and July, and the bias of BC was reduced using the posterior emission inventory compared to the prior emission inventory (Zhao et al., 2019).

4.3 Drivers of black carbon variations in China

Using the weather-normalization method, the contributions of emissions and meteorology to the inter-annual variations of air pollutants can be quantified (Chen et al., 2019; Zheng et al., 2020, 2023b). As shown in Fig. S12, compared to the observations, the emission-related (weather-normalized) concentrations of BC were higher in each year. For instance, the national mean concentrations of emission-related and observed eBC in 2008 were 3.59 and 3.47 µg m⁻³, respectively, yielding meteorology-related concentrations of −0.12 µg m⁻³. The negative value of the meteorology-related concentration indicates favorable meteorological conditions for reducing BC concentrations. Similarly, meteorology also contributed to lowering the concentrations of eBC_lf and eBC_sf at different types of stations (Fig. S12). To quantify the contributions of emission and meteorology to the reduction of BC concentrations, the slope between weather-normalized and observational trends was used. Results from Fig. 10a indicate that the reductions of eBC, eBC_lf, and eBC_sf were dominated by emission-related reductions, with relative contributions of 89 %, 87 %, and 62 %, respectively, implying the effectiveness of human efforts in BC emission reduction. These results are consistent with previous studies showing that the inter-annual variations of air pollutants were dominated by emission changes (Chen et al., 2019; Cheng et al., 2019; Zhang et al., 2019a; Zheng et al., 2020, 2024). Compared to eBC and eBC_lf, eBC_sf was more strongly impacted by variations in meteorological conditions, e.g., more coal and biomass were consumed for heating during colder months. To better understand which sources dominated the reduction of BC, the results from the Aethalometer model were used. As shown in Fig. 10b, the slopes of the fittings between eBC_sf and eBC_lf in baseline and rural stations were greater than 1, suggesting the dominant role of BC emission reduction from solid fuel in the decreases of ambient eBC concentrations in these two types of regions. In contrast, the slope between eBC_sf and eBC_lf in urban stations was less than 1, suggesting the dominant role of BC emission reductions from liquid fuel (e.g., vehicle emissions) in the decrease of eBC concentration in urban areas. The control measures related to BC emission reductions include stricter vehicle emission and industrial standards, eliminating outdated industrial capacity, prohibiting agricultural residue burning and scattered coal combustion, and replacing coal with natural gas and electricity (Zheng et al., 2018). The focus of BC mitigation measures, however, differs between rural and urban areas. In urban areas, BC emission reductions are related to stricter vehicle emission standards, improved fuel quality, and eliminating outdated industrial capacity (Xu et al., 2018; Zheng et al., 2018; Wang et al., 2020; Zhang et al., 2021), which result in less BC emissions from liquid fuel. In rural areas, BC emission reductions result from prohibiting agricultural residue burning, household stove upgrading, and energy switching, e.g., from coal to natural gas and electricity (Zhu et al., 2019; Meng et al., 2021; Shen et al., 2022). These measures lead to BC emission reductions from biomass burning and coal combustion. The results here are consistent with the differences in BC emission reduction measures between cities and rural areas, with rapid decreases of eBC_sf in rural areas and eBC_lf in urban regions. Due to the limitation of the Aethalometer model, the BC sources are apportioned into eBC_sf and eBC_lf, while the contributions of solid fossil fuel (e.g., coal) and biomass cannot be further separated from eBC_sf (Zheng et al., 2021). To further understand which types of fuels (e.g., biomass, coal, liquid fuel) dominated the BC reductions in China, we used the estimated BC emissions from different fuels derived from MEIC. As shown in Fig. 10c, BC from coal combustion decreased by 0.38 Tg from 2008 to 2020 in China, followed by liquid fuel (0.16 Tg), and biomass burning (0.09 Tg). Therefore, emission reduction from coal combustion dominated the BC reduction in China, with a contribution of 60.3 %.

Figure 10 The drivers of BC reduction including contributions from emission variations (a), different sources (b), and fuel types including liquid fossil fuel (LF), coal combustion, and biomass burning (BB) (c) in China.

In line with previous long-term observations of BC conducted in urban areas (Chen et al., 2016; Zheng et al., 2020; Sun et al., 2022; Xie et al., 2025; Abulimiti et al., 2025), these studies identified the predominant driving factor for BC reduction was the reduced emissions from fossil fuel combustion (e.g., coal and petroleum). Representative examples include studies in Wuhan (Zheng et al., 2020), Beijing (Xie et al., 2025), and Nanjing (Abulimiti et al., 2025). Unlike previous studies, this research also analyzed the driving factors behind BC reduction at rural and background sites. The results demonstrate that the decrease in BC concentrations is mainly attributed to reduced emissions from both coal and biomass combustion in rural and baseline areas. Evidence from BC emission data further support these findings. As shown in Fig. S13, BC emissions derived from the MEIC inventory (Geng et al., 2024) for biomass burning, solid fuel combustion, and liquid fuel combustion exhibited widespread reductions across most regions of eastern China between 2008 and 2020. Temporally, emissions from these sources showed consistent decreasing trends from 2008 to 2020 at rural, baseline, and urban stations. For example, BC emissions from rural residential energy consumption in northern China decreased from 0.85 Tg in 2010 to 0.55 Tg in 2020 (Zhang et al., 2023). These findings present novel and previously unreported conclusions based on observational data.

5 Conclusions

In this study, 13 years of continuous measurements of black carbon aerosol were analyzed from 48 stations in China. Using station-specific AAE_lf and AAE_sf, the sources of eBC were apportioned. The levels, spatial-temporal characteristics, and trends of black carbon aerosol were reported, and the key findings of this unique dataset are listed as follows.

Observations of black carbon aerosols from 2015–2017 showed averages of 2.05 ± 2.85, 1.08 ± 1.73, 0.97 ± 1.52, and 1.33 ± 0.29 for eBC, eBC_lf, eBC_sf, and AAE_370–950, respectively, in China. Long-term trends of eBC, eBC_sf, and eBC_lf indicated reductions from 2008 to 2020, with mean slopes of −0.17 ± 0.20 µg m⁻³ yr⁻¹ for eBC, −0.12 ± 0.14 µg m⁻³ yr⁻¹ for eBC_lf, and −0.06 ± 0.09 µg m⁻³ yr⁻¹ for eBC_sf on a national scale. Spatial distributions of eBC and eBC_sf were higher in eastern and northern China compared to western and southern regions. Among different types of stations, urban stations exhibited the highest eBC values, while baseline stations had the highest fractions of eBC_sf and AAE_370–950. The seasonal variations of eBC, eBC_sf, and AAE_370–950 typically showed the lowest values in summer and the highest in winter. However, some stations recorded abnormally high eBC levels during summer, which were associated with biomass burning and dusty weather events. The spatial-temporal patterns of BC in China can be explained by its emission variations and meteorological conditions. The weather normalized BC concentrations showed that emission reductions were the primary driver, accounting for 89 % of the decrease in eBC, 87 % in eBC_lf, and 62 % in eBC_sf. Urban stations showed the most rapid declines in eBC_lf, with the ratio of eBC_lf to eBC_sf trends exceeding 1, indicating that reductions in liquid fuel combustion emissions were the main contributor to BC decreases in urban areas. In contrast, rural and baseline stations exhibited the fastest declines in eBC_sf, suggesting that emission reductions from solid fuel combustion were the dominant factor in these regions.

This study provides valuable insights from China's long-term observations of BC aerosol, with significant implications for refining emission inventories and reducing model uncertainties. The comparison between observations and simulations revealed that most models failed to capture the inter-annual variation while the simulation with MEIC inventory was capable to reproducing the spatial and temporal variations of BC. Therefore, BC emissions from MEIC is recommended for modeling BC concentrations in China. Furthermore, a systematic model underestimation of BC was found particularly at rural stations, which emphasized the need to refine both emission source strengths and deposition processes to improve model accuracy. This study showed that the weather normalization is a promising technique for refining BC emission inventories, as the weather normalized BC concentrations showed a stronger correlation with the reported BC emissions than raw observations. Further analysis revealed a higher correlation between emissions and concentrations of BC from liquid fuel combustion. To reduce uncertainties in BC emissions, a target effort to improving accuracy of emission from the solid fuel combustion sub-sector is needed, including the collection of accurate BC emission factors from biomass and coal burning, as well as updating activity data in rural areas.