Modeling of polycyclic aromatic hydrocarbons (PAHs) from global to regional scales: model development (IAP-AACM_PAH v1.0) and investigation of health risks in 2013 and 2018 in China

Polycyclic aromatic hydrocarbons (PAHs) significantly impact human health due to their persistence, toxicity, and potential carcinogenicity. Their global distribution and regional changes caused by emission changes, especially over areas in developing countries, remain to be understood along with their health impacts. This study implemented a PAH module in the global–regional nested Atmospheric Aerosol and Chemistry Model of the Institute of Atmospheric Physics (IAP-AACM) to investigate the global distribution of PAHs and the change in their health risks from 2013 to 2018 in China. An evaluation against observations showed that the model could capture well the spatial distribution and seasonal variation in Benzo[a]pyrene (BaP), the typical indicator species of PAHs. On a global scale, the annual mean concentrations are the highest in China followed by Europe and India, with high values exceeding the target values of 1 ng m⁻³ over some areas. Compared with 2013, the concentration of BaP in China decreased in 2018 due to emission reductions, whereas it increased in India and southern Africa. However, the decline is much smaller than for PM_2.5 during the same period. The concentration of BaP decreased by 8.5 % in Beijing–Tianjin–Hebei (BTH) and 9.4 % in the Yangtze River Delta (YRD). It even increased over areas in the Sichuan Basin due to changes in meteorological conditions. The total incremental lifetime cancer risk (ILCR) posed by BaP only showed a slight decrease in 2018, and the population in east China still faced significant potential health risks. The results indicate that strict additional control measures should be taken to reduce the pollution and health risks of PAHs effectively. The study also highlights the importance of considering changes in meteorological conditions when evaluating emission changes from concentration monitoring.

1 Introduction

Polycyclic aromatic hydrocarbons (PAHs) are aromatic compounds with two or more aromatic rings. PAHs have been categorized as persistent organic pollutants (POPs) by the United Nations Economic Commission for Europe's (UNECE's) Convention on Long-Range Transboundary Air Pollution (CLRTAP) (Friedman and Selin, 2012), and they are widely distributed in the environment through atmospheric transport. PAHs have attracted significant attention in environmental research and risk assessment due to their persistence, toxicity, and potential carcinogenicity (Chen and Liao, 2006; Shen et al., 2014). These compounds are generated from both natural and anthropogenic sources (Haritash and Kaushik, 2009). Volcanic eruption, forest, and prairie fire are the major natural sources of atmospheric PAHs (Baek et al., 1991). Anthropogenic sources are the most important source of PAHs, including incomplete combustion of fossil fuels and biomass (Li et al., 2022; Ravindra et al., 2008).

Understanding the sources, distribution, and fate of PAHs is crucial for assessing their impacts on human health and the environment. Upon emission into the atmosphere, PAHs are redistributed by gas–particulate partitioning, gaseous-phase reactions, heterogeneous reactions, air–soil exchange, and wet/dry deposition during long-range transport (LRT; Inomata et al., 2013). Monitoring is the most commonly used way to investigate the concentration of PAHs in the atmosphere. Due to the high costs of observation and technical limitations, it is difficult to conduct a long-term and broad regional analysis through monitoring (Zhen, 2023). Up to now, there have been few continuous observations over the major continents at the same time (Dong et al., 2023). A transport model is an effective tool to simulate the distribution of PAHs and their LRT, which can greatly enhance our understanding of the distribution of PAHs on a regional and global scale (Byun and Schere, 2006; Wang et al., 2021).

As recently outlined by Galarneau et al. (2014), several numerical modeling studies have been reported in the literature. The models that can simulate PAHs include but are not limited to the following examples: GEOS-Chem (Friedman et al., 2014; Friedman and Selin, 2012), ECHAM5 (Lammel and Sehili, 2007; Lammel et al., 2009; Lammel et al., 2015; Octaviani et al., 2019), CAM5 (Lou et al., 2023; Shrivastava et al., 2017), and MOZART-4 (Shen et al., 2014). The horizontal resolutions of these reported models are primarily at 4° × 5° and 2.8° × 2.8°. Shen et al. (2014) simulated the transport of Benzo[a]pyrene (BaP), one of the most toxic and highly carcinogenic PAHs, in the global troposphere based on MOZART-4, and they showed that the model resolution was crucial for the health risk assessment. Lammel et al. (2015) demonstrated the significant impact of gas–particle partitioning mechanisms on the atmospheric lifetime, compartment distributions, and LRT of PAHs. The regional modeling studies focusing on Europe, east Asia, and North America have also been reported, with horizontal resolutions ranging mainly from 54 km × 54 km to 24 km × 24 km – CMAQ (Aulinger et al., 2009; Aulinger et al., 2007; Bieser et al., 2012; San José et al., 2013; Efstathiou et al., 2016), WRF-Chem (Mu et al., 2018), AURAMS (Galarneau et al., 2014), and CanMETO (Zhang et al., 2011a, b, 2009). Efstathiou et al. (2016) showed that considering absorption and adsorption processes can better capture the concentration levels and seasonal variations in BaP. In recent years, the effect of the heterogeneous reaction process of PAHs on transportation has also been studied. Mu et al. (2018) developed a new kinetic scheme describing the effects of temperature and humidity on the organic aerosol coating of BaP and BaP reaction rate. They found that low temperature and low humidity can significantly increase the lifetime of BaP and enhance its LRT capacity.

However, the resolutions and spatial range differed greatly between these models. Most of the models are either global or regional. There is a lack of simulation studies focusing on both global and key regions, making it difficult to investigate a specific focus in a global background in a consistent manner. Additionally, the resolution of most global models is low, which will further affect the health risk assessment of PAHs. Furthermore, the up-to-date mechanisms (gas–particle partitioning, heterogeneous reaction, and air–soil exchange) established for PAHs simulations are not considered in earlier modeling studies.

China is one of the largest PAH-emitting countries in the world (Inomata et al., 2012; Zhang and Tao, 2009). High concentrations of BaP have been reported (Bieser et al., 2012; Liu et al., 2014; Shrivastava et al., 2017; Su et al., 2023). Over the polluted regions in eastern China, annual concentrations of BaP exceeded 1 ng m⁻³, the target values proposed in the European Union (EU) and China. To improve air quality and protect public health, the State Council of China promulgated the Action Plan on Air Pollution Prevention and Control (hereafter the Action Plan) in 2013. This Action Plan established many effective emission reduction and energy-saving policies, such as strengthening industrial emission standards, eliminating outdated polluting industries, upgrading industrial boilers, and developing clean fuels in the residential sector. Since then, many studies have investigated the changes in concentration levels and health risks of conventional pollutants, such as PM_2.5 (Feng et al., 2019; Wang et al., 2018; Zhang et al., 2019; Zhu et al., 2021; Wang et al., 2019). Wang et al. (2019) pointed out that the annual average concentrations of PM_2.5 in the Beijing–Tianjin–Hebei (BTH) region, the Yangtze River Delta (YRD), and the Pearl River Delta (PRD) all decreased by more than 27 % in 2017, indicating that the control measures have achieved remarkable effects and the air quality has been significantly improved. However, for non-conventional pollutants, such as BaP and other PAHs, their concentration changes due to emission reduction in China after implementing of policies have not been quantified. The changes in health risks and the benefits from control measures have not yet been assessed.

Considering the aforementioned, we simulated PAHs from global to regional scales by coupling the key physical and chemical modules associated with PAHs in a global–regional nested atmospheric transport model. In particular, newly established parameterizations of gas–particle partitioning and heterogeneous reaction were incorporated into the model. Then the changes in global concentration and health risks of BaP over China were quantified based on model evaluation against a collected observation dataset. The study can advance our understanding of global PAH distribution and regional health risks and their responses to emission change. The paper is arranged as follows: Sect. 2 briefly describes the host model (the Atmospheric Aerosol and Chemistry Model of the Institute of Atmospheric Physics, IAP-AACM), the physical and chemical modules related to PAHs, and the method of assessing health risks. Section 3 presents the configuration of the model and the observations used in the evaluation. Section 4 shows the global and regional distributions of BaP concentrations and analyzes the health risks associated with BaP in China. Section 5 discusses the uncertainty of the model. In Sect. 6, the main conclusions are summarized.

2 Model description and development

2.1 Description of host model

The model used in this study is the Atmospheric Aerosol and Chemistry Model developed by the Institute of Atmospheric Physics, Chinese Academy of Sciences (IAP-AACM) (Wei et al., 2019), which was developed based on the Global Nested Air Quality Prediction Modeling System (GNAQPMS; Chen et al., 2015; Wang et al., 2001). IAP-AACM is a 3D Eulerian transport model that uses a multi-scale domain-nesting approach to simulate atmospheric chemistry and aerosol processes from global to regional scales. As recently described by Chen et al. (2015), compared with the traditional multi-scale modeling methods (Seigneur et al., 2001), the online nesting method uses the same parameters in the global and regional domains, which avoids uncertainties caused by different boundary conditions, and it also provides boundary conditions at higher time resolution (Zhang et al., 2012b; Chen et al., 2015), thus improving the performance of the model on the regional scale.

This model includes emission, horizontal and vertical advection (Walcek and Aleksic, 1998), diffusion (Byun and Dennis, 1995), dry deposition (Zhang et al., 2003), gaseous chemistry (CBM-Z, Carbon Bond Mechanism version Z, Zaveri and Peters, 1999), heterogeneous chemistry (Li et al., 2012), aqueous reactions in clouds, and wet scavenging (Stockwell et al., 1990). It has been successfully and widely applied to simulate the spatial–temporal distribution characteristics of gaseous pollutants, aerosol components, and long-distance transportation of mercury (Chen et al., 2015; Chen et al., 2014; Wei et al., 2019; Ye et al., 2021; Du et al., 2019). In addition, advanced particle microphysics (APM) has been incorporated to simulate new particle formation processes and predict the particle number concentrations on global and regional scales (Chen et al., 2021).

2.2 Development of the PAH module

The PAH processes in the IAP-AACM model include gaseous-phase reaction, heterogeneous reaction, gas–particle partitioning, air–soil exchange, dry deposition, and wet scavenging. The simulated species include BaP, Benzo[b]fluoranthene (BbF), Benzo[k]fluoranthene (BkF), and Indeno[1,2,3-cd]pyrene (IcdP) in the gas and particulate phases (Wu et al., 2024a). In this study, we mainly focus on BaP due to its highly carcinogenic nature and the relatively rich observations.

2.2.1 Gaseous-phase reactions

PAHs are degraded through reactions with various atmospheric oxidants such as hydroxyl radical (OH), nitrate radical (NO₃), and ozone (O₃) in the troposphere (Lammel and Sehili, 2007). Among these oxidants, the reactions with OH are considered to be the most important pathway for the removal of PAHs. The nighttime reaction of PAHs with NO₃ is also important in the atmosphere (Keyte et al., 2013). Therefore, reactions of gaseous-phase BaP with OH, NO₃, and O₃ are all considered in the model. The second-order rate coefficients are 5.0 × 10⁻¹¹, 5.4 × 10⁻¹¹, and 2.6 × 10⁻¹⁷ cm³ molecules⁻¹ s⁻¹, respectively (Inomata et al., 2013; Finlayson-Pitts and Pitts, 2000; Klöpffer et al., 2007).

2.2.2 Heterogeneous reaction

In the case of BaP, the heterogeneous reaction with O₃ is considered to be the dominant loss (Finlayson-Pitts and Pitts, 2000; Efstathiou et al., 2016). Studies have shown that the process of heterogeneous reaction can be well described by the Langmuir–Hinshelwood mechanism (Kahan et al., 2006; Kwamena et al., 2007), in which BaP is adsorbed to the surface while the O₃ is in phase equilibrium. The first-order reaction rate coefficient k (s⁻¹) of the Langmuir–Hinshelwood mechanism is as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }k={\displaystyle \frac{k_{\max }{K}_{{\mathrm{O}}_{\mathrm{3}}}\left[{\mathrm{O}}_{\mathrm{3}}\right]}{\mathrm{1}+K_{{\mathrm{O}}_{\mathrm{3}}}\left[{\mathrm{O}}_{\mathrm{3}}\right]}},\text{(2)}& {\displaystyle \frac{\partial C}{\partial t}}=-K_{{\mathrm{O}}_{\mathrm{3}}}\left[{\mathrm{O}}_{\mathrm{3}}\right],\end{array}$$

where k_max is the maximum rate coefficient and the value is 0.060 ± 0.018 s⁻¹. [O₃] is the concentration of O₃ (mol cm⁻³). $K_{{\mathrm{O}}_{\mathrm{3}}}$ is the O₃-to-surface equilibrium constant (0.028 ± 0.014 × 10⁻¹³ cm³).

In addition, we incorporated a more detailed parameterization (ROI-T) developed by Mu et al. (2018) based on the Langmuir–Hinshelwood mechanism. The scheme emphasizes the importance of representing the dependence of degradation on temperature and humidity when coated by organic aerosols The first-order reaction rate coefficient k (s⁻¹) is given by Eq. (3):

$$\begin{array}{}\text{(3)}& k=\mathrm{base}+{\displaystyle \frac{\max -\mathrm{base}}{\mathrm{1}+(\frac{\mathrm{xhalf}}{\left[{\mathrm{O}}_{\mathrm{3}}\right]}{)}^{\mathrm{rate}}}},\end{array}$$

where base, max, rate, and xhalf are all the parameterizations of the heterogeneous reaction, with specific values shown in Mu et al. (2018). In our study, we coupled these two parameterizations as two options for O₃ degradation by heterogeneous reaction in IAP-AACM. The model results using these two schemes were compared to analyze the influence of heterogeneous reaction schemes on BaP concentration. The ROI-T scheme was used as the default in this study.

2.2.3 Gas–particle partitioning

The partition of compounds between the gas and particulate phases is parameterized with the gas–particle partitioning coefficient (K_P; m³ µg⁻¹) (Harner and Bidleman, 1998):

$$\begin{array}{}\text{(4)}& K_{\mathrm{P}}=\left({\displaystyle \frac{{\left[\mathrm{PAH}\right]}_{\mathrm{p}}}{\left[\mathrm{TSP}\right]}}\right)/{\left[\mathrm{PAH}\right]}_{\mathrm{g}},\end{array}$$

where [PAH]_g and [PAH]_p are the concentrations of PAHs in the gas and particulate phase (µg m⁻³) and [TSP] is the concentration of total suspended particles (TSP; µg m⁻³) in the atmosphere (µg m⁻³).

Adsorption onto black carbon (BC) and absorption into aerosol organic matter (OM) are two important mechanisms of gas–particle partitioning (Odabasi et al., 2006). Therefore, we use the gas–particle partition coefficient equation, which was derived by Dachs and Eisenreich (2000), to represent these two mechanisms:

$$\begin{array}{}\text{(5)}& K_{\mathrm{P}}=\left[{\displaystyle \frac{\left(f_{\mathrm{OM}}{\mathrm{MW}}_{\mathrm{OCT}}{\mathit{\delta }}_{\mathrm{OCT}}\right)K_{\mathrm{OA}}}{\left({\mathit{\rho }}_{\mathrm{OCT}}{\mathrm{MW}}_{\mathrm{OM}}{\mathit{\delta }}_{\mathrm{OM}}{\mathrm{10}}^{\mathrm{12}}\right)}}\right]+\left[\left({\displaystyle \frac{f_{\mathrm{BC}}{a}_{\mathrm{BC}}{K}_{\mathrm{SA}}}{a_{\mathrm{AC}}{\mathrm{10}}^{\mathrm{12}}}}\right)\right],\end{array}$$

where MW_OCT and MW_OM are the mean molecular weights of octanol and the OM phase (g mol⁻¹) and δ_OCT and δ_OM are the activity coefficient of the absorbing compound in octanol and the OM phase, respectively. f_OM and f_BC are the mass fractions of the OM phase on TSP and the BC in the aerosol. ρ_OCT is the density of octanol (0.820 kg L⁻¹). a_BC and a_AC are the specific surface areas of BC (62.7 m² g⁻¹, Jonker and Koelmans, 2002) and activated carbon (AC), respectively. In this study, we use the same assumption as Odabasi et al. (2006) (${\mathrm{MW}}_{\mathrm{OCT}}/{\mathrm{MW}}_{\mathrm{OM}}=\mathrm{1}$, ${\mathit{\delta }}_{\mathrm{OCT}}/{\mathit{\delta }}_{\mathrm{OM}}=\mathrm{1}$, and $a_{\mathrm{BC}}/a_{\mathrm{AC}}=\mathrm{1}$).

$$\begin{array}{}\text{(6)}& \log K_{\mathrm{OA}}=A+B/\left(T\right),\end{array}$$

where K_OA is the octanol–air partitioning coefficient (temperature dependent). T is the temperature (K). The values of A and B are 5382 and −6.5, respectively (Odabasi et al., 2006).

$$\begin{array}{}\text{(7)}& {\displaystyle }\log P_{\mathrm{L}}=m_{\mathrm{L}}(T{)}^{-\mathrm{1}}+b_{\mathrm{L}},\text{(8)}& {\displaystyle }\log K_{\mathrm{SA}}=-\mathrm{0.85}\log P_{\mathrm{L}}+\mathrm{8.94}-\log \left({\displaystyle \frac{\mathrm{998}}{a_{\mathrm{BC}}}}\right),\end{array}$$

where P_L is the supercooled liquid vapor pressure (Pa). The values of b_L and m_L are 12.59 and −5252, respectively (Dachs and Eisenreich, 2000). K_SA is the soot–air partitioning coefficient (L kg⁻¹), which is a function of P_L and a_BC (van Noort, 2003).

2.2.4 Air–soil exchange

The semi-volatility and persistence of PAHs allow them to dynamically exchange between the atmosphere and soil by deposition and re-volatilization from ground surfaces (Semeena and Lammel, 2005). These processes can affect the distribution and long-distance transport of PAHs in the environment. As described by Hansen et al. (2004), air–soil exchange is parameterized following Strand and Hov (1996), which is based on Jury et al. (1983). Here, soil (standard soil) is considered to be a homogeneous layer of thickness z_s=0.15 m, and standard values and chemical properties are provided by Jury et al. (1983) (Table S1). The differential equation for the change in concentrations in soil and air can be expressed by Eqs. (9) and (10):

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{\partial c_{\mathrm{s}}}{\partial t}}={\displaystyle \frac{\mathrm{1}}{z_{\mathrm{s}}}}\left(F_{\mathrm{exc},\mathrm{soil}}+F_{\mathrm{wet}}\right)-k_{\mathrm{soil}}{c}_{\mathrm{s}},\text{(10)}& {\displaystyle \frac{\partial c_{\mathrm{a}}}{\partial t}}=-{\displaystyle \frac{\mathrm{1}}{z_{\mathrm{a}}}}{F}_{\mathrm{exc},\mathrm{soil}},\end{array}$$

where c_a and c_s are the concentrations of PAHs in the atmosphere and soil, respectively. z_a is the lowest atmospheric layer depth (m), and F_wet is the wet deposition flux (mol s⁻¹ m⁻²). k_soil is the degradation rate in soil, which is estimated to be 2.2 × 10⁻⁸ s⁻¹ (Finlayson-Pitts and Pitts, 2000; Klöpffer et al., 2007; Lammel et al., 2009). The air–soil exchange flux (F_exc,soil) is given by Eq. (11):

$$\begin{array}{}\text{(11)}& F_{\mathrm{exc},\mathrm{soil}}=K_{\mathrm{a}/\mathrm{s}}\left(c_{\mathrm{a}}-{\displaystyle \frac{C_{\mathrm{s}}}{K_{\text{soil-air}}}}\right).\end{array}$$

K_soil-air is the partitioning coefficient between soil and air, which is given by Karickhoff (1981):

$$\begin{array}{}\text{(12)}& K_{\text{soil-air}}=\mathrm{4.11}\times {\mathrm{10}}^{-\mathrm{4}}\times {\mathit{\rho }}_{\mathrm{s}}{f}_{\mathrm{oc}}{K}_{\mathrm{OA}},\end{array}$$

where f_OC is the fraction of OC in soil and 4.11 × 10⁻⁴ is a constant in units of m³ kg⁻¹. ρ_s is the density of soil. $K_{\mathrm{a}/\mathrm{s}}$ is the overall exchange velocity (m s⁻¹), which can be estimated by Eq. (13) (Strand and Hov, 1996) as follows:

$$\begin{array}{}\text{(13)}& \begin{array}{rl}& K_{\mathrm{a}/\mathrm{s}}\\ & ={\displaystyle \frac{D_{\mathrm{G}}^{\mathrm{air}}{a}^{\mathrm{10}/\mathrm{3}}(\mathrm{1}-l-a{)}^{-\mathrm{2}}+D_{\mathrm{L}}^{\mathrm{water}}{l}^{\mathrm{10}/\mathrm{3}}{K}_{\mathrm{WA}}(\mathrm{1}-l-a{)}^{-\mathrm{2}}}{z_{\mathrm{s}}/\mathrm{2}}},\end{array}\end{array}$$

where $D_{\mathrm{G}}^{\mathrm{air}}$ and $D_{\mathrm{L}}^{\mathrm{water}}$ are the air and liquid diffusion coefficient (m² s⁻¹), respectively. K_WA is the water–air partitioning coefficient. The differential equation is solved in ODEPACK (https://github.com/jacobwilliams/odepack, last access: 20 January 2024).

2.2.5 Dry and wet deposition

PAHs can be removed from the atmosphere and enter terrestrial ecosystems through dry and wet deposition (Cao et al., 2021). Dry deposition and wet scavenging have been included in IAP-AACM. For the gaseous species of PAHs, their wet scavenging is assumed to be the same as xylene in the CBM-Z mechanism, which is also an aromatic hydrocarbon like BaP; for the PAHs in the particle phase, these two processes are treated similarly to that of organic aerosol.

2.3 Risk assessment

The incremental lifetime cancer risk (ILCR) is widely used to calculate the risk of human exposure to PAHs (Nam et al., 2021). The carcinogenic risk of PAHs to humans through different exposure routes was calculated based on the health risk evaluation model proposed by the US Environmental Protection Agency (EPA) (Smith et al., 2007).

The national population data at 1 km × 1 km (at the Equator) resolution in 2013 and 2018 were obtained from LandScan (Oak Ridge National Laboratory; database can be accessed via https://landscan.ornl.gov, last access: 20 January 2024) and re-gridded to 1° × 1° and 0.33° × 0.33° to match the model resolution.

2.3.1 Daily exposure dose

Dermal contact and inhalation are regarded as the major routes of human exposure to BaP (Li et al., 2010; Ma et al., 2020; Zhang et al., 2016). In this study, the health risk for the entire population and three groups (adult women, adult men, and children) are calculated. The average daily exposure dose (ADD) to PAHs through the two exposure routes is calculated as follows:

$$\begin{array}{}\text{(14)}& {\displaystyle }{\mathrm{ADD}}_{\mathrm{der}}={\displaystyle \frac{C\times \mathrm{SA}\times \mathrm{ABS}\times \mathrm{AF}\times \mathrm{EF}\times \mathrm{ED}}{\mathrm{AT}\times \mathrm{BW}}},\text{(15)}& {\displaystyle }{\mathrm{ADD}}_{\mathrm{inh}}={\displaystyle \frac{C\times \mathrm{IR}\times \mathrm{EF}\times \mathrm{ED}}{\mathrm{AT}\times \mathrm{BW}}},\end{array}$$

where ADD_der and ADD_inh are the average daily exposure doses that enter the body through dermal contact and inhalation, respectively (ng kg⁻¹ d⁻¹), C is the concentration of PAHs (ng m⁻³). IR is the inhalation rate (m³ d⁻¹). AF is the dermal adherence rate (mg cm⁻² d⁻¹). EF and ED are the exposure duration (d yr⁻¹) and period (year), respectively. BW is the body weight (kg). SA is the skin exposed surface area (cm²). ABS is the skin absorption factor. AT is the average exposure time (d). The values are shown in Table S2.

2.3.2 Incremental lifetime cancer risk (ILCR)

The ILCR was calculated based on the ADD:

$$\begin{array}{}\text{(16)}& {\displaystyle }{\mathrm{ILCR}}_{\mathrm{der}}={\mathrm{ADD}}_{\mathrm{der}}\times {\mathrm{SFO}}_{\mathrm{der}}\times \mathrm{CF},\text{(17)}& {\displaystyle }{\mathrm{ILCR}}_{\mathrm{inh}}={\mathrm{ADD}}_{\mathrm{inh}}\times {\mathrm{SFO}}_{\mathrm{inh}}\times \mathrm{CF},\text{(18)}& {\displaystyle }\mathrm{TILCR}={\mathrm{ILCR}}_{\mathrm{der}}+{\mathrm{ILCR}}_{\mathrm{inh}},\end{array}$$

where ILCR_der and ILCR_inh are lifetime cancer risks through dermal contact and inhalation, respectively. TILCR is the total lifetime cancer risk of exposure through the two pathways. SFO is the cancer slope factor (kg d mg⁻¹) and CF is the conversion factor. The values are shown in Table S2. For a carcinogen, an ILCR of less than 1 × 10⁻⁶ indicates negligible cancer risk, an ILCR between 1 × 10⁻⁶ and 1 × 10⁻⁴ indicates potential cancer risk, and an ILCR larger than 1 × 10⁻⁴ indicates high potential cancer risk.

3 Experiment setup and observation data

3.1 Experiment setup

In this study, we used two nested domains covering the whole globe and east Asia as shown in Fig. S1. The horizontal resolutions are 1° × 1° and 0.33° × 0.33°, respectively. A total of 20 vertical layers are used in IAP-AACM. The first layer of the model is approximately 50 m deep, and the top layer extends to 20 km. The simulation results from 1 January to 31 December 2013 and from 1 January to 31 December 2018 were used for analysis. Each simulation had a 1-month spin-up before 1 January to reduce the influence of initial conditions. The global version of the Weather Research and Forecasting model (WRF, version 3.7.1) (Zhang et al., 2012a; Skamarock et al., 2008) provides the meteorological fields to drive the IAP-AACM. The initial and boundary conditions of the global WRF were produced by Final Analysis (FNL) data from the National Centers for Environmental Prediction (NCEP).

The emission inventory of BaP in 2013 and 2018 was derived from the Emissions Database for Global Atmospheric Research (EDGAR; Crippa et al., 2020, available from https://edgar.jrc.ec.europa.eu/dataset_pop60#sources, last access: 15 December 2023). We mainly analyzed the results using EDGAR emissions, which mainly includes anthropogenic sources such as power, transportation, industrial, agricultural, and energy for buildings. An additional simulation for 2013 using the emission inventory developed by the research group of Peking University (PKU) (http://inventory.pku.edu.cn, last access: 10 February 2023) was used to investigate the uncertainties from emissions. The resolution of both emission inventories is 0.1° × 0.1°. Therefore, we re-gridded the emission inventories to match the model grids at 1° × 1° and 0.33° × 0.33° resolution.

Figure 1 Spatial distributions of total emissions of BaP in (a) 2013 and (b) 2018 based on the EDGAR inventory. (c) BaP emissions for 16 regions (except oceans, Arctic, and Antarctic) in 2013 and 2018.

Figure 2 Comparison of simulated (sim_BaP) and observed (obs_BaP) annual mean concentrations of BaP in Asia (red), Europe (blue), United States (brown), and Canada (green) in 2013. The solid black line shows a ratio of 1:1 and the dashed gray lines show ratios of 5:1, 2:1, 1:2, and 1:5.

Figure 3 Comparison of the BaP annual mean simulated (red) concentrations with observed (black) values at European sites in (a) 2013 and (b) 2018.

The global total emissions of BaP in 2013 and 2018 are shown in Fig. 1a and b, respectively. The annual emissions in different regions (Fig. S1) were also calculated (Fig. 1c). The global emissions of BaP were 7166.9 t in 2013 and 7109.5 t in 2018, respectively. Figure S3 shows emissions and changes for different sectors in China, Africa, south Asia, Europe, North America, and South America in 2013 and 2018, where the residential sector, industry, and agriculture are the main sources of BaP. China is one of the largest BaP-emitting countries in the world. Its emissions were 1952.2 t in 2013 and 1750.2 t in 2018, respectively, accounting for about 27.2 % and 24.6 % of the world, which is generally consistent with the results of Shen et al. (2013) (20.3 %). The residential sector is the largest emission source of BaP in China, accounting for 63.2 % of the total emissions, followed by the industrial sector (35.6 %) in 2013. Africa and south Asia had the second- and third-largest emissions, with residential combustion accounting for 87.3 % and 84.0 % of total emissions, respectively. This is related to the widespread use of biomass fuels for heating and cooking in developing countries (Han et al., 2020). Emissions from China, southern Africa, and south Asia accounted for 62.1 % and 61.4 % of the world. China, Australia, south Asia, Europe, North America, South Korea, Japan, and North Korea displayed a declining trend from 2013 to 2018. China experienced the largest decline (10.4 %), due to the active emission control measures taken under the “Air Pollution Prevention and Control Action Plan” implemented in 2013. The industrial sector contributed the most to the decline in BaP emissions followed by residential combustion, and which decreased by 18.9 % and 5.1 % from 2013 to 2018, respectively, which are mainly related to the strengthened industrial emission standards, upgraded industrial boilers, and development of clean fuels. The above results are generally consistent with the conclusions of Wang et al. (2021) in which PAHs emissions decreased by 11.36 % from 2013 to 2017, with the industrial and residential sectors decreasing by 17.32 % and 10.58 %, respectively. The emissions increased in Africa (10.3 %) and South America (5.9 %), which was mainly caused by the emission increase in residential (10.6 %) and agricultural (9.5 %) sectors, respectively.

To understand the change in BaP concentrations, we conducted five experiments: the first and second experiments simulated the BaP concentration using the emissions in 2013 and 2018 driven by the corresponding meteorological fields. The third experiment used the emission in 2018 but kept the meteorological conditions in 2013 to investigate the effects of meteorological condition changes on the concentration of BaP. Studies neglecting the heterogeneous loss of BaP and using two different heterogeneous schemes were also performed to explore the impacts of heterogeneous reactions on BaP concentrations in the fourth and fifth experiments.

3.2 Observational data

To evaluate the model performance, we collected the PAH observations from several available datasets and more than 50 published papers. The observational data are summarized as follows: (1) European Monitoring and Evaluation Program (EMEP; available from https://projects.nilu.no/ccc/reports.html, last access: 15 December 2023), which includes annual and monthly averages of BaP concentrations at 36 European sites in Spain, Finland, France, Germany, Norway, Poland, and other countries in Europe; (2) National Air Pollution Surveillance network (NAPS; https://data-donnees.az.ec.gc.ca/data/air/monitor/, last access: 30 January 2024), which includes daily averages (autumn and winter) of BaP concentrations at Canadian stations; (3) Integrated Atmospheric Deposition Network (IADN; https://iadnviz.iu.edu/datasets/index.html, last access: 20 January 2023), which includes monthly mean concentrations of BaP at six sites in the United States and Canada from 1990 to 2021; (4) Chinese Persistent Organic Pollutants (POPs) Soil and Air Monitoring Program Phase II (SAMP-II; Ma et al., 2018), which is carried out by the International Joint Research Center for Persistent Toxic Substances (IJRC-PTS), focusing on 11 urban centers in China (Beijing, Xi'an, Nanchang, Kunming, Lanzhou, Chengdu, Harbin, Dalian, Lhasa, Guangzhou, and Shihezi), one suburb and three background/rural areas, with observational data only covering the period from August 2008 to July 2010; and (5) observational data collected from published papers (these sources are listed in the Supplement).

PAHs measurements data are very sparse compared to conventional pollutants (e.g., PM_2.5). Since most of the data are not continuous in time, we selected data covering at least 10 d in years as close as possible to the simulation year (2013) and used the mean values for comparison. The comparison of the monthly variation was conducted only for sites in Europe where observations were continuous and available. The locations of the BaP observation sites are shown in Fig. S2. The site information is listed in Tables S5 and S6.

4 Results

4.1 Global distribution of BaP and health risks

To evaluate the performance of the IAP-AACM model, annual mean simulated concentrations in Asia, Europe, the United States, and Canada were compared with observations (Fig. 2). The results show that the model can reproduce nearly half of the observation samples within a factor of 2 and most observations within a factor of 5 at sites in Asia, Europe, the United States, and Canada. The number of sites where BaP was underestimated was greater than the number where it was overestimated due to the averaging effect of subgrid emissions. Considering that some of the comparisons are not in the same year, a certain discrepancy between the model and observation is expected. Further, a specific comparison was performed using the data from about 30 stations in Europe (Fig. 3a and b). High concentrations were mainly found in polluted areas of central Europe, which is consistent with the simulation of Gusev et al. (2019), such as Poland (P05) and the Czech Republic (C03), with observed values of 0.70 and 0.62 ng m⁻³, respectively, and simulated concentrations of 0.66 and 0.75 ng m⁻³ in 2013. The model successfully reproduced the observed concentrations and differences between sites.

Figure 4 Spatial distributions of seasonal mean concentrations in (a) Europe, (b) contiguous United States, and (c) east Asia in 2013.

Figure 4 shows the seasonal mean concentrations of BaP in Europe, contiguous United States, and east Asia in four seasons: March–April–May (MAM; representing spring), June–July–August (JJA; representing summer), September–October–November (SON; representing autumn), and December–January–February (DJF; representing winter). Generally, BaP had the highest concentration in winter and the lowest in summer. This is caused by the larger emission and poorer atmospheric diffusion conditions in winter than in summer. In the contiguous United States, the concentrations were lower than 1 ng m⁻³ in all four seasons, which is consistent with the simulation of Galarneau et al. (2014). In east China, large areas have a concentration of >1 and even >5 ng m⁻³ in BTH in winter. Europe shows a distribution of high values in central areas and low values in remote areas. In central Europe (such as Poland and the Czech Republic), large areas have concentrations between 1 and 5 ng m⁻³ in winter. High concentrations were also reported by Bieser et al. (2012). The observation clearly shows higher concentrations at sites in Poland and the Czech Republic than at other sites in Europe (Fig. 5). The model successfully reproduces the seasonal variation in BaP at sites in Europe. The simulation had a good agreement with observations at C03 and P05 with correlation coefficient (R²) of 0.91 and 0.91, and normalized mean bias (NMB) of −0.04 and 0.14, respectively. The R² was higher than 0.8 at B13 and S08, and the NMB was 1.01 and −0.55. In summary, IAP-AACM can reasonably simulate the spatial distribution and seasonal variation in BaP.

Figure 5 Comparison of the BaP month mean simulated concentrations (red) with observed values (black) at six stations in Europe in 2013.

Figure 6 Spatial distributions of annual mean BaP concentrations based on the EDGAR in (a) 2013 and (b) 2018. The (c, e) absolute and (d, f) relative concentration changes from 2013 to 2018 are shown considering (c, d) both emissions and meteorological conditions or (e, f) only emissions, respectively.

The spatial distribution of annual mean BaP concentrations based on the EDGAR inventory in 2013 and 2018 is shown in Fig. 6a and b. The spatial distribution of BaP concentrations in 2018 was similar to that in 2013. The spatial pattern was consistent with the emission distribution in the EDGAR inventory. High concentrations of BaP were found in northern and eastern China and central Europe, even exceeding the European Union target value for BaP (1 ng m⁻³), indicating an urgent need to control BaP and other PAHs. The absolute and relative concentration changes from 2013 to 2018 are shown in Fig. 6c and d. The most significant decreases were seen in Russia, the United States, and eastern and northern China. By contrast, the concentration in India, Europe, southeast Asia, and south Africa shows an increase, with the average annual concentration increasing by 19.4 %, 1.2 %, 11.2 %, and 18.3 %, respectively. When only considering the impact of emission change (Fig. 6e and f), the decrease in the eastern United States is larger, and the increase in central Europe is larger. In particular, there is an obvious decline (about 8.0 %) across China, which demonstrates the effect of the emission reduction measures. These results clearly show the large influence of meteorological changes. It is crucial to consider meteorological factors when evaluating emission changes and reduction measures through monitoring concentrations in the atmosphere.

Figure 7 The distribution of health risks grade in (a) 2013 and (b) 2018, the distribution of TILCR in (c) 2013 and (d) 2018, and the absolute from 2013 to 2018 when considering the change in (e) emissions only and (f) both emissions and meteorological conditions.

The global distribution of health risks grade and TILCR (the sum of ILCR values of the two exposure routes after averaging the parameters of the different groups) in 2013 and 2018 are shown in Fig. 7. It can be seen that most of the countries just face negligible cancer risk and that no regions face high potential cancer risk when evaluation is based on annual mean concentration. However, eastern China and central Europe are in potential cancer risk, with the highest TILCR values of 1.54 × 10⁻⁵ (1.47 × 10⁻⁵) and 4.62 × 10⁻⁶ (4.23 × 10⁻⁶) in 2013 (2018), respectively. It should be noted that in other developing countries such as Africa, India, and southeast Asia, TILCR values increased in 2018 compared to 2013, with the highest (mean) values increasing from 1.47 × 10⁻⁶ (4.15 × 10⁻⁸) to 2.26 × 10⁻⁶ (4.78 × 10⁻⁸), from 1.41 × 10⁻⁶ (2.03 × 10⁻⁷) to 1.73 × 10⁻⁶ (2.43 × 10⁻⁷), and from 8.02 × 10⁻⁷ (6.07 × 10⁻⁸) to 1.14 × 10⁻⁶ (6.75 × 10⁻⁸), respectively. These countries will be likely to face potential or even high potential cancer risks, especially in winter. In addition, Lou et al. (2023) have indicated that health risks in south Asia and Africa will increase in the future with the increase in residential fuel use along with the population growth in these regions. It is clear that clean development is a necessary consideration for developing countries to avoid the health risks posed by PAHs.

4.2 Distribution of PAHs and their change in China

Figure 8 shows the annual mean distribution of BaP in China in 2013. The concentrations ranged from 0.02 to 6.14 ng m⁻³. Overall, high concentrations simulated in the North China Plain, east China, and northeast China were significantly higher than in northwest and southwest China, which is consistent with previous studies (Ma et al., 2020; Yan et al., 2019). This can be largely attributed to the industrial and residential coal combustion in these regions. Among the different provinces in China, there are 14 provinces with concentrations higher than the ambient air quality standards of China (1 ng m⁻³, GB 3095–2012; http://www.mee.gov.cn/, last access: 6 April 2014). Shanghai had the highest concentration of 6.14 ng m⁻³ followed by Tianjin (4.56 ng m⁻³), Beijing (3.41 ng m⁻³), and Shandong (3.10 ng m⁻³). The concentrations in the northwest and southwest regions were lower, with Tibet having the lowest concentration of only 0.02 ng m⁻³. This is due to lower levels of industrial activities and population density in these regions compared to eastern regions. In addition, the high topography of northwest regions has good air circulation and is conducive to the diffusion and dilution of atmospheric pollutants. In 2013, Beijing had the highest BaP concentration in winter (14.03 ng m⁻³), possibly due to the high population density, high number of vehicles, and frequent industrial activities in Beijing. Moreover, Beijing lies on the North China Plain, where the meteorological conditions make it easier for air pollutants to stay and accumulate, resulting in a high concentration of BaP.

Figure 8 Spatial distributions of annual mean concentrations in China. Comparison of the BaP month mean simulated concentrations (orange) with observed values (blue) at eight cities in 2013.

To reveal the seasonal variation in BaP concentrations in key regions, we analyzed the concentration in eight major cities, i.e., Beijing, Tianjin, Shijiazhuang, Xinxiang, Wuhan, Chengdu, Guangzhou, and Harbin (seen in Fig. 8). It can be seen that the seasonal variations in BaP in these cities are similar, with the highest values in winter and the lowest in summer. The seasonal difference in northern cities was significantly greater than that in southern cities. In Beijing, Xinxiang, Tianjin, Harbin, and Shijiazhuang, the differences in concentration between winter and summer were as high as 15.06, 11.76, 11.14, 9.45, and 12.42 ng m⁻³, respectively. This is caused by the fact that coal-fired heating is very common in northern China, which can significantly increase PAH emissions in winter (Yan et al., 2019). In addition, the meteorological conditions also affect the seasonal variation in PAHs, as lower temperature, less rainfall, and weaker solar radiation during the winter support the formation of a stable inversion layer, greatly limiting the diffusion of BaP in the air (Lin et al., 2015; Quan et al., 2014).

By comparing the simulated concentrations with the observed concentrations, we find that the model can capture the BaP concentrations and the seasonal pattern in different cities. For example, the observed and simulated concentrations show good consistency in the spring, summer, and autumn in Chengdu and in the summer, autumn, and winter in Beijing, with a deviation of only 0.04 to 1.1 ng m⁻³. However, there were some deviations between the simulated and observed concentrations. The most obvious underestimation is seen in Shijiazhuang. This is probably due to the underestimation of emissions and the model resolution that may not fully resolve the pollution in cities with urban areas smaller than the model grid. The model performance could be improved by more precise emissions and increasing grid resolution. Nonetheless, the model can capture the magnitude and seasonal variation in BaP concentration well in China and in other countries around the world and can therefore be used to evaluate the health effects of BaP exposure.

Figure 9 The (a, c) absolute and (b, d) relative concentration changes from 2013 to 2018 in mean annual BaP concentrations in China are shown considering (a, b) both emissions and meteorological conditions or (c, d) only emissions, respectively.

The changes from 2013 to 2018 are shown in Fig. 9. The trend in and magnitude of changes differ greatly across different regions. The largest decrease (>20 %) was seen in northeastern and southeastern China, which is consistent with the other studies in Beijing (Lin et al., 2024), Shanghai (Yang et al., 2021), Shenyang (Zhang et al., 2023), and Tianjin (Zhang et al., 2022). The concentration also decreased in the North China Plain, benefitting from emission reduction policies, such as the development of cleaner energy sources and the control of industrial emissions. The concentration decrease was larger than the emission reduction over these areas. By contrast, as shown in Fig. 9a, the concentration in the Sichuan Basin showed an inverse trend although the emission decreased. This phenomenon demonstrates the large impact of unfavorable meteorological conditions with decreasing temperatures and planetary boundary layer height compared to 2013 (Ding et al., 2019). When only considering the emissions changes, the concentration shows a decrease over most regions consistent with the emission change. It should be noted that the decrease in BaP in the two experiments is significantly lower than that of PM_2.5. As shown in Fig. 9b (Fig. 9d), the BaP concentration decreased by 9.1 % and 6.7 % (8.5 % and 9.4 %) in the BTH and the YRD, respectively. Wang et al. (2019) showed that compared with 2013, the concentration of PM_2.5 in the BTH, the YRD, and the PRD in 2017 decreased by 39.6 %, 34.3 %, and 27.7 %, respectively. For cities in north and east China, the concentration still exceeds the national limit value (1 ng m⁻³) although the concentration of BaP decreased significantly in 2018. For example, the BaP concentrations in Shanghai, Beijing, and Tianjin considering changes in both emissions and meteorology were 5.32, 3.31, and 3.38 ng m⁻³, respectively, and those with emission changes alone were 5.58, 3.11, and 4.17 ng m⁻³, indicating that the concentrations are mainly affected by the emission sources. The results in the central and western cities differed greatly between the two experiments, especially in Chongqing, Sichuan, and Guizhou, indicating that changes were mainly related to meteorological conditions. Therefore, when formulating emission reduction policies, it is necessary to take into account the effects of changes in meteorological conditions as well as emission sources.

4.3 Health risks of PAHs in China

In this section, the health risks of BaP are assessed based on the simulation over the nested domain (domain 2) covering China. As shown in Fig. 10a and b, the health risks in western China are negligible, while there is a potential cancer risk in eastern China. Figure 10c and d show the TILCR values in 2013 and 2018, and Fig. 10e and f show the change from 2013 to 2018, when only emissions or both emissions and meteorological conditions are considered. It can be seen that the spatial distribution of TILCR (Fig. 10c) is consistent with the spatial distribution of BaP annual concentrations (Fig. 6), showing a higher risk in eastern regions than in the western regions (Han et al., 2020). The values of the TILCR in China ranged from 3.0 × 10⁻⁹ to 3.5 × 10⁻⁵, with an average value of 1.4 × 10⁻⁶. Compared with 2013, the average TILCR in 2018 decreased by 8.0 × 10⁻⁸, which is mainly directly related to the decrease in concentration. From the perspective of two exposure routes, ILCR_inh and ILCR_der values ranged from 2.8 × 10⁻¹⁰–3.3 × 10⁻⁶ and 2.7 × 10⁻⁹–3.2 × 10⁻⁵, with an average value of 1.3 × 10⁻² and 1.3 × 10⁻⁶, respectively. The values of ILCR_der were 1 order of magnitude higher than ILCR_inh.

Figure 10 The distribution of health risks grade in China in (a) 2013 and (b) 2018, the distribution of TILCR in (c) 2013 and (d) 2018, and the TILCR changes from 2013 to 2018 when considering the change in (c) emissions only and (f) both emissions and meteorological conditions.

Figure 11 The TILCR values for the three age groups (children, women, and men) in different provinces of China in 2013.

The TILCR values of the three groups in 2013 and 2018 calculated using Eqs. (14)–(18), which ranged from 1.55 × 10⁻⁹ to 3.78 × 10⁻⁵ (1.60 × 10⁻⁹ to 3.41 × 10⁻⁵), are shown in Figs. 11 and S4, respectively. The order of TILCR was women (1.46 × 10⁻⁶) > men (1.31 × 10⁻⁶) > children (7.03 × 10⁻⁷), which was similar to that of dermal contact exposure. Overall, 29.2 % of TILCR was higher than 1.0 × 10⁻⁶, and 1.2 % of TILCR were higher than 1.0 × 10⁻⁵ in 2013. There was a slight decrease in TILCR values in 2018 due to the lower concentrations of BaP, with 27.9 % and 0.7 % of TILCR being higher than 1.0 × 10⁻⁶ and 1.0 × 10⁻⁵, respectively. There is no high cancer risk in China, but there are potential cancer risks in some areas, which should be paid attention to.

Figure 12 The ILCR_inh values for the three age groups (children, women, and men) in different provinces of China in 2013.

Figure 13 The ILCR_der values for the three age groups (children, women, and men) in different provinces of China in 2013.

The ILCR values of the three groups in China through inhalation and dermal exposure routes are shown in Figs. 12 and 13, and the ILCR in 2018 are shown in Figs. S5 and S6. For the inhalation pathway, the average ILCR_inh was 1.22 × 10⁻⁷ (<1.0 × 10⁻⁶). The order of ILCR_inh was men (1.53 × 10⁻⁷) > women (1.43 × 10⁻⁷) > children (6.95 × 10⁻⁸), and the risk for men was about twice that of children but was lower in women than in men. This may be caused by the fact that the inhalation and metabolic rate of women are weaker than those of men (Bai et al., 2020). The highest average value was found in Shanghai, where the average ILCR_ing for the three groups was 1.72 × 10⁻⁶, 1.62 × 10⁻⁶, and 7.84 × 10⁻⁷, respectively. Han et al. (2020) found cases of excess cancer due to exposure to PAHs in large cities such as Shanghai. Only 1.6 % of the three groups had ILCR_ing higher than 1 × 10⁻⁶, indicating that the health risks from inhalation exposure were low. A similar conclusion was mentioned in an earlier review (Yan et al., 2019).

For dermal contact exposure, the average ILCR_der was 1.04 × 10⁻⁶ (>1.0 × 10⁻⁶). Compared to ILCR_inh, the health risk to adults was slightly higher than that to children, but women had higher risk values than men. This may be caused by the fact that the body weight of women is lower than that of men. The order of ILCR_der was women (1.32 × 10⁻⁶) > men (1.15 × 10⁻⁶) > children (6.33 × 10⁻⁶), which is similar to the results of previous studies (Bai et al., 2020). Among the three groups, 27.4 % of the ILCR_der values was higher than 1 × 10⁻⁶, and 0.7 % was higher than 1 × 10⁻⁵. This shows that there is a greater potential carcinogenic risk through dermal contact exposure.

5 Discussion

It should be noted that model results have some uncertainties even though our model simulated the main features of PAH concentrations reasonably well. Firstly, we simulated lower BaP concentrations when using the PKU inventory than when using the EDGAR inventory over most continental areas, except for Inner Mongolia, eastern Russia, and north China (Fig. S7). The difference can be as high as 0.5 ng m⁻³ over some areas in wintertime although the spatial and temporal distributions are consistent. The emission inventory remains constrained by more observations. Current observations are too sparse to conduct detailed evaluation in areas where long-term measurements are not available. Secondly, we tested the influence of heterogeneous reaction schemes on simulation. When heterogeneous reactions were not considered, the model significantly overestimated the concentration of BaP (Fig. S8), suggesting the importance of heterogeneous loss of BaP. Using the Langmuir–Hinshelwood mechanism, we simulated lower concentrations in most regions of the Northern Hemisphere, especially in winter (Fig. S9). However, the difference between the simulated results of the two mechanisms is significantly reduced in summer due to high temperatures and high humidity. This is consistent with the results by Mu et al. (2018); i.e., low temperature and low humidity can significantly increase the lifetime of BaP. Comparison of model results using different schemes and model intercomparison would further help identify the uncertainties and improve model performance.

6 Conclusion

In this study, the PAHs modules were coupled into the IAP-AACM model to investigate the global and regional distribution of PAHs. The model has the state-of-the-art heterogeneous mechanism and allows us to consistently examine the multi-scale distribution of PAHs. Comparison with observations shows that the model can reproduce the different concentrations of BaP at the stations in Asia, Europe, the United States, and Canada. The model can capture the seasonal variation in BaP, with lower concentrations in summer and higher concentrations in winter over the continents in the Northern Hemisphere. The global distributions of BaP in 2013 and 2018 were very similar, with high concentrations centered in eastern China and central Europe that even exceed EU limits (1 ng m⁻³). Compared with 2013, BaP concentration in 2018 showed a decrease in the United States, Poland, France, Czechia, and some regions in China. By contrast, the concentrations increased by >10 % in India and south Africa. Populations in these regions are facing increased health risks posed by PAHs.

In China, the decline in the BTH (8.5 %) and the YRD (9.4 %) benefitted from the Action Plan. However, the decline was significantly lower than that of conventional pollutants, such as PM_2.5. Changes in meteorological conditions had a significant influence on changes in BaP concentration, which increased in the Sichuan Basin and central China even though the emissions over these areas decreased due to the control measures. There was a slight decrease in total ILCR (TILCR) values in 2018 compared to 2013. For the different exposure routes, dermal contact was an order of magnitude higher than the inhalation route. The TILCR for adults was greater than that for children. In total, 29.2 % of TILCR were higher than 1.0 × 10⁻⁶, indicating that there are still potential cancer risks in China. More attention must be paid to the non-traditional pollutants, and strict but different control measures are necessary to reduce PAHs' concentration and health risks. Especially in the fall and winter seasons, the concentration of BaP and the associated health risk are significantly higher than in other seasons. Management efforts in key sectors (e.g., industrial and residential sources) should be further strengthened. In heavily polluted cities, using clean energy to replace coal combustion, adjusting the energy structure of factory production, and developing innovative technologies with lower and even no emissions would be helpful to reduce PAH pollution.

In summary, our study developed an effective tool for simulating the global and regional concentrations of BaP and other PAHs and quantified the health risks in China from 2013 to 2018. Model analysis indicated that emission inventories and heterogeneous reactions can significantly affect the simulated BaP concentrations. Accurate emissions and reasonable representation of heterogeneous reactions would greatly reduce the gap between model results and observations. However, current observations are insufficient to fully evaluate and constrain the model. Long-term observations are particularly needed in Asia, India, and Africa. These regions are still facing significant health risks. In addition, monitoring in the background and remote regions (such as the Arctic) is necessary to quantify the long-range transport of PAHs.