Source-explicit estimation of brown carbon in the polluted atmosphere over the North China Plain: implications for distribution, absorption, and the direct radiative effect

Brown carbon (BrC) plays a significant role in altering atmospheric radiation. Beyond biomass and biofuel combustion, recent studies identify fossil fuel sources – especially residential coal burning and vehicle exhaust – as major contributors to BrC. This underscores a gap in climate models, which often assume fossil fuel organic aerosols (OAs) are non-absorbing or treat all OA as light-scattering. In this study, we simulate BrC over the North China Plain (NCP) during a winter pollution event using the WRF-Chem model, incorporating explicit BrC absorption properties. The model aligns well with observed pollutant and aerosol levels, revealing an average near-surface BrC concentration of 5.2 µg m⁻³, contributing 16.4 % to aerosol absorption at 365 nm. Using a diagnostic adjoint approach, we estimate that BrC exerts a direct radiative effect (DRE) averaging −0.09 W m⁻² at the top of the atmosphere, reducing the cooling effect of organic carbon by 28.0 % and producing a local warming effect of up to +0.40 W m⁻². Coal combustion is the largest BrC source in the NCP in 2014, though secondary BrC also significantly impacts the regional radiation balance.

1 Introduction

Brown carbon (BrC) is a collective component for those colored organic compounds with wavelength-dependent light absorption properties (Mukai and Ambe, 1986; Kirchstetter et al., 2004; Andreae and Gelencsér, 2006). BrC has been recognized as an important short-lived climate forcer contributing considerably to climate change by warming of the atmosphere (IPCC, 2013; Feng et al., 2013; Jacobson, 2014; Jo et al., 2016; Brown et al., 2018). A study suggested that the light absorption induced by BrC can be equal to or even higher than that of black carbon (Pokhrel et al., 2017) and substantially influences atmospheric radiative forcing. Recent studies have shown that BrC accounts for 30 %–50 % of the total absorption of aerosols in Atlanta, USA; Brazil; and Hebei, China (Hoffer et al., 2006; Yang et al., 2009; Liu et al., 2013). The direct radiation effect (DRE) caused by BrC is greater than +1 W m⁻² in some regions, such as South Asia, Africa and Southeast Asia, much higher than the global average (Park et al., 2010; Feng et al., 2013; Lin et al., 2014; Saleh et al., 2014; Jo et al., 2016; Wang et al., 2018; Yan et al., 2018). However, the modeled DRE associated with BrC remains highly uncertain, with variations spanning an order of magnitude difference. In particular, the estimated global DRE of BrC is in the range between +0.03 and +0.57 W m⁻² (Hammer et al., 2016), which is caused by the limited observations of BrC mass and absorption properties (Tuccella et al., 2020; Saleh, 2020).

It has been well established that BrC is not a single substance but a general term for light-absorbing organic aerosols. A series of laboratory measurements and observations in the earlier years demonstrate that BrC is mainly associated with smoldering biomass burning (BB) or biofuel (BF) combustion (Chakrabarty et al., 2010; Chen and Bond, 2010; Lack et al., 2012; Washenfelder et al., 2015; Kumar et al., 2018). On the other hand, OAs from fossil fuel combustion are generally assumed to be non-absorbing as the combustion conditions for fossil fuels (FFs) are typically not conducive to BrC formation (Hecobian et al., 2010; Shapiro et al., 2009; Bond et al., 2013). Therefore, earlier climate model studies have assumed that primary OA from BB and BF combustion is the main or sole BrC source (Feng et al., 2013; Jacobson, 2014; Saleh et al., 2014; Hammer et al., 2016; Brown et al., 2018). Recent studies have also incorporated the aging of secondary organic aerosol (SOA) (Jo et al., 2016; Wang et al., 2018; Zhang et al., 2020). However, more recent exceptions are being found in low-efficiency residential coal combustion (RCC) (Bond, 2001; Yan et al., 2017; Xie et al., 2019; Tian et al., 2019; Zhang et al., 2022a) and fuel oil combustion in vehicle and ship engines (Xie et al., 2017; Corbin et al., 2019; Tang et al., 2020; Huang et al., 2022). It is now generally accepted that the formation of BrC is not exclusively linked to the chemical make-up of biomass fuels but is most critically determined by the combustion conditions (Saleh et al., 2018; Cheng et al., 2020; Saleh, 2020; Wang et al., 2022a). The key factor contributing to the high levels of BrC observed from biomass fuels is their combustion under relatively low-temperature and fuel-rich conditions, which are highly favorable for BrC formation. In contrast, fossil fuels, such as those burned in internal combustion engines, typically undergo combustion at higher temperatures and under more fuel-lean conditions, which are less conducive to BrC production (Saleh, 2020). In China, a developing country, coal is commonly used for residential heating in the cold season, causing massive emissions of organic particles (Yan et al., 2017; Li et al., 2018). According to the National Bureau of Statistics of China (https://data.stats.gov.cn, last access: 10 March 2025), the coal consumption in 2014 was about 4000 Tg, accounting for 65.8 % of the total primary energy use of China. Of this, around 93 Tg is used as household fuel. The poor burning conditions and limited emission control facilities in this region could lead to substantial emissions of BrC. This could explain why, to date, all reported instances of coal-derived BrC have originated from China. Both Yan et al. (2017) and Mo et al. (2021) used dual carbon isotope-based source apportionment methods and reported that fossil fuel, especially coal combustion from the residential sector, is an important source in northern China and even the largest contributor in some regions.

These recent findings indicate a critical gap in the treatment of BrC in chemical transport models, atmospheric chemistry models and climate models as well, which present an even greater concern, as they typically do not consider BrC at all (Ma et al., 2021; Jo et al., 2023; Gao et al., 2025; Ge et al., 2025). Addressing this gap requires expanding the scope of BrC sources and assigning distinct optical properties for each source and incorporating those that have been underrepresented or overlooked in past assessments into numerical models. In this study, we include the main primary emission sources (RCC, BB, BF and FF emissions from RCC and on-road vehicles (FF-TRA)) of BrC and secondary-derived BrC in a regional model, the Weather Research and Forecasting model coupled with Chemistry (WRF-Chem). A representative region, the North China Plain (NCP), is chosen as the study domain with high anthropogenic carbonaceous aerosols due to the widespread use of coal and biomass burning for heating during winter and the increasing number of motor vehicles. We performed a month-long simulation to evaluate the surface distribution, absorption and the DRE of BrC in the NCP, by updating BrC optical properties of different sources. Sensitivity experiments have also been devised to assess the contribution of BrC from major sources.

2 Model and method

2.1 WRF-Chem model and configurations

The WRF-Chem model (Grell et al., 2005; Fast et al., 2006) modified by Li et al. (2010, 2011a, b, 2012) is used to quantitatively estimate the BrC in the NCP. A heavily polluted month from 1 to 30 January 2014 is selected for the simulation period. The anthropogenic emissions are developed by Zhang et al. (2009) and Li et al. (2017), including five sources, namely agriculture, industry, power generation, residential, and transportation. The biogenic emissions are calculated online using the MEGAN (Model of Emissions of Gases and Aerosol from Nature) model developed by Guenther et al. (2006). Additionally, the grid-based RCC, BB and BF combustion emissions are used to update the BrC sources in this study, described later. The model simulation domain is shown in Fig. 1. The detailed model description and configurations can be found in Sect. S1 and Table S1 in the Supplement, respectively.

Figure 1 WRF-Chem simulation domain with topography. The square denotes the field sites for simulation and observation comparison.

2.2 Aerosol radiative module

The aerosol radiative module developed by Li et al. (2011a) has been incorporated into the WRF-Chem model to calculate the aerosol optical depth (AOD or τ_a), single scattering albedo (SSA or ω_a), and asymmetry factor (g_a). In the aerosol module, aerosols are represented by a three-moment approach with a lognormal size distribution:

$$\begin{array}{}\text{(1)}& n\left(\ln D\right)={\displaystyle \frac{N}{\sqrt{\mathrm{2}\mathit{\pi }}\ln {\mathit{\sigma }}_{\mathrm{g}}}}\exp \left[-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\displaystyle \frac{\ln D-\ln D_{\mathrm{g}}}{\ln {\mathit{\sigma }}_{\mathrm{g}}}}\right)}^{\mathrm{2}}\right],\end{array}$$

where D is the particle diameter, N is the number distribution of all particles in the distribution, D_g is the geometric mean diameter, and σ_g is the geometric standard deviation. To calculate the aerosol optical properties, the aerosol spectrum is divided into 48 bins from 0.002 to 20.0 µm, with radius r_i. The aerosols are classified into four types: (1) internally mixed sulfate, nitrate, ammonium, hydrophilic organics and black carbon (BC) and water; (2) hydrophobic organics; (3) hydrophobic BC; and (4) other unidentified aerosols (generally dust-like aerosols). These four kinds of aerosols are assumed to be mixed externally. For the internally mixed aerosols, the complex refractive index at a specific wavelength (λ) is calculated based on the volume-weighted average of the individual refractive index. Given the particle size and complex refractive index and the extinction efficiency (Q_e), ω_a and g_a are calculated using the Mie theory at a certain wavelength (λ). The look-up tables of Q_e, ω_a and g_a are established according to particle sizes and refractive indices to avoid multiple Mie scattering calculation. The aerosol optical parameters are interpolated linearly from the look-up tables with the calculated refractive index and particle size in the module. The τ_a at a certain λ in a given atmospheric layer k is determined by the summation over all types of aerosols and all bins:

$$\begin{array}{}\text{(2)}& {\mathit{\tau }}_{\mathrm{a}}\left(\mathit{\lambda }k\right)={\sum }_{i=\mathrm{1}}^{\mathrm{48}}{\sum }_j^{\mathrm{4}}{Q}_{\mathrm{e}}(\mathit{\lambda },r_i,j,k)\mathit{\pi }{r}_i^{\mathrm{2}}n(r_i,j,k)\mathrm{\Delta }{Z}_k,\end{array}$$

where $n\left(r_i,j,k\right)$ is the number concentration of the jth kind of aerosols in the ith bin. ΔZ_k is the depth of an atmospheric layer. The weighted-mean values of ω_a and g_a are then calculated using D'Almeida et al. (1991):

$$\begin{array}{}\text{(3)}& {\displaystyle }{\mathit{\omega }}_{\mathrm{a}}(\mathit{\lambda },k)={\displaystyle \frac{\begin{array}{c}{\sum }_{i=\mathrm{1}}^{\mathrm{48}}{\sum }_j^{\mathrm{4}}{Q}_{\mathrm{e}}(\mathit{\lambda },r_i,j,k)\mathit{\pi }{r}_i^{\mathrm{2}}n(r_i,j,k)\\ \times {\mathit{\omega }}_{\mathrm{a}}(\mathit{\lambda },r_i,j,k)\mathrm{\Delta }{Z}_k\end{array}}{\begin{array}{c}{\sum }_{i=\mathrm{1}}^{\mathrm{48}}{\sum }_j^{\mathrm{4}}{Q}_{\mathrm{e}}(\mathit{\lambda },r_i,j,k)\mathit{\pi }{r}_i^{\mathrm{2}}\\ \times n(r_i,j,k)\mathrm{\Delta }{Z}_k\end{array}}}\text{(4)}& {\displaystyle }{g}_{\mathrm{a}}(\mathit{\lambda },k)={\displaystyle \frac{\begin{array}{c}{\sum }_{i=\mathrm{1}}^{\mathrm{48}}{\sum }_j^{\mathrm{4}}{Q}_{\mathrm{e}}(\mathit{\lambda },r_i,j,k)\mathit{\pi }{r}_i^{\mathrm{2}}n(r_i,j,k)\\ \times {\mathit{\omega }}_{\mathrm{a}}(\mathit{\lambda },r_i,j,k)g_{\mathrm{a}}(\mathit{\lambda },r_i,j,k)\mathrm{\Delta }{Z}_k\end{array}}{\begin{array}{c}{\sum }_{i=\mathrm{1}}^{\mathrm{48}}{\sum }_j^{\mathrm{4}}{Q}_{\mathrm{e}}(\mathit{\lambda },r_i,j,k)\mathit{\pi }{r}_i^{\mathrm{2}}n(r_i,j,k)\\ \times {\mathit{\omega }}_{\mathrm{a}}(r_i,j,k)\mathrm{\Delta }{Z}_k\end{array}}}.\end{array}$$

When the wavelength-dependent τ_a, ω_a, and g_a are calculated, they can be used in the Goddard shortwave module.

It is worth noting that the aerosol liquid water content in the study is predicted with the inorganic aerosols using a computationally efficient thermodynamic equilibrium model, ISORROPIA (version 1.7, Nenes et al., 1998; Fountoukis and Nenes, 2007). In this study, ISORROPIA is mainly used to predict the thermodynamic equilibrium between the ammonium–sulfate–nitrate–chloride–water aerosols and their gas-phase precursors H₂SO₄–HNO₃–NH₃–HCL–water vapor, and water uptake of aerosols is calculated using the Zdanovskii–Stokes–Robinson (ZSR) correlation (Stokes and Robinson, 1966):

$$\begin{array}{}\text{(5)}& W={\sum }_i{\displaystyle \frac{M_i}{m_{oi}\left(a_w\right)}},\end{array}$$

where W is the mass concentration of aerosol liquid water (kg m⁻³ air), M_i is the molar concentration of species i (mol m⁻³ air) and m_oi(a_w) is the molality of an aqueous binary solution of the ith electrolyte with the same a_w (i.e. relative humidity) as in the multicomponent solution.

The BrC in the model has an effective density of 1.2 g cm⁻³ for primary BrC (Turpin and Lim, 2001) and of 1.0 g cm⁻³ for secondary BrC (Hurley et al., 2001). The imaginary refractive index of BrC used in this study is discussed in Sect. 2.3.2.

2.3 Model modifications

2.3.1 Source separation of BrC

The definition of BrC in the model is dependent on its sources. Due to the lack of BrC emission inventories, according to the characteristics of energy structure in China, assumptions and code modifications of the WRF-Chem model have been made to consider the BrC from different sources. These involve three separate primary BrC sources, including BB emissions, FF emissions from RCC and on-road vehicles (FF–TRA), and SOA with light-absorbing properties, whereas other types of primary OA (POA) and SOA are assumed to be purely scattering. In this study, the BB source corresponds to open fire, household biomass burning and biofuel consumption emissions.

For the primary emissions, previous BrC simulations have substituted it with a proportion of POA directly (Feng et al., 2013; Lin et al., 2014; Wang et al., 2014; Tuccella et al., 2020; Xu et al., 2024), derived it from the relationship between the burning efficiency and the observed aerosol light absorption (Jo et al., 2016; Zhu et al., 2021), or determined it through parameterization where BrC absorption is a function of the BC-to-OA emission ratio (Zhang et al., 2020). In the present work, we calculated the primary BrC emissions based on the bottom-up OA emission inventory combined with reported annual BrC emissions from various primary sources, as shown in Table 1. Firstly, we collected the reported annual emissions of BrC from RCC, BB and FF-TRA by using bottom-up inventory method. It should be noted that given the proximity of the study period (January 2014) to 2013, we use the emissions of BrC from RCC and BB in 2013 provided by Sun et al. (2017, 2021), which is 592 and 712 Gg, respectively. The emissions of FF-TRA-derived BrC are 76 Gg, which is calculated based on the value of 2017 (Wang et al., 2022a) and scaled by a factor of 0.70 to reflect the change of annual civilian-owned motor vehicles. We assume that the spatial and seasonal variation of BrC is similar to OA. Then, the monthly BrC emissions from the bottom-up emissions inventory in the NCP in January 2014 are the annual BrC emissions multiplied by the ratio of OA emissions in the NCP vs China and the ratio of OA emissions in January 2014 vs the whole year, resulting in a value of 65.5, 56.8 and 4.4 Gg for RCC, BB and FF-TRA, respectively. Finally, the proportion of the three primary types of emissions of BrC used in the model is 36.3 %, 100.8 % and 15.8 %, respectively. Figure 2 shows the contributing regions and burdens of the three separate primary sources of BrC.

Table 1 The data for primary BrC emissions calculation.

^a The BrC emissions from China's RCC in 2013 were reported by Sun et al. (2017) based on experiments involving seven types of coal being burned in four typical stoves in both chunk and briquette style. ^b The calculated BrC emissions from China's household biomass burning in 2013 were reported by Sun et al. (2021) using 11 widely used biomass types in China burned in a typical stove. ^c The estimated BrC emissions from vehicle exhaust in 2017 were 109 Gg reported by Wang et al. (2022a). In this study, the emissions of FF-TRA-derived BrC is 76.0 Gg, with a yearly scale factor of 0.70, which is derived by the annual level of civilian-owned motor vehicles between 2014 and 2017. ^d The value of BrC emissions in the NCP in January 2014 is additionally added with OA emitted from open biomass burning (6 Gg), which is assumed to be entirely light-absorbing. ^e These values were derived from the OA emission inventory described in Sect. 2.1.

Figure 2 Monthly BrC emissions burdens in January 2014 in the NCP from RCC, BB and FF-TRA.

SOA has also shown light absorption in the atmosphere (Lin et al., 2014). Laboratory experiments have revealed that most of the light-absorbing SOA is associated with aromatic SOA (Jacobson, 1999; Laskin et al., 2015; Li et al., 2020), and the absorption from biogenic SOA in the field has been found to be negligible (Washenfelder et al., 2015). Therefore, here we assume aromatic-derived SOA as secondary BrC in the model following previous studies (Jo et al., 2016; Wang et al., 2018).

Moreover, it is worth noting that both POA light absorption and SOA light absorption were shown to be dynamic, where BrC can be bleached when they undergo photodissociation (Forrister et al., 2015; Wong et al., 2019) or be darkened by cloud and fog processing of aerosols (Moise et al., 2015; Lin et al., 2017; Cheng et al., 2020). These processes are not considered in this study yet. More detailed parameterization of the chemical aging of BrC is needed in future BrC models.

2.3.2 BrC optical properties

The refractive indices of BrC as a function of wavelength are used for radiative transfer calculations. The complex refractive index ($m=n+ik$) of aerosol components enables practical implementation in the model, where n is the real part primarily associated with the scattering efficiency, and k is the imaginary part primarily associated with the absorption efficiency (Bohren and Huffman, 1998). The real part of BrC refractive index is the same of non-absorbing OA, which is fairly constrained, with reported values typically ranging from 1.5 to 1.7 (Saleh et al., 2014; Browne et al., 2019; Li et al., 2020). This study follows the study by Li et al. (2011a) and Wu et al. (2020). The imaginary part (k) of the BrC refractive index exhibits strong wavelength dependence, and the values range over several orders of magnitude (Saleh et al., 2018; Sengupta et al., 2018). Limited studies on the optical properties of BrC from fossil fuel combustions are reported at present. The average mass absorption efficiency (MAE) of RCC at 365 nm ranges from 0.80 to 2.47 m² g⁻¹ (Yan et al., 2017; Li et al., 2019; Tang et al., 2020; Wang et al., 2020; Ni et al., 2021; Song et al., 2021; Wang et al., 2021). In this study, as shown in Table 2, two sets of MAE are used for the sensitivity experiments of BrC. We choose a field optical measurement of BrC from all sources made by Zhang et al. (2022b) as the high absorption case (HI-BRC-ABS). The optical properties of BB and FF-TRA obtained in laboratory by Xie et al. (2017), as well as MAE of RCC and secondary BrC obtained in laboratory by Ni et al. (2021), are adopted as the low absorption case (LOW-BRC-ABS) in the study. The imaginary part of the two cases has shown wavelength-dependent light absorption properties and that the changes in anthropogenic emissions affect the optical properties of BrC. The imaginary part of both two cases is interpolated to 11 wavelengths to match the aerosol radiation calculation of the Goddard module in WRF-Chem. The value of k in this work is derived from the measured MAE using the following equation (Liu et al., 2013; Lu et al., 2015), as shown in Table 2:

$$\begin{array}{}\text{(6)}& k_{\mathrm{BrC},\mathit{\lambda }}={\displaystyle \frac{\mathit{\rho }\times \mathit{\lambda }\times {\mathrm{MAE}}_{\mathit{\lambda }}}{\mathrm{4}\mathit{\pi }}},\end{array}$$

where MAE_λ (m² g⁻¹) is the bulk mass absorption efficiency of BrC at the corresponding wavelength λ. ρ (g cm⁻³) is the density of organic aerosols, which is calculated to be 1.2 g cm⁻³ (Turpin and Lim, 2001) in this study.

Table 2 The refractive index of BrC used in the model.

“–” means not available.

2.3.3 Shortwave direct radiative effect calculation and experimental design

The shortwave DRE calculations of BrC follow the method reported by Chen et al. (2021) as shown in Eq. (7). The DRE of BrC is calculated by the difference between the net radiant flux with and without BrC, where the net radiant flux is the difference between the downward (F_↓) and upward radiant flux (F_↑).

$$\begin{array}{}\text{(7)}& {\mathrm{DRE}}_{\mathrm{TOA}}=\left(F_{\downarrow \mathrm{TOA}}^{\mathrm{a}}-F_{\uparrow \mathrm{TOA}}^{\mathrm{a}}\right)-\left(F_{\downarrow \mathrm{TOA}}^{\mathrm{0}}-F_{\uparrow \mathrm{TOA}}^{\mathrm{0}}\right),\end{array}$$

where DRE_TOA represent the shortwave DRE at the top of the atmosphere (TOA). F^a and F⁰ are the radiant flux with and without BrC aerosols, respectively.

An adjoint methodology proposed by Zhao et al. (2013) and Huang et al. (2015) has been used to diagnose the optical depth and DRE of BrC aerosols. Optical properties and radiative transfer of different sources BrC are calculated multiple times with one or a group of aerosol mass removed or without BrC absorption from each calculation following Eqs. (8) and (9). In addition, the model also takes into account the reduced aerosol masses along with the change in aerosol number concentration and size distribution.

$$\begin{array}{}\text{(8)}& {\displaystyle }\begin{array}{rl}& {\mathrm{AOD}}_{\left[\mathrm{species}i\right]}={\mathrm{AOD}}_{\left[\text{all species}\right]}\\ & -{\mathrm{AOD}}_{[\text{without species}i/\text{without species}i\text{absorption}]}\end{array}\text{(9)}& {\displaystyle }\begin{array}{rl}& {\text{DRE Forcing}}_{\left[\mathrm{species}i\right]}={\text{DRE Forcing}}_{\left[\text{all species}\right]}\\ & -{\text{DRE Forcing}}_{[\text{without species}i/\text{without species}i\text{absorption}]}\end{array}\end{array}$$

This method is more efficient than the traditional approach of running the model multiple times with the exclusion of a specific aerosol component. It not only saves computational time but also provides a more accurate estimation focused solely on the direct radiative effect of aerosols.

3 Results and discussions

3.1 Model performance

Before evaluating the DRE of BrC, results from the standard simulation are used to validate the model performance. Using available measurements, we first validate the spatial distribution and temporal variation of air pollutants (PM_2.5, O₃, NO₂, SO₂) in the NCP; the temporal variation of downward shortwave flux at the surface (SWDOWN) in Beijing, Tianjin, Zhengzhou, Hefei and Ji'nan; and the temporal variation of aerosol species (OA, elemental carbon, ammonium, sulfate and nitrite) in Beijing and Tianjin and of primary OA from BB, RCC, motor vehicles and SOA in Beijing in January 2014. Detailed data descriptions and quantitative statements of model biases can be found in Supplement Sects. S2.1 and S3. In general, the model simulates the air pollutants, SWDOWN, and aerosol species against measurements reasonably well.

SSA determines the strength of aerosols in absorbing solar radiation. Here we conduct three sensitivity experiments to evaluate the effect of BrC with different k values on the simulated aerosol absorption. The first experiment is the control simulation in which all organic aerosols are treated as purely scattering particles with no absorption contribution of BrC, which is referred to as NOBRC. The high-absorption scenario (HI-BRC-ABS) and low-absorption scenario (LOW-BRC-ABS) characterize BrC light absorption using the higher and lower imaginary refractive index derived from Sect. 2.3.2, respectively. Figure 3 shows the comparisons of simulated versus observed SSA at 440 nm (SSA₄₄₀) at Sun-sky radiometer Observation NETwork (SONET) sites in Beijing, Songshan, Xi'an, Hefei and Nanjing in January 2014. Due to the influence of clouds, the observational data from SONET are not continuous, resulting in a total of 237 valid data points being available for comparison. Moreover, SSA retrieval typically has larger uncertainties at low AOD values (Dubovik et al., 2002). Therefore, we have excluded the SSA data when AOD is less than 0.5, which leaves 206 valid points in each case. We find that the inclusion of BrC in the model reduces the bias of simulated SSA. The HI-BRC-ABS case demonstrated the largest improvement, with the correlation coefficient increasing to 0.54, making it the best simulation in the study. This suggests that the stronger BrC absorption case, as prescribed in HI-BRC-ABS, better captures the aerosol optical properties observed in northern China during winter. Consequently, the HI-BRC-ABS case can serve as the base simulation for further investigation of radiative effects of BrC in this study. Overall, the model tends to underestimate SSA₄₄₀. The underestimation might be partly caused by the overestimation of absorbing aerosols like BC or dust. Meanwhile, the uncertainties of the simulated SSA can be caused by other factors, such as mixing state of aerosols, particle shape, wavelength, and mass ratio of non-black carbon to BC (Liu et al., 2017; Jeong et al., 2020).

Figure 3 Scatter plot and linear fitting of modeled and observed column integrated SSA at 440 nm in case HI-BRC-ABS (red), LOW-BRC-ABS (blue) and NOBRC (black) case.

The daily AOD at 550 nm (AOD₅₅₀) from the dataset of Long-term Gap-free High-resolution Air Pollutant (LGHAP), derived via tensor-flow-based multimodal data fusion method (Bai et al., 2022), is compared with the simulation. This gap-free daily AOD dataset at 1 km resolution for 2000–2020 in China was generated by integrating multimodal data from satellites, numerical models, and in situ measurements. Data gaps in Moderate Resolution Imaging Spectroradiometer (MODIS) AOD are reconstructed through spatial pattern recognition and statistical knowledge transfer. Validation against Aerosol Robotic Network (AERONET) observations showed strong agreement, with an R of 0.91 and an RMSE of 0.21. Figure 4a and b show the pattern comparison of the monthly simulated and retrieved AOD₅₅₀. The model reproduces the retrieved AOD distribution in the NCP reasonably. The monthly average simulated and retrieved AOD₅₅₀ is 0.45 and 0.48 on average in the NCP, respectively. Figure 3c shows the scatter plot of the daily simulated and retrieved AOD₅₅₀ averaged in the NCP during the simulation period. The simulated daily average AOD₅₅₀ correlates quite well with the retrieved value, with a regression slope of 1.08 and correlation coefficient of 0.83. Generally, the retrieved and simulated AOD increases with deterioration of the particulate pollution. Figure 5 provides the pattern comparison of the simulated and Ozone Monitoring Instrument (OMI) retrieved AOD at 440 nm (AOD₄₄₀) averaged during the simulated episode. OMI aboard NASA's Aura satellite offers global atmospheric measurements at a spatial resolution of 0.25° × 0.25°, with Beijing's overpass occurring at approximately 13:45 local time. The averaged AOD₄₄₀ from the model simulation at 14:00 local time shows generally agreement with the OMI retrieval. The average simulated and retrieved AOD₄₄₀ is 0.51 and 0.53 in the NCP, respectively. Overall, the model generally performs well in simulating the AOD distribution. It is worth noting that the simulated AOD is not only dependent on the column aerosol content and composition but is also substantially influenced by relative humidity (RH), which determines the aerosol hygroscopic growth. Additionally, the satellite-retrieved AOD is subject to contamination by the presence of clouds, and considering the high occurrence frequency of clouds during haze days, the retrieved AOD might be overestimated (Satheesh et al., 2010; Chand et al., 2012).

Figure 4 (a) Monthly simulated AOD₅₅₀ of WRF-Chem, (b) monthly retrieved AOD₅₅₀ of reanalysis dataset LGHAP, and (c) scatter plot of the daily simulated and retrieved AOD₅₅₀ averaged in the NCP from 1 to 30 January 2014.

Figure 5 (a) Monthly simulated AOD₄₄₀ of WRF-Chem and (b) monthly retrieved AOD₄₄₀ of OMI in the NCP from 1 to 30 January 2014.

3.2 Surface mass concentrations of BrC in the NCP

The simulated distribution of average near-surface BrC concentrations and each source contribution in January 2014 is shown in Fig. 6a. In January, the monthly mean concentrations of BrC in the NCP vary from 0.05 to 42.3 µg m⁻³, with an average of 5.2 µg m⁻³. The spatial distribution of near-surface BrC concentrations is like that of PM_2.5 in the NCP, with the highest concentration areas located in Hebei Province, with an average concentration of 14.9 µg m⁻³. The simulated BrC concentrations are higher than those reported by Zhu et al. (2021) in 2018, which is perhaps caused by the more severe particulate pollution in January 2014.

Figure 6 Simulated mean surface concentrations of BrC (a) and contributions from RCC (b), BB (c), FF-TRA (d) and secondary BrC (e) in January 2014 in the NCP.

Figure 6b–e present the average near-surface BrC concentrations from different anthropogenic sources and secondary formation in January 2014 in the NCP. The BrC in the NCP predominantly originates from RCC and BB, with an average contribution of 2.3 and 2.4 µg m⁻³ and a maximum of 33.8 and 31.6 µg m⁻³, respectively. A relatively small proportion of BrC is related to FF-TRA and secondary transformations with an average concentration of 0.2 and 0.4 µg m⁻³, respectively. The BrC from RCC accounts for 55.8 % of total BrC concentrations in the NCP, which is highest in Beijing, Hebei and Tianjin, reaching 67.7 %, 54.4 % and 53.3 %, respectively. The BrC from BB constitutes 36.9 % of total BrC concentrations, with a contribution of about 40 % in most provinces of the NCP but only 26.5 % in Beijing. This result shows that the RCC is one of the major sources of BrC in the NCP due to the wide use of coal for heating and cooking in winter with low combustion efficiency and little emission control. The BrC emitted by RCC is mainly concentrated in the Beijing–Tianjin–Hebei (BTH) region in January 2014, while the BrC emitted by BB is distributed in the whole NCP. The Fen-Wei Plain exhibits notably high contributions from BB, which is consistent with the emission distribution and with previous studies (Cao and Cui, 2021; Zhang et al., 2021). The Fen-Wei Plain is one of the most densely populated and heavily polluted areas in northern China where biomass is usually used for heating during winter.

3.3 BrC absorption in the NCP

BrC absorbs visible to near-ultraviolet light with its absorption capabilities extending prominently at shorter wavelengths. Therefore, we calculate the AAOD (aerosol absorption optical depth) to evaluate the absorption contribution of BrC versus bulk aerosols, each anthropogenic source and SOA versus BrC at 365 nm (Fig. 7), by differentiating the AAOD between the model runs with and without the contribution of BrC. The average contribution of BrC to the total AAOD of aerosols at 365 nm is 16.4 %, and the maximum is 39.5 % in the NCP in January 2014. The BrC to BC and OC ratios in surface air in the study is 2.2 and 0.31, respectively. They are higher than the surface ratio used in the global models (Jo et al., 2016; Park et al., 2010) but lower than the ratio of Feng et al. (2013). Although the concentration of BrC is relatively high and the absorption of BrC in the ultraviolet band is comparable to that of BC, the imaginary index of BrC (about 0.1) is still much lower than that of BC (about 0.76). As a result, the light absorption contribution of BrC at 365 nm is not as significant during the study period.

Figure 7 Simulated monthly BrC AAOD versus aerosol at 365 nm (a) and each anthropogenic source of RCC (b), BB (c), FF-TRA (d) and secondary BrC (e) AAOD versus BrC in the NCP.

Figure 8 Estimated DRE_TOA of BrC (a) and DRE_TOA increase of OA owing to the absorption of BrC (b). The DRE_TOA of total OA with absorbing BrC (c) and DRE_TOA of total OA with BrC is assumed to be scattering (d). The averages of DRE are shown in the upper left of each panel.

Figure 9 The estimation DRE_TOA of BrC from each anthropogenic source RCC (a), BB (b), FF-TRA (c) and secondary BrC (d) in the NCP in January 2014.

Furthermore, the light absorption properties of BrC during the winter season are predominantly attributed to RCC, followed by BB, SOA and FF-TRA in the NCP. The average contribution of RCC and BB to the AAOD of BrC is 59.3 % and 26.3 %, respectively. Although the concentration of BrC from RCC is comparable to that from BB, the much higher light-absorbing property of the BrC from RCC makes it a dominant contributor to the overall light absorption caused by BrC. Despite lower surface concentrations compared to primary BrC, secondary BrC contributes significantly to AAOD of BrC, averaging ∼ 10.0 % with elevated contributions in the sea and remote regions, which is likely due to the highly oxidized character of organic aerosols and its chemical aging in aging air masses, leading to the formation of BrC (Gouw et al., 2005; Kawamura et al., 2005; Tsigaridis and Kanakidou, 2018). While AOD represents column-integrated concentrations, the ratio of secondary BrC to primary BrC increases from 8.9 % at the surface to 12.0 % of atmospheric burden. It reaches 14.3 % at an altitude of 500 m as shown in Fig. S6, which could lead to a higher absorption contribution of secondary BrC (Wang et al., 2022b). Moreover, the observations indicate that a substantial amount of SOA is water-soluble (Maria et al., 2003; Peng et al., 2021) and therefore treated as hygroscopic components in the model, and its absorption could be magnified. The AAOD contribution of BrC from vehicles is the lowest, about 1.8 % on average, but its contribution in southern Jiangsu and southeast Anhui ranges from 4 % to 8 %, higher than in other regions. This may suggest that motor vehicle emissions account for a significant proportion of pollution in this area.

3.4 Direct radiative effect of BrC

Figure 8a shows the distribution of the average DRE_TOA caused by BrC at the top of the atmosphere during the episode. The DRE_TOA of BrC in the NCP is −0.09 W m⁻² on average, with a maximum of +0.46 W m⁻² and a minimum of −2.74 W m⁻². Compared to the average DRE_TOA of BC +3.9 W m⁻² and a maximum of +21.6 W m⁻², the average DRE_TOA of BrC in the NCP is close to zero. However, in terms of the spatial distribution, the DRE_TOA of total BrC in the NCP is predominantly negative, especially in those areas with high BrC concentrations including the BTH area and Fen-Wei Plain. The largest negative DRE_TOA is in Shanxi Province, where the highest contributor to BrC is BB. The overall scattering effect of BrC is greater than its absorption effect, so BrC has a net cooling effect. The solar irradiance in the UV band contributes only 10 % of the total solar irradiation, and the imaginary refractive indices of BrC are much lower than those strongly absorbing BC, especially in the visible band. On the other hand, although this leads to a small heating effect by BrC, the increased DRE_TOA induced by BrC, which is usually considered its scattering effect, is up to an average of +0.40 W m⁻² and a maximum of +1.83 W m⁻² as shown in Fig. 8b. The DRE_TOA of OA without BrC is −2.00 W m⁻² over the NCP (Fig. 8d) and is increased to −1.60 W m⁻² (Fig. 8c) after considering BrC. Consequently, the cooling effect of OA is reduced by BrC by an average of 28.0 %. So far, almost all estimates of the radiation effects of BrC are on a global basis. The study of Wang et al. (2014) reported a global DRE_TOA of BrC of −0.02 W m⁻², resulting in a DRE_TOA increase of +0.07 W m⁻². Brown et al. (2018) have also reported a global annual increased DRE_TOA of +0.13 W m⁻², with the maximum forcing ($\sim +\mathrm{1.75}$ W m⁻²) occurring on the west coast of southern Africa. BrC reduces the cooling effect caused by organic aerosols by approximately 16 % (Jo et al., 2016). However, all these estimates only considered BB, biofuels or SOA and the significance of BrC in the estimates for these studies typically based on the assumption that OA primarily scatters sunlight. Our results indicate that the solar radiation changes caused by BrC in the NCP are notable. The absorption effect of BrC should be considered in climate models to accurately assess the aerosol impact on atmospheric heating and climate change.

Figure 9a–d show the estimated DRE_TOA of BrC from RCC, BB, FF-TRA and secondary BrC in the NCP during the episode. Similar to the contribution of BrC sources to the AAOD at 365 nm, the most important source contributing to DRE_TOA of BrC is RCC (+0.26 W m⁻²), followed by BB (+0.11 W m⁻²), secondary BrC (+0.02 W m⁻²) and FF-TRA (+0.01 W m⁻²) in the NCP, as shown in Fig. 9. In addition, the DRE_TOA of BrC from various sources exhibits distinct spatial distribution characteristics in the NCP. The highest DRE_TOA of BrC from fossil sources, which include RCC and FF-TRA, is predominantly concentrated in Hebei. The highest positive DRE_TOA of BrC from BB is found in the Fen-Wei Plain. Meanwhile, the secondary BrC contributes most to the DRE_TOA in the south part of the NCP. The persistent high values over southern China might stem from the model's representation of secondary BrC as hygroscopic components (Peng et al., 2021), whose light-absorbing capacity is amplified in the region with high ambient humidity.

4 Conclusions

A source-explicit BrC simulation in January 2014 in the NCP is conducted by using the WRF-Chem model. We define the BrC based on varying proportions of RCC, BB, FF-TRA and SOA sources, assigning each a distinct imaginary refractive index to represent their differing light absorption characteristics. Model simulations are evaluated with various data sets. Besides the good reproduction of temporal and spatial variations of aerosol components and SWDOWN in the model, AOD and SSA are used to evaluate the aerosol optical properties.

Near-surface mass concentrations of simulated BrC in the NCP range from 0.05 to 42.3 µg m⁻³ with an average of 5.2 µg m⁻³, which is mainly attributed to RCC and BB, especially in the BTH region and Fen-Wei Plain. Estimation of the BrC contribution to AAOD shows that the BrC accounts for an average of 16.4 % and up to 39.5 % of the total aerosol absorption at 365 nm. The largest contributor to the absorption of BrC is RCC-derived BrC, reaching 59.3 %. BrC generally has a net cooling effect in the NCP if we consider both the absorption and scattering properties, with DRE_TOA of BrC of around −0.09 W m⁻² on average and ranging between −2.74 and +0.46 W m⁻². However, the absorption of BrC increases the DRE_TOA of OA by 28.0 % with an average of +0.40 W m⁻² and a maximum of +1.83 W m⁻². The average increase in DRE_TOA of BrC from RCC, BB, secondary formation and FF-TRA is +0.26, +0.11, +0.02 and +0.01 W m⁻², respectively. Our results indicate that BrC derived from RCC may have significant implications for regions relying heavily on coal as their primary energy source, such as northern China. Climate models should not only incorporate the absorption of BrC but also account for residential coal burning as a potentially important BrC emission source. Additionally, although we conducted simulations with a relatively conservative secondary BrC, the impact of secondary BrC on radiative processes should not be overlooked. More field observation and model experiments should be carried out in the future for better understanding of its role in atmospheric radiation balance.

It should be noted that China has started to switch from coal to cleaner and more efficient energy such as natural gas or liquid petroleum gas in recent years. According to the latest report of National Bureau of Statistics of China, the total coal consumption for residential use was 55.5 Gg in 2022 (https://data.stats.gov.cn, last access: 10 March 2025) with a 40.3 % decrease compared to 2014. Therefore, our diagnosis of the sources of BrC and their radiative effects is specifically targeted at the winter season in 2014. Moreover, future simulations should strengthen the parameterization for the evolution of BrC, such as bleaching or darkening processes.

The simulation of BrC in climate models is fraught with uncertainties due to its diverse sources which result in a wide range of optical properties. The absorption characteristics of BrC can change significantly as it undergoes atmospheric aging, impacting its radiative forcing estimates. Additionally, interactions of BrC with other atmospheric particles and its effects on cloud microphysics and albedo introduce further complexities in modeling its climate impact. These uncertainties necessitate enhanced observational data and better integration of BrC properties in climate models to improve the accuracy of climate predictions and assessments. Although the simulations have been evaluated with extensive aerosol mass and optical measurements, more field measurements and lab experiments are needed, especially for the inventory development of BrC, the vertical profiles from aircraft and physicochemical properties, which would be useful for further evaluating and improving model performance.