Expanding number of Western US urban centers face declining summertime air quality due to enhanced wildland fire activity

Combining multiple sources of information on atmospheric composition, wildland fire emissions, and fire area burned, we link decadal air quality trends in Western US urban centers with wildland fire activity during the months of August and September for the years 2000–2019. We find spatially consistent trends in extreme levels (upper quantile) of fine particulate matter (PM2.5), organic carbon, and absorption aerosol optical depth centered on the US Pacific Northwest during the month of August. Emerging trends were also found across the Pacific Northwest, western Montana, and Wyoming in September. Furthermore, we identify potential wildfire emission ‘hotspots’ from trends in wildfire derived PM2.5 emissions and burned area. The spatial correspondence between wildfire emissions hotspots and extreme air quality trends, as well as their concomitant spatial shift from August to September supports the hypothesis that wildfires are driving extreme air quality trends across the Western US. We derive further evidence of the influence of wildland fires on air quality in Western US urban centers from smoke induced PM2.5 enhancements calculated through statistical modeling of the PM2.5-meteorology relationship at 18 Western US cities. Our results highlight the significant risk of increased human exposure to wildfire smoke in August at these Western US population centers, while also pointing to the potential danger of emerging trends in Western US population growth, wildfire emissions, and extreme air quality in September.


Introduction
Atmospheric aerosols with an aerodynamic diameter of less than 2.5 µms (particulate matter (PM 2.5 )) are known to cause adverse health outcomes, such as increased emergency room visits for respiratory illness [1][2][3][4] and enhanced daily mortality [5,6]. In light of these adverse health outcomes, PM 2.5 is treated as a 'criteria pollutant' by the US Environmental Protection agency (EPA) and regulated under the Clean Air Act, with levels set by the National Ambient Air Quality Standards (NAAQS). Currently, NAAQS attainment requires annual mean and 98th percentile PM 2.5 concentrations of less than 12 and 35 µg m −3 , respectively, when averaged over three years.
Despite a 43% decline in national scale average PM 2.5 concentrations over the years 2000-2019 [7], trends toward increasingly extreme air quality have been identified across much of the Western US [8]. Efforts aimed at attribution of these increases in aerosol loading have indicated biomass burning as a potential culprit. Emerging trends in PM 2.5 associated with smoke have been identified in the Pacific Northwest and portions of the surrounding Western states [9]. Positive PM 2.5 trends have been associated with elevated total carbon at sites across the Western US [8], and enhanced surface PM 2.5 and aerosol optical depth [10] have been linked to fire smoke within the Colorado Front Range. Urbanization as the driver of these aerosol trends is unlikely given evidence of decreasing trends in the annual mean of PM 2.5 and organic aerosol concentrations throughout the Western US [7,8,11].
Investigations of wildfire activity point to Western US aerosol trends as the result of a shifting climate, by which increasing temperatures, aridity [12], and a declining snowpack [13] promote more frequent and intense wildfires [14][15][16]. Linkages between aridity, fire area burned (FAB), and aerosol loading in the Intermountain West have been documented [17]. Positive trends in FAB and the number of large wildfires from 1984 to 2011 have been identified [14], and modeling under a climate change scenario of 2 • C warming by 2050 suggests that the increasing trend will persist into the future [18]. The climate change-wildfire-aerosol linkage in the Western US will expose millions of additional people to elevated PM 2.5 annually by mid-century [19], producing a corresponding increase in premature deaths attributable to fire specific PM 2.5 [20]. Furthermore, the effects of climate change will be compounded by the legacy of decades of fire suppression, further amplifying the impact of large fires [21]. The extraordinary fire season of 2020 in the Western US [22] provides vivid testament to the societal costs of this trajectory.
With the potential nexus of increasing aerosol load, enhanced fire activity, and a growing Western US population [23] in mind, we carry out analyses to build upon existing work by: (a) improving current understanding of aerosol trends by combining a multitude of atmospheric composition data-PM 2.5 , organic carbon, and absorption aerosol optical depth (AAOD); organic carbon and AAOD serve as markers for biomass burning given the ability of Western US wildfire activity to explain organic carbon variability [24] and the sensitivity of AAOD to the presence of brown carbon via differential ultraviolet absorption [25]; (b) refining the temporal resolution of aerosol trends to monthly; (c) emphasizing human exposure to degraded air quality in urban centers, where the bulk of the population resides; and (d) synthesizing wildfire emissions and fire burned area with the atmospheric composition data to illustrate potential hotspot regions that are responsible for frequent human exposure to degraded air quality.

Data
The spatial information, temporal information, and usage of datasets included in this study are presented in table 1. Individual dataset descriptions and maps of aerosol observation sites can be found in supplementary information (SI) section S1 (available online at stacks.iop.org/ERL/16/054036/mmedia).

City clustering
As a means of emphasizing human exposure to degraded air quality in our analysis, Ozone Monitoring Instrument (OMI) swath AAOD trends and generalized additive modeling (GAM) were applied to observations of large 'cities' in the Western US.
Given that official city boundaries do not consistently subscribe to population characteristics, these cities were defined via a city clustering procedure applied to population density and population count data from NASA's Socioeconomic Data and Applications Center [29]. This procedure allows a 'city' to be defined beyond its administrative boundaries. When defining cities, the city clustering algorithm iteratively clusters grid cells meeting a threshold population density if they are within some distance (ℓ) of each other [30]. An ℓ value of 5 km was used here, as the literature suggests values between 2 and 6 km provide the best alignment between metropolitan statistical areas and city clusters [31]. A population density threshold of 1000 km −2 was selected via qualitative observations of returned city clusters using population density thresholds between 500 and 1200 km −2 . The combination of a 5 km clustering length and a threshold population density of 1000 km −2 (e.g. figure S3) seemed to provide the best distribution of large cities in the West when requiring the total population of a cluster exceed 200 000 people. In total, 33 large cities within the Western US were defined using this procedure.
The city clustering performed here can be recreated by applying the cca() function from the 'osc' package in R [32] to the referenced datasets. We have found that doing so recreates the clusters developed here with a maximum discrepancy in total city population of approximately 3%. City clusters will be identical for the majority of the 33 cities defined.
When applied to OMI swath AAOD data, these city clusters are used to look for intersections of the OMI swath with city polygons produced from the gridded clusters. Each OMI pixel that intersects a given city is treated as part of that cities AAOD dataset.

Quantile regression (QR)
To investigate the trajectory of extreme aerosol events in the Western US, QR was applied to the EPA PM 2.5 [33], Interagency Monitoring of Protected Visual Environments (IMPROVE) [34-37] PM 2.5 and organic carbon, EPA Chemical Speciation Network (CSN) organic carbon [33,38], and OMI swath AAOD [39] datasets on a month-by-month basis. QR was applied at upper quantiles (98th and 95th) to focus our analysis on episodic extreme aerosol events, consistent with the nature of biomass burning. A reduced quantile relative to surface aerosol datasets was used for AAOD (95th) to combat data limitations from cloud screening.
QR allows trends to be fitted to the extremes of a distribution rather than the mean (linear regression) by taking advantage of asymmetric weighting during regression [40,41]. This asymmetric weighting is in opposition to the symmetric ordinary least squares approach used for a linear regression. An added benefit of QR is that it allows for an understanding of trends in the extremes of a dataset without a need to subset the data. We perform QR using the rq() function from the 'quantreg' package in R [42]. Pvalues and 95% confidence intervals were developed via a bootstrapping procedure supplied by arguments within the summ() function of the 'jtools' package [43]. Essentially, each QR was reapplied to 100 000 resampled versions of the original dataset, with replacement, to see if the regression to the actual dataset was distinguishable at p < 0.05 from the scatter produced by the resampled datasets. Resampled datasets were also used to construct 95% confidence intervals. Specifics of QR application to each dataset are provided in SI section S2.

GAM
GAMs allow for the combination of many linear or non-linear predictor variables to estimate some response variable, and have the following general form [44]: where g is a link function relating the expected value of the response variable (µ i ) to a linear combination of ordinary linear model components (A i θ) and smooth nonlinear functions (f j ) of predictor variables (x k ).
GAMs have previously been used to investigate enhancements in urban ozone resulting from wildfire smoke inundation [45,46]. McClure and Jaffe [45] use a GAM framework to demonstrate significant increases in ozone on smoke impacted days at an EPA site in Boise, Idaho, when accounting for meteorological conditions and atmospheric transport. To achieve this, they fit a GAM to ozone concentrationmeteorology/back trajectory (x k ) relationships for smoke free days and then used the model to estimate ozone concentrations on smoke impacted days (µ i ). Model estimates could be considered the expected ozone concentrations for smoke impacted days based on meteorological and transport conditions alone, meaning the model residuals (observed values-µ i ) offered an estimate of the smoke enhancements to ozone.
Here we apply a GAM framework to 18 EPA PM 2.5 monitors in the Western US to investigate the role of wildfire smoke in anomalous PM 2.5 values. EPA monitors were utilized if they fell within one of the 33 large cities defined by our city clustering procedure and possessed daily data for 2008-2019. We rely on the meteorological datasets indicated in table 1 as our predictor variables (x k ) to estimate the daily averaged PM 2.5 on smoke impacted days (µ i ). We employ an identity link function, removing g() from the regression equation.
For each monitor, data was first divided into 'smoke' and 'non-smoke' days following the methods of McClure and Jaffe [45]. 'Smoke' days were required to meet two criteria: (a) the daily averaged PM 2.5 had to be in excess of one standard deviation above the mean across 2008-2019, and (b) the monitor had to fall within a NOAA Hazard Mapping System (HMS) daily smoke polygon [47]. While HMS smoke polygons have previously been used to demonstrate statistical differences in PM 2.5 concentrations when smoke is present [48], their subjective basis opens the door for false-negative 'non-smoke' designations. HMS smoke polygons are developed through manual analysis of satellite visible imagery, meaning they are subject to the obscuring effects of clouds and haze, and are best considered as a conservative estimate of smoke extent [49].
Once divided into 'smoke' and 'non-smoke' datasets, each GAM was trained on 'non-smoke' data, fitting the model based on PM 2.5 -meteorology relationships, and then applied to 'smoke' data. Training on 'non-smoke' data provided a baseline from which to predict 'smoke' day PM 2.5 from meteorology alone, allowing for an estimate of the PM 2.5 enhancement due to the presence of smoke on these days. Residuals for each dataset, 'smoke' and 'non-smoke' , were analyzed for statistical differences. It should be noted that the conservative nature of the 'smoke' day designation may provide for attribution of smoke driven PM 2.5 enhancements to meteorological variables within the models, resulting in a tempered estimate of typical PM 2.5 enhancements due to smoke for cities in the Western US. False negatives may also be responsible for reduced r-squared values in some GAM setups.
Specifics of GAM construction and crossvalidation within R may be found in SI section S3. Table S1 describes the implementation of individual variables in each of our 18 GAMs.

Identification of potential emissions hotspots
Recognition of potential wildfire emissions hotspots with relevance to human exposure to poor air quality was based on linear trend analyses of monthly aggregated FAB and wildfire emitted PM 2.5 in Western US and Western Canadian ecoregions/ecoprovinces [50,51]. Supporting statistics were generated via bootstrapping as outlined in the QR methods. Given the uncertainty in fire emissions inventories, numerous inventories were considered, including: the Monitoring Trends in Burn Severity (MTBS) [52] FAB dataset, Global Fire Emissions Database (GFED) [53][54][55], Fire Inventory from NCAR (FINN) [56], and Quick Fire Emissions Dataset (QFED) [57]. When possible, trends were calculated for time series ending in 2018 (a particularly active fire year) and 2019 (a rather quiet fire year) as a means of considering the sensitivity of results to emissions in the terminal year. Ecoregions/ecoprovinces demonstrating positive and statistically significant trends (p < 0.1, often p < 0.05) in FAB/PM 2.5 emissions across many inventories and time series were considered potential wildfire emissions hotspots.
FAB within the MTBS was calculated from the associated fire perimeter shapefiles using R. This calculation is to account for unburned islands within the maximum fire extent reported by MTBS. As a caveat, it should be noted that MTBS trends reflect FAB by ignition month, such that the associated emissions may occur in subsequent months. Given this, the relevant MTBS trend for the Central California Foothills and Coastal Mountains ecoregion reported here is for July (table S2).  1). The consistency of spatial patterns across datasets suggests wildfires are driving trends in extreme air quality across the Western US. These findings are particularly true in urban centers, as the EPA monitor sites used for PM 2.5 analyses tend to cluster in urban settings, and our city-clustered AAOD (described within section 2) focuses on locations with elevated population density (>1000 km −2 ) and total population (>200 000). The consistency between nearby EPA (primarily urban) and IMPROVE (primarily rural) PM 2.5 trends counters the hypothesis that these results are driven by urban development. Figures 1(a)-(c) depict the upper quantile trends for the month of August for PM 2.5 observations from the EPA and IMPROVE networks, organic carbon observations from the EPA CSN and IMPROVE networks, and satellite-based AAOD from the OMI swath-level data, respectively. From these panels, it is apparent that the Pacific Northwest and portions of the Central California Valley are experiencing intense degradation of air quality during August. Trends in the Pacific Northwest appear especially robust in terms of spatial consistency, with many of the strongest trends in each dataset found within the states of Washington, Oregon, and Idaho. Further, a linear fit to August mean daily PM 2.5 values (figure 2) indicates that the spatial coherence in trends persists beyond the extremes (98th quantile, figure 1(d)), underscoring the severity  Overlap of positive and statistically significant trends in organic carbon and AAOD datasets with positive and significant trends in PM 2.5 strongly suggests that the PM 2.5 trends may be fire driven. However, it should be noted that areas lacking agreement across datasets may partially arise from spatial under-sampling (figure S1) and/or the complexities associated with smoke transport in complex terrain and the wildfire plume rise. Plume rise may decouple surface (PM 2.5 and organic carbon) and column (AAOD) trends. Further evidence of a wildfire source is drawn from the maxima in trends of PM 2.5 and organic carbon located within or just east of the fire prone Klamath Mountains/California High North Coast Range. As detailed in later sections, the Klamath Mountains/California High North Coast Range ecoregion is characterized by increasing wildfire PM 2.5 emissions across several wildfire emissions inventories.

August
The presented decreasing trends in the Southwest demonstrate spatial consistency with identified trends toward reduced summertime organic aerosol concentrations [11]. Malm et al [11] indicate that significant decreasing summertime organic aerosol trends would prevail across the Western US were it not for outliers sourced from biomass burning.

September
For the month of September the 98th quantile daily PM 2.5 trends ( figure 1(d)) exhibit a shift to the north and east relative to patterns in August ( figure 1(a)), underscoring air quality degradation across portions of the Pacific Northwest, Western Montana, and Wyoming. Attribution of these trends to wildfires via aerosol observations is hindered by their reduced statistical significance and limited overlap with positive and significant organic carbon and AAOD trends, limited to the Seattle area. We suspect that the reduced spatial overlap of statistically significant trends, relative to August, is an artifact of the fire season maximizing in August (figure S11). For this reason, should the trends toward enhanced wildfire activity persist into the future, we expect September air quality trends to exhibit greater clarity given additional years of data. As an alternative hypothesis, increasing trends (2000-2014, p > 0.1) in average fine mode (PM 2.5 ) dust concentrations during the fall (September, October, November) suggest the potential role of dust in September aerosol trends at sites in Northwestern Montana and Western Oregon [58]. It should be noted that the reduced spatial consistency of PM 2.5 trends relative to August results, particularly in and around the Rocky Mountains, is largely the result of variable EPA/IMPROVE monitor site density (figure S1).

Smoke-driven PM 2.5 enhancements in urban centers
Clear enhancements in PM 2.5 levels on smoke impacted days in Western US urban centers are identified on the basis of GAM [37]. Mean PM 2.5 enhancements (figure 3) are calculated as the mean model residual (observed-GAM predicted PM 2.5 ) for this nonlinear statistical approach (see section 2). As such, model residuals reflect the GAMs' underprediction on 'smoke' days when only meteorological variables are used for prediction. Thus, residuals from the statistical model fitted through GAM (table S1) provide a way to quantify the additional PM 2.5 enhancement from wildfires after controlling for meteorology. The results of trend analyses for individual emissions inventories may be found in figures S12-S15 in SI section 6. The results of trend analyses for April-July and October-December are presented in figures S16-S28 in SI section 6.
Though applied to data for all months, our GAM results are generally in line with results from trend analyses of upper quantile air quality measures for the months of August and September. Much like the August trends, the largest mean PM 2.5 enhancements on smoke-impacted days can be found in the Pacific Northwest and California's Central Valley. Mean smoke-impacted PM 2.5 enhancements in these regions range from 18 µg m −3 to 33 µg m −3 at the Spokane site. Maximum enhancements similarly highlight the Pacific Northwest, peaking at 199 µg m −3 at Spokane.
While interpretation and comparison between GAM results are complicated by differing model forms (table S1), variable strength of model fits, and site-specific air quality issues (SI section S5), indicators of model performance highlight general patterns within the Pacific Northwest as being of greatest confidence. GAMs developed for urban centers in Washington and Oregon typically produced smaller mean squared normalized 'non-smoke' test residuals, possessed larger r-squared values, and were applied to a greater number of smoke impacted days, relative to GAMs developed for sites in other states. The mean and standard deviation of GAM r-squared values for sites in the Pacific Northwest were 0.61 and 0.10, respectively, as compared to 0.54 and 0.10 when considering all 18 sites. These attributes, paired with consistently larger mean enhancements relative to other regions, support the idea that urban centers in the Pacific Northwest are being subjected to especially pronounced degradation of air quality as a result of wildfire smoke.
Despite uncertainty in the exact magnitudes of GAM produced smoke driven PM 2.5 enhancements, the conservative nature, stemming from the 'smoke' day criterion (see section 2), and general sign of this analysis support the idea that wildfires are contributing to human exposure to extreme air quality in urban centers across West.

Potential wildfire emissions hotspots
Trends in monthly aggregated FAB and wildfire emitted PM 2.5 (section 2) across multiple fire emissions inventories produce a picture of wildfire activity (figure 4) that reflects the spatiotemporal evolution of trends in extreme air quality (figure 1), supporting the notion that wildfires are driving trends in extreme air quality across the Western US. Statistically significant and increasing trends in PM 2.5 and wildfire activity shift from the Pacific Northwest and the Central California Valley toward the Rocky Mountains from August to September. In August, the Pacific Northwest is uniquely situated between potential wildfire emissions hotspots in British Columbia and the Pacific Northwest/California, while also demonstrating the most extreme trends in upper quantile air quality. Progressing into September, trends in wildfire activity relax in British Columbia and Washington, while reflecting potential emissions hotspots in Northern California/Southern Oregon and the Rocky Mountains. Similarly, September air quality trends show an attenuation relative to August in the Pacific Northwest while increasing in extent near the Rockies.
In terms of human exposure to wildfire degraded air quality, this setup highlights the colocation of August trends and Western US population centers as being of particular concern (figure 4), while also pointing to the danger of concurrent trends in Western US population growth, wildfire emissions, and extreme air quality in September. Potential explanations for the spatial differences in the seasonality of wildfire trends arise from the complex interconnections between drought and aridity, surface and canopy fuels, temporal variability of wildfire-climate relationships, and variability of historic forest management practices, across local to regional scales [60].
In the context of the wildland-urban interface (WUI), where structures intermingle with wildland vegetation and human ignitions of wildfires are common [61], the risk posed by this setup is bolstered. Recent decades are characterized by an expansion of WUI across the Western US, with growth rates by area >75% in the Northern Rockies and portions of many Western states [61]. Further, the Western US counties with the greatest potential for WUI expansion (by area) are in Southwestern Oregon and Northern California, as well as the Northern Rockies, encompassing a suspected wildfire emissions hotspot and existing in relatively close proximity to Pacific Northwest urban centers, respectively [62]. The potential for WUI expansion to exacerbate the trends identified here is cause for concern with regards to the trajectory of regional air quality in the Western US. However, each region identified in figure 4(a) has a significant trend in one or more inventories through 2019. In light of the record-breaking 2020 fire season to date [22], it is suspected that extension of these analyses will produce a rebound in trend magnitudes and statistical significance.

September
Reminiscent of the spatiotemporal evolution of PM 2.5 , trends in FAB and fire-emitted PM 2.5 for the month of September depict the emergence of a possible hotspot for wildfire emissions in the Rocky Mountains while emission trends seem to flatten in the Pacific Northwest/California and British Columbia ( figure 4(b), table S2). September trends in FAB/wildfire emissions present overall weaker magnitudes with reduced statistical significance relative to August, much like the accompanying air quality trends.
A Rocky Mountains wildfire emission hotspot in September is supported by statistically significant trends in emitted PM 2.5 for timeseries ending in 2018 and 2019 in both the GFED and QFED inventories, as well as FAB in the Southern Rockies in the MTBS dataset. Dampening of potential hotspot activity in the Pacific Northwest/California is captured by statistically significant trends being limited to the Klamath Mountains/California High North Coast Range ecoregion, which garners statistical support from the FINN (2018 timeseries only), GFED, and QFED inventories. The magnitude of GFED trends through 2018 are of order 10 5 -10 6 kg PM 2.5 /year for each of these three ecoregions.
While Canadian contributions to wildfireemitted PM 2.5 trends in the month of September appear uncertain in our analysis, a possible explanation for the attenuation of emissions trends in British Columbia is provided by the southward progression of the polar jet/mid-latitude storm tracks during the transition from summer to fall. All trends pertaining to ecoregions/ecoprovinces depicted in figure 4 are provided in table S2.

Conclusions
We find that summertime Western US air quality is declining while wildland fire activity is increasing. In August, EPA and IMPROVE monitor sites across the Western US indicate positive trends in 98th quantile PM 2.5 , with the Pacific Northwest presenting particularly robust trends. Concurrent with August PM 2.5 trends, wildland fire emissions and FAB are increasing across portions of British Columbia, Washington, Oregon, and California, uniquely placing the Pacific Northwest between two potential emissions 'hotspots' . The notion that these August PM 2.5 trends are the result of spatially linked wildfire emissions trends is further supported by analyses of organic carbon, AAOD, and smoke-driven PM 2.5 enhancements, each of which points to the Pacific Northwest and portions of California as being impacted by wildfire smoke and extreme air quality.
In September, the spatial coverage of positive 98th quantile PM 2.5 trends shifts north and east relative to August, demonstrating clear enhancements in Western Montana and Wyoming. At the same time, the spatial distribution of wildland fire emissions trends is reconfigured such that positive trends are apparent in the Middle and Southern Rockies while trends relax in British Columbia and Washington. The correspondence in spatial shifts between wildfire emissions hotspots and extreme air quality trends from August to September provides further support for the hypothesis that wildfires are driving extreme air quality trends across the west. While we acknowledge that air quality and wildland fire trends for the month of September are less robust than those of August, we suspect the trends could emerge with greater clarity given additional data.
Thus far, the 2020 fire season has provided an exceptionally clear indication of the cost to human wellbeing should the trends we have identified persist over the coming decades. Extrapolation of the PM 2.5 trends we have identified 15 years into the future suggests that many cities in the Western US may struggle to meet NAAQS within the next few decades. Spokane, Washington provides a particularly striking example of this concern, as PM 2.5 trends at this city indicate a >16 µg m −3 (1.26-32.35 µg m −3 , 95% confidence interval) increase in the mean and >34 µg m −3 (15.01-53.85 µg m −3 , 95% confidence interval) increase in the 98th quantile of daily averaged PM 2.5 for the month of August by 2035. In terms of health outcomes, it has been documented that wildfire derived PM 2.5 enhancements in excess of 37 µg m −3 have been associated with a 7.2% increase in the risk of respiratory hospital admission across all ages [2].
While our statistical findings and previous work on the air quality-wildfire connection in the Western US [8][9][10]17], are highly suggestive, it is apparent that a more sophisticated, atmospheric model-based quantitative attribution of air quality trends to wildfire sources is needed. A quantitative attribution that identifies 'hotspots' for wildfire emissions with relevance to Western US population centers may provide a means for mitigation via targeted fuel treatments and related forest management practices.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).