We create four emission scenarios for emissions in 2017: actual emissions (2017 EFs) and three counterfactual scenarios in which county-level EFs for each vehicle type are the same as they were in 2014 (2014 EFs scenario), 2011 (2011 EFs), and 2008 (2008 EFs). Each scenario uses data for the NEI for the respective year (24⇓⇓–27), which is provided as total emissions (except refueling) and VMT by vehicle type, for each county. We calculate county-level EFs for each of the 13 vehicle types currently used by the EPA and apply them to their VMT in each county in 2017 (SI Appendix, Sections 1 and 2). These counterfactual scenarios therefore capture changes in fleet composition and VMT distribution. Refueling emissions for NEIs 2011, 2014, and 2017 are provided separately from vehicle emissions for each county, and we scale them up according to the total fleet VMT in each county. NEI emissions include only all on- and off-network processes and refueling, but not other nonroad life cycle emissions, which are therefore also not included in our study. For the state of California, we complement missing GHG emissions in NEI 2008 and missing N2O emissions in NEI 2014 with GHG emissions data from the California Air Resources Board (37).
We assess mortality attributable to chronic exposure to PM2.5 because it accounts for the vast majority of overall monetized air pollution impacts. Our analysis does not consider the following: 1) the effects of ozone, 2) the impacts of acute exposure, or 3) the many nonfatal effects of PM2.5 exposure. Mortality attributable to ozone has been estimated to be about an order of magnitude smaller than mortality attributable to PM2.5 for the transportation sector in the United States (1, 3, 12). Mortality attributable to acute exposure to PM2.5 has also been found to be much smaller than mortality attributable to chronic exposure. Even in China, where the haze episodes are severe, mortality attributable to acute exposure was an order of magnitude smaller (38). Finally, the incidence of nonfatal outcomes attributable to PM2.5 is higher, but an analysis of the benefits of the Clean Air Act from 1990 to 2020 attributes over 90% of the estimated monetized benefits in 2020 to PM2.5-attributable mortality (11).
We estimate the marginal impact on mortality associated with marginal changes in emissions of each of the five pollutants in each county. We follow the method used by Choma et al. (20), refined with county-level age-specific baseline mortality data and baseline ambient PM2.5 levels at a finer spatial resolution (1 km) (21). We also relax the assumption made by Choma et al. (20) that impacts occur in the same place as emissions to better capture spatial variability in baseline mortality rates and ambient PM2.5 levels. Baseline ambient PM2.5 levels affect marginal impacts if, as we assume, the relationship between ambient PM2.5 concentrations and mortality risk is nonlinear. The approach is based on three main components: 1) an estimate of the impacts of the ground-level emissions on fine particle exposure (i.e., changes in concentrations), 2) a CRF relating PM2.5 concentration to mortality risk, and 3) baseline mortality rates.
For the first component, we use the ISRM (2, 22, 36), which estimates changes in fine particle concentrations in a receptor cell as a function of emissions of each of the five air pollutants covered in this study in each cell. It has the advantage of fine spatial resolution with variable cell sizes that are as small as 1 × 1 km for densely populated areas. We map ISRM cells to US counties weighting by population as described in Choma et al. (20) but using more recent population estimates. We use 5-y estimates of population at the census block group level from the 2015 to 2019 American Community Survey (39) and geography boundaries from the US Census Bureau (40, 41).
For the second component, we use the GEMM CRF (16) because it incorporates the most recent epidemiological evidence coming from studies in North America, extending the range of exposures and risk estimates to concentrations as low as 2.4 μg/m3. The authors (a collaboration among research groups that includes those responsible for 15 of the largest individual cohorts to date) fit a unified model to individual-level data from these cohorts as well as published data for 26 other cohorts. Burnett et al. (16) provide different sets of coefficients, and we use the age-specific GEMM NCD+LRI coefficients, encompassing all nonaccidental deaths. We use the GEMM version that includes evidence from a recent Chinese male cohort study, although this inclusion did not substantially affect GEMM’s CRF and also did not result in substantial differences when estimating mortality in the United States in Choma et al. (20). The GEMM CRF is concave in ambient PM2.5 concentrations; therefore, risk ratios are a function of the baseline ambient PM2.5 concentration. This nonlinearity is supported by two other important recent syntheses of the epidemiological evidence, both of which estimated concave CRFs: GBD IER (17, 18) and the metaregression by Vodonos et al. (32). Nonlinearities in effects with respect to baseline PM2.5 concentrations are particularly relevant at the very low concentrations (<10 μg/m3) experienced by the vast majority of the US population (21).
We calculated baseline PM2.5 concentrations for each county by weighting concentration estimates at a 1-km resolution (21) by population at a census block level (n = 11,007,989) from the 2010 Decennial Census (42). We assigned an annual average PM2.5 level to each census block using the distances from block centroids to the centroids of the four closest 1-km model cells, weighting by inverse distance. The estimated concentrations result in national population-weighted averages of 10.5 μg/m3 in 2008 and 7.5 μg/m3 in 2016, similar to the decrease from 10.9 to 7.7 μg/m3 in those years shown by the EPA based on data from 406 trend sites (10). We used model outputs from Di et al. (21) for 2016, the last year available, as a proxy for 2017. The EPA data show a small increase from 7.7 μg/m3 in 2016 to 8.1 μg/m3 in 2017.
For the third component we used, for each county and GEMM age group, nonaccidental baseline mortality data from the Centers for Disease Control and Prevention (CDC) Wonder database (23) and population data from HHS (29). For each county and age group, we calculated a mortality rate from 5 y of deaths and population counts and multiplied it by the population of the year of interest (2008 and 2017) to estimate the number of deaths. For 2008, we used data from 2006 to 2010, whereas, for 2017, we used data from 2014 to 2018, the last 5 y available. For counties with <50 death counts in any given age group in the 5-y period, we used the state-level mortality rates for that age group and applied it to the county population for that age group. In our sensitivity analysis using other CRFs that estimate risks for all-cause mortality, we used the same procedure and data sources but collected data for all ages and causes of death.
Marginal impacts for emissions of each pollutant are assessed with Eq. 1. The ISRM provides increases in concentrations (ΔC) caused by an emission of 1 μg/s (0.03 kg/year); therefore, our marginal impacts are assessed for emissions of 0.03 kg of each pollutant. While GEMM’s CRF is nonlinear in ambient PM2.5 levels, even if we eliminated all 19,800 traffic-related PM2.5-attributable deaths we estimate for 2017, the change in GEMM’s slope and in the marginal effect would be smaller than 10%. This would cause an average error in our model that is smaller than 5% because we apply the same marginal impacts to all transportation emissions. Our model also does not include deaths occurring outside of the United States. We assume that any air pollution exports to other countries is very small—even exports across US states are responsible for just 24% of the total attributable deaths.Marginal Impacts,p=∑r∑aMr,a(ΔCs,r,p)=∑r∑a(RRa,d(Cr+ΔCs,r,p)−RRa,d(Cr) RRa,d(Cr)×Mr,a),where C is the PM2.5 concentration [μg/m3], ΔC is the increase in concentrations [μg/m3], RR is the risk ratio from the CRF [dimensionless], and M is the outcome measure [mortality in deaths/year], which is cause-specific (nonaccidental) for GEMM and all-cause for the other CRFs.
The index s represents the sources (52,411 InMAP cells, which we map to 3,108 US counties), r are the receptors (52,411 InMAP cells, which we map to 3,108 US counties), p = 1,2,…,5 are the five pollutants, and a are the age groups (12 for GEMM and 1 for all-cause mortality). In our analysis of impacts occurring in state and out of state, we sum receptors accordingly (r inside the state of s and r outside the state of s).
The computational implementation of Eq. 1, including the mapping of InMAP cells to US counties, is shown in Eq. 2. Eq. 2 is applied to each pollutant separately, as the ISRM is specific to each pollutant. The other matrices and parameters are the same for all pollutants.MI=P×ISRM×PT×M×k×VSL,where MI is a 3,108 × 3,108 matrix where MIi,j is the impact [2017 US dollars] occurring in county j (receptor) as a consequence of 1 metric ton of emissions in county i (source); P is a 3,108 × 52,411 matrix where element Pij is the percentage of the population of county i that is within InMAP cell j; ISRM is a 52,411 × 52,411 matrix where ISRMij is the increase in concentration (ΔC) [μg/m3] in InMAP cell j (receptor) as a consequence of emissions of 1 μg/s in InMAP cell i (source); M is a 3,108 × 3,108 diagonal matrix where Mij is the increase in mortality (deaths) for an increase in 1 μg/m3 in ambient concentration if i = j and Mij = 0 if i ≠ j; k = 1012/(24 × 3,600 × 365), representing the conversion from 1 μg/s to 1 metric ton/year; and VSL is the VSL [2017 US dollars] in the case of monetized damages.
M was calculated as M = Dg(diag(D × ST)), in which D is a 3,108 × Na matrix, where Dij is the number of deaths in county i and age group j; S is a 3108 × Na matrix, where Sij is the percent increase in baseline mortality in county i and age group j for an increase in ambient concentrations of 1 μg/m3; diag(X) denotes the vector containing the diagonal elements of matrix X; and Dg(x) denotes the square matrix where off-diagonal elements are 0 and diagonal elements are the elements of vector x.
For GEMM, D contains data only for nonaccidental mortality, and S only applies to those causes of death. For the other CRFs, all deaths are accounted for in D, and S applies to all-cause mortality. Na is the number of age groups, which equals 12 for GEMM and 1 for the other CRFs.
Matrix P maps InMAP cells to counties, allocating county death counts (in the case of receptors) and emissions (in the case of sources) to the InMAP cells within them, weighting by population and implicitly assuming that within-county spatial distributions of mortality and vehicle emissions follow that of population.
The shape of the CRF at very low concentrations is particularly uncertain, with a lack of evidence for exposures below the GEMM minimum (2.4 μg/m3). We assign zero marginal impacts for changes in concentrations in counties below that exposure level, which include 0.02% of the population in 2017 and 0% in 2008. Moreover, while we calculate marginal impacts for any level above 2.4 μg/m3, any policy that leads to meaningful public health improvements will lower concentrations by more than a marginal amount. If we consider that amount to be 0.5 μg/m3, it would include 0.1% of the population in 2017 and 0.006% in 2008 (i.e., that experienced baseline levels <2.9 μg/m3 or no more than 0.5 μg/m3 above the threshold). In these rare cases, while we estimate positive marginal benefits, meaningful policies of emission reductions might have little effect. We use the same minimum level of 2.4 μg/m3 for all CRFs, extending the risk estimates by Vodonos et al. (32) and Krewski et al. (34) to that level, even as the studies were not able to estimate risks for such low exposures.
We present our air pollution results in 2017 US attributable deaths as well as their monetized value in 2017 US dollars. The monetary value is calculated using a $10.2 million VSL, reflecting the $9.3 million in 2014 used by HHS (35) adjusted to 2017 for inflation [US gross domestic product (GDP) deflator from World Bank (43)] and income [median usual weekly earnings from US Department of Labor (44)], the latter assuming an income elasticity of the VSL of 1. There is a “cessation lag” between changes in emissions and changes in mortality, and we apply the cessation lag structure recommended by the EPA’s Advisory Council on Clean Air Compliance Analysis (45), discounting benefits using a 3% discount rate in our monetized results. In the EPA’s recommended structure, 30% of the benefits occur in the first year, another 50% uniformly in years 2 to 5, and the remaining 20% uniformly between years 6 and 20. This results in a net present value of 0.89 × $10.2 million, or $9.1 million, for each fatality. Our results presented as attributable deaths, on the other hand, are simply the undiscounted sum of attributable deaths that occur in years 1 through 20.