# Global trade drives transboundary transfer of the health impacts of polycyclic aromatic hydrocarbon emissions – Communications Earth & Environment

### An integrated framework

In this study, an integrated approach is applied to estimate the global PAH emissions and lifetime lung cancer deaths from the perspective of the whole supply chain in the international trades (primary inputs, production, final sales, final consumption) and to determine the contributions of different drivers to the changes of two indicators. The primary input is the investment in sectors in all regions, supporting the production process. Final sale refers to the total products sold from a certain sector in a region to meet the consumption of all regions. Final consumption represents the final demand of a certain region for all sectors’ products. A detailed explanation of different accountings has been provided in Supplementary Note 3 and Supplementary Fig. 1. The overall framework is shown in Supplementary Fig. 2.

First, the global PAH emission inventories were generated in this study based on the data of emission factors, fuel consumption, and industrial processes.

Second, the environmentally extended multi-regional input-output (EE-MRIO) model was employed to reallocate the responsibilities of emissions based on each region’s primary inputs, final sales, and final consumption. Within the model, the emission inventories were linked with the MRIO table through emission intensity. Then, due to the 13 aggregated worldwide regions and four supply chain perspectives, there are 52 (4 × 13) scenario-specific emission maps for each year. The details of all scenarios are presented in Supplementary Data 5.

Third, the GEOS-Chem model was applied to derive the surface BaP concentrations of each scenario by simulating the physical transport of BaP emissions in the atmosphere. Since BaP is the most toxic pollutant among the PAHs, its simulation could represent the health impacts of PAHs37. Through inputting the scenario-specific emission obtained from the MRIO model, the scenario-specific concentrations can be simulated.

Fourth, based on the simulated concentrations and the Global Burden of Disease (GBD) database, the lifetime lung cancer risk assessment was used to evaluate the PAH-related lung cancer deaths in each scenario.

Finally, the structural decomposition analysis (SDA) was employed to determine the contributions of socioeconomic determinates to the changes in PAH emissions from 2012 to 2015. Then, the health impacts driven by socioeconomic determinates and other factors (meteorological change, lung cancer rates) were obtained by combing the results of SDA, GEOS-Chem model, and lifetime lung cancer risk assessment. The uncertainty of the integrated model has been estimated by Monte Carlo simulation (Supplementary Fig. 9)20,38,39, and the details about the uncertainty analysis and limitations are provided in Supplementary Notes 7, 8.

### PAH emission inventory

The global PAH emission inventories with different sectors in different regions from 2012 to 2015 were generated in this study. Sixteen priority PAHs were included in the inventory, which are naphthalene (NAP), acenaphthene (ACE), acenaphthylene (ACY), fluorene (FLU), anthracene (ANT), phenanthrene (PHE), fluoranthene (FLUH), pyrene (PYR), benzo[a]anthracene (BaA), chrysene (CHR), benzo[b]fluoranthene (BbF), benzo[k]fluoranthene (BkF), benzo[a]pyrene (BaP), dibenzo[a,h]anthracene (DBA), indeno[1,2,3-cd]pyrene (IDP), and benzo[g,h,i]perylene (BgP). The production-based global emission inventory of PAHs was generated by multiplying fuel consumption and activity rates with different emission factors in various sectors in different countries:

$${{{{\boldsymbol{E}}}}}_{{{{\boldsymbol{i}}}},{{{\boldsymbol{m}}}}}=\mathop{\sum}\limits_{{{{\boldsymbol{k}}}}}{{{{\boldsymbol{FC}}}}}_{{{{\boldsymbol{i}}}},{{{\boldsymbol{m}}}}}^{{{{\boldsymbol{k}}}}}\times {{{{\boldsymbol{EF}}}}}_{{{{\boldsymbol{i}}}},{{{\boldsymbol{m}}}}}^{{{{\boldsymbol{k}}}}}+\mathop{\sum}\limits_{{{{\boldsymbol{p}}}}}{{{{\boldsymbol{AR}}}}}_{{{{\boldsymbol{i}}}},{{{\boldsymbol{m}}}}}^{{{{\boldsymbol{p}}}}}\times {{{{\boldsymbol{EF}}}}}_{{{{\boldsymbol{i}}}},{{{\boldsymbol{m}}}}}^{{{{\boldsymbol{p}}}}}$$

(1)

where Ei,m indicates PAH emissions released from consumed energy type k (including coal, natural gas, petroleum, and biomass) and industrial processes in sector i in country m, FC indicates the fuel consumption of different energy types, AR indicates the activity rate, representing produced products in different industrial process p, EF indicates the emission factor.

The global fuel consumption was obtained from the International Energy Agency (IEA), which includes 69 energy types in 103 flows (28 sectors) in 186 worldwide regions. The industrial processes contain primary aluminum production, pig iron production, crude steel production, coke production, petroleum catalytic cracking, and agricultural residue burning. The activity rates of these processes in all countries were derived from different data sources (Supplementary Data 1). The emission factors of PAHs for 69 sources were obtained from literature13. In addition, emission factors varied between countries due to the difference in technologies. Through the technology splitting approach, the source-specific emission factors were distinguished based on different emission mitigation measures in countries. The details about emission factors and the comparison with previous PAH inventories are presented in Supplementary Notes 1, 2 (Supplementary Data 2).

### Income, final sale, and consumption-based accounting

The EE-MRIO model has been widely applied in analyzing the influences of increasing international trade on pollution21,40,41. In this study, after obtaining the production-based emissions, the EE-MRIO model was employed to estimate the region-specific income-, final sale-, and consumption-based PAH emissions. The different emission accounting represents PAH pollutants produced throughout the supply chain caused by certain primary inputs, final sales, and final consumption. The EE-MRIO model can track PAH emissions from the region of primary input, final sales, and final consumption to the region of production through international trade. In the EE-MRIO framework, based on the detailed trade data about the import and export among all sectors of regions, all produced PAH emissions can be assigned to the primary suppliers, final sellers, and final consumers. In other words, a certain region can take the responsibility for PAH emissions emitted from the production of the goods and services due to this region’s primary input, final sale, and final consumption.

The latest MRIO tables for 2012, 2013, 2014, and 2015 were obtained from Eora (www.worldmrio.com, accessed Feb. 2021)42, which contains 190 regions and 26 sectors (Supplementary Data 6). Due to their time continuity, large coverage of countries and sectors, and sufficient data of final demand and primary input, the Eora MRIO tables were chosen. To make the temporal results comparable, the tables were converted into 2015 prices with a double deflation method based on the GDP deflator of different regions obtained from the World Bank (www.data.worldbank.org, assessed Feb. 2021)43. In addition, the MRIO tables and PAH emissions have been aggregated into 17 sectors to match each other (Supplementary Table 5). Furthermore, since deforestation and wildfire are not related to economic activities, the PAH emissions from natural sources were excluded from the MRIO model.

Assuming that there are m regions with i sectors in each region (superscripts and subscripts m and n are region identifiers; i and j are sector identifiers), the basic linear equation of the input-output model is:

$${{{\boldsymbol{X}}}}={({{{\boldsymbol{I}}}}-{{{\boldsymbol{A}}}})}^{-{{{\bf{1}}}}}{{{\boldsymbol{F}}}}$$

(2)

X represents the total output, whose element $${x}_{i}^{m}$$ indicates the sum of intermediate use in all sectors and the final demand of n regions for sector i in region m. I represents the identity matrix. The direct input coefficient matrix A can be expressed as:

$${{{\boldsymbol{A}}}}=\left[\begin{array}{c}\begin{array}{c}{{{{\boldsymbol{A}}}}}^{{{{\bf{11}}}}}\\ {{{{\boldsymbol{A}}}}}^{{{{\bf{21}}}}}\end{array}\\ \vdots \\ {{{{\boldsymbol{A}}}}}^{{{{\boldsymbol{m}}}}{{{\bf{1}}}}}\end{array}\begin{array}{ccc}\begin{array}{c}{{{{\boldsymbol{A}}}}}^{{{{\bf{12}}}}}\\ {{{{\boldsymbol{A}}}}}^{{{{\bf{22}}}}}\end{array} & \begin{array}{c}\cdots \\ \cdots \end{array} & \begin{array}{c}{{{{\boldsymbol{A}}}}}^{{{{\bf{1}}}}{{{\boldsymbol{n}}}}}\\ {{{{\boldsymbol{A}}}}}^{{{{\bf{2}}}}{{{\boldsymbol{n}}}}}\end{array}\\ \vdots & \ddots & \vdots \\ {{{{\boldsymbol{A}}}}}^{{{{\boldsymbol{m}}}}{{{\bf{2}}}}} & \cdots & {{{{\boldsymbol{A}}}}}^{{{{\boldsymbol{mn}}}}}\end{array}\right]$$

(4)

$${A}^{{mn}}=\left[\begin{array}{c}\begin{array}{c}{a}_{11}^{{mn}}\\ {a}_{21}^{{mn}}\end{array}\\ \vdots \\ {a}_{i1}^{{mn}}\end{array}\begin{array}{ccc}\begin{array}{c}{a}_{12}^{{mn}}\\ {a}_{22}^{{mn}}\end{array} & \begin{array}{c}\cdots \\ \cdots \end{array} & \begin{array}{c}{a}_{1j}^{{mn}}\\ {a}_{2j}^{{mn}}\end{array}\\ \vdots & \ddots & \vdots \\ {a}_{i2}^{{mn}} & \cdots & {a}_{{ij}}^{{mn}}\end{array}\right]$$

(5)

$${a}_{{ij}}^{{mn}}=\frac{{z}_{{ij}}^{{mn}}}{{x}_{j}^{n}}$$

(6)

The direct output coefficient matrix B can be expressed as:

$${{{\boldsymbol{B}}}}=\left[\begin{array}{c}\begin{array}{c}{{{{\boldsymbol{B}}}}}^{{{{\bf{11}}}}}\\ {{{{\boldsymbol{B}}}}}^{{{{\bf{21}}}}}\end{array}\\ \vdots \\ {{{{\boldsymbol{B}}}}}^{{{{\boldsymbol{m}}}}{{{\bf{1}}}}}\end{array}\begin{array}{ccc}\begin{array}{c}{{{{\boldsymbol{B}}}}}^{{{{\bf{12}}}}}\\ {{{{\boldsymbol{B}}}}}^{{{{\bf{22}}}}}\end{array} & \begin{array}{c}\cdots \\ \cdots \end{array} & \begin{array}{c}{{{{\boldsymbol{B}}}}}^{{{{\bf{1}}}}{{{\boldsymbol{m}}}}}\\ {{{{\boldsymbol{B}}}}}^{{{{\bf{2}}}}{{{\boldsymbol{m}}}}}\end{array}\\ \vdots & \ddots & \vdots \\ {{{{\boldsymbol{B}}}}}^{{{{\boldsymbol{m}}}}{{{\bf{2}}}}} & \cdots & {{{{\boldsymbol{B}}}}}^{{{{\boldsymbol{mm}}}}}\end{array}\right]$$

(7)

$${B}^{{mn}}=\left[\begin{array}{c}\begin{array}{c}{b}_{11}^{{mn}}\\ {b}_{21}^{{mn}}\end{array}\\ \vdots \\ {b}_{i1}^{{mn}}\end{array}\begin{array}{ccc}\begin{array}{c}{b}_{12}^{{mn}}\\ {b}_{22}^{{mn}}\end{array} & \begin{array}{c}\cdots \\ \cdots \end{array} & \begin{array}{c}{b}_{1j}^{{mn}}\\ {b}_{2j}^{{mn}}\end{array}\\ \vdots & \ddots & \vdots \\ {b}_{i2}^{{mn}} & \cdots & {b}_{{ij}}^{{mn}}\end{array}\right]$$

(8)

$${b}_{{ij}}^{{mn}}=\frac{{z}_{{ij}}^{{mn}}}{{x}_{i}^{m}}$$

(9)

where $${a}_{{ij}}^{{mn}}$$ indicates the direct consumption coefficient representing the direct intermediate use from sector i in region m to generate one monetary unit of production of sector j in the region n; $${b}_{{ij}}^{{mn}}$$ reflects the direct intermediate use from sector j in region n to generate one monetary unit of production of sector i in the region m; $${z}_{{ij}}^{{mn}}$$ indicates the intermediate monetary flows from sector i in region m to sector j in region n. (I-A)−1 indicates the Leontief inverse matrix, and (I-B)−1 indicates the Ghosh inverse matrix.

$${{{\boldsymbol{F}}}}=\left[\begin{array}{c}\begin{array}{c}{{{{\boldsymbol{f}}}}}_{{{{\bf{1}}}}}^{{{{\bf{11}}}}}\\ {{{{\boldsymbol{f}}}}}_{{{{\bf{2}}}}}^{{{{\bf{11}}}}}\end{array}\\ \begin{array}{c}\vdots \\ {{{{\boldsymbol{f}}}}}_{{{{\boldsymbol{i}}}}}^{{{{\bf{11}}}}}\\ \vdots \end{array}\\ {{{{\boldsymbol{f}}}}}_{{{{\boldsymbol{i}}}}}^{{{{\boldsymbol{m}}}}{{{\bf{1}}}}}\end{array}\begin{array}{ccc}\begin{array}{c}{{{{\boldsymbol{f}}}}}_{{{{\bf{1}}}}}^{{{{\bf{12}}}}}\\ {{{{\boldsymbol{f}}}}}_{{{{\bf{2}}}}}^{{{{\bf{12}}}}}\end{array} & \begin{array}{c}\cdots \\ \cdots \end{array} & \begin{array}{c}{{{{\boldsymbol{f}}}}}_{{{{\bf{1}}}}}^{{{{\bf{1}}}}{{{\boldsymbol{n}}}}}\\ {{{{\boldsymbol{f}}}}}_{{{{\bf{2}}}}}^{{{{\bf{1}}}}{{{\boldsymbol{n}}}}}\end{array}\\ \begin{array}{c}\vdots \\ {{{{\boldsymbol{f}}}}}_{{{{\boldsymbol{i}}}}}^{{{{\bf{12}}}}}\\ \vdots \end{array} & \ddots & \begin{array}{c}\vdots \\ {{{{\boldsymbol{f}}}}}_{{{{\boldsymbol{i}}}}}^{{{{\bf{1}}}}{{{\boldsymbol{n}}}}}\\ \vdots \end{array}\\ {{{{\boldsymbol{f}}}}}_{{{{\boldsymbol{i}}}}}^{{{{\boldsymbol{mn}}}}} & \cdots & {{{{\boldsymbol{f}}}}}_{{{{\boldsymbol{i}}}}}^{{{{\boldsymbol{mn}}}}}\end{array}\right]$$

(10)

$${{{\rm{y}}}}=\left[\begin{array}{c}\begin{array}{c}\mathop{\sum}\limits_{n}{f}_{1}^{1n}\\ \vdots \end{array}\\ \mathop{\sum}\limits_{n}{f}_{i}^{1n}\\ \begin{array}{c}\vdots \\ \mathop{\sum}\limits_{n}{f}_{i}^{{mn}}\end{array}\end{array}\right]$$

(11)

$${{{\rm{v}}}}=\left[\begin{array}{ccc}{v}_{1}^{1} & \cdots & \begin{array}{ccc}{v}_{j}^{1} & \cdots & {v}_{j}^{n}\end{array}\end{array}\right]$$

(12)

where F indicates the matrix of final demand, whose element $${f}_{i}^{{mn}}$$ represents the total final demand in region n for the products of sector i in region m; y indicates the vector of final sales, whose element $$\mathop{\sum}\limits_{n}{f}_{i}^{{mn}}$$ refers to the sum of total products sold from sector i in region m to all destination regions; v indicates the vector of primary input, whose element $${v}_{j}^{n}$$ represents the primary input of sector j in region n.

To link PAH emissions with the monetary flow, the emission intensity u, indicating the ratio of the total emissions of sector i to the total output of sector i, can be expressed as follows:

$${{{\boldsymbol{u}}}}=\left[\begin{array}{c}\begin{array}{c}{{{{\boldsymbol{u}}}}}^{{{{\bf{1}}}}}\\ {{{{\boldsymbol{u}}}}}^{{{{\bf{2}}}}}\end{array}\\ \vdots \\ {{{{\boldsymbol{u}}}}}^{{{{\boldsymbol{m}}}}}\end{array}\right],\,{{{{\boldsymbol{u}}}}}^{{{{\boldsymbol{m}}}}}=\left[\begin{array}{c}\begin{array}{c}{{{{\boldsymbol{u}}}}}_{{{{\bf{1}}}}}^{{{{\boldsymbol{m}}}}}\\ {{{{\boldsymbol{u}}}}}_{{{{\bf{2}}}}}^{{{{\boldsymbol{m}}}}}\end{array}\\ \vdots \\ {{{{\boldsymbol{u}}}}}_{{{{\boldsymbol{i}}}}}^{{{{\boldsymbol{m}}}}}\end{array}\right],{{{{\boldsymbol{u}}}}}_{{{{\boldsymbol{i}}}}}^{{{{\boldsymbol{m}}}}}=\frac{{{{{\boldsymbol{E}}}}}_{{{{\boldsymbol{i}}}}}^{{{{\boldsymbol{m}}}}}}{{{{{\boldsymbol{X}}}}}_{{{{\boldsymbol{i}}}}}^{{{{\boldsymbol{m}}}}}}$$

(13)

where u represents a vector of the emission intensity for all sectors in all regions. Then, the PAH emissions attributed to primary input, final sale, and consumption can be expressed mathematically as:

$${{{{\boldsymbol{Q}}}}}_{{{{\boldsymbol{I}}}}}={{{\boldsymbol{V}}}}{({{{\boldsymbol{I}}}}-{{{\boldsymbol{B}}}})}^{-{{{\bf{1}}}}}{{{\boldsymbol{U}}}}$$

(14)

$${Q}_{S}=U{(I-A)}^{-1}Y$$

(15)

$${Q}_{C}=U{(I-A)}^{-1}F$$

(16)

QI, QS, QC indicate the total income-based, final sale-based, and consumption-based emission matrixes, whose element Qmn represents the emissions released in region m for the primary input, final sale, and final consumption in region n. U, V, and Y are the diagonal matrix of u, v, and y.

Geographical aggregation was required due to the computational constraints and consumed time of each GEOS-Chem model run. According to the roles in international trade, pollution levels, and levels of economic development, the world was classified into 13 regions which is consistent with the previous research3: India, China, rest of east Asia, rest of Asia, Canada, USA, Latin America, eastern Europe, western Europe, Russia, Middle East and north Africa, Sub-Saharan Africa, and rest of the world (Supplementary Fig. 10). The result matrixes of QI, QS, QC were aggregated into these 13 regions so that the income-, production-, final sale-, and consumption-based emissions of 13 regions can be obtained.

### GEOS-Chem chemical transport model

GEOS-Chem is a global 3-D model of atmospheric chemistry driven by assimilated meteorological observations from the Goddard Earth Observing System (GEOS) of the NASA Global Modeling Assimilation Office (GMAO)44. This study employed the GEOS-Chem model (version 13.0.2) to simulate the global near-surface BaP concentrations of different emission scenarios.

The simulation for the persistent organic pollutants was chosen in the model run. This simulation can generate the concentrations of PAH in the gas phase, PAH partitioning in/onto organic carbon and black carbon aerosols through introducing temperature-dependent gas-particle partitioning14. Furthermore, the PAH loss through oxidation reaction with hydroxyl radical and O3 was considered in the simulation by using an empirically derived rate constant. In addition, wet depositions including rainout and washout from scavenging in convective updraft and precipitation for both gas and particulate PAH were compatible in the GEOS-Chem model by applying the temperature-dependent air-water partition coefficient and scavenging ratio. Dry deposition of PAH was included by calculating the dry deposition velocities with a resistance-in-series scheme. Moreover, the simulation was updated by calculating the soil-air and vegetation-air exchanges with a level-III fugacity model29. The incorporation of re-emissions of PAH increased the accuracy of the model results.

Before the running of the GEOS-Chem model, the input files need to be prepared. The income-, production-, final sale-, and consumption-based BaP emissions of 13 regions have been obtained through the construction of the PAH emission inventories and EE-MRIO model. The PAH emission distribution was based on 0.1° × 0.1°-resolution PKU-PAH inventory (www.inventory.pku.edu.cn, accessed Feb. 2021)13. Supposing that the spatial distribution patterns of the emissions did not change greatly, the BaP emission maps with different scenarios from 2012 to 2015 in this study were generated based on calculated country-specific income-based, production-based, final sale-based, and consumption-based emissions and PKU-PAH emission distributions.

Due to the high performance and convenience, the Amazon Web Services cloud computing platform was employed to run the GEOS-Chem model. In all runs, the simulations were conducted for the entire year with a 12-month spin-up in a 1-min time-frequency. The results of the simulation are 12 monthly average concentrations with 4° × 5° resolution and 47 vertical levels, which were further remapped into annual average concentrations with 1° × 1° resolution. The sum of BaP concentrations of the gas phase, particulate in/onto organic carbon, and black carbon aerosols at the bottom level of the model was considered as the ground-level concentrations (Supplementary Fig. 11).

The model validation was conducted by comparing the simulated BaP concentrations with the ground measurements from European Monitoring and Evaluation Programme (EMEP) (www.projects.nilu.no/ccc/emepdata, assessed Mar. 2021), and the National Air Pollution Surveillance Network of Canada (NAPS) (www.data.ec.gc.ca, assessed Mar. 2021), and observations in the world from literature (Supplementary Data 7). The locations of the observations from EMEP, NAPS, and previous studies are shown in Supplementary Fig. 12. The comparison between the simulated values and observation is presented in Supplementary Fig. 13. In general, the model-simulated concentrations agree well with the ground observations (R2 = 0.84), representing that the simulations are dependable.

### Lifetime lung cancer assessment

The lung cancer risk assessment based on the dose-response concept has been widely applied to estimate the lung cancer risks associated with the general exposure to ambient atmospheric PAH within a lifetime9,45,46. The population attributable fraction (PAF) is the percentage of lung cancer deaths caused by exposure to certain BaP concentrations in a lifetime. For each grid cell, PAF can be calculated as follows:

$${{{\bf{rr}}}}\left({{{{\bf{c}}}}}_{{{{\bf{BaP}}}}}\right)={\left[{{{{\bf{URR}}}}}_{{{{\bf{cum}}}}.{{{\mathbf{exp }}}}={{{\bf{100}}}}}\right]}^{({{{{\bf{c}}}}}_{{{{\bf{BaP}}}}}\times \frac{{{{\bf{70}}}}}{{{{\bf{100}}}}})}$$

(17)

$${{{\rm{PAF}}}}=\frac{{{{\rm{rr}}}}\left({{{{\rm{c}}}}}_{{{{\rm{BaP}}}}}\right)-1}{{{{\rm{rr}}}}\left({{{{\rm{c}}}}}_{{{{\rm{BaP}}}}}\right)}$$

(18)

where cBaP indicates the near-surface BaP concentrations (μg∙m−3), URR is the unit relative risk and rr(cBaP) is the relative risk associated with a certain BaP concentration. URRcum.exp=100 refers to the unit relative risk at 100 μg∙m−3 years BaP, which is different between regions (Asia: 1.3, Europe: 1.13, north America: 1.16, the average value of the world: 1.2)9.

In this study, a lifetime exposure of 70 years was assumed, and the grids were at 1° × 1° resolution. The lifetime lung cancer deaths due to exposure to BaP in each grid cell can be calculated as:

$${{{{\bf{D}}}}}_{{{{\bf{BaP}}}}}={{{\bf{PAF}}}}\times \mathop{\sum }\limits_{{{{\bf{y}}}}={{{\bf{1}}}}}^{{{{\bf{70}}}}}{{{{\bf{LD}}}}}_{{{{\bf{y}}}},{{{\bf{m}}}}}$$

(19)

$${{{{\rm{LD}}}}}_{{{{\rm{y}}}},{{{\rm{m}}}}}={{{{\rm{LR}}}}}_{{{{\rm{y}}}},{{{\rm{m}}}}}\times {{{{\rm{P}}}}}_{{{{\rm{y}}}},{{{\rm{m}}}}}$$

(20)

$${{{{\rm{D}}}}}_{{{{\rm{BaP}}}},{{{\rm{m}}}}}={\sum }_{g}^{m}{{{{\rm{D}}}}}_{{{{\rm{BaP}}}}}$$

(21)

where DBaP represents the BaP-related lung cancer deaths in the next 70 years in each grid cell, assuming people are exposed to certain BaP concentrations for the whole lifetime; DBaP, m is the sum of DBaP in all grid cells of the region m; $$\mathop{\sum }\nolimits_{{{{\rm{y}}}}=1}^{70}{{{{\rm{LD}}}}}_{{{{\rm{y}}}},{{{\rm{m}}}}}$$ refers to the sum of lung cancer deaths in countries m in the lifetime (70 years); LRy,m and Py,m are the lung cancer rate (deaths per 100,000 people) and population of certain country m in a certain year y. The lung cancer rates of countries were obtained from the Global Burden of Disease (GBD) database (www.healthdata.org, assessed Mar. 2021). The population for the next 70 years was obtained from World Population Prospects 2019 (www.population.un.org, assessed Mar. 2021)47. The population distributions for different years were obtained from the LandScan global population database (www.landscan.ornl.gov, assessed Mar. 2021) at a 1 km resolution, which was further aggregated to the same resolution of the BaP concentration (1° × 1°). The lifetime lung cancer deaths for base cases in each year can be estimated based on the above method, and those for each scenario can be calculated as follows. For example, a scenario is about the emissions due to the consumption of China in 2015. Firstly, the emissions caused by the consumption of China in 2015 were imported into the GEOS-Chem model to obtain the concentrations of this scenario. Then, the deaths attributed to the Chinese consumption-based emissions in each grid can be calculated by multiplying the base deaths and the fractional contribution of scenario concentration to base concentration3. Since the results cover all the grids, the deaths in each region caused by the Chinese consumption-based emissions can be obtained by summing the values in the region’s grid cells (Eq. 21). The sum of the scenario’s deaths in all grids is the Chinese consumption-based lifetime lung cancer deaths. Through the lifetime lung cancer assessment, the lifetime lung cancer deaths due to each region’s primary input, production, final sales, and final consumption can be estimated. The emissions and lifetime lung cancer deaths that occurred in each region from 2012 to 2015 are provided in Supplementary Table 6.

### Structural deposition analysis (SDA) for emissions

Input-output SDA has been widely employed for quantifying the influence of socioeconomic drivers on pollutant emissions48,49. To estimate the influence of drivers on the changes in PAH emissions in the world from 2012 to 2015, the income-, final sale-, and consumption-based emissions can be decomposed into six factors, as follows:

$${{{{\boldsymbol{Q}}}}}_{{{{\boldsymbol{I}}}}}={{{\boldsymbol{V}}}}\cdot {{{\boldsymbol{G}}}}\cdot {{{\boldsymbol{U}}}}={{{\boldsymbol{P}}}}\cdot {{{\boldsymbol{W}}}}\cdot {{{\boldsymbol{N}}}}\cdot {{{\boldsymbol{G}}}}\cdot {{{\boldsymbol{T}}}}\cdot {{{\boldsymbol{K}}}}$$

(22)

$${Q}_{S}=U \cdot L \cdot Y=K \cdot T \cdot L \cdot Z \cdot S \cdot P$$

(23)

$${Q}_{C}=U \cdot L \cdot F=K \cdot T \cdot L \cdot M \cdot R \cdot P$$

(24)

where L represents the Leontief inverse matrix indicating the production input structure, G represents the Ghosh inverse matrix indicating the production output structure, and U is the diagonal emission intensity matrix, which can be decomposed into two factors: K indicates the emission factors (PAH emissions per unit of energy consumption) and T indicates the energy efficiency (energy consumption per unit of total output and fuel mix of sectors); and V, Y, F are the primary input, final sale and final demand matrixes, which can be further decomposed into three explicit determinants: N, Z, and M are the primary input, final sale and final consumption structures (the share of each variable of total values); W, S, and R are the primary input, final sale and final consumption levels (variable per capita); and P is the population. Each of the six factors denotes the contributions to PAH emission change triggered by one driving factor when the other factors are kept constant. Theoretically, in this six-factor SDA model, the non-uniqueness problem will result in 6! types of first-order decomposition forms, and different results can be obtained from different procedures. Following previous studies50,51, two-polar decomposition using the arithmetic average of all possible first-order decomposition forms was applied to address this problem. Thus, the changes in consumption-based PAH emissions across different years can be expressed as follows:

$$\triangle {{{{\boldsymbol{Q}}}}}_{{{{\boldsymbol{C}}}}}= \,\triangle {{{{\boldsymbol{Q}}}}}_{{{{\boldsymbol{K}}}}}+\triangle {{{{\boldsymbol{Q}}}}}_{{{{\boldsymbol{T}}}}}+\triangle {{{{\boldsymbol{Q}}}}}_{{{{\boldsymbol{L}}}}}+\triangle {{{{\boldsymbol{Q}}}}}_{{{{\boldsymbol{M}}}}}+\triangle {{{{\boldsymbol{Q}}}}}_{{{{\boldsymbol{R}}}}}+\triangle {{{{\boldsymbol{Q}}}}}_{{{{\boldsymbol{P}}}}}\\ = \frac{\triangle {{{\boldsymbol{K}}}}{{{{{\boldsymbol{T}}}}}_{{{{\boldsymbol{t}}}}}{{{\boldsymbol{L}}}}}_{{{{\boldsymbol{t}}}}}{{{{\boldsymbol{M}}}}}_{{{{\boldsymbol{t}}}}}{{{{\boldsymbol{R}}}}}_{{{{\boldsymbol{t}}}}}{{{{\boldsymbol{P}}}}}_{{{{\boldsymbol{t}}}}}+\triangle {{{\boldsymbol{K}}}}{{{{{\boldsymbol{T}}}}}_{{{{\bf{0}}}}}{{{\boldsymbol{L}}}}}_{{{{\bf{0}}}}}{{{{\boldsymbol{M}}}}}_{{{{\bf{0}}}}}{{{{\boldsymbol{R}}}}}_{{{{\bf{0}}}}}{{{{\boldsymbol{P}}}}}_{{{{\bf{0}}}}}}{{{{\bf{2}}}}}\\ +\frac{{{{{\boldsymbol{K}}}}}_{{{{\bf{0}}}}}\triangle {{{\boldsymbol{T}}}}{{{{\boldsymbol{L}}}}}_{{{{\boldsymbol{t}}}}}{{{{\boldsymbol{M}}}}}_{{{{\boldsymbol{t}}}}}{{{{\boldsymbol{R}}}}}_{{{{\boldsymbol{t}}}}}{{{{\boldsymbol{P}}}}}_{{{{\boldsymbol{t}}}}}+{{{{\boldsymbol{K}}}}}_{{{{\boldsymbol{t}}}}}\triangle {{{\boldsymbol{T}}}}{{{{\boldsymbol{L}}}}}_{{{{\bf{0}}}}}{{{{\boldsymbol{M}}}}}_{{{{\bf{0}}}}}{{{{\boldsymbol{R}}}}}_{{{{\bf{0}}}}}{{{{\boldsymbol{P}}}}}_{{{{\bf{0}}}}}}{{{{\bf{2}}}}}\\ +\frac{{{{{\boldsymbol{K}}}}}_{{{{\bf{0}}}}}{{{{\boldsymbol{T}}}}}_{{{{\bf{0}}}}}{\triangle {{{\boldsymbol{LM}}}}}_{{{{\boldsymbol{t}}}}}{{{{\boldsymbol{R}}}}}_{{{{\boldsymbol{t}}}}}{{{{\boldsymbol{P}}}}}_{{{{\boldsymbol{t}}}}}+{{{{\boldsymbol{K}}}}}_{{{{\boldsymbol{t}}}}}{{{{\boldsymbol{T}}}}}_{{{{\boldsymbol{t}}}}}\triangle {{{\boldsymbol{L}}}}{{{{\boldsymbol{M}}}}}_{{{{\bf{0}}}}}{{{{\boldsymbol{R}}}}}_{{{{\bf{0}}}}}{{{{\boldsymbol{P}}}}}_{{{{\bf{0}}}}}}{{{{\bf{2}}}}}\\ +\frac{{{{{\boldsymbol{K}}}}}_{{{{\bf{0}}}}}{{{{\boldsymbol{T}}}}}_{{{{\bf{0}}}}}{{{{\boldsymbol{L}}}}}_{{{{\bf{0}}}}}{\triangle {{{\boldsymbol{MR}}}}}_{{{{\boldsymbol{t}}}}}{{{{\boldsymbol{P}}}}}_{{{{\boldsymbol{t}}}}}+{{{{\boldsymbol{K}}}}}_{{{{\boldsymbol{t}}}}}{{{{\boldsymbol{T}}}}}_{{{{\boldsymbol{t}}}}}{{{{\boldsymbol{L}}}}}_{{{{\boldsymbol{t}}}}}{\triangle {{{\boldsymbol{MR}}}}}_{{{{\bf{0}}}}}{{{{\boldsymbol{P}}}}}_{{{{\bf{0}}}}}}{{{{\bf{2}}}}}\\ +\frac{{{{{\boldsymbol{K}}}}}_{{{{\bf{0}}}}}{{{{\boldsymbol{T}}}}}_{{{{\bf{0}}}}}{{{{\boldsymbol{L}}}}}_{{{{\bf{0}}}}}{{{{\boldsymbol{M}}}}}_{{{{\bf{0}}}}}\triangle {{{\boldsymbol{R}}}}{{{{\boldsymbol{P}}}}}_{{{{\boldsymbol{t}}}}}+{{{{\boldsymbol{K}}}}}_{{{{\boldsymbol{t}}}}}{{{{\boldsymbol{T}}}}}_{{{{\boldsymbol{t}}}}}{{{{\boldsymbol{L}}}}}_{{{{\boldsymbol{t}}}}}{{{{\boldsymbol{M}}}}}_{{{{\boldsymbol{t}}}}}\triangle {{{\boldsymbol{R}}}}{{{{\boldsymbol{P}}}}}_{{{{\bf{0}}}}}}{{{{\bf{2}}}}}\\ +\frac{{{{{\boldsymbol{K}}}}}_{{{{\bf{0}}}}}{{{{\boldsymbol{T}}}}}_{{{{\bf{0}}}}}{{{{\boldsymbol{L}}}}}_{{{{\bf{0}}}}}{{{{\boldsymbol{M}}}}}_{{{{\bf{0}}}}}{{{{\boldsymbol{R}}}}}_{{{{\bf{0}}}}}\triangle {{{\boldsymbol{P}}}}+{{{{\boldsymbol{K}}}}}_{{{{\boldsymbol{t}}}}}{{{{\boldsymbol{T}}}}}_{{{{\boldsymbol{t}}}}}{{{{\boldsymbol{L}}}}}_{{{{\boldsymbol{t}}}}}{{{{\boldsymbol{M}}}}}_{{{{\boldsymbol{t}}}}}{{{{\boldsymbol{R}}}}}_{{{{\boldsymbol{t}}}}}\triangle {{{\boldsymbol{P}}}}}{{{{\bf{2}}}}}$$

(25)

where ∆ represents the change in a factor, and the subscripts t and 0 indicate two specific years. The changes in income-based and final sale-based emissions can be decomposed into six socioeconomic drivers in a similar method.

### Drivers in changes of health impacts

In this study, the contributions of determinates to the health risks and lung cancer deaths were determined. Since lifetime lung cancer deaths were obtained through the GEOS-Chem model and lifetime lung cancer risk assessment, not only the socio-economic drivers but also the meteorological change and country-specific lung cancer death between different years can affect the change of PAH-related environmental health risks. The changes in deaths from year 0 to year t (∆D0~t) can be decomposed as:

$$\triangle {{{{\boldsymbol{D}}}}}_{{{{\bf{0}}}} \sim {{{\boldsymbol{t}}}}}=\triangle {{{{\boldsymbol{D}}}}}_{{{{\boldsymbol{LD}}}}}+\triangle {{{{\boldsymbol{D}}}}}_{{{{\boldsymbol{MC}}}}}+\mathop{\sum}\limits_{{{{\boldsymbol{\theta }}}}}\left({\triangle {{{\boldsymbol{D}}}}}_{{{{{\boldsymbol{Q}}}}}_{{{{\boldsymbol{\theta }}}}}}\right)$$

(26)

$${\triangle D}_{{Q}_{\theta }}=\frac{1}{2}\left[\left({D}_{{0,Q}_{0}+\triangle {Q}_{\theta }}-{D}_{{0,Q}_{0}}\right)+\left({D}_{t,{Q}_{t}}-{D}_{{t,Q}_{t}-\triangle {Q}_{\theta }}\right)\right]$$

(27)

where ∆DQθ is the change in deaths caused by the socioeconomic driver of θ from year 0 to year t; $${D}_{{0,Q}_{0}}$$ is the lifetime lung cancer deaths determined by original emissions (Q0) simulated in the previous year of 0; $${D}_{{0,{Q}}_{0}+\triangle {Q}_{\theta }}$$ refers to the deaths determined with the combination of original emissions (Q0) and the emissions changed by the socioeconomic driver of θ ($$\triangle {Q}_{\theta }$$) simulated in the year 0; $${D}_{{t,{Q}}_{t}-\triangle {Q}_{\theta }}$$ is the deaths determined with the subtraction of the emissions ($${Q}_{t}$$) and the emissions changed by the socioeconomic driver of θ (∆Qθ) simulated in the year of t. To increase the accuracy of the results, forward and backward assumptions were employed through controlling variates. The forward assumption was that only the driver i affected the emissions, and the backward assumption was that all the drivers affected the emissions except the driver θ. In other words, for each assumption, two simulations were conducted in the same situation (year of 0 or t) with different emission maps, and the difference in maps was the addition or removal of emissions caused by the change in a certain socioeconomic driver (determined by SDA). Through controlling variates from two sides, the average of the results represents the influences of the changes in socioeconomic drivers on health impacts.

DMC is the change in deaths due to meteorological changes, which can be determined by taking the average of the results of importing the same emissions into different years (different meteorological data) in the GEOS-Chem model; ∆DLD is the changes in lifetime lung cancer deaths caused by country-specific lung cancer deaths in different years, which can be obtained by taking the average of the results of importing the same concentrations into different years in the lifetime lung cancer risk assessment. The contributions of meteorological changes and lung cancer deaths can be calculated as:

$$\triangle {{{{\boldsymbol{D}}}}}_{{{{\boldsymbol{MC}}}}}=\frac{{{{\bf{1}}}}}{{{{\bf{2}}}}}\left[\left({{{{\boldsymbol{D}}}}}_{{{{\boldsymbol{t}}}},{{{{\boldsymbol{Q}}}}}_{{{{\bf{0}}}}}}-{{{{\boldsymbol{D}}}}}_{{{{\bf{0}}}},{{{{\boldsymbol{Q}}}}}_{{{{\bf{0}}}}}}\right)+\left({{{{\boldsymbol{D}}}}}_{{{{\boldsymbol{t}}}},{{{{\boldsymbol{Q}}}}}_{{{{\boldsymbol{t}}}}}}-{{{{\boldsymbol{D}}}}}_{{{{\bf{0}}}},{{{{\boldsymbol{Q}}}}}_{{{{\boldsymbol{t}}}}}}\right)\right]$$

(28)

$$\triangle {D}_{{LD}}=\frac{1}{2}\left[\left({D}_{{{LD}}_{t},{{PAF}}_{0}}-{D}_{{{LD}}_{0},{{PAF}}_{0}}\right)+\left({D}_{{{LD}}_{t},{{PAF}}_{t}}-{D}_{{{LD}}_{0},{{PAF}}_{t}}\right)\right]$$

(29)

where $${D}_{t,{Q}_{0}}$$ is the lung cancer deaths when importing the emission map of previous year 0 to the GEOS-Chem model with the year of t; $${D}_{{{LD}}_{t},{{PAF}}_{0}}$$ is the deaths when multiplying PAF results of the previous year 0 and the country-specific lung cancer deaths at the year of t. Other factors in the equations have similar expressions.

#Global #trade #drives #transboundary #transfer #health #impacts #polycyclic #aromatic #hydrocarbon #emissions #Communications #Earth #Environment