© The Author(s) 2023. Published by Science Press and EDP Sciences.

INTRODUCTION

Global climate change and its consequential health effects are among the most concerning environmental and public health issues [1–4]. Climate change, manifesting as changes in temperature, precipitation rate, and wind speed, can directly or indirectly affect the environmental fate of pollutants, and further affect human exposure and other related human health risks. The ubiquity of hydrophobic organic contaminants (HOCs) throughout the environment has put the ecosystem and human health at risk. HOCs are able to distribute in multiple media and accumulate in food webs due to their lipophilic characteristics. Climate change can alter the transport [5–7], volatilization [8–11], degradation [12,13], and bioaccumulation [14,15] of HOCs, thus affecting their distribution in air, freshwater, sediment, soil, and organisms. These changes can further affect human exposure to HOCs through inhalation, dermal contact, and food ingestion [16–20].

Existing investigations have primarily focused on the impact of climate change on the environmental levels and the fate of HOCs or human health risks induced by exposure to HOCs via a single exposure pathway [21–23]. Few studies have explored how climate change drives multimedia distribution, multi-pathway human exposure of HOCs, and the associated health risks at a continental scale in areas with large populations and across different climate zones, such as China. In addition to inhalation and dermal contact, dietary exposure is a significant exposure pathway for many HOCs. However, variations in food chain accumulation and dietary exposure to HOCs are seldom quantified under climate change scenarios. This makes it challenging to understand the climate change-related health burdens of HOCs in environments, particularly those attributed to dietary exposure [4,24]. Therefore, a critical research question for the scientific community is how climate change may impact the environmental fate of HOCs and the associated health risks in the future. Our goal is to link modules of climate change, multimedia environmental concentrations, and human exposure to better address this issue.

Polycyclic aromatic hydrocarbons (PAHs), a group of HOCs, are known for their carcinogenic, mutagenic, and teratogenic properties, posing a serious threat to human health [25]. Among the carcinogenic effects of toxic chemicals in the air, PAHs contribute the largest portion to lung cancer, mainly through respiratory exposure [26–28]. Exposure levels of PAHs and associated health risks are potentially high for the Chinese population, especially in areas with extensive emissions.

In this study, we comprehensively evaluate the human health risks induced by multi-pathway exposure to PAHs under two different climate change scenarios (representative concentration pathway (RCP) 4.5 and RCP8.5) with 1.5°C–4°C warming targets across the Chinese mainland by comparison with the base year 2009. We couple a well-verified spatially explicit high-resolution multimedia -SESAMe v3.4 model (Sino Evaluative Simplebox-MAMI model) with a bioaccumulation model to predict PAH multimedia concentrations in air, freshwater, sediment, soil, vegetation, and fish and further the resulting human intake exposed by inhalation, dermal contact, and oral intake. We discuss the driving mechanism of environmental processes to PAH multi-pathway exposure and associated health risks for people at different ages and across different climate zones in China under climate change scenarios. Finally, we provide emission reduction targets relative to the baseline year to reduce cancer risk based on the result of PAH exposure and human health risk with climate change.

RESULTS

Multimedia partitioning and spatial distribution of PAHs in the base year

In this study, we set 2009 as the base year, and validate the model by comparing a large set of observations with a sampling year around 2009 (1999–2020) collated from the literature with model predictions for 2009 (see Methods). The large observation dataset has ensured the more complete spatial coverage of the Chinese mainland targeted in this research. We find that the SESAMe v3.4 model is well verified. Most points cluster around the 1꞉1 line of prediction against observation, and the average of the root mean squared logarithmic error (RMSLE) and correlation coefficient is 0.51 and 0.49, respectively, for different PAHs and environmental compartments (Figure S1). Driven by meteorological data in 2009 taken from regional climate models (RCMs) under the RCP4.5 scenario (well verified, see Figure S2), PAHs are predicted to be mainly transported from atmosphere to soils and freshwater through deposition and sorption which are dominated transport processes at the air-land surface interface (Figures S3 and S4). Soil is the primary sink of PAHs. At the steady state, 66% on average of total PAHs are found in soil and 31% in sediment across China (Figure S5A). High-molecular-weight (HMW) PAHs, i.e., PAHs with four or more benzene rings, account for over 60% of total PAHs in all environmental compartments except the air, and even higher than 90% in sediment and soil (Figure S5B).

The model predicts total concentration ranges (5th–95th percentiles, plus median) of 16 PAHs in different media as follows: air, 0.3–392 ng m⁻³ (median, 27.9 ng m⁻³); soil, 0–626 ng g⁻¹(44.3 ng g⁻¹); freshwater, 0–2.24 × 10³ ng/L(80.5 ng/L); freshwater sediment, 0–1.28 × 10³ ng g⁻¹(42.1 ng g⁻¹); vegetation, 0–43.0 ng g⁻¹(1.9 ng g⁻¹); fish, 0–4.27 × 10³ ng g⁻¹(143 ng g⁻¹). The geographic distribution pattern of PAHs is generally similar for each compartment across the whole country, which is mainly driven by PAH emissions (Figure 1). The PAH concentrations are the highest in North China mainly because of the extensive use of biomass fuels for cooking and heating in the residential/commercial sector, the mass coke production in the industry sector, and the surge in motor vehicles in the transportation sector [29–31], especially for the cities and provinces with the highest emissions, including Beijing, Tianjin, Hebei, Shanxi, Shandong, and Henan, where the “2 + 26” cities are located (Figure S6). The “2 + 26” cities are Beijing, Tianjin, and 26 surrounding cities in North China proposed and focused in the Air Pollution Prevention and Control Action Plan enforced, where the air is the most contaminated in China [32]. Due to their significantly higher population densities and urbanization than the national average and vast contribution to coke production (51% in 2010 and 44% in 2020) [33], these regions produce most residential, vehicle and industrial emission density of PAHs [30,31]. The second highest concentrations are in South China, followed by Northwest China and the Qingzang Gaoyuan, which have the lowest PAH emissions mainly caused by lower human population density and social-economic activities. The PAH concentrations in Northwest China and the Qingzang Gaoyuan are mainly affected by the surrounding atmospheric transport [6], except in some areas with local PAH emissions such as northwest Xinjiang and the southernmost part of Xizang. The regional distribution pattern of PAH concentration is in accordance with previous observations in the air [34,35], the soil [36,37], and the freshwater [38].

Figure 1 Emission and modelled spatial distribution of PAHs for the base year. (A) The average PAH emissions (Mg year−1) in 2008–2010. The spatial distribution patterns of PAH concentrations in (B) air (ng m−3), (C) freshwater (ng L−1), and (D) soil (ng g−1) in the base year (2009) under the RCP4.5 scenario. Dark blue lines indicate four geographical regions of China, including North China, South China, Northwest China, and the Qingzang Gaoyuan. Light green lines indicate the provinces and cities where the “2 + 26” cities are located, including Beijing, Tianjin, Hebei, Shanxi, Shandong, and Henan. Data from Hainan and Taiwan are unavailable so that Hainan and Taiwan are uncolored in all the maps.

Effects of future climate change on the multimedia environmental fate of PAHs

We evaluated the influence of climate change under the RCP4.5 and RCP8.5 scenarios on the multimedia fate of PAHs. The RCP4.5 scenario incorporates climate change mitigation measures, while RCP8.5 is the baseline scenario that includes no climate change mitigation measures [39]. Future changes in meteorological data in China under the RCP4.5 and RCP8.5 scenarios are projected using the ensemble mean of an RCM downscaling simulation driven by four global climate models (GCMs) from CMIP5 (see Methods). Specifically, the global annual average surface air temperature (ASAT) rises by 1.5°C,2°C, and 3°C in 2028, 2043 and 2091, respectively, under the RCP4.5 scenario, and the global ASAT reaches an increase of 1.5°C,2°C,3°C, and 4°C in 2024, 2037, 2055, and 2073, respectively, under the RCP8.5 scenario (Figure S7). As shown in Table S1, the RCP8.5 scenario predicts an increase of 4.5°C in air temperature and an increase of 10.3% in precipitation rate in the 2090s compared with the 2009 base year, and the increase in temperature and precipitation rate under the RCP4.5 scenario is half of that under the RCP8.5 scenario.

The gridded environmental fate models forced under the RCP4.5 and 8.5 climate change scenarios for different warming years of 1.5°C–4°C are intended to represent future climate change years compared with the base year. To project the impact of future climate change on PAH fate in China, we assume a constant emission of PAHs after 2009 to remove the influence of emission changes. The projection indicates that changes in temperatures, wind speeds and precipitation rates have varying effects on the environmental fate of PAHs under different climate scenarios. More PAHs migrate from freshwater systems and soils to the air due to the greater volatilization with increasing temperatures compared with the base year, especially under the RCP8.5 scenario, although still with a dominant transport direction from air to land surfaces. The net flux at the air-land surface processes would decrease by 0–24% in different warming years across China (Figure S8). Climate change will affect PAH distribution among different environmental compartments. With a rising temperature from 1.5°C to 4°C, the proportion of PAHs would increase respectively by 6.9%–35.5% and 1.6%–8.2% on average in the atmosphere and vegetation, and would decrease respectively by 1.2%–5.9%, 0.1%–0.3%, and 0–0.2% in the water, soil, and sediment (Figure 2). Soils would still be the primary sink of PAHs in China under future climate change in both scenarios. Proportions of low-molecular-weight (LMW) PAHs would decrease while HMW PAHs would increase in all environmental compartments under future climate change (Figure S9).

Figure 2 The percentage change (%) of modelled PAH proportions in different media under different scenarios. (A) RCP4.5. (B) RCP8.5. Among them, soil includes natural soil, agricultural soil, and urban soil; vegetation includes natural vegetation and agricultural vegetation.

Across China, the total concentration of 16 PAHs would decrease in all environmental compartments on average when the temperature increases by 1.5°C–4°C. This decrease is mainly caused by accelerated degradation of PAHs in environments [16,20], with the degradation rate of 16 PAHs rising by 0.4%–33.0% in different environmental media under 1.5°C–4°C warming. Compared with the base year, the concentration of PAHs would averagely decrease by 40.0–161 ng L⁻¹ (4.3%–19.7%) in freshwater, 17.8–77.7 ng g⁻¹ (3.8%–17.7%) in freshwater sediment, 7.1–34.2 ng g⁻¹ (4.4%–20.4%) in soils, 0.3–1.7 ng g⁻¹ (3.0%–14.5%) in vegetation, 58.3–246 ng g⁻¹ (3.6%–16.8%) in fish, and 0.1–0.5 ng m⁻³ (0.1%–0.5%) in air under different warming thresholds (see Figure 3 and Figure S10). The percentage of concentration changes in the air would be 1–2 orders of magnitude lower than that in other environmental media (Figure 3C).

Figure 3 Changes in modelled PAH concentrations in different media. The changes in modelled concentrations of PAH (%) in (A) freshwater, (B) soil, and (C) atmosphere on average in China under the RCP4.5 and RCP8.5 scenarios. (D) The changes in modelled concentrations of LMW PAHs and HMW PAHs (%) in the atmosphere for different regions in China under the RCP4.5 and RCP8.5 scenarios.

However, in the atmosphere, it is worth noting that the concentrations of HMW PAHs would increase (0–17.5 ng m⁻³, 0–18.5% for different regions), while those of LMW PAHs would decrease (0–22.3 ng m⁻³, 0–3.9% for different regions) (Figure 3D). For example, under the RCP8.5 scenario, when the temperature increases by 4°C, the atmospheric concentration of phenanthrene, one of the LMW PAHs, would decrease by 0–6.4 ng m⁻³ (0–3.4%) in 99.8% of the country area. In contrast, the atmospheric concentration of benzo(a)pyrene, one of the HMW PAHs, would increase by 0–1.7 ng m⁻³ (0–47.2%) in 95.5% of the country area. The reversed changing pattern of phenanthrene and benzo(a)pyrene in the atmosphere should be a result of their different physical-chemical properties and response to rising temperature and precipitation in most regions as the main meteorological drivers. Phenanthrene has a vapor pressure nearly five orders of magnitude higher than benzo(a)pyrene and an air degradation rate of ca. 3.5 times that of benzo(a)pyrene. It is more likely to partition in the gaseous phase than the particulate phase, and degrades more rapidly in the air than benzo(a)pyrene. Meanwhile, the washout rate of gaseous-phase chemicals is over two orders of magnitude higher than particulate-bound chemicals. Therefore, although higher temperature promotes the volatilization rate of phenanthrene (14.0%–34.5%) more than that of benzo(a)pyrene (4.6%–29.6%) under the RCP8.5 scenario (with the rise of temperature by 4°C), the projected elevation in precipitation levels is expected to intensify the scavenging of phenanthrene to a greater extent than benzo(a)pyrene in the air. The greater reduction of the LMW PAHs (0–22.3 ng m⁻³) than the increase of HMW PAHs leads to the decreases of the total PAH concentrations in the air.

As shown in Figure 4 and Figures S11–S26, the spatial distribution of changes in PAH concentrations is consistent with the distribution of PAH emissions. For example, the concentration of 16 PAHs in the provinces and cities with the highest emissions where the “2 + 26” cities are located would have reduced the most under future climate change. For the change percentage of PAH concentration, its spatial distribution is controlled by climate change. The Qingzang Gaoyuan, where the polar tundra climate is the dominant climate type (60%) [40], is an area sensitive to climate change, with the largest changes in temperatures, precipitation rates, and wind speeds in China. As a result, the concentrations of PAHs in various environmental media in this area respond the strongest to climate change. When the temperature increases by 4°C under the RCP8.5 scenario (in 2073), the change percentage of PAH concentration in different environmental media is 3.9%–220% higher than that in other regions. The other three regions, North China, South China, and Northwest China, are dominated by cold climate (84%), temperate climate (95%), and arid climate (86%), respectively, defined by the temperature and precipitation level according to the climate classification map of Beck et al. [40]. Compared with the Qingzang Gaoyuan, the change percentage of PAH concentration in these three areas is lower due to its relatively weak sensitivity to climate change. Similarly, the spatial distribution of the increase in HMW PAHs and the decrease in LMW PAHs in the atmosphere is affected by both emissions and climate change.

Figure 4 Spatial distribution of change in modelled atmosphere PAHs concentration. The spatial distribution patterns of change in atmospheric concentrations (ng m−3) of (A) 16 PAHs, (B) LMW PAHs, and (C) HMW PAHs, when the temperature rises 4°C under the RCP8.5 scenario. Green lines indicate the provinces and cities where the “2 + 26” cities are located, including Beijing, Tianjin, Hebei, Shanxi, Shandong, and Henan. The spatial distribution patterns of percentage change of atmosphere concentrations (%) of (D) 16 PAHs, (E) LMW PAHs, and (F) HMW PAHs, when the temperature rises 4°C under the RCP8.5 scenario. Data from Hainan and Taiwan are unavailable so that Hainan and Taiwan are uncolored in all the maps.

Effects of future climate change on human exposure and health risks of PAHs

We quantify the impact of climate change on PAH intake through various human exposure pathways, including respiratory (inhalation of air and soil particles), dermal contact (bathing, swimming, and soil contact), and oral intake (drinking water, accidental soil ingestion, and dietary ingestion (agricultural vegetation and fish)). Benzo(a)pyrene-equivalent concentrations are used to calculate the cumulative lifetime average daily dose (LADD) and incremental lifetime cancer risk (ILCR) of PAH mixtures via multi-pathway exposure. Moreover, population-weighted ILCRs are calculated in each 0.5° grid cell to consider the potentially higher health impact of PAH exposure in densely populated areas (see Methods).

We find that oral intake is the most critical exposure pathway for PAHs in 99.6% of the country, contributing more than 50% (Q1–Q3, 97.4%–99.4%) of the multi-pathway human health risks in the base year 2009. Among the various oral intake pathways, fish ingestion (median, 61.9%; Q1–Q3, 30.7%–84.3%) contributes the most, followed by vegetation ingestion (38.0%; Q1–Q3, 15.7%–69.2%), while drinking water and accidental soil ingestion contribute the least (Table S2). The percentage of the risk induced by respiratory intake (1.4%; Q1–Q3, 0.7%–2.5%) is greater than that induced by dermal contact (0.1%; 0.1%–0.2%) in 94.7% of the area. Therefore, the order of human health risks exposed to PAHs through different pathways is oral intake > respiratory intake > dermal contact.

When considering all exposure pathways, the total cumulative ILCR of human health in the base year varies extensively with the three quartiles of 3.23 × 10⁻⁷, 4.20 × 10⁻⁶, and 1.56 × 10⁻⁴ across China (Table S3). Without considering individual susceptibility, people aged 30–60 have the greatest PAH exposure risks due to their more environmental exposure behaviours than other age groups. People aged 3–10 and 10–30 followed, with those aged less than 3 and more than 60 at the lowest risk. The cumulative health risks induced by all three exposure pathways exhibit a similar spatial pattern to health risks induced by individual exposure pathways. The risks are the highest in North China, particularly in the provinces and cities where the “2 + 26” cities are located, followed by South China, and are the lowest in the Qingzang Gaoyuan and Northwest China. In addition, ILCRs are generally higher on the east side of the Hu line (Figure S6), which is the result of its higher PAH concentration and population density, reflecting the higher risks in areas with higher PAH emissions and denser populations (Figure S27).

Under future climate change, the order of the main exposure pathways remains unchanged: oral intake > respiratory intake > dermal contact, and oral intake remains dominant in most areas (99.4% of the country), accounting for more than 48% (96.8%–99.2%) of the multi-pathway human health risks. However, the contribution of oral intake would decrease by 0.3% (Q1–Q3, 0.2%–0.6%), while respiratory exposure would increase by 24.2% (16.3%–31.7%) compared with the base year for a temperature increase of 4°C under the RCP8.5 scenario (in 2073).

The cumulative overall health risks of multi-pathway exposure would decrease with climate warming (Figure 5). Specifically, the overall health risks would decrease by 2.8% (median; Q1–Q3, 1.7%–4.5%), 5.4% (3.3%–7.5%) and 10.1% (6.3%–14.1%) when the temperature increases by 1.5°C,2°C and 3°C in 2028, 2043 and 2091, respectively, under the RCP4.5 scenario. The overall health risks would decrease by 3.2% (2.0%–4.7%), 5.5% (3.4%–7.6%), 10.8% (6.9%–14.7%) and 16.4% (10.3%–20.5%) when the temperature increases by 1.5°C,2°C,3°C and 4°C in 2024, 2037, 2055 and 2073, respectively, under the RCP8.5 scenario. The ILCR under the RCP8.5 (high CO₂ emissions) scenario decreases much more at the same warming threshold; the reduction in the ILCR under the RCP8.5 scenario is 1.2%–11.8% higher than that under the RCP4.5 scenario in the same warming threshold year. However, the cumulative overall health risks would still be higher than 1.00 × 10⁻⁶ with potential cancer risk in 66%–67% of China and higher than 1.00 × 10⁻⁴ with high cancer risk in 27%–28% of China under future climate change. These regions with high cancer risk are mainly located in the “2 + 26” urban and rural areas with high PAH emissions.

Figure 5 Changes in modelled ILCR through different exposure pathways. The changes in the modelled ILCR (%) through different exposure pathways on average under scenarios (A) RCP4.5 and (B) RCP8.5. Oral intake, including drinking water, dietary ingestion, and accidental soil ingestion; dermal contact, including bathing, swimming, and soil contact; respiratory intake, including inhalation of air and soil particles.

However, as a result of the increase in HMW PAHs in the atmosphere, the health risk induced by respiratory exposure would increase by 0.8% (0.2%–1.5%), 1.4% (0.8%–2.1%) and 2.8% (1.3%–4.0%) when the temperature increases by 1.5°C,2°C and 3°C in 2028, 2043 and 2091, respectively, under the RCP4.5 scenario. The health risk induced by respiratory exposure would increase by 0.9% (0.4%–1.5%), 1.5% (0.8%–2.2%), 2.9% (1.5%–4.2%) and 4.1% (2.9%–5.8%) when the temperature increases by 1.5°C,2°C,3°C and 4°C in 2024, 2037, 2055 and 2073, respectively, under the RCP8.5 scenario (Figure 5). The health risk induced by respiratory exposure in 2073 with the rise of temperature by 4°C under the RCP8.5 scenario varies extensively with the three quartiles of 9.18 × 10⁻¹⁰, 6.05 × 10⁻⁸, and 2.30 × 10⁻⁶, respectively across China (Table S4). The increase in the ILCR under the RCP8.5 scenario is 1.3%–16.2% higher than that under the RCP4.5 scenario in the same warming threshold year. In addition, the change in human health risks induced by inhalation is similar among different age groups under climate change scenarios.

Regarding the spatial distribution of risk change, the response of PAH-exposed human health risks to future climate change in the Qingzang Gaoyuan and Northwest China is stronger than that in other regions (see Figure 6 and Figure S28). The risk of respiratory exposure in the Qingzang Gaoyuan would increase by 7.0% (5.0%–9.4%) when the temperature increases by 4°C under the RCP8.5 scenario, which is 1.6 times higher than the national average level (Figure 6A). The multi-pathway exposure risk in the Qingzang Gaoyuan would decrease by 23% (18.2%–27.6%) when the temperature increases by 4°C under the RCP8.5 scenario, which is 1.3 times higher than the national average level. For the provinces and cities where the “2 + 26” cities are located with the highest PAH emissions, although the respiratory risks and multi-pathway exposure risks in this area are higher than those of other regions, the change percentages in respiratory risk and multi-pathway exposure risk are lower than the national average level with future climate warming.

Figure 6 Spatial distribution of modelled ILCR change. The spatial distribution patterns of change in the modelled ILCR (%) exposed to (A) respiratory intake, (B) dermal contact, and (C) oral intake, when the temperature rises to 4°C under the RCP8.5 scenario. (D) The spatial distribution change pattern in the total cumulative ILCR (%) induced by the above three exposure pathways when the temperature rises 4°C under the RCP8.5 scenario. Data from Hainan and Taiwan are unavailable so that Hainan and Taiwan are uncolored in all the maps.

DISCUSSION

This study represents the first comprehensive assessment of the potential effect of future climate change on the environmental fate of PAHs and their associated human health risks induced by three primary exposure pathways in China with different climate zones and high PAH emissions and population density. We find that the total concentrations of 16 PAHs would decrease in all environmental compartments, and the multi-pathway cumulative exposure ILCR would decrease with the warming climate. However, due to the increase in HMW PAHs in the atmosphere, projected ILCR through respiratory intake would increase, particularly under the RCP8.5 scenario. Our projection shows that dietary ingestion would be predominant among all exposure pathways, including inhalation, dermal contact, and oral intake, regardless of climate change. However, past research failed to include it in the human health risk assessment, which may have led to an underestimation of the cumulative human health risk, and caused the knowledge gap regarding how climate change will impact the dietary exposure of PAHs and associated human health risks. Therefore, human health risk assessment has been made with more comprehensive exposure pathways in our study. In addition, China covers a variety of climate types [40], and the responses of PAH exposure and associated human health risks to climate change are distinct in different regions at the continental scale. For example, our results indicate that in the Qingzang Gaoyuan, an area more sensitive to climate change, the PAH concentrations and associated health risks are likely to vary more significantly with climate change, compared with other regions. Our study is conducted at a continental scale and covered a variety of climate types in China, making it difficult to compare our findings with previous studies on different regions and populations [4,24,41] which focused on local and regional scales.

Our study finds a significant impact of temperatures on the fate and health effects of PAHs under future climate change. Temperatures, precipitation rates and wind speeds jointly affect the environmental fate of PAHs and related human health risks, among which temperature has the greatest impact. Considering only the change in meteorological factors, the variance contribution rate of temperatures to the multi-pathway cumulative risk would be 71.2%, with precipitation rate and wind speed only being 28.4% and 0.5%, respectively. Moreover, temperatures exhibit a negative correlation with the multi-pathway cumulative health risk, while the precipitation rates and wind speeds show a positive correlation with the health risk (Figure S29A). Under future climate change scenarios, the projected increase in the precipitation rate could promote deposition, and thus increase the concentration of PAHs in freshwaters and soils, which would lead to an increase in health risks caused by oral intake and dermal contact. However, the increase in temperatures would cause a decrease in the total PAH concentrations in all media to varying degrees by promoting the degradation, resulting in a reduction in health risks induced by all exposure pathways. According to the variance contribution rates mentioned above, the effect of temperatures is greater than that of precipitation rates, leading to a final decrease in health risk from all exposure pathways under climate change scenarios. In addition, among the environmental parameters, population density, local emissions, and water flow rates explain more than 85% of the variability in multi-pathway health risks (Figure S29B).

Uncertainties and limitations in this study are analysed as follows. As we aim to investigate the potential effects of climate change without the interaction of emissions and future PAH emission data over such a long period is not available, we assume constant PAH emissions in the 21st century in this study. However, the change in PAH emissions in the future could significantly affect human health risks. For instance, Nam et al. [4] expected that the overall cancer risks in Korea would increase by more than 50% nationwide in 2050 under the RCP8.5 scenario compared with 2018 based on the assumption that future PAH emissions would increase continuously at the country level. Taking the current declining emissions of PAHs in China into consideration [29], our results provide a conservative health risk assessment. More discussions about the uncertainty of the results caused by the variation of PAH emissions can be found in the Supplementary information C. In addition, the lack of temporal and spatial information on some exposure parameters, such as intake rates, may cause regional uncertainties in risk assessment, which deserves further investigation. Furthermore, to accommodate available spatial resolution of different environmental parameters for multiple compartments, the resolution of SESAMe v3.4 is 0.5° currently, which may make capturing regional pollution hotspots challenging compared with some atmospheric models. However, the resolution of SESAMe v3.4 is the highest among various large-scale spatial explicit multimedia models, which means it is an optimal choice as it enables the simulation of the multiple environmental compartments and exposure pathways in a national scale. Finally, our study was conducted under RCP climate change scenarios from CMIP5 rather than Shared Socioeconomic Pathway (SSPs) from CMIP6. However, both the CMIP6 multi-model ensemble and the CMIP5 multi-model ensemble offer good performance in simulating temperature and precipitation indices, and present similar trends and spatial patterns of climatic factors as well as the uncertainties in China over the 21st century [42,43]. Therefore, the coupling models with CMIP5 exhibit similar trends and spatial patterns of PAH concentrations and related human health risks compared with CMIP6, and our results under RCP scenarios from CMIP5 can reflect how future climate change could affect the multimedia environmental fate of PAHs, human exposure, and associated health risks.

Although our future climate change assessment is associated with uncertainties, our comparative analysis clearly shows an increasing tendency of inhalation-induced risk, and this may increase the likelihood of causing lung cancer and other respiratory diseases due to the significant contribution of PAH to lung cancer through respiratory exposure [26‒28]. In addition, although the cumulative overall health risk of multi-pathway exposure to PAHs shows a decreasing trend in the context of future climate change, the value of health risk is still very high, resulting in potential cancer risk. We need to start from at least two aspects to reduce potential cancer risk. First, our findings indicate that when the temperature increases by 1.5°C‒4°C under the RCP8.5 scenario, the risk of cancer would possibly increase by 0.2%‒5.8% through respiratory intake exposure compared with the baseline year, which is 1.3%‒16.2% higher than that at the same warming threshold under the RCP4.5 scenario. Restricting the global mean temperature rise to 1.5°C will lower this increasing proportion to less than 1%, regardless of RCP4.5 and 8.5.

Second, our findings suggest that it is necessary to fundamentally cut PAH emissions to reduce PAH-induced cancer risk in the future considering climate change. If we are committed to reducing the median multi-pathway exposure cumulative overall ILCR value of China (3.48 × 10⁻⁶ when the temperature rises by 4°C under the RCP8.5 scenario) to be at a safe level (less than 1.00 × 10⁻⁶), the total PAH emissions of China need to be reduced by 69.2% compared with the 2009 base year, and PAH emissions should be less than 170 kg/(km² year) for each grid. In addition, as revealed in our study, the PAH associated human health risks in different regions followed a descending order as North China > South China > Northwest China and Qingzang Gaoyuan. The PAHs emission in these regions should be reduced by 75.1%, 75.0%, 69.3%, and 0% (no need to reduce), respectively, compared with the 2009 base year to achieve a safe level of multi-pathway exposure cumulative overall ILCR value under future climate change scenarios (see Figure 7A). Furthermore, for 30 cities and provinces in the Chinese Mainland, the emission reduction percentages compared with the base year range from 0% (no need to reduce) to 75.3% (see Figure 7B). The cities and provinces including Beijing, Tianjin, Hebei, Shanxi, Shandong, and Henan where “2 + 26” cities are located as well as Shanghai, Jiangsu and Anhui should make the greatest efforts to reduce PAH emissions, with the reduction percentage higher than 75.2%. However, Xizang and Qinghai do not need emission reduction measures due to the low human health risks under future climate change based on our results.

Figure 7 PAH emission reduction percentages of regions and provinces under future climate change. The PAH emission reduction targets of reducing the median multi-pathway exposure cumulative overall ILCR values to be at a safe level (less than 1.00 × 10−6) compared with the emission in the base year when the temperature rises to 4°C under the RCP8.5 scenario for (A) four geographical regions (including North China, South China, Northwest China, and the Qingzang Gaoyuan), and (B) 30 provinces (excluding Hong Kong, Macao, Hainan, and Taiwan) of China. Data from Hainan and Taiwan are unavailable so that Hainan and Taiwan are uncolored in all the maps.

In fact, some measures have been taken to control PAH emission sources in above mentioned areas, such as the residential/commercial sector. Policies such as the Winter Clean Heating Plan in North China (available at: http://www.gov.cn/xinwen/2017-12/20/content_5248855.htm) aim to reduce indoor biomass burning for daily cooking and heating which contribute more than half of PAH emissions in China in 2007 [29]. Clean energy has been replacing biomass fuels in growing rural areas with local residents being subsidized. While the phase-out of solid fuels for heating and cooking in North China has great potential for reducing PAH emissions, further evaluation is needed to determine its effectiveness in reducing health risks of multi-pathway exposure to PAHs.

Our study based on the coupling models has provided a conservative assessment of PAH-related human health risks under future climate change scenarios. However, climate change would be accompanied and affected by changes in socioeconomic factors, for instance, land use types, population, energy structure, and pollutant emissions. Therefore, the future investigation on the use of CMIP6 models under SSP, which have additionally considered varieties of more socioeconomic factors, is proposed for a more comprehensive assessment in the next-step research. Besides, considering the geographically different impact patterns of climatic factors, we also suggest taking account of other social variables, such as human population and its alteration, human movement (e.g., travel or migration), human behavior patterns in relation to exposure (e.g., hand-to-mouth behaviors, frequencies of washing etc.) and so on, in predicting the future health burden of HOC exposure quantitatively. More studies are required to appropriately assess the complex combination of climate change and dynamic human activities in the future to promote a better understanding of HOCs fate and associated human health risks, thereby providing support for making measures and policy decisions.

METHODS

Prediction of PAHs in multiple compartments using SESAMe v3.4 combined with a bioaccumulation model

To analyse the effect of climate change on the environmental fate of PAHs in multiple compartments, we used an improved version of our previously published model, the SESAMe v3.4 model [15,44,45]. It is a steady-state spatially explicit environmental fate model for the Chinese mainland that utilizes (1) emissions of PAHs, (2) physical-chemical properties of PAHs, and (3) environmental parameters, to estimate the concentration of chemicals in different environmental compartments with a spatial resolution of 0.5°. It includes four basic processes of PAH behaviours: emissions to the environment, advective transport, diffusive transport, and degradation. The model has been well verified for a wide range of trace organic chemicals, including PAHs [16,44–47]. However, the aquatic organism phase is not included in the current version of the SESAMe model. In this study, a bioaccumulation model is loosely coupled with the SESAMe v3.4 model to investigate the concentrations of PAHs in aquatic organisms such as fish, which is one of the predominant sources of dietary exposure pathway for humans. We used bioaccumulation factors (mid trophic, including biotransformation rate estimates) estimated by the Arnot-Gobas model (obtained from EPI Suite, a chemical property estimation tool [48]) and each grid-mean freshwater concentration of PAHs to calculate the concentrations of PAHs in fish in each grid cell, which provides general estimates in absence of site-specific measurements or estimates (not for a particular fish species or size) [49–52].

In these calibrated and validated models, we substituted present-day climate with future climate projections to simulate environmental PAHs with new climate data forcing. It can be efficiently coupled with climate models and is useful for evaluating the impacts of climate change projections on the fate of PAHs and the associated human health risks.

Input data

Emission data from 2008–2010 were taken from the PKU-PAHs-v2 emission inventory at a spatial resolution of 0.1°, where PKU stands for Peking University (available at http://inventory.pku.edu.cn/ [29]). The sixteen parent PAHs listed among U.S. EPA priority pollutants [53], including naphthalene (NAP), acenaphthylene (ACY), acenaphthene (ACE), fluorene (FLO), 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(l,2,3-cd)pyrene (IcdP), benzo(g,h,i)perylene (BghiP), were chosen as the case HOCs for this study (Table S5). PAHs with four or more benzene rings are generally considered HMW PAHs, showing potentially genotoxic and carcinogenic effects. They are more recalcitrant than LMW ones with three or fewer benzene rings, persisting in the environment longer [54].

In this study, the driving data of regional surface meteorological elements in 2008–2010, including 2-meter air temperature, 10-meter wind speed, and precipitation rate, were downloaded from the China Meteorological Forcing Dataset (1979–2018) with a spatial resolution of 0.1° in China [55]. Future changes in meteorological factors in China under global warming in this century for the RCP4.5 and RCP8.5 scenarios were generated by an ensemble mean of RegCM4 (the latest version of RCM) downscaling simulations, which were driven by four GCMs, including CSIRO Mk3.6.0, EC-EARTH, HadGEM2-ES and MPI-ESM-MR from CMIP5. Among them, RCP8.5 is a baseline scenario that does not include any climate mitigation policies, while RCP4.5 represents a medium mitigation scenario [39]. In brief, first, the annual average data of each grid of each GCM were calculated on the determined time scale, and then the global average data were calculated on the spatial scale to obtain the annual average data of the global grid of each GCM. Second, we computed anomalies relative to the base period (1860–2005) and its 11-year smoothed global mean surface air temperature derived from the ensemble mean of the above four driven GCMs for the RCP4.5 and RCP8.5 scenarios from 1860 to 2099, which showed that the ensemble global ASAT of GCMs would increase by 1.5°C,2°C, and 3°C in 2028, 2043, and 2091 for the RCP4.5 scenario and increase by 1.5°C,2°C,3°C and 4°C in 2024, 2037, 2055, and 2073 for the RCP8.5 scenario, respectively (Figure S7). Similarly, the calculation of RCMs data driven by GCMs was the same, except that the base period was 1986–2005, and the same method was used to calculate the wind speed and precipitation. Finally, we extracted 2009 for the RCP4.5 scenario as “the present year” and year data corresponding to the above temperature increase threshold (1.5°C–4°C) for the RCP4.5 and RCP8.5 scenarios as “the future years” input data from processed RCMs data. The meteorological data of the present and future years were used as inputs to drive each SESAMe v3.4 model after bias adjustment at the yearly timescale from the China Meteorological Forcing Dataset (1979–2018). The bias adjustment for precipitation rate and wind speed used ratio method, and the temperature fits well with the observed data without adjustment. It should be noted that these datasets were regridded to match the spatial resolution of the SESAMe v3.4 model, which is 0.5° × 0.5°, yielding a representative year of climatic conditions for each scenario that could be used to drive dynamic model calculations.

Projection scenarios

To evaluate the response of PAH environmental concentration levels to future climate change, in this study, we ran three group simulations with the SESAMe v3.4 model with constant PAH emissions: (1) one with the China Meteorological Forcing Dataset for temperature, wind speed and precipitation rate in the yearly mean of 2008–2010 to validate the model, where each of the three years is input into the model separately, (2) one with the RCP4.5 for temperature, wind speed and precipitation rate in the yearly mean of 2009 after being corrected with the observed meteorological data from Simulation (1) to present “the present period” or the base year, (3) the other with the RCP4.5 and RCP8.5 for future meteorological data of 11-yearly smoothed ensemble mean in the year under 1.5°C‒4°C to present “the future period”. Differences between projected future period concentrations and baseline concentrations produced from model simulation data of the RCP were used to describe changes in PAH concentrations.

Model validation

Annual emissions and values of meteorological variables during 2008‒2010 in addition to other environmental data were input into SESAMe v3.4 and bioaccumulation models to predict multimedia environmental concentrations of the 16 PAHs at the steady state in the individual years. The predictions were then compared with the measurements across China, collated from the literature. We conduct a complete literature search using Web of Science and China National Knowledge Infrastructure (available at https://www.cnki.net/) to collect measurements of the 16 PAHs in air, soil, freshwater, freshwater sediment, vegetation and fish with a sampling year ranging from 1999 to 2020 (ten years before and after 2009). A total of approximately 7200 sets of available data for all compartments were collected from around 1300 pieces of literature for the model validation. The RMSLE and correlation coefficient were used to evaluate the relative error and the significance of correlation between the predicted and the observed values. In this study, we insist on using 2009 as the base year, which can enable us to use the extensive database covering 1999‒2020 collected from the literature for more reliable model validation and projection on a long-time scale in the future.

Health risk assessment

We applied the ILCR method to estimate gridded PAH-induced human health risks for each year of warming under the RCP4.5 and RCP8.5 scenarios. To determine human exposure to all the 16 PAHs through multiple exposure pathways, concentrations of individual PAHs in multimedia compartments as previously mentioned are predicted by SESAMe v3.4 and bioaccumulation models, and then converted to benzo(a)pyrene equivalent concentrations using toxic equivalency factors [56]. The population-weighted ILCR was calculated as the product of the ILCR and population density in each 0.5° × 0.5° grid cell divided by the nationally averaged population density over the land area. The population density data used in this study were obtained from the China Population Spatial Distribution Dataset with a grid resolution of 1 km in 2010 (https://www.resdc.cn/DOI/DOI.aspx?DOIid=32) [57]. We assumed that the baseline population density data remain unchanged across the warming scenarios when applied to conduct the population-weighted health risk assessment under climate change scenarios. The population-weighted ILCR of a region with more than one grid cell was calculated as Σ(P_i·C_i)/P_t, where P_i and C_i are population and ILCR at the ith grid cell, and P_t is the total population in the region, which represents the average exposure situation [58,59].

The ILCR was calculated by multiplying the cancer slope factor (CSF) of PAHs and LADD to conservatively estimate the potential lifetime cancer risk associated with PAH exposure [60] as follows:

where a and g are subscripts representing age and gender: the population was divided into 14 groups according to gender and age: all the lifetimes (all), infants (<3 years), children (3–10 years), teens (10–18 years), the young (18–30 years), the middle-aged (30–60 years) and the old (>60 years);C_eq (ppm) is modelled benzo(a)pyrene-toxic equivalent concentration in different environmental compartments for 16 USEPA PAHs; IR (volume/day) is the intake rate; ED (year) is exposure duration; EF (350 day/year [61]) is exposure frequency; ET (h/day) is exposure time; BW (kg) is the body weight; AT (day) is the average exposure duration, here refers to the average life expectancy of the Chinese population (74.7 years for males and 80.5 years for females based on the WHO reports [62] (available at https://apps.who.int/gho/data/node.main.688?lang=en).

Three main exposure pathways of PAHs to the human body were identified, including oral intake, dermal contact and respiratory intake, and they have slight differences in the calculation of human health risk. Among them, the IR of inhalation was calculated on the basis of oxygen consumption associated with energy expenditures [21], and other exposure parameters are shown in Table S6. We calculated the cumulative risks of multiple exposure pathways by adding them together for a comprehensive evaluation of changes in future health risks. According to the USEPA guideline [63], one in a million chance of human cancer risk (ILCR = 1.00 × 10⁻⁶) is considered an acceptable threshold, 1.00 × 10⁻⁶–1.00 × 10⁻⁴ means a potential cancer risk, and ILCR ≥ 1.00 × 10⁻⁴ indicates high risk of cancer and needs to be of particular concern.

Uncertainty analysis

We analysed the uncertainty of the coupled model of the high-resolution climate model, the environmental chemical fate model, and the health risk model, introduced by the variance of input parameters using the Monte Carlo simulation. The input parameters, including emission vectors of the 16 PAHs and 20 environmental variable vectors, were evaluated by the Kolmogorov-Smirnov test and found to conform to normal or log-normal distributions (Figure S30). Random values for these input parameters were generated based on their probability distributions for 10,000 runs to obtain the total cumulative ILCRs induced by all exposure pathways for individual PAHs. As shown in Figure S31, The total cumulative ILCRs conform to log-normal distributions, and the interquartile range of the total cumulative ILCRs is almost within one order of magnitude for each PAH.

Note: Data associated with this study, including source data files for Figures 1–7, are available within the main text and Supplementary Information files or from corresponding authors upon reasonable request. The dataset referenced and analysed during this study is available (with some institutional limitations) from the corresponding authors and institutions, which are clearly described in the Methods section. Data processing was conducted using MATLAB R2017b. All figures were created using Origin 2018, and ArcMap 10.2 was used to extract gridded data and create maps in Figures 1, 4, 6 and 7, and Figures S6 and S11–S28. Uncertainty analysis and sensitivity analysis were performed using the Oracle Crystal Ball software.

Acknowledgments

We gratefully acknowledge Prof. Shu Tao from Peking University for sharing PAH emissions data and his suggestions to improve the paper.

Funding

This work was supported by the National Key Research and Development Program of China (2017YFA0605001) and the National Natural Science Foundation of China (52039001, 92047303 and 41977359).

Author contributions

X.X. and Y.Z. designed the research. J.B. built and ran the models. J.B., B.Z., E.G., Y.L., Z.Z., N.X. compiled literature data. J.B. and X.Y. handled the emission data. X.X., Y.Z., J.B., B.Z., Z.Z., N.X., Y.G., Y.X., Z.Y., D.M. performed data analysis and wrote the manuscript with significant contributions from all authors.

Conflict of interest

The authors declare no conflict of interest.

Supplementary information

The supporting information is available online at https://doi.org/10.1360/nso/20230010.

References

All Figures

Figure 1 Emission and modelled spatial distribution of PAHs for the base year. (A) The average PAH emissions (Mg year−1) in 2008–2010. The spatial distribution patterns of PAH concentrations in (B) air (ng m−3), (C) freshwater (ng L−1), and (D) soil (ng g−1) in the base year (2009) under the RCP4.5 scenario. Dark blue lines indicate four geographical regions of China, including North China, South China, Northwest China, and the Qingzang Gaoyuan. Light green lines indicate the provinces and cities where the “2 + 26” cities are located, including Beijing, Tianjin, Hebei, Shanxi, Shandong, and Henan. Data from Hainan and Taiwan are unavailable so that Hainan and Taiwan are uncolored in all the maps.

In the text

	Figure 2 The percentage change (%) of modelled PAH proportions in different media under different scenarios. (A) RCP4.5. (B) RCP8.5. Among them, soil includes natural soil, agricultural soil, and urban soil; vegetation includes natural vegetation and agricultural vegetation.
In the text

	Figure 3 Changes in modelled PAH concentrations in different media. The changes in modelled concentrations of PAH (%) in (A) freshwater, (B) soil, and (C) atmosphere on average in China under the RCP4.5 and RCP8.5 scenarios. (D) The changes in modelled concentrations of LMW PAHs and HMW PAHs (%) in the atmosphere for different regions in China under the RCP4.5 and RCP8.5 scenarios.
In the text

Figure 4 Spatial distribution of change in modelled atmosphere PAHs concentration. The spatial distribution patterns of change in atmospheric concentrations (ng m−3) of (A) 16 PAHs, (B) LMW PAHs, and (C) HMW PAHs, when the temperature rises 4°C under the RCP8.5 scenario. Green lines indicate the provinces and cities where the “2 + 26” cities are located, including Beijing, Tianjin, Hebei, Shanxi, Shandong, and Henan. The spatial distribution patterns of percentage change of atmosphere concentrations (%) of (D) 16 PAHs, (E) LMW PAHs, and (F) HMW PAHs, when the temperature rises 4°C under the RCP8.5 scenario. Data from Hainan and Taiwan are unavailable so that Hainan and Taiwan are uncolored in all the maps.

In the text

Figure 5 Changes in modelled ILCR through different exposure pathways. The changes in the modelled ILCR (%) through different exposure pathways on average under scenarios (A) RCP4.5 and (B) RCP8.5. Oral intake, including drinking water, dietary ingestion, and accidental soil ingestion; dermal contact, including bathing, swimming, and soil contact; respiratory intake, including inhalation of air and soil particles.

In the text

Figure 6 Spatial distribution of modelled ILCR change. The spatial distribution patterns of change in the modelled ILCR (%) exposed to (A) respiratory intake, (B) dermal contact, and (C) oral intake, when the temperature rises to 4°C under the RCP8.5 scenario. (D) The spatial distribution change pattern in the total cumulative ILCR (%) induced by the above three exposure pathways when the temperature rises 4°C under the RCP8.5 scenario. Data from Hainan and Taiwan are unavailable so that Hainan and Taiwan are uncolored in all the maps.

In the text

Figure 7 PAH emission reduction percentages of regions and provinces under future climate change. The PAH emission reduction targets of reducing the median multi-pathway exposure cumulative overall ILCR values to be at a safe level (less than 1.00 × 10−6) compared with the emission in the base year when the temperature rises to 4°C under the RCP8.5 scenario for (A) four geographical regions (including North China, South China, Northwest China, and the Qingzang Gaoyuan), and (B) 30 provinces (excluding Hong Kong, Macao, Hainan, and Taiwan) of China. Data from Hainan and Taiwan are unavailable so that Hainan and Taiwan are uncolored in all the maps.

In the text