A coupled global and regional climate modeling system was developed to provide 2001–2004 and 2057–2059 air pollution levels in the continental U.S. at 12 km resolution. The configuration and evaluation of this system are briefly described here (More details are addressed in our previous study, Gao et al. 2013). The state-of-the-art earth system model Community Earth System Model version 1.0 (CESM1.0), developed at the National Center for Atmospheric Research (NCAR), can simultaneously simulate the Earth’s atmosphere, ocean, land surface, and sea-ice (Gent et al. 2011). In this study, the atmospheric chemistry was integrated with the atmospheric component, Community Atmosphere Model, referred to as the CAM-Chem (Lamarque et al. 2012). CESM/CAM-Chem was used in the long-term global climate and chemistry simulations, from 1850 to 2005 as present climate and 2005 to 2100 for RCP4.5 and RCP8.5, with a spatial resolution of 0.9° (latitude) ×1.25° (longitude). The CESM-projected meteorological fields were used as the initial and boundary conditions of the regional climate model, Weather Research and Forecasting model (WRF) 3.2.1 (Skamarock and Klemp 2008). Chemical species simulated by the atmospheric chemistry component of the CESM/CAM-Chem were mapped to the regional air quality model, Community Multi-scale Air Quality Model version 5.0 (CMAQ 5.0; Wong et al. 2012).

We selected RCP4.5 and RCP8.5 to evaluate and compare future O₃ level and its related health impact. RCP4.5 corresponds to the stabilization of radiative forcing at 4.5 W/m² in the year 2100 without ever exceeding that value, and it assumes a lower population growth rate, high income, and rapid technology development (Thomson et al. 2011). RCP8.5 is a high-GHG emission pathway corresponding to the stabilization of radiative forcing at 8.5 W/m² in the year 2100, and it assumes high population growth rate, lower income, and lower technology development rate in developing countries (Riahi et al. 2011). The emissions trends in the RCPs are mainly controlled by the driving forces (population, income, energy use, etc.), climate policy, and air pollution control policy (Van Vuuren et al. 2011). In terms of climate policy, RCP4.5 has more stringent GHG emissions control, while RCP8.5 assumes no climate policy, even though both scenarios are projected to increase the use of renewable resources. Unlike the SRES scenarios, where anthropogenic emissions related to air pollutants have a modest decrease or even increase (under A1FI), all RCPs include stringent air pollution control policies with any increase in income. Thus, under the RCPs, O₃ precursors including carbon monoxide (CO), nitrogen oxides (NOx), and non-methane volatile organic compounds (NMVOCs), are projected to decrease in the U.S. (See Table 2 in Gao et al. (2013) for projection factor for anthropogenic emissions). However, the GHG emissions increase, particularly the methane concentrations under RCP8.5, will negatively affect O₃ air quality. As a result, the projected O₃ levels in 2057–2059 in this study reflect the influence of both climate change and emissions control on O₃ precursors.

Using CMAQ-simulated hourly O₃, we computed annual means of maximum daily 8-hour average O₃ (MDA8 O₃) and maximum daily 1-hour average O₃ (MDA1 O₃) for each 12 km grid cell. We then aggregated the 12×12 km grid of O₃ concentration to the county level (3,109 counties in the continental U.S.). We determined the population-weighted centroid of each county and then averaged the data in the nine grid cells closest to the centroid. We calibrated the simulated O₃ level for bias reduction. To compare the CMAQ-simulated values with observations at the county level, we used monitored data from the U.S. Environmental Protection Agency (USEPA) Air Quality System (USEPA-AQS) and the Environmental Benefits Mapping and Analysis Program Community Edition 1.0.8 (BenMAP-CE) developed by the USEPA (USEPA 2012). Using BenMAP-CE, we calculated county-level O₃, applying the Voronoi neighbor averaging algorithm, which interpolates monitored O₃ data to unmonitored locations. BenMAP-CE first identifies the set of monitors that best “surround” the population center of the grid cell and then takes an inverse-distance weighted average of the monitoring values (Fann et al. 2012; USEPA 2012). We matched the CMAQ-simulated values with the interpolated, county-level observations for 2001–2004. After calculating the county-level ratios of CMAQ simulations over observations in the base years, we calibrated the CMAQ-simulated MDA8 O₃ and MDA1 O₃ concentrations for both 2001–2004 and 2057–2059 using these calibration ratios for each county (Figure S1, Supplementary Material).

The county-level population projections for 2055 were taken from the Integrated Climate and Land-Use Scenarios (ICLUS) project (USEPA 2009). The ICLUS population projections were based on a standard cohort-component method, in which the initial U.S. population in 2005 was divided into age-, gender-, and race/ethnicity-specific cohorts. Cohort size was estimated over time separately for each cohort, with various scenarios regarding fertility, mortality, and migration (Preston et al. 2001). In ICLUS, four scenarios were considered: A1, low fertility, high net domestic migration, and high net international migration; B1, low fertility, low net domestic migration, and high net international migration; A2, high fertility, high net domestic migration, and medium net international migration; and B2, medium fertility, low net domestic migration, and medium net international migration, with the categories (i.e., low, medium, and high) defined in the IPCC-SRES. Mortality rate was kept constant across all scenarios because the U.S. Census Bureau did not release alternative scenarios for mortality rates. Although the four ICLUS scenarios did not cover the widest possible range of population size, they represented combinations of potential fertility, mortality, and migrations in the future. A more detailed discussion on ICLUS population projection was presented elsewhere (USEPA 2009).

To calculate baseline future mortality incidence without the impact of climate and O₃ precursors changes, we used predicted county-specific mortality rates adopted in BenMAP-CE, which is derived from a projected, age-specific ratio of the 2050 mortality rate to the 2005 mortality rate. The mortality rates do not, however, explicitly account for climate change. The basis for the future mortality rate projection is described in detail elsewhere (USEPA 2012). We used non-accidental, all-cause mortality to be consistent with the CRF used in this study.

To calculate EM due to O₃ exposure, we based the CRFs on the association between non-accidental, all-cause mortality and short-term exposure to MDA8 O₃ (Bell et al. 2004) and MDA1 O₃ (Bell et al. 2004; Levy et al. 2005). Bell et al. (2004) analyzed the relationship between O₃ and non-accidental, all-cause mortality using the National Morbidity, Mortality, and Air Pollution Study (NMMAPS) dataset, which covers 95 U.S. cities. Levy et al. (2005) derived the relative risk from meta-analyses of 27 studies for North America. There might be publication bias in the meta-analyses results. Nevertheless, the cities chosen in these studies have a wide range of O₃ concentrations and represent both rural and urban populations. We therefore applied the CRF coefficients in Bell et al. (2004) and Levy et al. (2005) (Table S1, Supplementary Material).

Finally, we estimated the change in county-level EM as follows (Fann et al. 2012; Post et al. 2012): (1)Δyi=POPi×MRi×eβ×ΔCi-1 where Δy is the expected number of deaths per year attributable to changing O₃ level at county i; POP_i is population in county i; MR_i is population mortality rate; POP_i×MR_i indicates baseline mortality incidence (i.e., assuming no climate change and O₃ precursors emission change); β is the coefficient of the concentration-response function for O₃; ΔC_i is the concentration difference of O₃ between 2057–2059 and 2001–2004.

We estimated EM for non-accidental, all-cause mortality associated with O₃ for both year round and warm season (May through September). Lastly, we analyzed the spatial variation of EM by county and the nine climate regions defined by the National Climatic Data Center (http://www.ncdc.noaa.gov/temp-and-precip/us-climate-regions.php) (Fig. 1).

To determine the sensitivity of the EM to various factors (i.e., two RCPs, four future population projections, and three CRF coefficients), we conducted an analysis of variance (ANOVA) for EM estimates to decompose the total variability in estimated EM into the contribution of each factor. In contrast to the approach in a previous study (Post et al. 2012) which used national-level EM estimates, we conducted the ANOVA by county and then calculated the county-level percentage of sum of squares (SS) for RCPs, population projections, and CRF coefficients. We also plotted the county-level percentage of SS to visually analyze the spatial distribution of the sensitivity to the factors.

The ANOVA apportions the variability in the mean estimated air pollution EM in each county to the mean levels of each factor in Equation 1 without considering their uncertainties (i.e., each factor has a distribution of possible values). We used Monte Carlo simulations to evaluate the uncertainty of EM estimates attributable to the ranges of CRF coefficients, mortality rates, and concentration changes of O₃. CRF coefficients and mortality rates are assumed to be normally distributed with independent, county-specific means and standard errors. As there is no standard error for the projected 2050 mortality rate, we used the standard error for the year 2009 at the county level provided by the Centers for Disease Control and Prevention (CDC 2012). For concentration changes of O₃, triangular-distributed random variables were generated using the minimum, maximum, and mean of three annual mean values. Random sampling and EM calculations were repeated 1,000 times for each county. Then, we computed the Monte Carlo means and standard errors of the 1,000 EM estimates by county. By summing all the county-level EMs derived from the Monte Carlo simulations, we estimated national-level EM estimates and their 95 % confidence intervals (CIs). We also calculated the probability distributions of all the national EM due to future O₃ change.

The spatial distributions of county-level O₃ change between the future (2057–2059) and the base years (2001–2004) are shown in Figure 2. Projected O₃ concentrations by climate region are presented in Supplementary Material, Table S2. Changes in estimated county-level MDA8 O₃ ranged from −8.8 to 7.5 ppb under RCP4.5 and from −3.6 to 11.5 ppb under RCP8.5 across 3,109 counties. Higher O₃ increases in the Northwest, West, and West North Central region are likely caused by higher background O₃ in the spring and winter, as well as an increase of methane emissions in RCP8.5 (Gao et al. 2013). Nationally, O₃ was projected to decrease under RCP4.5 and to increase under RCP8.5. However, O₃ levels of most counties in the Northwest, West North Central, and West region were projected to increase under RCP8.5 (Fig. 2(a)). Even under RCP4.5, some counties were expected to have higher O₃ levels (Fig. 2(b)).

Table 1 shows the estimated national-level O₃-related EM changes between 2001–2004 and 2057–2059 under the two RCPs. We derived national EM estimates in Table 1 from the summation of county-level EMs for each ICLUS population scenario and the average of the four ICLUS-based national EMs. Our 95% CIs in Table 1 are based on 4,000 national-level EMs (1,000 times of Monte Carlo simulations for four population projections) for each RCP. O₃-related EMs vary by emission pathway, population projection, and CRF. Under RCP8.5, increased MDA8 O₃ resulted in ~1,900 (95% CI: 1,700 to 2,100) premature deaths per year nationwide, when county-level EMs derived from the Monte Carlo simulation were summed using the four population projections in 2057–2059 and CRF from Bell et al. (2004). In contrast, decreased MDA8 O₃ under the RCP4.5 resulted in the avoidance of ~1,600 (95% CI: 1,300 to 1,800) deaths per year. Projected county-level baseline mortality incidence for 2055 is shown in Supplementary Material, Figure S2. During the warm season, lower future MDA8 O₃ levels could result in ~2,400 (95% CI: 2,100 to 2,700) avoided premature deaths under RCP4.5 (Table 1). Most of the counties would see a decrease in all-cause mortality under both RCPs, mainly as a result of lower O₃ levels, while some counties in the West and Northwest would see an increase due to higher O₃ levels (Supplementary Material, Figure S3(a) and (b)). Decreases in the EMs under RCP8.5 during the warm season, unlike year-round mortality, are related to decreases in O₃ in summer. In other seasons, increases in EMs are due to higher O₃ under RCP8.5 (Gao et al. 2013).

County-level, O₃-related EM estimates had high spatial variation (Fig. 3). Under RCP8.5, most counties in the Northwest, West North Central, and West were estimated to have an increased O₃-related EM. Counties with a high population, such as Los Angeles County in California and Cook County in Illinois, showed higher increases in O₃-related EM. The EMs in these counties appeared to increase even under RCP4.5 (Fig. 3(a)), despite the overall reduction in O₃-related EM (Table 1). The spatial distribution for O₃-related EM per 100,000 persons was different from that of the absolute values showing greater spatial heterogeneity. All counties in the Northwest showed higher increases in O₃-related EM per 100,000 persons whereas many counties in the South and Southeast had decreased rates under RCP8.5 (Fig. 3(c)). The standard errors of county-level EM estimates also had high spatial variations as presented in Fig. 3(d) and (e) with Southern California showing the greatest standard errors.

Table 2 shows the percentage contribution of county-level SS from factors in Equation 1 (i.e., CRF coefficient, population, and RCP) to the total SS based on the ANOVA analysis. A greater percentage reflects higher sensitivity of the EM to a specific factor. Overall, O₃-related EM estimates were the most sensitive to the RCP, which account for almost 80% of the total SS.

However, factor contributions also vary widely in space (Fig. 4). Counties in the Northwest, West North Central, West, and Southwest had higher sensitivity to the RCP scenarios (Fig. 4(a)), since the projected O₃ levels in these regions switched from a projected decrease under RCP4.5 to an increase under RCP8.5 (Fig. 2(a) and (b)). In contrast, counties in the Central, South, and Southeast had a higher sensitivity to the population projections.

Figure 5 shows the probability distribution of the estimated, national-level O₃-related EMs (the national EM is the sum of county-level EMs and each scenario has 1,000 national EM estimates based on the county-level Monte Carlo simulations), accounting for MDA8 O₃ and MDA1 O₃, two different CRF coefficients, and four population scenarios. O₃-related EMs ranged from 427 to 2,198 (95% CI) deaths per year, with a mean of 1,312 deaths per year under RCP8.5; and from −3,021 to −1,216 (95% CI) deaths per year, with a mean of −2,118 deaths per year under RCP4.5. The RCPs played a critical role in estimated, O₃-related EM because O₃ levels were found to increase under RCP8.5 in 2057–2059, but decrease under RCP4.5 compared with 2001–2004. The probability distribution of the estimated national EM is not smooth because of discrete population projections, averaging hours, and CRF ranges, but its two distinct modes follow the RCPs (Fig. 5(a)). The combined probability distribution function from the two RCPs predicts a ~50% chance that national O₃-related EM will be greater than the base years, and a ~10% chance that O₃-related EM in the future will be over 1,900 deaths per year more than the base years (Fig. 5(b)).