Global Cancer Risk From Unregulated Polycyclic Aromatic Hydrocarbons.

In assessments of cancer risk from atmospheric polycyclic aromatic hydrocarbons (PAHs), scientists and regulators rarely consider the complex mixture of emitted compounds and degradation products, and they often represent the entire mixture using a single emitted compound-benzo[a]pyrene. Here, we show that benzo[a]pyrene is a poor indicator of PAH risk distribution and management: nearly 90% of cancer risk worldwide results from other PAHs, including unregulated degradation products of emitted PAHs. We develop and apply a global-scale atmospheric model and conduct health impact analyses to estimate human cancer risk from 16 PAHs and several of their N-PAH degradation products. We find that benzo[a]pyrene is a minor contributor to the total cancer risks of PAHs (11%); the remaining risk comes from other directly emitted PAHs (72%) and N-PAHs (17%). We show that assessment and policy-making that relies solely on benzo[a]pyrene exposure provides misleading estimates of risk distribution, the importance of chemical processes, and the prospects for risk mitigation. We conclude that researchers and decision-makers should consider additional PAHs as well as degradation products.

Introduction

Polycyclic aromatic hydrocarbons (PAHs) are a class of chemicals that contain multiple fused aromatic rings, and are emitted into the atmosphere as byproducts of burning organic matter (Keyte et al., 2013). Several PAHs have been identified as mutagenic or carcinogenic (Bostrom et al., 2002) and therefore have the potential to harm the health of humans (Hansen et al., 2007; Park & Park, 2009) and ecosystems (Gray, 2002). In the atmosphere, PAHs are present as a complex mixture, with different people inhaling different combinations of these carcinogens (Dixon et al., 2019). However, scientific research and environmental guidelines often represent this complex PAH mixture using a single surrogate compound, benzo[a]pyrene (BAP).

Epidemiological (Armstrong et al., 2004; Moolgavkar et al., 1998) and animal (Collins et al., 1991; Heinrich et al., 1994; Thyssen et al., 1981) studies have been used to estimate the cancer risk of human exposure to BAP, even though humans are exposed to many different kinds of PAHs. Risk estimates derived from epidemiological studies imply that exposure of 1 ng/m3 of BAP to a population of 1 million people will induce 213–850 cancer cases over their lifetime (∼70 years). This is equivalent to an epidemiologically derived unit risk (URE) that range from 213–850 × 10−6 per (ng/m3) (Armstrong et al., 2004; Moolgavkar et al., 1998). Existing epidemiological studies have not accounted for the confounding exposure to other PAHs, or exposure to other pollutants such as heavy metals, sulfur dioxide, and nitrogen oxide; applying URE thus implicitly reflects the impact of exposure to a mixture of PAHs. Cancer risks derived from studies in which animals were exposed to BAP alone are lower than those derived from epidemiological studies (URA = 0.3–1.7 × 10−6 per [ng/m3]) (Collins et al., 1991; Thyssen et al., 1981) suggesting an increased risk associated with other components of PAH mixtures.

In global‐scale modeling and air quality management, scientists and policy makers frequently use BAP as an indicator for calculating risk from the entire PAH mixture. Shen et al. (2014) and Shrivastava et al. (2017) estimate global‐average human lung cancer risks of 20–31 × 10−6 by combining global‐scale models of BAP in the atmosphere with an epidemiologically derived BAP unit risk. However, although estimates of URE vary by a factor of ∼4, both Shen et al. (2014) and Shrivastava et al. (2017) test only the upper limit of this value. By using BAP as the sole indicator, these studies also assume that variations in BAP concentrations reflect proportional variations in risk. While PAHs are regulated as a class of substances, national and international governing bodies also use BAP as an indicator for all species: BAP is the only PAH to have a guideline concentration.

In regional‐scale modeling and field campaigns, scientists are able to capture a greater range of PAHs, but still neglect several highly toxic PAHs, and may double count the health effects of PAHs in cancer risk assessments. Toxic equivalent quotients (TEQs) are the toxicity of a compound relative to a reference compound, and in the case of PAHs, the reference compound is BAP. Several studies estimate the cancer risk of individual components of PAHs mixtures by first using the TEQ to convert the PAH concentration to an equivalent concentration of BAP, and then combining it with estimates of URA, through which the total cancer risk from PAHs can be estimated. In those studies, BAP is estimated as the major contributor to the cancer risk of PAHs (40%–80%) (Delgado‐Saborit et al., 2011; Nielsen et al., 1996; Norramit et al., 2005; J. Zhang et al., 2016; Y. X. Zhang et al., 2009), but many of those studies do not include highly toxic emitted PAHs (e.g., dibenz[a,h]anthracene), and none include degradation products. An additional complication is that several such studies, however, use URE together with TEQs to assess combined risks from the mixture, which leads to a potential double‐counting (Elzein et al., 2019; Zhang et al., 2016).

Recent work suggests that the atmospheric degradation products of PAHs, such as those containing a nitro group (‐NO2) that we refer to here as N‐PAHs, are highly toxic, but their impact on human health remains uncertain. N‐PAHs, including both nitro‐PAHs (one ‐NO2 group) and dinitro‐PAHs (two ‐NO2 groups), can be up to 1,000 times more toxic than their respective parent compound (Wislocki et al., 1986). Laboratory studies have shown that N‐PAHs are formed under several different oxidation reactions in the atmosphere (Keyte et al., 2013). N‐PAHs have been detected in a variety of environments, from urban (Albinet et al., 2007; Elzein et al., 2019) to remote (Drotikova et al., 2020; Lammel et al., 2017). The chemical formation of N‐PAHs in the atmosphere has been simulated in regional‐scale (Mulder et al., 2019) and global‐scale (Wilson, 2020) atmospheric modeling studies; those studies, however, did not include a key particle‐phase reaction with the nitrate radical (NO3) that laboratory studies suggest could be an extremely efficient source of N‐PAHs (Zelenov et al., 2018; P. Zhang et al., 2014). Furthermore, as prior studies did not perform human health impact assessments, the importance of N‐PAHs in the context of human health, as well as in relation to other advancements in this field, are both unknown.

Researchers in previous studies have argued that uncertainties in heterogeneous oxidation kinetics (Poschl et al., 2001; S. Zhou et al., 2019; S. M. Zhou et al., 2013) and gas‐particle partitioning (Dachs & Eisenreich, 2000; Shahpoury et al., 2016) have a large effect on exposure and human health impacts of BAP, but the importance of these uncertain processes in the context of other emitted PAHs and degradation products remains unknown. Shrivastava et al. (2017) found that reducing the rate of BAP heterogeneous oxidation resulted in a five‐fold increase in estimated human exposure to BAP, and a three‐fold increase in PAH human cancer risk when using BAP as an indicator of health risks. However, reductions in the oxidation rate will also diminish human exposure to degradation products, which may themselves be toxic—this effect was not considered by Shrivastava et al. (2017), who used BAP concentrations as a proxy for overall PAH exposure and risk. Additional studies (Friedman, Pierce, & Selin, 2014; Friedman & Selin, 2012; Friedman, Zhang, & Selin, 2014; Mu et al., 2018; Thackray et al., 2015) have advanced understanding of PAH chemistry in the atmosphere through simulations of 1–3 PAH species. The sensitivity of process‐based conclusions to the inclusion of additional PAHs and degradation products remains unassessed.

Here, we reevaluate the importance of BAP and its suitability as an indicator compound for global‐scale cancer risk of PAHs, and we perform a bounding exercise to assess how uncertainties in concentrations and atmospheric processes affect conclusions drawn in previous studies which were based on BAP alone. To do this, we develop and use a global‐scale atmospheric chemistry model to estimate concentrations of PAHs in the atmosphere, and use human health impact analyses to estimate the human cancer risk associated with atmospheric PAHs, using traditional epidemiologically based functions as well as a novel framework based on toxicity data from animal studies which allows us to explicitly estimate the risks of individual components of the PAH mixture. In contrast to previous global‐scale modeling studies which typically considered only BAP (Friedman, Pierce, & Selin, 2014; Friedman & Selin, 2012; Friedman, Zhang, & Selin, 2014; Lammel et al., 2009; Mu et al., 2018; Octaviani et al., 2019; Shen et al., 2014; Shrivastava et al., 2017; Thackray et al., 2015), we account for the 16 PAHs identified as priority pollutants by the United Stated Environmental Protection Agency (“USEPA16”) as well as N‐PAH degradation products. We quantify the sensitivity of estimates of global cancer risk to (a) inclusion of additional PAHs, (b) gas‐particle partitioning, (c) heterogeneous oxidation kinetics, and (d) model resolution. We conclude that BAP accounts for only a small (11%) fraction of the human cancer risk of PAHs globally, while the N‐PAHs, which are unregulated and commonly omitted by measurement and modeling studies focused on the atmosphere, are a potentially large (17%) source of carcinogenic risk. We also find that atmospheric kinetic and partitioning uncertainties have a much lower impact on risk magnitudes than was identified in previous studies. We suggest that future research and regulatory guidelines explicitly consider a broader range of PAHs and their degradation products in assessments of cancer risk from these compounds.

Materials and Methods

Experimental Design

We provide a total of 14 estimates of global human cancer risk from PAHs by combining seven different PAH concentration distributions from a global‐scale atmospheric chemistry model, with two distinct methods for estimating human cancer risk. Together, this analysis allows us to quantify the suitability of BAP as an indicator of risk of PAH mixtures, and to bound the importance of contributors to uncertainty.

Description of Atmospheric Chemistry Model (GEOS‐Chem)

We use a numerical, three‐dimensional, atmospheric chemistry model, GEOS‐Chem (Bey et al., 2001) version 11. We perform model simulations at two different horizontal resolutions, of 4° × 5° and 2° × 2.5°—both with 47 vertical levels, extending from the surface of the Earth to ∼80 km altitude. Meteorological fields are driven by MERRA‐2 reanalysis from the NASA Global Modeling and Assimilation Office (Global Modeling and Assimilation Office [GMAO]). We use the “tropchem” chemical mechanism (Eastham et al., 2014) for this work, which includes a detailed and fully coupled treatment of HOx–NOX–VOC–O3 and a bulk aerosol scheme with fixed log‐normal modes. Aerosol components considered include sulphate, ammonium, nitrate, sea salt, black carbon, mineral dust, and organic carbon. Inorganic aerosol thermodynamics are calculated using ISORROPIA (Fountoukis & Nenes, 2007). Both OC and BC are further separated into hydrophilic and hydrophobic components. The tropchem mechanism is expanded to also consider PAHs (described below Section 2.3). Dry deposition of both gases and aerosol are parameterized in a scheme which applies a resistance‐in‐series approach (Wesely, 1989; L. M. Zhang et al., 2001). Wet deposition occurs both within and below clouds, and is dependent on the species‐specific effective Henry's Law constant (Amos et al., 2012). These model simulations use a variety of different global and regional‐scale emission inventories for non‐PAH species. Global‐scale emission inventories used here for non‐PAH species include EDGAR (Crippa et al., 2018) and RETRO (Hu et al., 2015). Where necessary, these global‐scale emission inventories are overwritten with regional‐scale emission inventories—e.g., NEI over USA (Travis et al., 2016), EMEP over Europe (van Donkelaar et al., 2008). For biomass burning and biogenic emissions, we use GFED4 (Giglio et al., 2013) and MEGAN (Guenther et al., 2012), respectively. Emissions of PAHs are discussed in greater detail in the following section. In the next sections, we describe how PAHs are simulated within the GEOS‐Chem model, highlighting new developments.

Treatment of PAHs and N‐PAHs in GEOS‐Chem

The GEOS‐Chem model has previously been used to examine PAH chemistry and transport in the atmosphere, both for the present day (Friedman, Pierce, & Selin, 2014; Friedman & Selin, 2012; Thackray et al., 2015) and future climate scenarios (Friedman, Zhang, & Selin, 2014). We extend the model (Friedman & Selin, 2012) such that PAHs are now fully interactive with other atmospheric species. This allows a two‐way chemical feedback between PAHs and all other gaseous and aerosol species. Previous simulations used an offline version of the model, whereby gas and aerosol concentrations from the “full” chemistry simulation were archived and used as input for the PAH simulation.

Whereas previous modeling studies only consider between 1 and 3 PAHs, we extend the GEOS‐Chem model to provide global‐scale concentration information for a total of 48 PAH species. This consists of 16 emitted PAHs; the US EPA's list of priority PAHs. These include naphthalene (NAP), acenaphthylene (ACY), acenaphthene (ACE), fluorene (FLU), phenanthrene (PHEN), anthracene (ANT), fluoranthene (FLA), pyrene (PYR), benzo[a]anthracene (BAA), chrysene (CHR), benzo[b]fluoranthene (BBF), benzo[k]fluoranthene (BKF), benzo[a]pyrene (BAP), benzo[g,h,i]perylene (BGHIP), indeno[1,2,3‐c,d]pyrene (ICDP), and dibenz[a,h]anthracene (DAHA). The remaining 32 PAH species are the corresponding nitro‐PAHs (×16) and dinitro‐PAHs (×16).

Although commonly neglected from atmospheric modeling studies, we account for N‐PAH formation in the atmosphere by building a chemical mechanism. We build a degradation mechanism of pyrene (PYR) that accounts for the formation of nitropyrene (nitro‐PYR) and dinitropyrene (dinitro‐PYR) (Table 1). This mechanism is generalized, as it is based on the findings from laboratory studies which are not all on pyrene. We chose pyrene for three reasons. First, pyrene is the only species where N‐PAH formation yields have been determined in laboratory studies for each reaction pathway. Second, N‐PYR is included in multiple field campaigns, allowing us to evaluate our predicted concentrations for this species. Many other N‐PAHs have not been measured in the atmosphere. Third, the toxicity of N‐PYR and DN‐PYR are known, allowing us to quantify the human cancer risk of the oxidation products. The mechanism that describes the source and sinks for pyrene and its N‐PAHs is displayed in Table 1. For this species, we account for the formation of nitropyrene via (a) gas‐phase photooxidation, (b) gas‐phase direct nitration, and (c) particle‐phase direct nitration. Photolysis of the N‐PYRs is also accounted for. For the remaining 15 emitted PAHs, the chemical mechanism only accounts for the chemical removal of the parent compound, and not the production of the N‐PAHs. For these PAHs, we account for the major sinks (gas‐phase photooxidation and particle‐phase ozonolysis), without directly tracking the products of chemical mechanism online in GEOS‐Chem. This way of representing PAHs (i.e., the processes considered, and the neglect of oxidation products), is consistent with previous global and regional‐scale modeling studies (e.g., Friedman & Selin, 2012). N‐PAH concentrations from these 15 emitted PAHs are estimated in our bounding exercise and uncertainty analyses by applying the spatial distributions in the N‐PYR/PYR and DN‐PYR/PYR ratios (Figure S2) to concentrations of the remaining PAHs as a proxy for spatial patterns of their N‐PAH products.

Table 1

Reaction Kinetics for N‐PAHs

	Reaction		Second order rate coefficient	Reference
Pyrene (PYR)	Gas‐phase
	R1.	PYR(g) + NO3 → NPYR	1.6 × 10−27 × [NO2]	Atkinson et al. (1990) and Keyte et al. (2013)
	R2.	PYR(g) + OH → PYR‐OH	5.0 × 10−11	Atkinson et al. (1990)
	R3.	PYR‐OH + O2 → products	1.0 × 10–17	Koch et al. (2007) (from benzene)
	R4.	PYR‐OH + NO2 → NPYR	3.6 × 10–11	Feilberg et al. (1999) (from naphthalene)
	Particle‐phase
	R6.	PYR(p) + NO3 → NPYR	6.4 × 10–12	Liu et al. (2012)
	R7.	PYR(p) + O3 → products	4.27 × 10–17	Mean = Liu et al. (2012) and Perraudin et al. (2007)
Nitropyrene (NPYR)	Gas‐phase
	R8.	NPYR(g) + NO3 → DNPYR	1.6 × 10−27 × [NO2]	Identical to PYR
	R9.	NPYR(g) + OH → NPYR‐OH	5.0 × 10–11	Identical to PYR
	R10.	NPYR‐OH + O2 → products	1.0 × 10–17	Koch et al. (2007) (from benzene)
	R11.	NPYR‐OH + NO2 → DNPYR	3.6e × 10−11	Feilberg et al. (1999) (from naphthalene)
	Particle‐phase
	R12.	NPYR(p) + NO3 → DNPYR	1.3 × 10–12	Liu et al. (2012)
	R13.	NPYR(p) + O3 → products	2.2 × 10–17	Miet et al. (2009)
	R14.	NPYR(p) + hν → products	1.3–5.0 × 10−4
Dinitropyrene (DNPYR)	Gas‐phase
	R15.	DNPYR(g) + NO3 → products	5.0 × 10–11	Identical to PYR
	R16.	DNPYR(g) + OH → products	1.6 × 10−27 × [NO2]	Identical to PYR
	Particle‐phase
	R17.	DNPYR(p) + NO3 → products	1.3 × 10–12	Same as nitro‐PYR
	R18.	DPYR(p) + O3 → products	2.2 × 10–17	Same as nitro‐PYR
	R19.	DPYR(p) + hv → products	1.3–5.0 × 10−4	Calculated based on photolysis cross sections of 2‐methyl‐1‐nitronaphthalene

Note. Chemical and photolytic reactions of pyrene (PYR) included in the updated version of the chemical‐transport model (GEOS‐Chem).

Figure 1

Spatial distributions of annual‐average surface PAH concentrations (ng/m3) in the GEOS‐Chem model and overlaid with observed values. The left column shows three emitted PAHs: (a) benzo[a]pyrene (BAP), (b) dibenzo[a,h]anthracene (DAHA), and (c) f fluoranthene (FLA). The right column shows a parent compound and its N‐PAH degradation products: (d) pyrene (PYR), (e) nitropyrene (nitro‐PYR), and (f) dinitropyrene (dinitro‐PYR). Circles represent observed concentrations. Green circles correspond to urban environments, and blue circles correspond to non‐urban (background/remote) environments.

We provide two descriptions for heterogeneous oxidation kinetics which differ only by their reaction rate coefficients, allowing us to quantify how uncertainty in the rate of this process contributes to uncertainty in PAH distributions and human cancer risk. Within the GEOS‐Chem model, particle‐phase ozonolysis kinetics follow the Arrhenius equation using a second‐order rate coefficient (k) from Perraudin et al. (2007). Although alternative laboratory studies show that heterogeneous PAH oxidation follows a Langmuir–Hinshelwood type reaction mechanism, implying that is variable (dependent on ozone), the parameters required to account for these more realistic descriptions of heterogeneous PAH oxidation kinetics have only been developed for a limited number of PAHs, and are not implemented in this study. However, as we know this process could be much slower, we conduct a sensitivity simulation where k is reduced to 10% of its laboratory‐derived value.

We implement two widely used approaches to estimate gas‐particle partitioning, allowing us to bound the importance of uncertainties associated with this process. PAHs are semi‐volatile, meaning they can partition between gas and particle phases. In the particle phase, PAHs are observed to either be absorbed within organic aerosol (OA), or adsorbed onto the surface of black carbon (BC). We chose to implement a poly‐parameter linear free energy (Shahpoury et al., 2016) (ppLFER) scheme, but for comparison, we also conduct simulations using a single parameter scheme following Dachs and Eisenreich (Dachs & Eisenreich, 2000) (D&E).

PAH emissions for the year 2014 are from the global‐scale emission inventory developed by Shen et al. (2013). This emission inventory is used widely across global‐scale atmospheric modeling studies (Friedman, Pierce, & Selin, 2014; Friedman & Selin, 2012; Mu et al., 2018; Shrivastava et al., 2017; Thackray et al., 2015). The combined USEPA16 global‐total annual‐total emission rate is 504 Gg/a, with an interquartile range of 331–818 Gg/a (Shen et al., 2013). The sectors included in this inventory are residential and commercial, industry, transportation, deforestation and wildfire, agriculture, and energy production.

Gas‐phase PAHs undergo dry and wet deposition in a similar fashion to other gases. For all gas‐phase species, we assign a Henry's Law solubility constant of 3.1 × 10−5 m3 atm−1 mol−1, taken from Sander (2015). Particle‐phase PAHs are assumed to undergo dry and wet deposition according to the aerosol that it is bound to—that is, dry and wet deposition parameters describing the aerosol particle are used to describe deposition of the particle‐phase PAH.

GEOS‐Chem Model Simulations Performed in Study

We perform four global‐scale model simulations, which are presented in Table 2. For all simulations, we discard the first month of simulation as spin up, and base our analysis on the remaining 12 months: January 2014–December 2014. For the base simulation, we use 4° × 5° horizontal resolution, second order rate coefficients describing heterogeneous oxidation kinetics are taken directly from the laboratory studies, and gas‐particle partitioning follows the ppLFER scheme. We then perform three sensitivity simulations, where we change one parameter at a time, allowing us to isolate the importance of uncertainties in each of these processes. In the first sensitivity simulation, we reduce the heterogeneous oxidation rate coefficient to 10% of its original value (“Het_0.1”; Table 2). In the second sensitivity simulation, we change the gas‐particle partitioning scheme from the ppLFER to the D&E scheme (“D&E”; Table 2). In the third sensitivity simulation, we change the model horizontal resolution from 4° × 5° to 2° × 2.5° horizontal (“2 × 2.5”; Table 2).

Table 2

Overview of GEOS‐Chem Model Simulations Performed in This Study

Simulations	Gas‐particle partitioning	Gas‐particle partitioning scheme	Resolution	Bias‐correction
Base	Laboratory‐derived	ppLFER	4° × 5°	None
Het_0.1	10% of laboratory value	ppLFER	4° × 5°	None
D&E	Laboratory‐derived	D&E	4° × 5°	None
2 × 2.5	10% of laboratory value	ppLFER	2° × 2.5°	None
PAH_Corr	Laboratory‐derived	ppLFER	4° × 5°	Corrected to PAH concentrations from
N‐PAH_Min	Laboratory‐derived	ppLFER	4° × 5°	Corrected to minimum N‐PAH yields
N‐PAH_Max	Laboratory‐derived	ppLFER	4° × 5°	Corrected to maximum N‐PAH yields

Note. Whereas the first four simulations in table are unique model simulations, the final three simulations are based on the Base simulation, but with various bias‐correction techniques applied.

Observations Used to Evaluate Simulated PAHs and Provide Simple Bias‐Corrections

We use observations to evaluate performance of the four aforementioned GEOS‐Chem model simulations, and also to provide three additional “bias‐corrected” PAH concentration distributions (from the Base simulation). Observations for Europe and North America were taken from continuously monitoring air quality networks (US Environmental Protection Agency Air Toxics Program; European Monitoring and Evaluation Program) and those for continental Asia (Saha et al., 2017), Asia‐Arctic ship cruise transect (Ma et al., 2013), and Africa (Klánová et al., 2008) are from field campaigns. We apply three different “bias corrections” to the PAH distributions from the base simulation. In “PAH_Corr” (Table 2), we test for the effects of a bias in the differences in concentrations across the different PAHs, by multiplying the simulated PAH concentrations by the average bias between the simulated and observed mean PAH concentration across the non‐urban sites (these bias‐correction factors are shown in Table S2). In “N‐PAH_Max” and “N‐PAH_Min” (Table 2), we test for potential biases in simulated N‐PAH/PAH yields by scaling the model‐derived N‐PAH/PAH ratio by the maximum/minimum bias between the simulated and observed N‐PAH/PAH (these scaling factors are shown in Table S2). This provided upper and lower bound estimates for N‐PAH formation potential.

Calculation of Incremental Lifetime Cancer Risk

We estimate the incremental lifetime cancer risk (ILCR) using two different methods. We combine these methods with PAH distributions from the 4 model simulations and 3 bias corrected model simulations, leading to a total of 14 different estimates of ILCR.

In the epidemiologically based method, we estimate ILCR (unitless) following. where URE is the epidemiologically derived BAP unit‐risk (unit = per [ng/m3]), which is estimated at 21.3 (Moolgavkar et al., 1998), 32.7 (Armstrong et al., 2004), 85.0 (Armstrong et al., 2004) ×10−6 per (ng/m3), and [BAP] is the atmospheric concentration (ng/m3) of BAP, which is derived from the model. Under this method, overall ILCR is assumed to scale directly with BAP concentrations, and the impacts of the entire PAH mixture are accounted for (but the mixture is assumed to be fixed across the world).

We also develop an animal‐based method for estimating ILCR, which allows us to compare the human cancer risk of individual PAH species without double counting. This is calculated as following where URA is the animal‐derived BAP unit risk (unit = per [ng/m3]), which is estimated by Collins et al. (1991) at 0.37 1.0, or 1.7 × 10−6 per (ng/m3), [PAH] is the atmospheric concentration (ng/m3) of PAHs, and TEQ are their toxic equivalent quotients (unitless). Where possible, we use estimates of TEQ from the primary literature, which have not been rounded to the nearest significant figure or order of magnitude, as is the case in many literature reviews; this rounding would introduce additional uncertainty in the relative importance of different PAHs. However, where the primary literature is not available, we use the recommended values from the literature reviews. TEQ used in this study from the literature are shown in Table 3. We use the terms TEQ and relative potency factor (RPF) interchangeably. This method assumes that the cancer risk of individual PAHs combines linearly, as there is no conclusive evidence to suggest otherwise. This animal‐based method only includes the cancer risk of PAHs for which both exposure concentrations and toxicity information were available (28 of the 48 species: all 16 emitted of the emitted PAHs, 6 out of the 16 nitro‐PAHs, and 6 out of the 16 dinitro‐PAHs). When attributing ILCR to different PAHs in the animal‐based method, we also account for possible biases in the simulated distribution of PAHs (i.e., differences in concentrations among different PAH species). To account for the effect of possible biases, we used the “bias‐corrected” spatial distributions of PAHs concentrations (see Table S1 for scaling factors). We also test for any biases in our predicted N‐PAH/PAH ratios by performing sensitivity calculations, where these ratios are scaled using observed values (see Table S2 for scaling factors).

Table 3

Relative Toxicity of PAHs (TEQ)

	PAH	Nitro‐PAH	Dinitro‐PAH
NAP	0.001a
ACY	0.001a
ACE	0.001a
FLO	0.00075a
PHE	0.00075a
ANT	0.155a
FLA	0.052b	0.13b	0.13 (assumed)
PYR	0.065c	0.1c	5.1c
BAA	0.35c	0.1c	0.1 (assumed)
CHR	0.011b	10.8c	10.8 (assumed)
BBF	0.210000
BKF	0.03d
BAP	‐	0.47c	0.47 (assumed)
ICDP	0.08d
DAHA	3.0e	1.6e	1.6e
BGHIP	0.01a, f

Note. These values are used in the animal‐based method for estimating ILCR. Note, we use the terms TEQ and RPF interchangeably.

Nisbet and Lagoy (1992);

Busby et al. (1989);

Wislocki et al. (1986);

Deutschwenzel, Brune, Grimmer, Dettbarn, and Misfeld (1983);

Fu et al. (1998);

EPA (2009).

In the two equations above, the unit for the ILCR is cancer risk, and therefore dimensionless. We can combine these risks with gridded human population density and assume a human lifetime of 70 years to express the ILCR in the form of cancer rates per year. Gridded human population density is taken from the Socioeconomic Data and Applications Center (SEDAC). Both of these methods for calculating ILCR are applied identically to all regions of the world. This means that while we account for regional variability in PAH exposure, we do not account for variability in human lifetime expectancy and cancer susceptibility, which are known to depend on a wide range of biological, environmental, and socioeconomic factors. However, using this relatively simple method for calculating ILCR ensure the results of this scoping study are easy to interpret and reproducible.

Across both methods for estimating ILCR, we test minimum, median, and maximum values of BAP toxicity (UR) from the literature. We also used both the epidemiological and animal‐based methods to evaluate how uncertainties in heterogeneous oxidation kinetics affect estimates of PAH cancer risk, by applying them to two model simulations which differ in the reactivity of particle‐phase PAHs. In the discussion, we provide a more detailed evaluation of advantages and disadvantages in the epidemiological‐ and animal‐based methods for estimating ILCR.

Results

Global and Regional Concentrations of Emitted PAHs

We evaluate model performance of simulating 16 emitted PAHs by comparing simulated and observed PAH concentrations in a variety of environments. Simulated PAH concentrations were compared to a variety of different measurements and are displayed in Figure 1. In Figure 1, the simulated annual average surface concentrations for three illustrative PAHs are shown, and observed values are overlaid (circles). The left column shows three emitted PAHs, which are discussed in this section. The right column shows an emitted PAH and the N‐PAHs degradation products, which are discussed in the next section.

Global PAH simulations are poorly constrained by available data, and many uncertainties exist in their emissions and atmospheric chemistry that affect the ability to model them accurately (see Global Model Performance in SI for more detail). As noted in Section 2.6, observations for Europe and North America were taken from continuously monitoring air quality networks (US EPA and EMEP) and those for continental Asia (Saha et al., 2017), Asia‐Arctic ship cruise transect (Ma et al., 2013), and Africa (Klánová et al., 2008) are from field campaigns. Simulated and observed data at the location of measurements in Figure 1 are also represented in the form of box and whisker plots (Figure 2), where red reflects observed data, green reflects our base simulation, and blue shows a sensitivity simulation to test the influence of chemical uncertainties (described below). In Figure 2, PAH species on the x‐axis are ordered from lowest molecular weight (left) to highest molecular weight (right). A summary of statistics are also shown in the Supporting Information (Table S1).

Figure 2

Box and whisker plot of PAH concentrations for all measurement sites displayed in Figure 1. Panel a shows concentrations (ng/m3), with red representing the observed data, green representing the base model simulation, and blue representing simulated PAH concentrations under the sensitivity simulation (a 90% reduction in the second order rate coefficients describing heterogeneous oxidation). The boxes denote the 25th and 75th percentiles, the whiskers denote the 5th and 95th percentiles, the horizontal line denotes the median, and dots denote outliers. Asterisks (*) indicate where p‐value is less than 0.05. Panel b shows the ratio of simulated to observed PAH concentrations (unitless), with green representing the base, and blue representing the sensitivity simulation. The 16 directly emitted PAHs considered are naphthalene (NAP), acenaphthylene (ACY), acenaphthene (ACE), fluorene (FLO), phenanthrene (PHEN), anthracene (ANT), fluoranthene (FLA), pyrene (PYR), benzo[a]anthracene (BAA), chrysene (CHR), benzo[b]fluoranthene (BBF), benzo[k]fluoranthene (BKF), benzo[a]pyrene (BAP), benzo[g,h,i]perylene (BGHIP), indeno[1,2,3‐c,d]pyrene (ICDP), and dibenz[a,h]anthracene (DAHA). Observations are described in greater detail in Section 2, where full citations are provided.

The model captures average PAH concentrations for most PAHs, but with some low biases, especially in urban areas. For 11 out of 16 emitted species, the p‐value is less than 0.05 (indicated by * in Figure 2). Overall, simulated PAH concentrations are lower than observed (Figure 2b). For 13 of the emitted species, the model underpredicts the observed global average PAH concentration (normalized mean bias (NMB) ranges from −97% to −42%; Figure 2b). For the remaining three emitted species (ACY, BKF, and DAHA) the model overestimates the observed global‐average PAH concentration (NMB = 104%–464%; Figure 2b). For non‐urban sites (i.e., outside of cities), simulated PAH concentrations are lower than observed over the United States, and higher than observed over Europe. The bias over the United States is likely a result of our choice of emission inventory: higher‐resolution regional‐scale emission inventories predict 2–3 times higher PAH emissions over the United States (J. Zhang et al., 2017) compared to the global‐scale emission inventory used in this study (Shen et al., 2013). Similar to other global‐scale models (Friedman & Selin, 2012), our simulation underpredicts the high observed PAH concentrations typical of urban environments (Figures 1a–1c; green circles), especially in cities across Asia and Africa. In remote regions, the model reproduces heavier molecular weight species measured across a ship cruise from Beijing (CH) to the Arctic (Ma et al., 2013), but underpredicts concentrations of lighter molecular weight species. In Svalbard (Norway), however, the model captures annual average PAH concentrations for the lighter molecular weight species, but underpredicts the heavier molecular weight species by up to several orders of magnitude. Our model performance at capturing observed BAP concentrations is comparable to that of previous global‐scale modeling studies (see Discussion section in the SI—Global Model Performance).

We further evaluated model performance at capturing differences in concentrations between the different PAHs, both globally and in different regions, which are important to capture the relative cancer risk of different species (not shown). The lighter PAH species are much more abundant than the heavier species, both globally (Figure 2a) and regionally. The model generally captures these relative differences in PAH concentrations. Nevertheless, we applied a simple “bias‐correction” to the spatial distributions in PAH concentrations to test for the effects of a bias in the differences in concentrations across the different PAHs. Bias correction was conducted by multiplying the simulated PAH concentrations by the average bias between the simulated and observed mean PAH concentration across the non‐urban sites (these bias‐correction factors are shown in Table S1). These bias‐corrected PAH concentration distributions were used in the cancer risk assessment in a sensitivity calculation.

To test the ability of the model to simulate the atmospheric lifetime of different PAHs, which affects the composition of PAH mixtures when comparing source and receptor environments, we examined concentration gradients between Central Europe and the Arctic (Figure S1). The model captures the observed PAH concentration gradient between Kocetice (Czech Republic) and Svalbard (Norway) for most of the lighter weight PAH species (ACY, ACE, PHE, ANT, FLA, and PYR; Figure S1). The model underpredicts this gradient for heavier PAHs (Figure S1). While uncertainties in lifetimes and emissions combine to influence concentrations at remote sites (Thackray et al., 2015), these biases in lifetime are large enough to offset the likely overestimates in European emissions. To test the importance of these biases, we used the model sensitivity simulation that reduces the second order rate coefficient describing heterogeneous oxidation to 10% of its original value (“Het_0.1”). Lighter PAHs are insensitive to changes in heterogeneous oxidation kinetics, as these species exist mostly in the gas‐phase. The heavier molecular weight PAHs, which mostly exist in the particle‐phase, are extremely sensitive to this sensitivity simulation. Under this simulation where heterogeneous oxidation kinetics are reduced to 10% of the value in the base simulation, the simulated PAH concentration gradients between Kocetice and Svalbard agree with observed values, and biases in PAH concentrations over Svalbard are minimized (Figure S1). We use this sensitivity simulation below to test the influence of uncertainties in heterogeneous oxidation on estimates of human cancer risk (Figure 3a).

Figure 3

Global and regional impacts of ambient PAH concentrations on human cancer risk. Panel a shows a breakdown in global human cancer risk from different PAHs (%). Panel b shows global annual cancer rates (cancer cases a−1) induced by PAHs, as estimated by two different methods (epidemiologically based, and animal‐based), and under three different model simulations (Base, Het_0.1, and 2 × 2.5). Estimates from the literature are also shown. Panel c shows the spatial pattern in PAH‐induced human cancer rates under the base simulation, applying the animal‐based method to estimate ILCR (but the pattern is similar for other formulations, which are not). Note that an ILCR of 1 × 10−6, applied to the global population (∼7 × 109) is equal to 7 × 103 lifetime cancer cases, and assuming a life expectancy of 70 years, equates to an annual cancer rate of 100.

For almost all continental regions, PAH concentrations are insensitive to uncertainties in gas‐particle partitioning. The poly‐parameter linear free energy relationship (ppLFER) gas‐particle partitioning scheme used in the base model simulation captures observed particle fractions better than the D&E scheme used in the sensitivity simulation. While particle‐phase fractions differed greatly between the two gas‐particle partitioning schemes, annual‐average PAH concentrations are within 5% of each other over most continental environments under the two schemes. This is because of the very small differences between simulated gas‐ and particle‐phase lifetimes.

Global and Regional Concentrations of PAH Degradation Products

We estimated the atmospheric concentrations of N‐PAHs (nitro‐PAHs and dinitro‐PAHs) for all 16 of the emitted PAH discussed above (Section 2.4). Within the GEOS‐Chem model, the chemical mechanism for pyrene (PYR) accounts for the formation of N‐PYR (nitro‐PYR and dinitro‐PYR) (Table 1). Oxidative processes considered include gas‐phase photooxidation (+OH) and direct nitration (+NO3), and heterogeneous ozonolysis (+O3) and direct nitration (+NO3), all of which contribute to N‐PYR formation, except heterogeneous ozonolysis. For the remaining N‐PAHs, which are not incorporated into the online chemical mechanism, concentrations of the N‐PAHs were predicted offline by scaling concentrations of the remaining parent PAHs by spatially resolved ratios of nitro‐PYR/PYR and dinitro‐PYR/PYR (Figure S2). This approach recognizes that detailed atmospheric degradation data for each individual PAH are not available, and thus provides a bounding estimate of the magnitude of N‐PAH impacts under chemically relevant background conditions. We used sensitivity calculations (described below) to further test the uncertainties introduced by this approach.

Simulated global mean atmospheric concentrations of nitro‐PYR in the atmosphere are not statistically different from observations (p < 0.05, Figure 2b). Spatial patterns in the N‐PAHs (nitro‐PYR and dinitro‐PYR) and the parent PAH, PYR, are very similar (Figures 1d–1f). Across non‐urban sites, the model underpredicts the observed nitro‐PYR concentration (0.016 ng/m3) by a factor of 2.5. This bias is due to a combination of two factors. First, the model underestimates pyrene, the parent compound, by a factor of 1.4, likely due to underestimates of emissions, as described above. Second, our simulated nitro‐PYR/PYR ratio of 0.021, is lower (×0.6) than the observed value (0.036) across the aggregated data set of measurements, implying that we underestimate nitro‐PAH formation on average. However, this is not consistent across all environments.

We compared the simulated N‐PAH/PAH ratio to observed values from field campaigns, which measure the two species simultaneously (Table S2). The simulated nitro‐PYR/PYR ratio lies within the range of estimates from field campaigns; our value is ×0.2, ×0.7, and ×3 times the value observed across China, France, Hungary, respectively. These biases vary across the different PAH species, but there were no systematic patterns (Table S2). The simulated nitro‐PAH/PAH ratios range from 2.0‐ to 4.0‐fold of the observed ratios for FLA, 0.2–3.0 fold of the observed ratios for PYR, 0.06–10 fold of the observed ratios for CHR, and 0.02–5.0 fold of the observed ratios for BAP. To account for the impact of potential biases in estimating the cancer risk of PAHs and N‐PAHs, we conducted sensitivity calculations by scaling the model‐derived N‐PAH/PAH ratio by the maximum/minimum bias between the simulated and observed N‐PAH/PAH (these scaling factors are shown in Table S2). This provides upper and lower bound estimates for N‐PAH formation potential.

Heterogeneous direct nitration (+NO3) is the major source of N‐PAHs in the atmosphere, as discussed in greater detail in the SI (Global and Regional Concentrations of PAH Degradation Products). Globally, heterogeneous direct nitration accounts for 99% of nitro‐PYR production in the model. In laboratory studies, this process is a combination of multiple elementary reaction steps. Because the exact mechanism is unknown, we simulated it here using a single‐step reaction. Laboratory studies find that the yield of N‐PAHs from this process ranges 0.04% (Zelenov et al., 2018)—100% (Ringuet et al., 2012). In our model, we assumed a fixed N‐PAH formation yield, and chose 100% in order to bound this reaction pathway; despite this maximal assumption, our model still underestimates the nitro‐PYR/PYR ratio as discussed above.

The base simulation provides a better representation of the relative importance of parent PAHs and N‐PAHs compared to the sensitivity simulation in which heterogeneous oxidation kinetics are reduced to 10% of their laboratory values. PYR is insensitive to assumptions in heterogeneous oxidation kinetics, as this chemical reaction represents a minor removal term for this species compared with gas‐phase oxidation. However, heterogeneous oxidation is the major source of N‐PAHs, so as this process slows down, concentrations of N‐PAHs reduce. Reductions in the heterogeneous oxidation rate reduce the model's ability to capture the observed N‐PYR concentration and N‐PYR/PYR ratio (Figure 2b). Hence, while slower oxidation improves reaction kinetics for heavier molecular weight PAHs as discussed above, it decreases the model's ability to capture the observed concentration of N‐PYR.

Human Cancer Risk of Ambient PAH Mixtures

We calculated the incremental lifetime cancer risk (ILCR) of PAHs (Figure 3) using two different methods (Section 2.6). Briefly, the epidemiologically based method uses BAP as proxy for risk of the whole PAH mixture, and the animal‐based method accounts for regional variations in the PAH mixtures.

Using the animal‐based method allowed us to compare the relative importance of different PAHs to ILCR (Figure 3a). BAP accounts for just 11% of the calculated global human cancer risk of the entire pollutant mixture, with the remaining emitted PAHs accounting for 72%, and the 12 N‐PAHs (6 nitro‐PAHs and 6 dinitro‐PAHs) for which toxicity information is available account for the remaining 17% of global ILCR (Figure 3a). The finding that BAP was of low importance to global human cancer risk (11%) contradicts the findings of previous studies, but was robust across all sensitivity calculations conducted here (5%–19%; Figure 3a). Across our sensitivity calculations designed to capture the entire range of observed nitro‐PAH/PAH ratios of the 6 PAHs considered, the contribution of N‐PAHs to global ILCR ranges from 3% to 62% (Figure 3a). There are very few constraints on the nitro‐PAH/PAH ratio; thus, we tested the limits of these values. While these sensitivity calculations gave some indication of the uncertainty related to the importance of N‐PAHs, because each calculation used a single field study (i.e., single point) to constrain the entire global distribution of N‐PAHs, they should be considered as extreme estimates for the importance of N‐PAHs, and only applicable to environments close to where the field study is located. Despite these limitations, our model results suggest that N‐PAHs could contribute substantially to global human health, even though they are routinely neglected in risk assessments.

Under the animal‐based method, PAHs induce a global annual total human lung cancer rate of 231/year, which is ∼3 times lower than the rate estimated under the epidemiologically based method (759/year, Figure 3b). As noted above, epidemiologically based methods implicitly address the impacts of mixtures of PAHs, because they are derived from studies in which people were exposed to multiple compounds simultaneously. Our animal‐based method results in lower estimates because it includes only a small sample of all known PAHs—the 16 emitted species and the N‐PAHs for which toxicity information is available. Our animal‐based method also does not include other groups of PAHs which may be extremely toxic, such as oxygenated, halogenated, and alkylated species (Andersson & Achten, 2015), and assumes that effects of individual PAHs add linearly. Thus, the overall magnitude of ILCR associated with PAHs is likely more realistic using the epidemiologically based method. However, the epidemiologically based method does not provide insight as to the toxicity and cancer risk of different chemical species and their respective degradation products. Further, the magnitude and spatial distribution of risk predicted by the epidemiologically based method are only accurate to the extent that global PAH concentrations reflect the mixtures to which people were exposed in the original epidemiological studies, and the degree that the overall risk scales with changes in BAP concentration. Under both methods for estimating ILCR, global ILCR exceeds the commonly applied threshold level of 1 in 1,000,000 (1 × 10−6); this is comparable in magnitude to estimates from previous studies (Figure 3a). Under the epidemiologically based method for estimating human cancer risk of PAHs, 70 (63–82)% of the global population breathe air which exceeds this safe threshold level (with the ranging representing the upper and lower bounds in BAP toxicity).

Under the animal‐based method, human cancer risk is much less sensitive to uncertainties in PAH heterogeneous oxidation kinetics than that reported in previous studies using BAP as an indicator species (Shrivastava et al., 2017). This is because when particle‐phase PAHs are assumed to be less reactive, concentrations of the parent compounds increase, while concentrations of the oxidation products decrease. Previous studies, which have used the epidemiologically based method, only accounted for the former (increased human exposure to the parent compound), whereas the animal‐based method used here also accounts for the latter (decreased concentrations of the degradation products). Compared to the base simulation (Base = 0.19 ng/m3), the global‐average population‐weighted BAP concentration is 90% higher in the sensitivity simulation that tests the impact of oxidation kinetics (Het_0.1 = 0.36 ng/m3). Under the epidemiology‐based method for estimating ILCR, where the cancer risk of PAHs scales closely with BAP, global ILCR is 76% higher in the sensitivity simulation (1,335/year) compared to the standard version of the model (759/year). However, the animal‐based method shows a much weaker sensitivity in global ILCR to particle‐phase reactivity. For the same increase in BAP exposure (+90%), the global ILCR increases by only 18%, from 231 to 273/year (Figure 3b). Under this method, while slower particle‐phase reactivity increases human exposure to the parent compounds, this is partially offset by reductions in exposure to the oxidation products. The animal‐based method thus provides a more realistic estimate of the impact of heterogeneous oxidation uncertainty on cancer risks. Human cancer risk is also insensitive to uncertainties in gas‐particle partitioning. Both globally and regionally, the ILCR changes by less than 1% when the gas‐particle partitioning scheme is changed from the ppLFER scheme used in the base model simulation, to the D&E scheme used in the sensitivity simulation. Similarly, when the model resolution is increased from 4° × 5° in the base simulation to 2° × 2.5° in a sensitivity simulation (“2 × 2.5”), global ILCR increases by 29%—however, the relative importance of each PAH species to global ILCR remains unchanged.

Although omitted from previous global‐scale assessments, uncertainties in BAP toxicity also play a substantial role in influencing the magnitude in global ILCR from PAHs. In Figure 3b, for our estimates of global ILCR, the length of the bars represents the uncertainty bounds in BAP toxicity, but for Shen et al. (2014) and Shrivastava et al. (2017), they represent uncertainties associated with genetic susceptibility and heterogeneous oxidation kinetics, respectively. From Figure 3b, under the epidemiologically based method, uncertainties associated with BAP toxicity (length of red and orange bars in Figure 3b) have a larger influence on global ILCR than uncertainties in genetic susceptibility (lengths of yellow colored bar in Figure 3b) and uncertainties in heterogeneous oxidation kinetics (length of peach colored bar in Figure 3b), highlighting the importance of future research on this parameter.

Under the epidemiologically based and animal‐based methods, human cancer risk associated with PAHs is highest in urban and industrial regions, but the two methods differ in their assessments of the spatial variability of these risks. Both methods predict PAH exposure leads to the highest human cancer risk over regions such as China, India, Central, and Eastern Europe (Figure 3c). Under the epidemiologically based method, however, differences in ILCR are driven solely by BAP, whereas under the animal‐based method, the ILCR varies spatially with BAP and many other PAHs. When assessing the impact on human cancer risk of reducing BAP emissions (not shown), the epidemiologically based method would estimate a proportional reduction in the human cancer risk, but the animal‐based method would estimate a reduction in human cancer risk a factor of three or more lower (which would vary based on the regional variations in the contribution of BAP to total ILCR).

Conclusions

We developed and evaluated a new, global‐scale model that accounts for 16 emitted PAHs in addition to their degradation products. We used this model to calculate the human cancer risk of exposure to these PAH mixtures using two different methods: an epidemiologically based method based on BAP concentrations that quantifies overall risk, and an animal‐based method that allowed us to attribute risk to individual components of the pollutant mixture without double counting. We then evaluated the relative importance of BAP to global risk, and assessed the utility of using BAP as an indicator compound.

We found that BAP is only a small contributor to the global human cancer risk of PAHs (11%), suggesting it is an inadequate indicator of human cancer risk from this pollutant mixture. Atmospheric modeling studies typically only consider a single PAH species (BAP), and our work suggests that conclusions based on modeling this single compound can be misleading or erroneous. In previous studies, BAP accounted for 40%–80% of the cancer risk of PAHs (Delgado‐Saborit et al., 2011; Nielsen et al., 1996; Norramit et al., 2005; Zhang et al., 2009, 2016). However, Zhang et al. (2016) only considered eight PAHs, which partially explain their high BAP contribution to ILCR over the United States (9–154/year) of 40%–60%. In our study, where we considered 16 emitted PAHs and 12 N‐PAHs, BAP accounted for 5%–15% of ILCR over the same region. The assumed toxicity of DAHA (ranging from 1 to 10 times that of BAP), and whether it is even included, varies from study to study. Nielsen et al. (1996), who estimated that BAP accounts for 70% of cancer risk in a field study, did not include DAHA in their analysis. In addition, several highly toxic PAHs were not included in this study, due to lack of data to constrain their emissions and chemistry. Anderson and Achten (2015) argue that, in addition to the USEPA16 and N‐PAHs, scientific research should be expanded further still, to include other highly toxic parent PAHs and degradation products with an oxy group (O‐PAHs). However, atmospheric emission inventories are available only for the USEPA16, and current understanding of O‐PAHs is insufficient to build chemical mechanisms within atmospheric models. However, including further PAHs and degradation products would almost certainly further diminish the importance of BAP, strengthening our main conclusions.

In addition to the emitted PAHs, we also considered N‐PAHs, and we found them to be an important contributor to human cancer risk, but unlike BAP and the other USEPA16, they are not regulated or routinely monitored. Previous assessments of the impact of N‐PAHs were limited to a small number of field campaigns and a single box‐model study. Our model simulations showed that wherever PAHs are emitted, there is sufficient NO3 to allow the formation of N‐PAHs. Accounting for 15%–20% of the carcinogenic potential of PAH mixtures, we estimated that N‐PAHs are comparably dangerous for human cancer risk to BAP (11%). In our sensitivity calculations, the uncertainty in N‐PAH/PAH concentration ratios led to considerable variance in the contribution of N‐PAH to human cancer risk. Increased confidence in this class of chemical would be provided by (a) a deeper understanding of the formation processes (yields and mechanisms), (b) a wider understanding of the toxicity of N‐PAHs, and (c) greater geographical coverage and density of observations. Furthermore, in addition to being formed during the oxidation of parent PAHs, as simulated here, N‐PAHs can also be directly emitted into the atmosphere. We do not consider direct emissions of N‐PAHs in our atmospheric model. Our simulations, however, are constrained by observed values of N‐PAHs, thus we believe our result that N‐PAHs contribute 15%–20% to global ILCR is a robust bounding estimate. Nevertheless, providing better constraints on the source of N‐PAHs in the atmosphere should be a future research priority, especially if mitigation measures are to be considered.

PAHs and their degradation products, including N‐PAHs, can have a complex impact on human cells, altering transcriptional profiles, signaling networks, and in many cases causing DNA adducts that eventually can progress to DNA mutations. While most research has been focused on understanding the response of cells to single PAH species such as BAP, cancer risk data highlight the need to consider more complex, real‐world mixtures of PAHs in order to better define the interactions between different compounds and develop more accurate predictive models.

PAHs pose a substantial threat to global human cancer risk across all of our simulations and the two methods for estimating ILCR. We estimated an overall annual cancer risk of 231–759/year from ambient exposure to PAHs globally. Across each of the model simulations and methods in this study, as well as in previous studies, global ILCR exceeded the commonly applied threshold level of 1 in 1,000,000 (1 × 10−6) (Figure 3b). Our epidemiologically based estimate using our base simulation nevertheless calculated a lower global cancer risk than previous studies. Results from Shen et al. (2014) and Shrivastava et al. (2017) are shown for comparison in Figure 3b. Our estimates were lower for three reasons. First, Shen et al. (2014) and Shrivastava et al. (2017) choose to “downscale” their simulated BAP concentrations to reduce bias in urban environments, whereas we did not. Downscaling introduces additional uncertain parameters; we chose instead to apply the best available physical, process‐based and explore the importance of simulation biases more directly through supplementary calculations. Second, we used median estimates of the toxicity of BAP, whereas Shen et al. (2014) and Shrivastava et al. (2017) used maximum estimates from the literature. Third, to facilitate comparison of the variability driven by the concentration of different PAHs, we did not account for variability in cancer susceptibility in either of our methods for estimating ILCR, which has been shown to double global ILCR.

Our animal‐based method provides a more realistic description for spatial differences in the human cancer risk associated with PAHs, as it captures regional differences in PAH mixtures. For example, simulated annual‐average BAP concentrations were 3.5 times higher over Hong Kong compared to southern India. Using the epidemiologically based method, the calculated difference in ILCR between these two locations also differed by the same amount (×3.5), but the animal‐based method suggested that cancer risk in Hong Kong is 12 times higher than over southern India. Hong Kong had a particularly high contribution of DAHA, which the epidemiologically based method did not account for. This suggests that variations in BAP should not be viewed as indicators of variation in human cancer risk due to PAH mixtures, which is a common conclusion drawn from atmospheric models of BAP alone.

While toxicity information used in the animal‐based method may not be exactly representative of humans, it does allow us to compare the individual risk of different species in the PAH mixture. PAHs impart their toxic effects on cells through complex pathways that include responses to DNA damage and protein‐mediated cellular signaling pathways that alter gene expression of several cytochrome P450s and other enzymes. The expression of PAH‐responsive enzymes can vary widely between animals and humans, resulting in differences in susceptibility to these compounds. However, in this study we used literature values of the relative toxicities of PAHs compared to BAP, where mechanisms of action for different PAH species are often similar between animals and humans. There are no human data on the relative toxicities of these PAHs, making us reliant on toxicity data derived from animals where multiple different PAHs have been tested individually in animals. Mechanisms for differences in relative toxicities of PAHs and their nitro derivatives are not yet well understood; however, in most cases, the differences have been attributed to differences in absorption, transport and solubility of compounds in the body (Fu et al., 1990). More comprehensive understanding of how these pharmacokinetic parameters differ between animal and humans and among different PAH species could further extend the applicability of our animal‐based method.

In contrast to previous studies, we found that the cancer risks associated with PAHs are not sensitive to uncertainties in heterogeneous oxidation kinetics. As discussed above, when particle‐phase PAHs were assumed to be less reactive, concentrations of the parent compounds increase at the expense of concentrations of the oxidation products. Hence, in the animal‐based method for estimating human cancer risk, where both emitted PAHs and the degradation products are considered, global ILCR increased by 18%. Contrastingly, in the epidemiologically based method for estimating ILCR, which only considers concentrations of a single parent compound (BAP), global ILCR increased by 76%. This reduced sensitivity contradicts previous results by Shrivastava et al. (2017), who estimated, using an epidemiologically based method, that reductions in reactivity (due to a hypothesized mechanism including particle shielding) increased global ILCR by a factor of 3, corresponding with five‐fold estimates of global‐average population‐weighted BAP concentration. Hence, our holistic view of PAHs, considering both parent compounds and oxidative derivatives, weakens the sensitivity of PAH human cancer risk to uncertainties in heterogeneous oxidation kinetics. In addition, the human cancer risk associated with PAHs was insensitive to uncertainties in gas‐particle partitioning. When then the gas‐particle partitioning scheme was changed from the ppLFER to the D&E scheme, global cancer risk changed by less than 1%.

Overall, we conclude that BAP is a poor indicator of human health risks, and that other emitted PAHs and N‐PAHs are the dominant contributors to the human cancer risk of PAHs. Researchers and governing bodies should consider extending assessment and monitoring beyond BAP in order to better capture who is affected, and how the health impacts could be mitigated. Further understanding of PAHs and their impacts will require additional measurements of PAHs in the atmosphere, especially outside North America and Europe, increased understanding of the mechanisms and degradation products of PAH oxidation, better quantification of the toxicity of PAHs through toxicological and epidemiological studies, and new policy frameworks that adequately capture complex pollutant mixtures.

Conflict of Interest

The authors declare no conflicts of interest relevant to this study.

Supporting Information S1

Click here for additional data file.