Variations of the atmospheric polycyclic aromatic hydrocarbon concentrations, sources, and health risk and the direct medical costs of lung cancer around the Bohai Sea against a background of pollution prevention and control in China

The Bohai Sea (BS) region, an area of China that is severely polluted with atmospheric polycyclic aromatic hydrocarbons (PAHs), has received wide attention in recent decades. To characterize the variations of the concentrations and sources of PAHs from June 2014 to May 2019, 15 PAH congeners (Σ₁₅PAHs) were measured in atmospheric samples (N = 228) collected at 12 sites around the BS, and the health risk and direct medical costs associated with lung cancer caused by exposure to PAHs were also estimated. The annual daily average concentration of Σ₁₅PAHs was 56.78 ± 4.75 ng m⁻³, which was dominated by low molecular weight PAHs (LMW-PAHs) containing three rings (58.7 ± 7.8 %). During the 5-year sampling period, the atmospheric Σ₁₅PAHs concentration reduced by 17.5 % for the whole BS region, with an especially large drop of 51.7 % observed in the tightly controlled area of Tianjin (TJ), which was mainly due to a decrease in the concentration of high molecular weight PAHs (HMW-PAHs) containing five or six rings. Generally, the concentration of Σ₁₅PAHs was highest in winter and lowest in summer, which can mainly be attributed to the change in LMW-PAH concentration. Based on the positive matrix factorization (PMF) model, PAHs in the BS region were mainly ascribed to coal combustion and biomass burning. The PAH contributions of coal combustion and motor vehicles showed different trends for the BS region (PAHs from coal combustion rose by 7.2 % and PAHs from motor vehicles fell by 22.4 %) and for TJ (PAHs from coal combustion fell by 12.6 % and PAHs from motor vehicles rose by 6.9 %). The incidence of lung cancer (ILCR) caused by exposure to atmospheric PAHs decreased by 74.1 % and 91.6 % from 2014 to 2018 in the BS region and in TJ, respectively. This was mainly due to the decrease in the concentration of highly toxic HMW-PAHs and was reflected in USD 10.7 million of savings in direct medical costs of lung cancer caused by exposure to PAHs, which represent a decrease of 46.1 % compared to the corresponding costs before air pollution prevention and control was implemented around the BS. There was an even higher reduction in medical costs of 54.5 % in TJ. Hence, this study proves that implementing pollution prevention and control not only effectively reduced the concentration of pollutants and the risks caused by them, but it also significantly reduced the medical costs of diseases caused by the corresponding exposure.

1 Introduction

Polycyclic aromatic hydrocarbons (PAHs), a class of classical organic compounds with at least two benzene rings, have long received attention because of their cytotoxic, teratogenic, mutagenic, or carcinogenic properties (Colvin et al., 2020; Marvin et al., 2020). The United States Environmental Protection Agency (USEPA) identified 16 PAH congeners as priority pollutants (Lv et al., 2020). Previous studies have shown that PAHs in the atmosphere of heavily polluted areas such as factories and urban environments pose a threat to human health, especially that of the respiratory system (Agudelo-Castañeda et al., 2017; Ramírez et al., 2011). Because of their relatively high concentrations, strong toxic potency, and capacity to be transported over long distances, PAH congeners in the atmosphere are considered as a major factor in the lung cancer risk to the public (Ma et al., 2010; Gong et al., 2011; Ma et al., 2013; Hong et al., 2016). In terms of its incidence and the mortality resulting from it, lung cancer is ranked higher than any other cancer around the world, and so the lung cancer risk owing to exposure to PAHs is of particular concern and has been widely assessed (Křůmal and Mikuška, 2020; Liao et al., 2011; Taghvaee et al., 2018; Zhang et al., 2023).

PAHs are emitted primarily via the incomplete combustion and pyrolysis of carbon-contained materials, such as fossil fuels and biomass (Biache et al., 2014). China was found to be the largest emitter of PAHs all over the world for the last two decades because of the rapid development of its economy and its increasing consumption of carbon-containing materials (Y. X. Zhang et al., 2007), particularly in the Bohai Economic Zone, as the third-fastest-developing economic region (Sun et al., 2022). PAH pollution in the atmosphere above the Bohai Sea (BS) is in a severe situation (Wang et al., 2018). The Bohai Economic Zone includes the Beijing–Tianjin–Hebei (BTH) region, the Liaodong Peninsula, and the Shandong Peninsula. The BTH region is the center of economic development in the Bohai Rim economic area (Liang et al., 2018; Zhang et al., 2016). Hence, the Beijing–Tianjin–Hebei (BTH) region is one of the regions with the highest PAH emission intensities and the heaviest atmospheric PAH concentrations in China (Zhang et al., 2007; Y. X. Zhang et al., 2007). Given the presence of such serious pollution, the health risk due to exposure to PAHs is the cause of great concern. The population attributable fraction (PAF) for lung cancer caused by the inhalation of PAHs in the atmosphere in the BTH area was more than twice higher than the mean value for the whole of China in 2009 (Zhang et al., 2009). The incremental lifetime cancer risk (ILCR) due to PAH exposure in Tianjin was in the range of 1 × 10⁻⁵ to 1 × 10⁻³ in 2008, which was much higher than the mean level of 4.56 × 10⁻⁶ for China (Lian et al., 2021). The annual lung cancer morbidity in Tianjin (within the BTH region), 6.99 × 10⁻⁶, was the highest for any city among 35 cancer registries in China (S. W. Zhang et al., 2007). Meanwhile, given the frequent occurrence of haze in the BTH region, more attention has been paid to the concentration levels and the health risk of fine particulate matter with an aerodynamic equivalent diameter of ≤ 2.5 µm (PM_2.5) since 2013 (Chen et al., 2020).

PM_2.5 pollution in China has clearly improved since the Air Pollution Prevention and Control Action Plan (2013–2017) and the Three-Year Action Plan for Winning the Blue-Sky Defense Battle (2018–2020) were proposed by the Chinese government in 2013 and 2018, respectively (Zhao et al., 2023). As it is one of the severely polluted areas in China, the improvement was more significant in the BTH region, which implemented the strictest pollution control policy (Li et al., 2020). As reported, the concentration of PM_2.5 in the BTH region dropped by 52 % from 106 µg m⁻³ in 2013 to 51 µg m⁻³ in 2020 (Bulletin of the State of China's Ecological Environment, 2021). Reducing emissions from coal combustion and motor vehicles is a major aim of policies aimed at the prevention and control of pollution (Guo et al., 2018; Li et al., 2019). These two sources have been recognized as primary contributors to the PAHs in the atmosphere as well (Lin et al., 2015; Han et al., 2018). As a result, imposing controls on these two sources reduced not only PM_2.5 emissions but also PAH emissions (Zhi et al., 2017). As part of the control processes, the variations in the concentrations and the health risk of PM_2.5 in the BTH region have been identified in detail (D. Fang et al., 2016; Yan et al., 2019), but our understanding of the PAHs in this region needs to be updated urgently. In particular, the statistical data on the lung cancer risk due to exposure to PAHs were established 10 years ago (Zhang et al., 2009; Lian et al., 2021).

To track changes in the concentrations and sources of atmospheric PAHs and to estimate the health risk and the direct medical costs associated with lung cancer caused by exposure to PAHs during the air pollution control actions, a field monitoring campaign was conducted at 12 sites around the BS for 5 years from June 2014 to May 2019. The air pollution control measures implemented were different for the BTH region, the Liaodong Peninsula, and the Shandong Peninsula (Huang et al., 2017). Thus, this study provided us with an opportunity to understand the variations in the environmental concentrations, source contributions, and health risk of PAHs. The main aims of this study were (1) to characterize the spatial and temporal changes in the concentrations and components of PAHs in the atmosphere around the BS, (2) to evaluate differences in PAH source contributions, and (3) to assess the changes in direct medical costs of treating lung cancer caused by inhalation exposure to PAHs against a background of atmospheric pollution prevention and control over 5 years.

2 Materials and methods

2.1 Sampling sites and sample collection

The sampling sites used for this study were reported in previous literature (Sun et al., 2021) and are briefly introduced here. Information on the sites is shown in Table S1 in the Supplement. The 12 air sampling sites were located at Beihuangcheng (BH), Dalian (DL), Donggang (DG), Dongying (DY), Gaizhou (GZ), Longkou (LK), Laoting (LT), Rongcheng (RC), Tianjin (TJ), Xingcheng (XC), Yantai (YT), and Zhuanghe (ZH). A passive air sampler with polyurethane foam (PUF, 14.00 cm diameter × 1.35 cm thickness) was used to collect atmospheric samples at each sampling site (Eng et al., 2014). The PUF disks were deployed around 1.5–2.0 m above the ground, and the sampling duration was about 3 months for one batch. 228 samples were collected from June 2014 to May 2019. The sampling rate of atmospheric PAHs was 3.5 m³ d⁻¹ (Jaward et al., 2005; Moeckel et al., 2009). Prior to sampling, the PUF disks were pre-cleaned with methanol, acetone, and hexane. The extracted PUF disks were placed in airtight containers and stored at −18 ^∘C before the sampling campaign. After sampling, the samples were prepared and then stored in a freezer at −18 ^∘C in the lab for further analyses.

2.2 Sample pretreatment and instrumental analysis

The five PAH surrogates (naphthalene-D₈, acenaphthene-D₁₀, phenanthrene-D₁₀, chrysene-D₁₂, and perylene-D₁₂) and the activated copper fragments were added in advance (Qu et al., 2022). The samples were extracted for 24 h with an eluent of acetone and hexane (200 mL, v:v = 1:1) through a Soxhlet apparatus. The extracted solution was concentrated to 1 mL with a rotary evaporator (Zhengzhou Greatwall Ltd., China). Then, a silica-alumina column was used to obtain the aromatic components. The targets were obtained with 40 mL of a mixed solution of dichloromethane and hexane (v:v = 1:1). Finally, the eluent was concentrated and reduced to 500 µL by a gentle nitrogen stream. 400 ng of hexamethylbenzene (Supelco, USA) were added as an internal standard substance to each sample solution before the instrumental analysis.

The targets were detected using a gas chromatograph equipped with a mass spectrometer (GC-MS, Agilent 5975C-7890A, USA), and the chromatographic column was a DB-5MS (Agilent Technologies, 30 m × 0.25 mm × 0.25 µm). Each 1 µL extract was injected in splitless mode. High-purity helium (purity ≥ 99.99 %) with a flow rate of 1.3 mL min⁻¹ was used as the carrier gas. The oven temperature was set and held at 80 ^∘C for 3 min; it was then raised at 10 ${}^{\circ }\mathrm{C}{\min }^{-\mathrm{1}}$ to 310 ^∘C, where it was held for 10 min. The temperatures of the inlet and the ion source were 290 and 230 ^∘C, respectively. Details regarding the targeted compounds are shown in Table S2. Seven gradients of mixed solutions were established for the quantitative calculation of PAHs. More details were reported in a previous study (Wang et al., 2018).

2.3 Quality assurance and quality control

The mean recovery values of naphthalene-D₈, acenaphthene-D₁₀, phenanthrene-D₁₀, chrysene-D₁₂, and perylene-D₁₂ were 77.3 %, 85.9 %, 87.5 %, 88.3 %, and 92.8 %, respectively, with the total range of 66.5 %–123.1 %. All the relative deviations were within 20 % except that for naphthalene-D₈. Naphthalene (Nap) was excluded from this investigation because of its low recovery, and another 15 PAHs were determined in this study. For each batch of 12 PUF samples, a field blank and a procedural blank were also analyzed following the same treatment process. In this study, the method detection limits (MDLs, defined as the mean blank value plus 3 times the standard deviation) for the 15 PAH congeners ranged from 0.02 to 0.13 ng sample⁻¹, as shown in Table S2. The final concentrations were not surrogate corrected. The glassware was all cleaned and burned for 8 h in a muffle oven at 450 ^∘C before the experiment. The solvents were chromatography pure or had been redistilled and purified before using.

2.4 PAH source apportionment

The model of positive matrix factorization (PMF) released by the USEPA (PMF 5.0) was used to apportion the emission sources of PAHs in this study. The basic calculation formula for the PMF method is given by Eq. (1):

$$\begin{array}{}\text{(1)}& x_{ij}=\sum _{k=\mathrm{1}}^p{g}_{ik}{f}_{kj}+e_{ij},\end{array}$$

where p represents the number of sources identified by the PMF model, x_ij represents the original concentration data for the ith chemical species and the jth sample, g_ik represents the source profile of the kth source and the jth chemical species, f_kj represents the contribution ratio of the kth source to the jth sample, and e_ij represents the simulated residual error of the ith chemical species and the jth sample. Source contributions and profiles are solved in the PMF model by minimizing the objective function Q as in Eq. (2):

$$\begin{array}{}\text{(2)}& Q_{min}=\sum _{i=\mathrm{1}}^n\sum _{j=\mathrm{1}}^m{\left({\displaystyle \frac{x_{ij}-{\sum }_{k=\mathrm{1}}^p{g}_{ik}{f}_{kj}}{u_{ij}}}\right)}^{\mathrm{2}},\end{array}$$

where x_ij, g_ik, and f_kj are the same as in Eq. (1). u_ij is the uncertainty in x_ij, and the calculation method for the uncertainty is shown in Text S2. More details have been previously documented (Sofowote et al., 2011; Paatero et al., 2014).

Before the source apportionment, principal component analysis (PCA) was applied to pre-estimate the minimum number of emission sources in this study because PCA was able to explain the overall variability with fewer variables with a minimum loss of information (Liu et al., 2021). SPSS Statistics 25.0 was used to perform the PCA analysis in this study.

2.5 Health risk assessment

The total toxicity equivalent (TEQ, ng m⁻³) of the 15 PAHs with benzo[a]pyrene (BaP) used as the reference is calculated as in Eq. (3):

$$\begin{array}{}\text{(3)}& \mathrm{TEQ}=\sum _{i=\mathrm{1}}^n(C_{\mathrm{i}}\times {\mathrm{TEF}}_i),\end{array}$$

where C_i is the concentration of the ith PAH compound (ng m⁻³) and TEF_i is the cancer potency of the ith PAH compound (dimensionless), as shown in Table S2.

In this study, ILCR refers to the cancer risk in a population due to exposure to a specific carcinogen (Zhuo et al., 2017). Its calculation formula is given by Eq. (4):

$$\begin{array}{}\text{(4)}& \mathrm{ILCR}={\mathrm{UR}}_{\mathrm{BaP}}\times \mathrm{TEQ}.\end{array}$$

In the above, UR_BaP represents the cancer risk when the concentration of BaP is 1 ng m⁻³. According to the regulations of the World Health Organization (WHO), UR_BaP can be 8.7 × 10⁻⁵ per ng m⁻³. That is, for a life span of 70 years, a lifetime of exposure to a BaP concentration of 1 ng m⁻³ results in a risk of cancer by inhalation of 8.7 × 10⁻⁵ (Luo et al., 2021).

2.6 Medical costs assessment

In this study, the medical costs were assessed by comparing the total direct medical costs for treating lung cancer caused by respiratory exposure to PAHs in the atmosphere under the assumption that there is no air pollution control with the corresponding costs when air pollution control has been implemented. The total direct medical costs for treating lung cancer (C_t) are calculated as in Eq. (5):

$$\begin{array}{}\text{(5)}& C_{\mathrm{t}}=C_{\text{pc}}\times P\times I_{\text{add}},\end{array}$$

where C_t is the total direct medical costs of lung cancer induced by PAH exposure and C_pc is the per capita direct medical costs of lung cancer, for which USD 8700 in China in 2014 was used in this study (Shi et al., 2017). P is the annual population and I_add is the additional incidence of lung cancer due to PAH inhalation exposure, which is calculated as in Eq. (6):

$$\begin{array}{}\text{(6)}& I_{\text{add}}=I\times \mathrm{PAF},\end{array}$$

where I is the incidence of lung cancer. The estimated value of I at Tianjin in 2012 was 87.37 × 10⁻⁵ (Cao et al., 2016), which was used in this study. PAF is the population attributable fraction, defined as the decrease in the incidence or mortality of a disease when a certain risk factor is completely removed or reduced to another lower reference level (Menzler et al., 2008). The PAF can be calculated as in Eq. (7):

$$\begin{array}{}\text{(7)}& \begin{array}{rl}\mathrm{PAF}& ={\displaystyle \frac{\mathrm{rr}\left(\mathrm{TEQ}\right)-\mathrm{1}}{\mathrm{rr}\left(\mathrm{TEQ}\right)}}\text{ and }\mathrm{rr}\left(\mathrm{TEQ}\right)\\ & =[{\mathrm{URR}}_{\mathrm{cum}},\exp =\mathrm{100}{]}^{(\mathrm{TEQ}\times \mathrm{70}/\mathrm{100})},\end{array}\end{array}$$

where rr is the relative risk; that is, the risk of exposure to a specific concentration relative to no exposure. URR is the unit relative risk, for which a reference value of 4.49 per 100 µg m⁻³ years of BaP exposure was adopted in this study (Zhang et al., 2009). This reference value is based on an epidemiological study of lung cancer conducted in Xuanwei, China (Menzler et al., 2008) (Gibbs, 1997). This study assumed that the mean life expectancy in China was 70 years and that the lifetime exposure was equivalent to 70 years.

Figure 1 Atmospheric concentrations of polycyclic aromatic hydrocarbons (PAHs) around the BS from June 2014 to May 2019. Acenaphthylene (Acy), Acenapthene (Ace), Fluorene (Flu), Phenanthrene (Phe), Anthracene (Ant), Fluoranthene (Fla), Pyrene (Pyr), Benzo[a]anthracene (BaA), Chrysene (Chr), Benzo[b]fluoranthene (BbF), Benzo[k]fluoranthene (BkF), Benzo[a]pyrene (BaP), Dibenzo[a,h]anthracene (DahA), Indeno[1,2,3-cd]pyrene (InD), Benzo[ghi] perylene (BghiP), low molecular weight PAHs (LMW-PAHs), middle molecular weight PAHs (MMW-PAHs), high molecular weight PAHs (HMW-PAHs), and total concentrations of 15 PAH congeners.

3 Results and discussions

3.1 Concentration and composition of PAHs

3.1.1 General information on PAHs

Figure 1 summarizes the annual daily average concentrations of the 15 PAHs in the atmosphere at the 12 sampling sites around the BS from June 2014 to May 2019. The annual daily average concentration of Σ₁₅PAHs around the BS was 56.78 ± 4.75 ng m⁻³, with a range of 51.39–63.55 ng m⁻³. The highest concentration was presented by the low molecular weight PAHs (LMW-PAHs, three rings), followed by the middle molecular weight PAHs (MMW-PAHs, four rings) and high molecular weight PAHs (HMW-PAHs, five or six rings), which accounted for 58.7 %, 34.8 %, and 6.7 % of the total concentration, respectively. In this study, the atmospheric PAH concentration was dominated by the LMW-PAHs, of which Phe, Fla, and Flu were the main compounds, accounting for 37.7 %, 19.8 %, and 12.6 % of the total. The atmospheric PAH concentrations around the BS were higher than those around the Yangtze River Delta and the Pearl River Delta, such as those in Ningbo (45 ng m⁻³) (Tong et al., 2019) and Guangzhou (9.72 ng m⁻³) (Yu et al., 2016). The atmospheric concentrations of PAHs around the BS were also much higher than those in the atmosphere above the Great Lakes (1.3 ng m⁻³) (Li et al., 2021) and southern European cities (3.1 ng m⁻³) (Alves et al., 2017). Overall, it was found that the pollution from atmospheric PAHs around the BS was still worrying.

3.1.2 Temporal variations of PAHs

To obtain a better understanding of the characteristics of the variation in PAHs in the atmosphere, the period from the summer of the previous year to the spring of the next year was taken as a statistical cycle. The concentrations of Σ₁₅PAHs around the BS during five annual cycles were 63.55 ± 58.43 ng m⁻³ (2014–2015), 55.50 ± 37.94 ng m⁻³ (2015–2016), 60.90 ± 31.13 ng m⁻³, (2016–2017), 51.39 ± 29.41 ng m⁻³ (2017–2018), and 52.50 ± 40.08 ng m⁻³ (2018–2019) (Table S3). Overall, the concentrations of Σ₁₅PAHs from June 2014 to May 2019 showed a slow downward trend, with a decrease of 17.5 %. The decrease in atmospheric PAH concentrations was mainly due to the decline in the HMW-PAH concentrations. The HMW-PAH composition ratio decreased from 11.3 % (2014–2015) to 3.4 % (2018–2019), while that of the MMW-PAHs rose from 35.5 % (2014–2015) to 41.2 % (2018–2019). The LMW-PAH composition ratio was stable, varying from 53.4 % (2014–2015) to 55.4 % (2018–2019). One factor that affects the concentrations of PAHs in the atmosphere after they have been discharged from the emission source is the meteorological conditions (Fan et al., 2021); the other is the amount of direct emission from the emission source (Ma et al., 2018). PAHs with different ring numbers arise from different sources (Li et al., 2021). LMW-PAHs are mainly produced in the combustion of non-petroleum sources, while HMW-PAHs are mainly from the products of processes based on the high-temperature combustion of fossil fuels, including some activities involving pyrolysis, such as vehicle emissions, industrial production, and other high-temperature source emissions (Zhang et al., 2018; Xing et al., 2020). The significant decrease in HMW-PAH concentrations in the BS region during the 5-year observation period might be related to the decrease in high-temperature emission sources. Due to the high toxicity of HMW-PAHs (Biache et al., 2014; Ma et al., 2020), the decrease in the concentration of HMW-PAHs might indicate a decrease in the environmental toxicity of PAHs.

The seasonal distribution of PAH concentrations in the atmosphere of the BS region showed that concentrations were high in the cold season and low in the warm season. The concentrations of Σ₁₅PAHs in the four seasons were as follows: winter (104.32 ± 9.50 ng m⁻³) > autumn (53.94 ± 9.10 ng m⁻³) > spring (43.89 ± 19.54 ng m⁻³) > summer (26.28 ± 13.42 ng m⁻³) (Table S5). The concentration of PAHs in winter was about 4 times higher than that in summer, and the maximum and minimum of the annual daily average concentration at the 12 sampling points mostly occurred in winter and summer. In addition, there were significant differences between the total PAH concentration and the concentrations of PAHs with different numbers of rings (p < 0.05; the difference level is shown in Table S6). The seasonal characteristics of PAH concentrations seen in this study were consistent with results reported for North China (Ma et al., 2018; Zhang et al., 2019). Interestingly, it was found that the differences between the PAH concentrations in the four seasons were mainly due to the LMW-PAHs. This indicated that there were nonnegligible pollution sources for LMW-PAHs, especially in winter in the BS region. Thus, identifying the source of LMW-PAHs was crucial for improving the environmental quality of the BS. Studies have shown that coal-burning emissions and biomass burning are the main sources of atmospheric PAHs in this region (Liu et al., 2019). In terms of the spatial distribution of the per capita fuel consumption, the consumption in North and West China is significantly higher than that in Southeast China, principally because of the difference in heating fuel consumption in winter. Therefore, there were significant seasonal variations of per capita fuel consumption, with peak consumption in the winter months being about twice as high as that in the summer months (Zhu et al., 2013). In addition, due to the migration characteristics of atmospheric PAHs, meteorological conditions such as the temperature and wind direction would also affect the observed concentrations in different seasons (Tan et al., 2006). A low temperature and inversion layer in winter are not conducive to atmospheric diffusion, resulting in a relatively high concentration of PAHs in the atmosphere near the surface (Wang et al., 2018).

Figure 2 The mean concentration distribution of Σ₁₅PAHs at 12 sites around the BS from June 2014 to May 2019.

3.1.3 Spatial characteristics of PAHs

Figure 2 and Table S7 display the distribution of the 5-year mean concentration of Σ₁₅PAHs from June 2014 to May 2019 at the 12 sampling sites around the BS. The concentrations of atmospheric Σ₁₅PAHs ranged from 25.92 ± 6.41 ng m⁻³ (RC) to 103.71 ± 39.11 ng m⁻³ (XC). The concentrations of PAHs on the north coast of the BS were twice higher than those on the south coast. PAHs are a class of pollutants that can undergo long-range transport in the atmosphere (Wang et al., 2018), and their spread is greatly affected by local meteorological conditions (Ding and Chan, 2005). The climate in North China and the adjacent oceanic area is greatly affected by the East Asian monsoon, and the characteristic weather phenomenon in the winter monsoon is the strong north and northwest winds (Tian et al., 2009). Due to the additional emissions from fuel consumption for domestic heating in the source areas, the atmospheric PAH concentrations significantly increased (Feng et al., 2007; Gao et al., 2016). Combined with the backward trajectory shown in Fig. S4, this suggests that the elevated PAH concentrations in winter in the north of the BS region were mainly attributable to the outflow of PAHs from the north and northwest source regions, carried by the winter monsoon winds. According to the distribution of atmospheric PAHs in some representative parts of northern China, it was found that the Beijing–Tianjin–Hebei region was greatly affected by nearby sources, while Shandong Province and other places were mainly affected by regional emissions (Zhang et al., 2016). However, the composition of PAHs in the north showed consistency with those in the south (there were no significant differences; Table S8). Overall, the compositions of PAHs at the 12 stations showed that the highest contribution came from LMW-PAHs (north: 60.0 %, south: 57.4 %), followed by MMW-PAHs (north: 32.7 %, south: 32.4 %), while HMW-PAHs provided the lowest (north: 7.3 %, south: 10.8 %). The above results indicate that the emission sources of PAHs in the atmosphere were the same all around the BS region.

However, for TJ, the study found that there was a more significant change in the concentration of atmospheric PAHs over time, as they decreased from 68.61 ng m⁻³ (2014–2015) to 33.14 ng m⁻³ (2018–2019). The reason for this was mainly that TJ is located in the Beijing–Tianjin–Hebei region, which has the strictest air pollution prevention and control measures as it is a key area in China's 12th Five Year Plan. To explore the potential differences between the source emissions at the 12 sampling points, Pearson correlation analysis was used to analyze the seasonal distribution of PAH concentrations, as shown in Table S9. Among the five stations LK, DY, TJ, LT, and XC in the western BS region centered on TJ, the correlation coefficients from comparisons of atmospheric PAH concentrations among the four stations except for TJ (0.72–0.89) were greater than the corresponding correlation coefficients between each site and TJ (0.50–0.68). The covariability in PAH concentrations between TJ and the other four stations was weaker. This indicated that there were certain differences in PAH emission sources between TJ and adjacent areas.

3.2 Source apportionment of PAHs

To further probe the causes of the variations of the concentration and composition of PAHs, the source apportionment of PAHs in the atmosphere around the BS region and TJ from 2014–2015 to 2018–2019 was investigated via PCA and PMF. The PCA analysis results showed that when four factors (eigenvalues > 1) were extracted from the data set, the total cumulative load accounted for more than 85 % of the variance (Table S10). This indicated that the atmospheric PAHs were best explained by at least four types of emission sources. For the PMF model, the key process was to determine the correct number of factors, and this study was based on the results of the PCA. Based on a random seed, 4–7 factors were used with the PMF model for source analytical simulation. The source analytical simulation using five factors yielded the most stable results and the most easily interpreted factors. The solution produced Q values (both robust and true) that were close to the theoretical Q values, which indicated that the PAH data set in the modeling input provided appropriate uncertainty (Sun et al., 2021). The data set used for PMF analysis included the concentrations of 228 samples of 15 PAHs and their uncertainties. The diagnostic regression R² value for the overall concentrations of the 15 PAH components was 0.986. The concentrations of the 15 PAHs predicted via the PMF model were almost the same as the actual concentrations of the 15 PAHs around the BS (Figs. S5 and S6 and Sect. S2). This meant that the model results were good and could be used as the basis for judging the sources of the target species, so these five factors explain the sources of the PAHs well. The sources identified by PCA and PMF were coal combustion, biomass burning, industrial processes, gasoline emissions, and diesel emissions. Detailed information on source identification is shown in Sect. S3.

Fossil fuel combustion emissions were the reason for the significant increase in atmospheric pollutants, and they were also responsible for the elevated PM_2.5 levels around the BS region (Yang et al., 2017). To explore the relationship between Σ₁₅PAHs and PM_2.5 concentrations, PM_2.5 data available online for eight cities that corresponded to sampling sites (DG, DL, DY, GZ, LT, TJ, XC, and YT) (Air quality historical data query, 2014–2019) were collected, and the concentrations were averaged according to the sampling periods in the study (Table S12). The Pearson correlation coefficients of the concentrations of atmospheric PAHs and PM_2.5 ranged from 0.485 to 0.868, and the significance levels were greater than 95 %, as listed in Table S13. During the 5-year observation period, the PM_2.5 concentration in the BS region decreased by 29.6 % (from 57 to 40 µg m⁻³), and that at TJ showed an even greater decrease of 33.8 % (from 78 to 51 µg m⁻³). From 2013, PM_2.5 has been strictly controlled year by year by the government, and the significant correlation indicates that this control affected the PAH concentrations in this region. To explore the factors that potentially influence the difference in atmospheric PAH composition between the BS area and TJ, the average annual contributions of various PAH emission sources from 2014–2015 to 2018–2019 were compared, as shown in Fig. 3. During the sampling period of the BS region, the main source of atmospheric PAH emissions in 2014–2015 was coal combustion (44.7 %) followed by biomass burning (24.4 %), but these switched to coal combustion (51.8 %) and industrial processes (24.4 %) in 2018–2019. For TJ, coal combustion was also the main source of atmospheric PAHs emissions (52.5 %) in 2014–2015, followed by biomass burning (20.1 %), which switched to coal combustion (40.0 %), industrial processes (17.8 %), and gasoline emissions (18.3 %) in 2018–2019. The contribution of coal combustion as a source of atmospheric PAHs around the BS increased by 7.2 %, while the corresponding contribution in TJ fell by 12.6 %. The absolute contribution (the total concentration of PAHs multiplied by the percentage contribution of the source) decreased, which indicated that the reduction in the coal contribution as a source significantly improved the atmospheric PAH pollution.

Figure 3 Concentrations of Σ₁₅PAHs and contributions of Σ₁₅PAHs sources around the BS (a) and TJ (b) from 2014–2015 to 2018–2019.

The main source of atmospheric PAHs around the BS was coal combustion (Liu et al., 2019; Qu et al., 2022). However, TJ, one of the key areas for air pollution control in China, took stricter measures to control emissions from coal combustion (Wu et al., 2015). For instance, the city took the lead in switching domestic fuel from coal to natural gas and electricity in 2017 to reduce emissions of air pollutants (Zhang et al., 2021). These targeted measures led to more forcefully controlled coal-combustion emissions of PAHs in TJ than in the other places around the BS region (Guo et al., 2018). Vehicle emissions (from gasoline and diesel exhausts) of atmospheric PAHs experienced a sharp drop of 22.4 % in the BS area, while those in TJ rose by 6.9 %. The same trend for vehicle emissions was found in studies of Beijing and Tianjin (Zhang et al., 2016; Chao et al., 2019). The decrease was mainly due to the elimination and scrapping of substandard vehicles by the Chinese government in 2015. Based on the China Vehicle Environmental Management Annual Report, car ownership around the BS increased by about 17.5 million but emissions of hydrocarbons, including PAHs, reduced by 95 000 t from 2014 to 2018 (Fig. S8). The source apportionment showed that the contribution of vehicle emissions to PAHs underwent a sharp decline from the spring of 2016 (Fig. S9), with a decrease of 38 % (19 % for gasoline and 19 % for diesel) around the BS (Huang et al., 2017). Although the contribution of vehicle emissions increased for TJ, the concentration of PAHs decreased. This indicates that these measures also controlled vehicle emissions and kept emissions of PAHs at a low level. Therefore, targeted control measures can effectively control PM_2.5 and PAH pollution in the atmosphere in the BS region and TJ. Moreover, PAHs are a kind of organic compound produced along with black carbon (BC), and, to some extent, the molecular characteristics of PAHs have also provided the basic data needed to analyze the sources of BC in the atmosphere of the BS region (Y. Fang et al., 2016). At the same time, the results of the PAH source analysis in this study revealed that the composition and source of atmospheric BC in the BS region also changed from 2014 to 2019. This problem needs our attention and confirmation.

Figure 4 shows the seasonal distribution of five sources of atmospheric PAHs in the BS region. Generally, the seasonal distribution of the five sources of atmospheric PAHs in TJ was consistent with that in the BS region, so they are not separately discussed here. The relevant information for TJ is shown in Fig. S10. Coal combustion was the main emission source across the four seasons, followed by biomass burning, while the contributions of the other sources (industrial processes, gasoline emission, and diesel emission) were similar. The contribution of coal combustion to atmospheric PAHs in the BS region was the highest in winter, followed by spring, while its lowest contribution occurred in summer. This was consistent with the seasonal distribution of the concentration of PAHs in the atmosphere in the BS region. Based on the seasonal distribution of concentration, we can conclude that the increased concentrations of atmospheric PAHs in winter were mainly caused by coal combustion. This might be due to people in northern China relying on coal combustion for heating in cold winters. The contribution from biomass combustion was higher in summer and autumn, which was related to straw burning after the harvest. Given all this, the seasonal distribution of PAH sources indicated that the pollution from atmospheric PAHs was mainly influenced by human activities.

Figure 4 The seasonal and average contributions of five sources of Σ₁₅PAHs, as derived from PMF; (a) the 5-year-average contributions of the five sources.

Figure 5 The annual mean TEQ values of the 15 PAHs, BaP, Dibenzo[a,h]Anthracene (DahA) and in the atmosphere at the BS (a) and in TJ (b) from June 2014 to May 2019.

3.3 Health risk due to exposure to PAHs

On the basis of Eq. (3), the annual mean TEQ value around the BS region was 1.37 ± 1.05 ng m⁻³ from June 2014 to May 2019, which is below the national standard (10 ng m⁻³) but slightly higher than the WHO standard (1 ng m⁻³). The HMW-PAHs contributed the most (76.4 %) to the total TEQ. However, the concentration of HMW-PAHs in the atmosphere accounted for 6.5 % of the total PAH concentration. The two major TEQ contributors were BaP (38.2 % ± 8.0 %) and DahA (16.6 % ± 9.0 %). For TJ, the annual mean TEQ value was 1.69 ± 1.50 ng m⁻³, which was slightly higher than that in the BS region. This indicates that a higher health risk was caused by PAH exposure in TJ than around the BS. The HMW-PAHs contributed most (90.9 %) of the total TEQ. However, the concentration of HMW-PAHs in the atmosphere accounted for 8 % of the total PAH concentration. Among the HMW-PAHs, the two major contributors were BaP (47.2 % ± 9.2 %) and DahA (19.7 % ± 16.2 %).

Information on the TEQ in the BS region and in TJ from June 2014 to May 2019 is shown in Fig. 5. The average value of TEQ in the BS region in each of the five annual cycles was 2.55 ± 1.49, 2.49 ± 1.63, 0.69 ± 0.76, 0.47 ± 0.66, and 0.67 ± 0.84 ng m⁻³, respectively. The value of TEQ in the BS region showed a downward trend year by year. The environmental health risk from PAHs in the fifth year was decreased by 3 times compared to the first year. It was found that the decrease in the HMW-PAH concentration was the main reason for the decrease in the toxicity of the PAHs. For example, the concentration of BaP in the atmosphere in the BS region decreased by 79.1 % in 5 years, and the concentration of DahA, a species with a carcinogenic toxicity equivalent to that of BaP, decreased by 96.1 %. In TJ, the average value of TEQ in each of the five annual cycles was 3.63 ± 0.14, 3.38 ± 0.72, 0.84 ± 0.38, 0.28 ± 0.10, and 0.31 ± 0.15 ng m⁻³, respectively. The TEQ value of the PAHs in the atmosphere decreased by 91.5 % in TJ during the 5 years. In TJ, BaP and DahA, the major factors contributing to the TEQ in the atmosphere, also showed more significant declines than they did around the BS. To sum up, the results showed that pollution control not only reduced the total concentration of PAHs in the atmosphere in the BS region, but it also affected the composition of the PAHs. It also mainly affected the concentration of HMW-PAH compounds, meaning that the total toxic equivalent of PAHs in the atmosphere in the BS region was remarkably reduced.

Simultaneously, the incremental lifetime cancer risk (ILCR) was used to assess the potential carcinogenic risk from PAHs in the atmosphere in the BS region. According to the USEPA, an ILCR value of less than 1 × 10⁻⁶ is an acceptable risk level. When the ILCR value is equal to or higher than 1 × 10⁻⁶ but less than 1 × 10⁻⁴, there is a serious risk of cancer, and the health issues should be taken seriously. When the ILCR value is equal to or greater than 1 × 10⁻⁴, it is considered life-threatening to humans. The specific calculation is shown in Eq. (4). It was found that the range of ILCR values for the atmospheric PAHs in the BS region across the 5 years was 4.1 × 10⁻⁵ to 2.2 × 10⁻⁴, with an average value of 1.2 × 10⁻⁴, which means that the risk of cancer in this region was at a serious level and more attention should be paid to the consequent health problems. Similarly to the TEQ, the ILCR values were dominated by the HMW-PAHs. The ILCR caused by PAHs is listed in Table S14. The ILCR in the BS region decreased significantly by 74.1 % from 2.2 × 10⁻⁴ in the first year to 5.7 × 10⁻⁵ in the fifth year. Compared with the BS region, the ILCR in TJ decreased more significantly, from 3.2 × 10⁻⁴ to 2.7 × 10⁻⁵ (i.e., by 91.6 %). As shown in Table S15, the study found that the concentration variations of the highly toxic BaP and DahA were basically synchronized with the changes in ILCR, which implied that the decrease in the concentrations of both was the main reason for the cancer risk reduction. The significant reduction in cancer risk in the BS region indicated that the emission of highly toxic HMW-PAHs in the atmosphere was effectively controlled, which also implied that the prevention and control of air pollution had effectively reduced the health risk. In particular, the reducing effect on the PAH exposure risk was more obvious in TJ, where the air pollution control was strict.

Figure 6 The medical costs of lung cancer caused by PAH exposure in TJ and around the BS from 2014 to 2018, before and after air pollution was controlled.

3.4 Direct medical costs of lung cancer caused by exposure to PAHs

This reduction in PAH-related health risk would lead to a reduction in the number of people who develop cancer, thus saving on the costs of cancer treatment. In this study, the direct medical costs of lung cancer caused by respiratory exposure to PAHs were estimated based on the additional incidence of lung cancer caused by PAH exposure, the population in the study area, and the direct medical costs per capita of lung cancer patients. The specific calculation is shown in Eq. (5). In addition to PAH exposure, there are many environmental risk factors that could induce lung cancer. To derive the lung cancer burden caused by respiratory exposure to atmospheric PAHs from the incidence of lung cancer, the percentage of population risk attribution (PAF) was performed in this study. The details are shown in Eqs. (6) and (7). Here, PAF represents the percentage reduction in lung cancer incidence that would occur if PAHs, an environmental factor, were to be completely removed or their concentration were to be reduced. According to the above introduction to the PAF and analysis of TEQ, the directly calculated PAF around the BS during the 5 years ranged from 0.5 ‰ to 2.7 ‰, with an average value of 1.4 ‰. The 5-year PAF in TJ ranged from 0.3 ‰ to 3.8 ‰, with an average value of 1.7 ‰. Remarkably, the PAF values around the BS region and in TJ decreased significantly during the 5 years, from 3.8 ‰ and 2.7 ‰ in the first year to 0.3 ‰ and 0.7 ‰ in the fifth year, respectively. The additional lung cancer incidence (I_add) due to respiratory exposure to PAHs was calculated using the product of lung cancer incidence and the PAF. Previous studies reported that the incidence of lung cancer in TJ in 2012 was 87.37 × 10⁻⁵ (Cao et al., 2016). In this study, 87.37 × 10⁻⁵ was used as the reference value of lung cancer incidence. The average I_add caused by respiratory exposure to PAHs around the BS region and in TJ was 1.26 × 10⁻⁶ and 1.55 × 10⁻⁶, respectively. During the 5-year observational period, I_add around the BS region and in TJ decreased from 2.34 × 10⁻⁶ and 3.33 × 10⁻⁶ in the first year to 6.15 × 10⁻⁷ and 2.87 × 10⁻⁷ in the fifth year, respectively. The population numbers in the study area were all sourced from public data in the statistical yearbook. The estimated results for the BS region and TJ are shown in Tables S16 and S17, respectively. It was reported that the direct cost of an average case of lung cancer in China in 2014 was USD 9042.79 (Shi et al., 2017; Huang et al., 2016). Since no reference data were available for the other years, we took the direct cost per case of lung cancer patients in 2014 as the baseline in this study, and the estimate assumed the same direct medical costs per capita for lung cancer across the 5 years.

Figure 6 shows the comparative results for the direct medical costs of lung cancer in the BS region and in TJ from 2014 to 2018, before and after pollution control. During the 5 years, under the implementation of air pollution control, the total direct medical costs of lung cancer caused by respiratory exposure to PAHs in the Bohai Rim region were USD 12.6 million. Assuming that no air pollution control was implemented, the total direct medical costs of lung cancer caused by PAH exposure did not change in the 5 years, and the total direct medical costs were USD 23.2 million. The effect of the actual implementation of air pollution control on the total direct medical costs of lung cancer was to save USD 10.7 million. In TJ, the total direct medical costs of lung cancer induced by respiratory exposure to PAHs under actual air pollution control were USD 1.0 million. Under the assumption that no air pollution control was implemented, the total direct medical costs of lung cancer caused by PAH exposure were USD 2.2 million, implying a saving of about USD 1.2 million in TJ. Compared to the costs incurred without air pollution control, the total direct medical costs of lung cancer caused by PAH exposure decreased by 46.1 % around the BS region and by an even greater 54.5 % in TJ. This illustrates that the implementation of air pollution control not only reduced the risk of lung cancer caused by PAH exposure around the BS region, but it also created a significant health benefit in the direct medical costs of lung cancer, especially in tightly controlled areas such as TJ. Therefore, the above results show that more precise pollution prevention and control could better reduce the emissions of pollutants and consequently reduce the health risk from human exposure.

4 Conclusions

Atmospheric PAH observations were conducted at 12 sites around the BS across 5 years from June 2014 to May 2019. The 5-year atmospheric concentration of Σ₁₅PAHs was 56.78 ± 4.75 ng m⁻³ and was characterized by a dominance of LMW-PAHs (58.7 ± 7.8 %). The maximum annual concentration and seasonal concentrations occurred in the first year and every winter, respectively. The concentrations of ∑₁₅ PAHs in the atmosphere around the BS reduced significantly during the sampling period, especially at the sampling site in TJ. The contributions of coal combustion and vehicle emissions to the PAHs in the atmosphere around the BS during the sampling period showed an increase and a decrease, respectively. However, the variations of coal combustion and vehicle emissions as source contributions in TJ showed just the opposite behavior. From 2014 to 2018, the additional lung cancer incidence caused by PAH exposure around the BS dropped by 74.1 %, and a higher drop of 91.6 % occurred in TJ. From a statistical standpoint, the drop in incidence saved about USD 10.7 million from the total direct medical costs of lung cancer caused by PAH exposure around the BS. Compared to the costs without air pollution control, the total direct medical costs of lung cancer caused by PAH exposure decreased by 46.1 % around the BS region and by an even greater 54.5 % in TJ. These results further certify that pollution reduction is beneficial to human health. In the fight against air pollution, more precise pollution prevention and control strategies are needed.