Trade-driven relocation of air pollution and health impacts in China

Recent studies show that international trade affects global distributions of air pollution and public health. Domestic interprovincial trade has similar effects within countries, but has not been comprehensively investigated previously. Here we link four models to evaluate the effects of both international exports and interprovincial trade on PM2.5 pollution and public health across China. We show that 50–60% of China’s air pollutant emissions in 2007 were associated with goods and services consumed outside of the provinces where they were produced. Of an estimated 1.10 million premature deaths caused by PM2.5 pollution throughout China, nearly 19% (208,500 deaths) are attributable to international exports. In contrast, interprovincial trade leads to improved air quality in developed coastal provinces with a net effect of 78,500 avoided deaths nationwide. However, both international export and interprovincial trade exacerbate the health burdens of air pollution in China’s less developed interior provinces. Our results reveal trade to be a critical but largely overlooked consideration in effective regional air quality planning for China.

a. These values appear too small, for unclear reasons.
b. No data, due to different classification of these sectors in Japan and/or the U.S. compared to China. 10

S1
Baseline scenario, production-based emissions from MEIC and EDGAR E_S1

S2
Exclude emissions embodied in international exports (E_GE) from each sector of all Chinese provinces based on S1 Substitute emissions embodied in domestic interprovince exports (E_DE) with those embodied in interprovincial imports (E_DI) for each sector of all provinces in S1 Exclude emissions embodied in the exports to US (E_US) based on S1 E_S4 = E_S1 -E_US

S6
Exclude emissions embodied in the exports to East Asia (E_EAsia) based on S1

Supplementary Methods
The linked multi-regional input output model Framework for the linked multi-regional input output model. We integrated the Chinese multiregional input-output (MRIO)  Table Establishment Method, we removed the processing and assembling trade of 30 provinces and other countries and then adjusted the charges for processing and assembling trade due to the international trade gap between the Chinese IO table and custom trade data 5 . The provincial sectoral import and export matrices are expressed as V Im and V Ex , respectively: where, v m is a column vector of sectoral import trade with other countries in province m, and v n is a column vector of provincial sectoral export trade to a country n.
We assumed that the roles of imported goods are similar to those of local goods in China's provinces, both of which are used to satisfy local industrial production processes or consumption directly as final demand. For example, the iron ore imported from Australia is mainly used as raw material in Chinese steel production, which is similar to the iron ore produced in China 6 . Therefore The intermediate input structure (i.e., production structure) and final demand structure of the sectoral import goods to the 30 provinces are assumed to be the same as the domestic proportion in each province in China. Therefore, we can use the "equal proportion disassembly" method, using the production structure and final demand structure of each province in the Chinese MRIO model to allocate the provincial sectoral import intermediate inputs and final demand sections to local production structures and final demand structures. Here, we introduced the production structure allocation coefficient matrix PIm and the final demand structure allocation coefficient matrix Q Im for Chinese provincial imports. The production structure matrix Z FC-CP and final demand matrix D FC-CP of provincial imports can be expressed as: Similarly, we assumed further that the allocation of other countries' intermediate input and final demand goods imported from China's provinces should be similar to their local production structure and final demand structure. Therefore, we adopted each country's production structure and final demand structure in the WIOD to disassemble the intermediate input and final demand goods from the individual Chinese provinces. Here, we introduce the production structure allocation coefficient matrix PEx and the final demand structure allocation coefficient matrix Q Ex of provincial exports. The production structure matrix Z CP-FC and final demand matrix D CP-FC for provincial exports can be expressed as: Compilation of trade data for construction and services. The international trade data for Chinese construction and services sectors are derived primarily from the Balance Table of International Payments. As there are no open official statistics of such data at the provincial level, we adopted the equal-proportion method to evaluate the provincial trade. We assumed that international imports in the construction and services sectors of each individual province have the same structures in terms of intermediate input and final demand as the domestically supplied sectors. Therefore, we disaggregated international imports for the Chinese construction and services sectors in the WIOD into 30 provinces according to the proportions of the export and import values of individual provinces in the Chinese provincial MRIO model (Supplementary Table 6).  Table 5), one as construction, and the remaining four as trade in services (Supplementary Table 6). The linked MRIO model can be expressed as:

Description
where the technical coefficient A * and final demand Y * are: where A rs is a 21  21 submatrix of A*, representing the technical coefficients of region s from region r, which can be calculated by dividing the total input row vector (X input ) s by the intermediate input matrix Z rs . Y rs is a 21  5 submatrix of Y*, representing the final demand of region s from region r. The total input matrix X input is the transpose of matrix X* (total output), which can be described as: As the consumption structure is distant from the objective of this paper, we had not specifically discussed the investment-related emissions although we included them in the consumptionbased emissions as a whole. Investment is described as being part of final demand in the IO table and accounts for a big share of the air pollutant emissions in China. According to our previous study 10 , for example, investment accounted for nearly 50% of the total consumptionbased SO 2 emissions in China, much higher than some developed countries (e.g. 23% in USA, 36% in Japan).
A previous study 7 calculated the sectoral international imports and exports of Chinese provinces by disassembling provincial international trade data according to the sectoral proportion of international trade at the national level. However, the availability of natural resources, economic development, industrial structure, and level of technological advance are significantly different among Chinese provinces, and thus international trade features are distinct in various provinces. The assumption of identical sectoral shares of international trade in various provinces will bias regional differences in international trade. To avoid such bias, we have accessed to the real trade matrix of Chinese provinces to other countries at the sector level 4 when developing our linked MRIO model. We adopted the WIOD instead of the Global Trade Analysis Project (GTAP) because the former provides more reliable information in the following respects. The original data of the WIOD originate from supply and use tables (SUTs) of countries, while those of the GTAP are obtained from their IO tables. The SUTs provide more detail on local industry, trade, production, socioeconomic conditions, environment, and other factors compared with IO tables. And, some IO tables are even compiled from the SUTs. In addition, the WIOD consists of full time-series data, while the GTAP only has data for specific years 2 .
The raw data used in linking the Chinese provincial MRIO model with WIOD are derived primarily from statistics from official sources (e.g., statistical yearbooks and customs data), authoritative research centers, and peer-reviewed materials. The individual province IO tables themselves have inherent uncertainties. For example, the quality of IO tables in developed provinces (e.g., Beijing and Shanghai) tends to be better developed than those in less developed regions 11 . The data decomposition, transformation, and deduction processes used to compile the Chinese provincial MRIO model comprising 30 provinces in 2007 1 introduced additional uncertainties. Different statistical methods used to collect these raw data can also introduce uncertainty. For example, data on goods trade from customs underwent some revisions to match the data needs of the MRIO model, which increased the uncertainty of the international trade data. However, national and local official statistics are the most reliable data sources available in research on climate change and macro-environmental policy 12 .
It should be noted that detailed data at a smaller scale usually have higher uncertainties compared to those at larger scales [13][14] . Therefore, uncertainties were introduced during the process of integrating, matching, and decomposing international trade data at the sectoral level when we linked the Chinese provincial MRIO model with WIOD. Fortunately, the uncertainty of small-scale data in IO analysis has less effect on the final result because the error can be largely counteracted by adding or multiplying the small-scale data using Monte Carlo methods in IO analysis 14,[15][16] . Therefore, the analysis of trade-related health impacts at China's macro provincial level, based on our linked MRIO model, is expected to be comparatively reliable.

Supplementary
where E* is the environmental impact, which in this analysis refers to air pollutant emissions (e.g., SO 2 , NOx, and PM 2.5 ) from anthropogenic activities; F is a row vector for per unit output of emissions (i.e., emission intensity) for each sector in the 30 Chinese provinces and the other 40 countries and regions; fj r is an element of a sub-vector F r in F, which can be calculated using  Because the sector classifications of the emission inventories and our linked MRIO model differ, we conducted a mapping process to relate them to each other.
The MEIC is a unit/technology-based, bottom-up air pollutant emission inventory for China, updated from the widely used INTEX-B data set 17 , which provides rich emission data from more than 700 emission sources and production categories 18 . This detailed information enables the direct allocations of most emission sources in MEIC into economic sectors according to the products or/and usage categories, such as the emission sources of power plants, industrial processes, and transportation. For the emissions from residential and service sectors and industrial boilers, data from the China Energy Statistical Yearbook 19 Supplementary Fig. 9). The major emission sources in EDGAR could be mapped to our linked MRIO model sectors according to their products or/and usage categories, such as the emission sources of public electricity and heat production, and transportation. For the emissions from other sources (the italic labeled sectors in Supplementary Fig. 9), e.g., the residential and other sectors, the corresponding sectoral energy consumptions from IEA (http://www.iea.org/statistics/) or the economic activities were used as proxies to allocate the aggregated emissions into the corresponding 21 economic sectors classified in our linked MRIO model.

Supplementary
where E d rr denotes emissions associated with the products both produced and consumed in region r, and

Simulations of impact of trade on ambient PM 2.5 in China
GEOS-Chem model. The impact of trade on Chinese air quality is simulated by the nested-grid chemical transport model, GEOS-Chem (v10-01, http://geos-chem.org/). Simulations are performed at a horizontal resolution of 1/2°latitude by 2/3°longitude over the East Asian region including China. This high-resolution model domain is embedded in a global chemical transport simulation of horizontal resolution 4°latitude by 5°longitude, with the latter providing initial and boundary conditions for all simulated chemical species. The model is driven by the meteorological data from the Goddard Earth Observing System (GEOS, version 5) of the NASA Global Modeling Assimilation Office (GMAO). It contains 47 vertical layers up to 0.01 hPa. GEOS-Chem uses the same advection algorithm as the GEOS general circulation model (http://gmao.gsfc.nasa.gov/GEOS/). Convective transport in GEOS-Chem is computed from the convective mass fluxes in the meteorological archive. Boundary layer mixing in GEOS-Chem is calculated by a non-local scheme 24 . Wet deposition by rain is considered for both water-soluble aerosols and gases, and scavenging by snow and cold/mixed precipitation is also considered for aerosols 25 . Dry deposition is calculated based on the resistance-in-series scheme for all the species, with gravitational settling for dust and coarse sea salt 26 .
The simulation contains a gas phase HOx-NOx-VOC-ozone-BrOx chemistry, which considers the production and loss of ozone through reactions with HOx, NOx, VOC and BrOx [27][28][29] . GEOS-Chem includes detailed sulfate-nitrate-ammonium-carbonaceous-dust-sea salt aerosol chemistry, which is coupled to gas phase chemistry. The model considers the thermodynamics of inorganic aerosols and the in-cloud sulfate formation based on cloud water pH. GEOS-Chem also simulates the dust and the sea salt aerosol in different size bins. Aerosols interact with gas-phase chemistry in GEOS-Chem through the effect of aerosol extinction on photolysis rates, heterogeneous chemistry, and gas-aerosol partitioning of semi-volatile compounds.
GEOS-Chem simulations require information on the spatial distributions of emissions. The impacts of trade on emissions, however, can only be quantified at the provincial level in China through the linked the WIOD and MRIO model. Following the approach adopted by earlier studies [30][31] , we assumed that trade-related emissions exhibit the same spatial distribution as total emissions specified in the MIX Asian emission inventory (http://mix.greenresource.cn), which has already integrated the MEIC as an emission inventory for China and has been applied in GEOS-Chem. In the current research, the emission inventories for Scenarios 2-6 (Supplementary Table 2) can be used as inputs for the GEOS-Chem model. The original emission inventory has a spatial resolution of 1/4°by 1/4°and was re-gridded into the GEOS-Chem grid by the HEMCO emission module 32 . Regarding species for which we have estimates from all anthropogenic sources, including SO2, NOX, NMVOCs, NH3, CO and PM2.5 (BC and OC are scaled for the same factors as PM 2.5 because of similar source contributions) [33][34][35][36][37][38] , the emissions from each grid cell E i,j GEOS-Chem are calculated based on the MIX inventory E i,j MIX adjusted by the ratio of provincial emissions from our scenarios E k MRIO and the MIX inventory E k MIX : where i and j are the grid box index in latitudinal and longitudinal directions, respectively, and k is the provincial index for this grid box.
The anthropogenic emissions outside of the East Asia domain, which are required by the coarser resolution global model, are obtained from the EDGAR global emission inventory, with various For all non-anthropogenic emissions, such as wild fires and lightning, we used the default emission inventories recommended by the GEOS-Chem model (http://wiki.seas.harvard.edu/geos-chem/index.php/Emissions_overview). All emissions are assigned to a province based on their geographical locations, and the trans-boundary transports of emissions are simulated by GEOS-Chem.
Modelling the effect of trade on PM 2.5 pollution. In the ideal case, exposure can be estimated using a perturbation in emissions with compared to a baseline case to evaluate the impacts of emission changes on human exposure 39 . This comparison can be based on measured concentrations before and after a natural "experiment", such as the strict air-quality controls enacted for the 2008 Beijing Olympics or the Clean Air Act in the US [40][41] .
Alternatively, studies can employ chemical transport models such as GEOS-Chem to simulate PM 2.5 concentrations based on different alternative emission scenarios relative to a baseline scenario to analyze the impact of policies on air pollution 39 . We also employ this sensitivity analysis to evaluate the impact of various trade activities on the ambient PM 2.5 pollution and related premature deaths across China. We simulate ambient PM2.5 concentrations across China by running GEOS-Chem using emissions of the above-defined 6 scenarios. The spatial change ratios of modeled PM 2.5 concentrations of Scenario 2-6 compared to Scenario 1 represent the impacts of various types of trade (e.g., international exports from China, exports to the US only, interprovincial trade) on PM2.5 concentrations across China (Supplementary Fig. 1). These spatial change ratios were then multiplied by satellite-retrieved PM 2.5 concentrations developed by van Donkelaar et al. 42 to calculate the ambient PM2.5 concentrations influenced by various types of trade.
Similar methods have already been widely applied in earlier studies to quantify the impact of policies on air quality or public health, such as the health effects attributable to policies addressing climate change [43][44] and the health benefit of controlling emissions from power plants [45][46] . We should note that the relationship between emissions and ambient PM 2.5 concentrations is nonlinear, which means that the removal of individual emission sources does not engender a linear response in the modeled contributions to ambient PM 2.5 concentrations (i.e., the sum of modeled contributions of various sources to ambient PM2.5 concentrations does not necessarily equal 100%). For example, Lelieveld et al. 47 evaluated the global premature deaths attributable to ambient PM 2.5 from various emission sources (e.g., power plants and transportation) and applied scaling corrections (about 10%) to ensure the contributions of various sources add up to 100% at the country level. Lin et al. 31 applied a similar method to evaluate the impact of Chinese emissions on O 3 pollution in the US, although the nonlinear chemistry in O 3 formation 48 is even more complicated than that of secondary PM 2.5 , which is dominated by local emission sources (e.g., over 60% of PM 2.5 concentrations in each Chinese province are contributed by emissions from the province itself) [49][50] . Although it is not perfect, in this sense, such methodology offers an effective means for the purpose of our analysis.
Model evaluation and uncertainties. The GEOS-Chem model adopted in this analysis has been widely used for studying air quality of different regions across the world including North America 36, 51-54 Europe 55-58 , and Asia [59][60][61][62][63] . Detailed model performance has been extensively evaluated against observed data obtained from satellite retrievals, ground-based observation sites and networks, and aircraft campaigns 60,[64][65][66][67][68] . To correct any potential model bias caused by the relative coarse resolution, imperfect meteorological data, or incomplete atmospheric chemical reaction schemes, we used the modeled relative changes in concentrations ( Supplementary Fig. 1) of different scenarios (as defined in Supplementary Table 2) and then applied them to surface concentrations derived from satellite observations. As satellites only measure the column density of aerosols through retrieval of aerosol optical depth (AOD), van Donkelaar et al. 42,69 applied GEOS-Chem-modeled aerosol vertical profiles and vertical profile information from the CALIOP (Cloud-Aerosol Lidar with Orthogonal Polarization) satellite instrument to allocate the retrieved AOD to near-ground concentrations. According to Boys et al. 70 and van Donkelaar et al. 42 , the uncertainty in satellite-derived PM 2.5 decreases with an increase in sampling days, and annual mean PM 2.5 satellite retrievals (1999-2011) are estimated by using the 3-year moving average from 1998 to 2012. The resultant ground-level PM 2.5 estimates are highly consistent (R 2 =0.81) with PM 2.5 concentrations from in-situ surface monitors (the Figure 5 in Donkelaar et al. 69 ). This satellite-derived PM 2.5 concentration dataset 42 has been applied broadly in the literature, notably including the Global Burden of Disease (GBD) assessments 71 , to represent the spatiotemporal distribution of PM 2.5 exposures globally.

Evaluation of the impact of trade on public health in China
Background of the health impact assessment. The air pollution to which individuals are exposed is multifaceted, typically including numerous individual gaseous compounds and particles of complex physicochemical composition 72 . Accordingly, indicator pollutants are often used to assess exposures for risk assessment and epidemiologic analysis. Evaluation of adverse health impacts associated with elevated ambient PM 2.5 concentrations has been supported by extensive epidemiological literature [73][74][75] , and PM 2.5 is the most robust indicator of adverse impacts (e.g., mortality) of long-term exposure to air pollution in epidemiologic cohort studies 76 .
Epidemiological associations between elevated ozone concentrations and premature mortality have also been observed, which are independent of associations between PM 2.5 and mortality [77][78] . Estimates of the global burden of disease attributable to outdoor air pollution can be enhanced by the inclusion of ozone in addition to PM 2.5 , as in Lelieveld et al. 47 and the GBD 71 . However, the premature deaths attributable to O 3 exposure are far fewer, around 5% globally compared to that of PM 2.5 in 2010 47,71 . Relatively scarce O 3 monitoring data in China make current evaluation of large-scale health impacts attributable to O 3 mainly dependent on atmospheric chemistry models (e.g., TM5, CEOS-Chem, EMAC, or CMAQ), which introduce large uncertainties due to the more complicated chemical reactions of O 3 production compared to secondary PM 2.5 48 .
For this reason, we focused on ambient PM 2.5 in this study and applied an integrated exposureresponse (IER) model 79 to evaluate the impact of trade on premature deaths attributable to air pollution. The IER model is a mass-based concentration-response (C-R) model that integrates the relative risks (RR) of death from a variety of illnesses, including ischemic heart disease (IHD), cerebrovascular disease (stroke), chronic obstructive pulmonary disease (COPD) and lung cancer (LC), attributable to exposure to PM2.5 over diverse concentrations (i.e., ambient air pollution, active smoking, secondhand smoking, and indoor burning of solid fuels). It thus covers the global exposure range of PM2.5, including the high concentrations of countries like China and India. The IER model has already been widely applied, including in the GBD, to evaluate the health impacts attributable to ambient PM 2.5 in various countries. Burnett et al. 79 compared the RRs observed in the limited Chinese cohort and those predicted by the IER model. The similarity for different health endpoints (e.g., IHD, Stroke, LC) suggests that the IER model can yield reasonable predictions of the change of risk over the range of concentrations that prevail in China.

Uncertainty and limitations.
It is known that analysis of the health impacts of air pollution generally has large uncertainties [80][81] . Although we have used state-of-the-art methods to estimate the impact of trade on mortality attributable to ambient PM2.5 pollution in China, there are a number of assumptions and limitations inherent in our estimates. First, the IER model itself has uncertainty. Besides the issues such as the shape of the concentration-response functions and the possible existence and specific levels of C-R thresholds, there is considerable evidence that the chemical composition, size distribution, and sources of PM 2.5 may influence the health effects 82 . However, the relative toxicity of various constituents of ambient PM 2.5 has not been well established 47 , and definitive estimates for various chemical forms or sources would go beyond the capacity of current scientific evidence required for accurate determination 83 . Therefore, current cohort epidemiologic results for adverse (e.g., mortality) impacts are still based on long-term exposure to PM 2.5 total mass. Second, uncertainty of trade-related emissions results from that of emission inventories (i.e., MEIC and EDGAR) and the economic input-output analysis. For example, the uncertainties of the MEIC were estimated to be 12% for SO2, 31% for NOX, 68% for NWVOC, and 107% for primary PM2.5 17,84 , while the uncertainty of input-out analysis might be less than 50% 31 . Third, GEOS-Chem simulations are unavoidably affected by errors in meteorological inputs and imperfect representation of tropospheric chemistry. The importance of GEOS-Chem model uncertainties, however, might be substantially reduced because we focus on change ratios between the modeled concentrations of the baseline and alternative scenarios (Supplementary Table 2), not the absolute concentration levels simulated by GEOS-Chem.
Epidemiologic studies suggest logarithmic relationships between ambient PM 2.5 concentrations and RR 79,81 , and uncertainties in the C-R functions linking ambient PM 2.5 concentrations with mortality are particularly important for health impact assessment 80 . Therefore, we evaluated the uncertainties of the premature mortality analysis associated with the C-R function for PM2.5 based on the methods of the GBD project 47,71 , and emphasize that the 95% confidence intervals of premature death described in this study reflect only the statistical uncertainty of the RR calculated through the IER model.