Characterizing emissions, chemistry, and health impacts of aged wildfire smoke in a western US city

We report hourly surface observations of PM_2.5, CO, NO_x, O₃, and 75 speciated VOCs in Missoula, Montana, during a strong smoke event in 2020. This study tests our current understanding of wildfire emissions, chemistry, and health effects as implemented in the GEOS-Chem chemical transport model. Three-or-more-day-old smoke transported from California and the Pacific Northwest increased CO, PM_2.5, and total measured VOCs by factors of 2–8, with hourly maxima of 800 ppb, 120 µg m⁻³, and 85 ppb, respectively. In contrast, NO_x levels were not elevated compared to the urban background. O₃ showed a non-monotonic response to wildfire smoke: MDA8  O₃ increased under light smoke but flattened or declined when PM_2.5 exceeded ∼ 30–40 µg m⁻³, a feature that GEOS-Chem failed to reproduce. A 2020-style wildfire season recurring annually would yield an excess lifetime cancer risk of 100-in-1 million or approximately 7 times the non-smoke baseline. The chronic non-cancer hazard index (HI) would reach 3.0, indicating appreciable potential for chronic non-cancer effects. About 90 % of cancer risks are from PM_2.5 whereas non-cancer risks are dominated by formaldehyde, benzene, acrolein, and acetaldehyde. GEOS-Chem captured major smoke intrusions but underestimated CO, PM_2.5, and VOCs by 30 %–90 %. These model biases propagate to health metrics, with GEOS-Chem underestimating smoke-attributable cancer risk by ∼ 40 % and chronic HI by ∼ 10 times. We attribute the model errors to underpredicted fire emissions and unrepresented VOC chemistry, which together led to an overestimation of OH and insufficient secondary production.

1 Introduction

Wildfires are the second largest source of volatile organic compounds (VOCs) in the western US (Hoesly et al., 2018; Jin et al., 2023). Some VOCs are designated hazardous air pollutants (HAPs) due to their carcinogenicity and potential for acute or chronic health effects (Naeher et al., 2007; Reid et al., 2016). Lengthening and intensified fire seasons in recent decades have posed direct air quality risks to rural populations and the inhabitants of ∼ 50 million US homes that already lie within the fire-prone wildland–urban interface (WUI), a figure expected to grow by one million every three years (Burke et al., 2021). Moreover, recent years have shown that western US and Canadian biomass burning (BB) smoke frequently impacts cities in western North America and can reach population centers thousands of kilometers downwind on the East Coast (Yu et al., 2024). Despite the growing threat of more frequent and intense wildfires, the health effects of HAPs emitted from wildfires are not well documented, and residents often lack chemically resolved data on the pollutants to which they are exposed, because of limited ground-level measurements. Here, we combine hourly ground-based observations of VOCs and criteria pollutants with the GEOS-Chem chemical transport model (CTM) to characterize the composition, chemical evolution, and health implications of aged wildfire smoke.

Past observational efforts have provided valuable yet incomplete insights into the health risks associated with wildfire smoke exposure (Gould et al., 2024). For example, some studies have targeted near-source fireline PM (Navarro et al., 2021), CO (Semmens et al., 2021), or high-altitude plumes sampled by aircraft (O'Dell et al., 2020), and others have examined speciated VOCs and/or aerosol composition (Akagi et al., 2014; Fiddler et al., 2024; Joo et al., 2024; Liang et al., 2022). Speciated, near-surface data where the majority of human exposure to fresh and aged smoke actually occurs are still seriously under sampled. The lack of sufficient measurements also limits the evaluation of CTMs, which are widely used for smoke forecasting and health-risk assessments. Recent aircraft observations show that GEOS-Chem (a commonly used CTM) underestimates BB CO and key HAPs such as formaldehyde by a factor of three or more (Jin et al., 2023). However, analogous ground-based validation remains absent. Because near-surface chemistry differs fundamentally from lofted plumes (e.g., ambient temperature, nighttime oxidants, boundary layer mixing) (Decker et al., 2019; Pagonis et al., 2023), ground observations are urgently needed to improve model validation and exposure assessments.

The western US fire season of 2020 produced some of the most extreme smoke levels on record in OR, WA, and CA (Albores et al., 2023; Reilly et al., 2022), with September standing out as the peak month (Abatzoglou et al., 2021; Mass et al., 2021). More than 1.6 million ha burned in California alone, releasing an estimated ∼ 127 Tg CO₂-equivalent, nearly twice the state's cumulative 2003–2019 greenhouse-gas reductions (Jerrett et al., 2022). These BB events provide an ideal natural laboratory for quantifying the composition of aged smoke, and testing CTM capabilities for representing regional smoke and assessing the continental-scale impacts of wildfires.

In this study, we present hourly ground-based measurements of 75 VOCs (including 15 HAPs), PM_2.5, CO, NO_x, and O₃ collected in aged (mostly > 3 d) regional wildfire smoke during September 2020 in Missoula, Montana – a representative northwestern US city frequently impacted by regional wildfire smoke. Leveraging this comprehensive dataset, we (i) characterize the temporal evolution of VOCs and criteria pollutants (Sect. 4), (ii) quantify species-specific BB enhancements (Sect. 5), (iii) investigate O₃–PM_2.5 relationships (Sect. 6), (iv) assess public health risks using regulatory exposure metrics for PM_2.5 and HAPs (Sect. 7), and (v) evaluate GEOS-Chem simulations against observations to diagnose model biases in BB emissions and chemistry, and their implications for health-risk estimates (Sect. 8).

2 Methods

2.1 Missoula Science Observatory

Missoula, MT (46.8721° N, 113.9940° W) lies in a mountain valley in the northern Rocky Mountains with a population of approximately 120 000. Missoula has been studied in the past for its frequent impacts from both local (Montana and Idaho) and regional wildfires, including long-range smoke transport from the West Coast and the Pacific Northwest, such as California, Oregon, Washington, and British Columbia (Selimovic et al., 2019, 2020). Its frequent exposure to surface smoke makes it a representative site for evaluating human exposure and air quality under real-world wildfire conditions.

Long-term air quality monitoring at the Missoula Science Observatory (MSO) began in 2017 with measurements including four criteria pollutants (PM_2.5, CO, NO₂, and O₃) and aerosol optical properties (Selimovic et al., 2019, 2020). In 2020, the monitoring scope expanded to include 75 individual VOCs, of which 15 are classified as HAPs by the U.S. EPA. The University of Montana (UM) campus serves as the primary site for CO, NO, NO₂, O₃, and VOC measurements. The NO, NO₂, and O₃ inlet was located 12.5 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$ at the Charles H. Clapp Building; the CO and VOCs inlet was initially located 10 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$ on the Chemistry Building, ∼ 70 m away.

Measurements of PM_2.5 were obtained from the Montana Department of Environmental Quality site at Boyd Park, roughly 3 km southwest of the UM campus. Despite the spatial separation between measurement sites, previous studies have validated the temporal and spatial agreement in pollutant concentrations between the two sites for 2017–2019 summers (Selimovic et al., 2019, 2020); we confirm similar agreement among PM_2.5 measured at the Boyd Park and other tracers measured at the UM campus for this study (R² > 0.9), supporting the use of combined datasets in this analysis.

Hourly meteorological data were obtained from the MesoWest database and accessed via the Synoptic Data API (https://synopticdata.com/, last access: 30 July 2026). Incoming shortwave radiation (W m⁻²) was obtained from the Blue Mountain site (station ID: BLMM8; 46.73° N, 114.09° W), while air temperature (°C), relative humidity (%), and precipitation accumulation (mm) were obtained from the Missoula Valley site (station ID: E0591; 46.86° N, 114.02° W). These stations were selected as the closest available meteorological observations to the MSO (10 and 5 km away, respectively). We used broadband shortwave (SW) as a first-order proxy for actinic flux and photolysis as they were not measured; we note that the SW (200–4000 nm) can diverge from the near-UV range relevant to J(NO₂) and J(O¹D) (∼ 300–420 nm) under smoke. Planetary boundary layer (PBL) mixing height (m) was obtained from the High-Resolution Rapid Refresh (HRRR, 3 km), an operational forecast and analysis system developed by NOAA with hourly updates at 3 km horizontal resolution (Dowell et al., 2022). We also obtained the mixing height data from the NASA Goddard Earth Observing System Forward Processing system (GEOS-FP, 0.25° × 0.3125°, hourly) as it serves as one of the inputs for GEOS-Chem.

Figure 1 illustrates the climatology of PM_2.5, defined here as the multi-year statistical distribution of daily-mean PM_2.5 during the wildfire season (June–October) over 2010–2024 in Missoula. Despite large interannual variability, a clear seasonal progression is evident. The multi-year daily average PM_2.5 concentrations remain below 15 µg m⁻³ throughout June, increase steadily during July, and peak from mid-August to occasionally mid-September, coincident with increases in wildfire smoke episodes. Peak hourly PM_2.5 even reached ∼ 470 µg m⁻³ in 2017. Reflecting this chronic smoke burden, in 2024, Missoula ranked 14th among 223 US metropolitan areas for 24 h particle pollution, measured by the weighted annual average number of unhealthy 24 h PM_2.5, and 29th for annual particle pollution, measured by the annual average PM_2.5 concentration (American Lung Association, 2024). In this study, we focus specifically on September 2020, when regional smoke transport led to one of the most prolonged and chemically distinct episodes of aged wildfire smoke observed at the surface (blue line in Fig. 1; Sect. 4).

Figure 1 Wildfire smoke over the western US and climatological PM_2.5 in Missoula, Montana. (A) Cropped GOES-17 GeoColor satellite image (15 September 2020) over the Western United States, with the location of Missoula, Montana indicated by the red dot. (B) Daily-mean PM_2.5 concentrations at Missoula during the 2020 fire season (blue) compared with the 2010–2024 climatology (black line; grey box plots give the interquartile range and whisker range (± 1.5×IQR; vertical lines) for each calendar day). Dashed horizontal lines indicate 15 and 35 µg m⁻³, as the WHO 24 h guideline and the U.S. EPA 24 h standard, respectively.

2.2 VOC measurements

Ambient mixing ratios of 75 VOCs were measured as 2 min averages using a custom-built Proton Transfer Reaction Time-of-Flight Mass Spectrometer (PTR-ToF-MS; PTR-ToF-4000, Ionicon Analytik GmbH, Innsbruck, Austria). Instrument operation followed protocols developed in our previous campaigns (Cope et al., 2024; Permar et al., 2021; Selimovic et al., 2022), with drift tube conditions set to $E/N$ = 130 Td, 3.00 mbar, 60 °C, and 800 V.

Air was continuously sampled from an inlet ∼ 10 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$ through ∼ 10 m of heated (60 °C) 6.35 mm ($\mathrm{1}/\mathrm{4}$ in.) outer diameter (O.D.) PTFE tubing. The PTR-ToF-MS then subsampled via ∼ 1 m of 1.59 mm ($\mathrm{1}/\mathrm{16}$ in.) O.D. PEEK tubing, also maintained at 60 °C. The sampling inlet featured a 2 µm PTFE filter to prevent particle intrusion. Filters were replaced every 2 weeks, or more frequently during high pollution periods. Instrument backgrounds were quantified every 2.5 h by sampling VOC-free air generated via a heated platinum catalyst (375 °C; Sigma-Aldrich, 1 wt %).

In this work, we included 75 VOC species measured by PTR-ToF-MS in the analysis. Calibrations were conducted every other day using two compressed gas standard cylinders for 25 individual VOCs (± 5 % at ∼ 1 ppmv; Apel-Riemer Environmental, Inc.). Among them, isomers calibrated at the same $m/z$ (i.e., methyl vinyl ketone and methacrolein at $m/z$ 71.049; ethylbenzene and o-xylene at $m/z$ 107.086; 1,2,4- and 1,3,5-trimethylbenzene at $m/z$ 121.101), which used a weighted average sensitivity based on the corresponding isomeric contributions (Permar et al., 2021). A third standard gas cylinder containing 10 species was calibrated after the campaign (± 5 % at ∼ 1 ppmv; Apel-Riemer Environmental, Inc.). The overall measurement uncertainty for these species is estimated to be less than 15 %. Formaldehyde was calibrated after the campaign using a gas standard cylinder (stated accuracy 5 % at ∼ 2 ppm; Airgas USA LLC, Plumsteadville, PA, USA). Gases were mixed in a Liquid Calibration Unit (LCU, Ionicon Analytik GmbH, Innsbruck, Austria) to derive the dependence of instrument sensitivity on changing humidity. Formic and acetic acids were calibrated after the campaign using LCU, which also applied humidity dependence. Uncertainty for these species is estimated at ∼ 30 % (Permar et al., 2021). For 40 uncalibrated species, instrument sensitivities were estimated using theoretical methods (Sekimoto et al., 2017) refined by field comparisons (Permar et al., 2021), yielding overall uncertainties of ∼ 50 %. Sensitivity of all compounds used in this study ranges from 0.9 to 14 ncps ppbv⁻¹, where ncps denotes normalized counts per second. The background signal of all compounds ranges from 0.05 to 10.1 ncps.

2.3 Measurements of CO, NO_x, O₃, and PM_2.5

CO was measured every ∼ 3 min using a Reducing Compound Photometer (Peak Performer 1, PEAK Laboratories, Mountain View, CA, USA), which separates CO by gas chromatography and detects it via photometric absorption of mercury vapor produced from CO reduction of HgO. The instrument was calibrated weekly via CO standard in a compressed standard gas cylinder, with a detection limit of 0.3 ppb. The instrument shared the inlet with PTR-ToF-MS and sampled ambient air at 40 sccm through 5 m of 3.175 mm ($\mathrm{1}/\mathrm{8}$ in.) O.D. PTFE tubing.

NO₂ and NO were measured every 1 min using a model 405 nm NO_x monitor (2B Technologies, Boulder, CO, USA), which directly quantifies NO₂ by measuring optical extinction at 405 nm, where NO₂ strongly absorbs. NO is measured by adding excess ozone into the optical cell, which quantitatively converts NO to NO₂. The instrument was calibrated by the manufacturer with a detection limit of 1 ppb for NO and NO₂. Ambient air was drawn at ∼ 1.5 L min⁻¹ through ∼ 1 m of 6.35 mm ($\mathrm{1}/\mathrm{4}$ in.) O.D. tubing and passed through a 47 mm, 5 µm PTFE filter (Savillex), which was replaced biweekly or upon visible loading.

PM_2.5 and O₃ were obtained from the hourly measurements at the Montana DEQ site at Boyd Park. PM_2.5 measurements were conducted using a Met One BAM-1020 with VSCC inlet (Volumetric Size-Selective Cyclone), which determines particle mass by measuring the attenuation of beta radiation through the filter tape as particles accumulate. The instrument has a detection limit of 5 µg m⁻³ for hourly measurements. O₃ was measured using an ultraviolet photometric analyzer (Thermo Scientific Model 49i UV photometric O₃ analyzer), which quantifies absorption of UV radiation at 254 nm, and with a detection limit of 5 ppb. We note that the Thermo 49i instrument can exhibit positive interferences in wildfire smoke (Bernays et al., 2022; Long et al., 2021). This interference would not affect our conclusion and would, if anything, make the inferred O₃ suppression and model ozone biases more conservative.

2.4 GEOS-Chem chemical transport model

We use the nested-grid GEOS-Chem CTM over North America (version 13.3.0; 10–70° N, and 140–60° W) to interpret Missoula ground-based measurements. The model is driven by the NASA GEOS-FP meteorological data, with 0.25° × 0.3125° horizontal resolution (∼ 25 km × 30 km; latitude × longitude) and 47 vertical layers extending up to 0.01 hPa (Kim et al., 2015; Wang et al., 2004). Boundary conditions were updated every 3 h from a 4° × 5° global simulation. Transport and convection time steps were 5 min; emission and chemistry time steps were 10 min. Model spin-up included a 1-year global simulation followed by a 1-week nested simulation to reduce the influence of initial conditions prior to 1 September 2020.

The chemical mechanism included detailed HO_x-NO_x-VOC-ozone-halogen-aerosol chemistry with fully coupled troposphere and stratosphere (Eastham et al., 2014; Mao et al., 2010; Park et al., 2004; Schmidt et al., 2016). Dry deposition used a resistance-in-series approach (Wesely, 1989), and wet deposition included scavenging of soluble tracers in convective updrafts, as well as rainout and washout of soluble tracers (Liu et al., 2001). GEOS-Chem uses the non-local scheme for PBL mixing (Lin and McElroy, 2010). As a sensitivity test, we also conducted simulations using the “full mixing” scheme, assuming instantaneous vertical mixing of tracers evenly through the mixing depth, which did not affect the conclusions we derived.

Emissions were computed using the HEMCO module (Keller et al., 2014). These include biogenic VOC emissions from MEGANv2.1 (Guenther et al., 2012) as implemented in GEOS-Chem (Hu et al., 2015). Anthropogenic emissions are from the CEDS global emission inventory, overwritten with the 2011 EPA NEI inventory for the US (Hoesly et al., 2018). Daily BB emissions were taken from the Global Fire Assimilation System (GFAS) version 1.2, chosen for its relatively better VOC performance when compared to other commonly used BB inventories over the western US (Jin et al., 2023; van der Werf et al., 2017). Following previous work (Jin et al., 2023), we expanded the default GFAS VOC speciation by adding lumped ≥ C₃ aldehydes (RCHO), MEK, formic acid, and acetic acid by scaling CO BB emissions with corresponding emission ratios reported from temperate forest burns (1.01, 0.73, 9.5, and 8.61 ppb ppm⁻¹, respectively) (Permar et al., 2021). All fire emissions followed a climatological diurnal profile, which allocates ∼ 85 % of daily fire emissions to the afternoon (local time) (Western Regional Air Partnership, 2005).

Biomass burning injection heights are prescribed with the daily 0.1° × 0.1° mean altitude of maximum injection (“mami”) product derived from a satellite-constrained plume rise model (Freitas et al., 2007; Latham, 1994; Rémy et al., 2017). For each GEOS-Chem grid cell, we computed emission-flux-weighted “mami” values to correct the grid-dependence inherent in the standard model (Jin et al., 2023). BB emissions were then distributed evenly from the surface up to this “mami”.

Overall, four nested simulations were conducted for September 2020: (i) a base run using default GFAS emissions and non-local PBL mixing scheme, (ii) a sensitivity run as the base but without fire emissions (noBB), (iii) a second sensitivity run as the base but with GFAS VOC and CO emissions tripled (3 × BB), and (iv) a third sensitivity run using GFAS and the full instantaneous mixing within the PBL mixing height.

3 Identification of smoke-impacted events

Wildfire smoke at ground level is often detected using particle-based diagnostics, for example: (i) satellite-derived aerosol products such as the overhead information from the NOAA Hazard Mapping System (HMS) (Jaffe et al., 2022; O'Dell et al., 2021), (ii) column-integrated aerosol optical depth from AERONET, (iii) sustained elevations in surface PM_2.5 relative to the general urban background (Selimovic et al., 2019, 2020), and (iv) combinations of these metrics (Kaulfus et al., 2017; McClure and Jaffe, 2018). Many studies also incorporate gas-phase tracers – initially CO, hydrogen cyanide (HCN), acetonitrile (ACN), and more recently, furan and maleic anhydride (MA) (Coggon et al., 2019; de Gouw et al., 2003, 2006; Li et al., 2000). Among these VOCs, furan is a sensitive indicator of fresh BB smoke due to its short lifetime and smaller emission amount from anthropogenic sources. In addition, MA is an oxidation product of furan and its derivatives, with a longer lifetime; thus, it has been proposed as a tracer of aged smoke (Coggon et al., 2019).

Here, we evaluate some current BB-impacted identification metrics in the context of our MSO case study and examine whether our VOC measurements provide additional value in identifying BB events in the urban setting. We first derive the hourly median diurnal urban background for BB tracers and other species using HMS “no-fire” days for September 2020. Any enhancement (ΔX) > 1.5 times the background is considered BB-impacted, and a provisional smoke day is defined when at least six such hours occur. We also tested the ${\mathrm{PM}}_{\mathrm{2.5}}/\mathrm{CO}$ criterion (> 30 $\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{ppm}}^{-\mathrm{1}}$) proposed for urban smoke identification (Jaffe et al., 2022). However, at MSO this threshold classified nearly all days as smoke-impacted, consistent with the limitation that ${\mathrm{PM}}_{\mathrm{2.5}}/\mathrm{CO}$ becomes non-discriminating when background PM_2.5 and CO are low (Jaffe et al., 2022). We therefore do not use $\mathrm{\Delta }{\mathrm{PM}}_{\mathrm{2.5}}/\mathrm{\Delta }\mathrm{CO}$ as a standalone classifier for this site.

Figure 2 summarizes the various daily BB classifications for September 2020 in Missoula. HMS alone identifies 20 smoke days (the greatest number among the diagnostics), but it can flag a day as smoke-impacted even when surface concentrations remain low if smoke is primarily aloft. For example, HMS triggers on 2 and 30 September (daily average PM_2.5 is 8.5 and 5.5 µg m⁻³, respectively), which likely reflect elevated plumes that were not mixed down to the surface (Liu et al., 2024). By comparison, the ΔPM_2.5 criterion identifies 16 d. Its diagnostic power depends strongly on aerosol loading at the surface.

Figure 2 Wildfire smoke identification in Missoula, Montana, during September 2020. Rows show the six individual smoke metrics evaluated in this study. A filled marker indicates that HMS smoke was detected, or that the enhancement of a biomass burning (BB) indicator (ΔX) exceeded 1.5× its diurnal background for at least 6 h on that calendar day. Gray bars (right axis) show the daily mean PM_2.5 concentration.

When daily mean PM_2.5 exceeded 35 µg m⁻³ (n = 7; 13–19 September), the ΔPM_2.5 approach aligned with HMS and with all gas-phase tracers considered here (i.e., CO, ACN, furan, and MA), indicating concentrated smoke that coupled large aerosol mass with gaseous enhancements. A k-means clustering analysis (k = 3: smoke, no-smoke, uncertain) applied to pairwise combinations of CO, PM_2.5, and each tracer corroborated this classification, reinforcing that at high aerosol loadings, routine ground-level PM_2.5 alone can reliably identify smoke during summer in Missoula. We therefore classify these seven days as a high-confidence BB period.

At intermediate daily mean PM_2.5 between 20 and 35 µg m⁻³ (n = 5; 5–6, 11–12, and 23 September), ΔPM_2.5 still overlaps at least two gas-phase tracers but not all, signaling moderately aged or dispersed smoke and therefore medium diagnostic confidence. When the daily-mean PM_2.5 falls below 20 µg m⁻³ (n = 4; 7, 21–22, and 24), PM_2.5 overlaps with at most one gas tracer, and HMS does not indicate smoke on 24 September, so these low-load cases warrant only low-to-medium confidence. One additional case (21 September) meets the PM_2.5 threshold but lacks corroboration from the gas-phase tracers, suggesting a potential false positive, although HMS indicates overhead smoke. Across September (a total of 30 d), the remaining 14 d do not meet the ΔPM_2.5 criterion and are classified as non-smoke by this approach.

Gas-phase tracers extend the detectability of wildfire smoke beyond what is captured by HMS and PM_2.5 alone, particularly during precipitation events. For example, on 20 September, CO and furan enhancements were observed in the absence of corroborating HMS smoke or elevated PM_2.5. A rain event from the preceding day into the afternoon of 20 September efficiently scavenged aerosol mass, resulting in a sharp decrease of ∼ 30 µg m⁻³ in PM_2.5, whereas gas-phase species were less affected and their enhancements persisted in the boundary layer.

Because no single tracer is definitive, we adopt a balanced multi-indicator rule: an hour is classified as smoke-impacted when at least three independent smoke indicators support smoke influence, based on HMS smoke detection and/or BB-associated species exceeding 1.5 × their background values for at least 6 h on that day. Applied to September 2020, this criterion identified 14 smoke days (5–7, 11–19, and 22–23 September). For subsequent sections, we describe the three smoke events, and we focus on the high-confidence BB period of 13–19 September to analyze its enhancement ratios, thereby minimizing uncertainties from marginal or non-BB sources.

4 Overview of BB smoke in September 2020

Wildfires across the Pacific Northwest and northern California were exceptionally active in 2020 (Albores et al., 2023; Higuera and Abatzoglou, 2021; Neyra-Nazarrett et al., 2025; Reilly et al., 2022), and September delivered strong smoke signatures in Missoula. Applying the three-tracer classification described in Sect. 3, we identified 16 smoke-impacted days across confidence levels.

Figures 3 and S1 in the Supplement show three distinct multi-day smoke episodes affecting Missoula identified using the criteria described in Sect. 3, including 5–7 September (Event 1), 11–20 September (Event 2), and 22–24 September (Event 3). Because each episode involved multiple overlapping plumes from different fires with varying transport times, we refer to them as “events” rather than discrete “plumes”. Analysis of GOES-17 and VIIRS satellite imagery, and synoptic meteorological conditions indicates that these events were predominantly driven by southwesterly transport of wildfire smoke from Oregon, Idaho, and California, with likely influence from local smaller fires.

Figure 3 Time series of selected VOCs and criteria pollutants in Missoula, Montana, during September 2020. Time is shown in local time (MDT, UTC−6). (A) Hourly measurements of PM_2.5 (black; left axis), carbon monoxide (CO; green, right axis), total measured VOCs (blue, right axis), acetonitrile (red, right axis), and maleic anhydride (orange, right axis). (B) Concurrent hourly NO₂ (black; left axis), furan (purple; right axis), and ozone (orange, right axis). Gray shading indicates smoke days based on the criteria in Sect. 3.

Event 1 (5–7 September) resulted from the California Creek Fire (United States Forest Service, 2020), which started on 4 September mixing with a localized, narrow plume moving east-northeast over the Bitterroot Mountains toward Missoula, as captured by GOES-17 imagery. Correspondingly, hourly PM_2.5 increased rapidly from < 10 to 30 µg m⁻³ by 06:00 Mountain Daylight Time (MDT, UTC−6) on 6 September and CO from 100 to ∼ 300 ppb by 09:00 MDT the same day, with brief spikes (≤ 50 ppt) of the fresh-smoke tracer furan at 01:00 MDT. Pollutant concentrations returned to the urban background within 24 h after a post-frontal north-westerly wind flushed the valley.

Event 2 (11–20 September) was the longest and strongest impact on Missoula from the Oregon–California “megafire” complex that ignited in August and erupted on 7 September (Abatzoglou et al., 2021). Intense smoke on 8 September was advected westward offshore, where it became entrained into a cut-off low over the NE Pacific (around 42° N, 135° W), and subsequently advected eastward inland, mixing with smoke freshly emitted by the megafire complex and fires in Idaho. By 11 September, dense smoke formed a clearly visible “smoke river” that stretched across Washington, Oregon, and Idaho into western Montana. Back-trajectory analyses indicate that air parcels reaching Missoula had aged for at least two days. Consistently, smoke-impacted hours during this event exhibited pollutant enhancements of four- to fifteen-fold for CO, PM_2.5, acetonitrile, and total measured VOCs (TVOCs), reaching peak values of 800 ppb CO, 120 µg m⁻³ PM_2.5, 2 ppb acetonitrile, and 85 ppb TVOCs in mid-September. Strong correlation (r²=0.9; Fig. S2) among PM_2.5 measured at the Boyd Park and other tracers measured at the UM campus indicates spatially coherent smoke influence across Missoula at the surface. Brief furan spikes (≤ 0.7 ppb) on 12–13 and 18–19 September suggest additional influence from local fires. Overall, this event accounts for ∼ 80 % of smoke hours in this month and is the focus of the enhancement ratio analysis (Sect. 5).

Event 3 (22–24 September) was a weaker resurgence of smoke from the same Oregon–California fires in Event 2. A preceding rain event on 20 September coincided with a rapid decrease in the 24 h mean PM_2.5 from 40 to 10 µg m⁻³, whereas most gas-phase species decreased less, with CO declining from 500 to 350 ppb and total VOCs from 50 to 30 ppb. This pattern suggests preferential particle removal via precipitation scavenging and may provide a physical explanation for “smoke-present but low-PM_2.5” days. It also highlights that smoke identification based on column indicators (e.g., HMS) may not always reflect surface PM exposure when precipitation occurs. Subsequent south-westerly flow reintroduced the residual, aged smoke into the region. As a result, 24 h PM_2.5 concentrations reached 40 µg m⁻³ and hourly CO peaked at 350 ppb before a front cleared the valley, which was especially evident for PM_2.5 but still leaving fire influences for trace gases.

Figure 4 compares pollutant concentrations during smoke-impacted and background periods in September 2020. On average, MDA8 ozone increased by 3 ppb (7 %) during smoke episodes, while NO_x showed no visible change. In contrast, CO, PM_2.5, and TVOCs were significantly elevated, with hourly concentrations increasing by a factor of 3–8 compared to urban background. Average hourly CO concentrations increased from 160 ± 80 ppb (mean ± SD of hourly data) to 420 ± 205 ppb on smoke days. PM_2.5 showed even larger enhancements, rising from 6.1 ± 4.5 to 43 ± 34 µg m⁻³. The background CO and PM_2.5 concentrations agree with typical western US urban levels (150–200 ppb and 5–10 µg m⁻³) (Cope et al., 2024; Lill et al., 2022; Lopez-Coto et al., 2020; Pfister et al., 2011).

Figure 4 Enhanced surface-level pollutant concentrations during smoke-impacted periods. Boxplots show (A) maximum daily 8 h average ozone (MDA8  O₃) and hourly concentrations of (B) nitrogen oxides (NO_x), (C) carbon monoxide (CO), (D) fine particulate matter (PM_2.5), and (E) total measured volatile organic compounds (TVOCs) during smoke-impacted (red) and background (blue) periods in Missoula, Montana, during September 2020. Boxes represent interquartile ranges (25th–75th percentiles) with medians shown as center lines; whiskers extend to 1.5× the interquartile range. Black diamonds and error bars indicate mean ± standard deviation for each group. The right-hand summary block lists mean enhancement factors (smoke-impacted $/$ background) for each pollutant.

Measured hourly TVOCs tripled from 23 ± 15 to 55 ± 24 ppb, with individual species enhancements summarized in Table S2 in the Supplement. Ten VOCs increased by factors of ∼ 4–7, five of which were furanoids or their derivatives. These species included furfural, methylfurfural, 2-furanmethanol, 2-furanone, 5-hydroxy-2-furfural/2-furoic acid, maleic anhydride, acrylonitrile, methyl benzoic acid, methyl methacrylate, and anisole. Forty-eight VOCs increased by a factor of 2–4, while 16 species increased by a factor of 1–2. Fifteen of the 75 measured VOCs are on the EPA's list of HAPs, and they collectively represented 51 % of TVOCs by average molar mixing ratio. During smoke events, total HAPs (THAPs) increased threefold relative to background periods (28 ± 13 ppb vs. 11 ± 7 ppb), with individual species exhibiting enhancements ranging from 40 % (dichlorobenzene) to nearly fivefold (acrylonitrile and maleic anhydride). Figure S3 shows that OH reactivity (OHR) from CO and TVOCs approximately doubled during smoke-impacted days compared to background (13.4 vs. 6.5 s⁻¹). During background periods, isoprene was the largest individual contributor to OHR (0.9 s⁻¹; 13.6 % of total OHR), followed by monoterpenes (10.4 %), formaldehyde (9.0 %), CO (8.8 %), and other individual VOCs. During smoke-impacted days, CO became the largest combustion-related contributor (1.4 s⁻¹; 10.7 %), and the relative contribution of monoterpenes fell to 5.9 %, consistent with previous urban smoke OHR analyses in Boise, ID (Permar et al., 2023). Other VOC groups show similar contributions to OHR in both periods. Together, these enhancements underscore the strong influence of regional wildfire smoke on western US air quality and motivate continued monitoring of smoke exposure in communities such as Missoula.

5 Smoke enhancement ratios (EnRs) and emission ratios (ERs)

We further apply a photochemical-age framework to determine whether VOC enhancement ratios (EnRs) measured in aged smoke reaching Missoula still preserve interpretable source information after multi-day transport and mixing (Fig. 5). Originally developed for urban plumes, this framework provides a useful diagnostic for separating direct emissions, chemical removal and formation (de Gouw et al., 2017). Here, photochemical age is an oxidation-based metric inferred from the VOC ratio clock and should not be interpreted as the exact physical time since emission. Negative values occur when the observed toluene-to-benzene enhancement ratio exceeds the assumed initial ratio and indicate that the clock assumptions are not fully satisfied for those observations, and large uncertainty exists (Sect. S1 in the Supplement). For primary VOCs, linear fits of ln(EnR) versus photochemical age yield back-extrapolated time-zero emission ratios (ERs) and effective OH rate constants (k_OH) under the assumption of plume-integrated aging. For oxygenated VOCs (OVOCs), the same relationships are interpreted qualitatively because secondary formation can offset chemical loss. We provide methodological details in Sect. S1.

Figure 5 Observed and modeled enhancement ratios (EnR) of key volatile organic compounds (VOCs) versus photochemical age. Panels show ln(EnR) versus photochemical age for (A) benzene, (B) toluene, (C) C₈ aromatics, (D) methanol, (E) formaldehyde, (F) acetaldehyde, where EnR is defined as $\mathrm{\Delta }\mathrm{VOC}/\mathrm{\Delta }\mathrm{CO}$ (ppb ppm⁻¹ CO). Observations are shown as black circles and results from GEOS-Chem driven by GFAS fire emissions are shown as red circles. Points represent consecutive 3 h binned data during the main aged-smoke event. Photochemical age is derived from the toluene-to-benzene ratio clock (Sect. S1).

For relatively long-lived primary aromatics, the framework-inferred ERs are physically reasonable, indicating that aged smoke arriving in Missoula still retains measurable source signatures for these species, though the high general urban background, the initial fire emission ratio, and the potential mixing of different sources during transport introduce uncertainties that are hard to resolve. Nevertheless, if we infer back-extrapolated time-zero ERs, this approach estimates 2.05 ± 0.02 ppb ppm⁻¹ for benzene and 1.42 ± 0.02 ppb ppm⁻¹ for toluene relative to CO. These values agree with near-source ranges reported from the WE-CAN and FIREX-AQ campaigns within ∼ 20 % (benzene: 1.8–2.3 ppb ppm⁻¹; toluene: 1.2–1.5 ppb ppm⁻¹) (Gkatzelis et al., 2024; Permar et al., 2021), and are the same as what we inputted for the photochemical age calculation (benzene/toluene fire emission ratio of 0.7 ppb ppb⁻¹; Sect. S1; Eq. S1). This agreement indicates that the framework can recover plausible source-like ERs for relatively unreactive primary species even after multi-day transport. However, inferred ERs for C₈–C₁₀ aromatics are higher than published near-source wildfire values by roughly a factor of 2–3, suggesting that these species were substantially influenced by non-fire sources during transport. Urban mixing is a likely contributor, consistent with city-based measurements showing elevated aromatics relative to CO during wildfire-smoke periods (Cope et al., 2024).

Beyond ER derivation, the ln(EnR)-age framework also provides insight into the net photochemical evolution of VOCs and is not sensitive to an accurate time-zero estimate, unlike ERs. For primary VOCs, EnRs decreased significantly with photochemical age (p≪ 0.01), and the relative decay rates broadly tracked their OH rate constants (k_OH). For example, benzene exhibited only modest decay over 1 week of photochemical aging, whereas toluene and C₈ aromatics decayed more rapidly, consistent with their higher reactivity. The inferred k_OH values (Sect. S1; Eq. S7) were (9.8 ± 1.6) × 10⁻¹³ and (5.6 ± 0.2) × 10⁻¹² ${\mathrm{cm}}^{\mathrm{3}}\mathrm{molec}{.}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$ for benzene and toluene, respectively. These values agree with the literature within ∼ 10 %–20 % and, again, agree with our photochemical age estimates. For C₈ aromatics, the inferred k_OH (7.9 × 10⁻¹² ${\mathrm{cm}}^{\mathrm{3}}\mathrm{molec}{.}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$) lies between those of ethylbenzene and the xylene isomers, as expected for a lumped mixture.

The framework becomes less robust for more reactive hydrocarbons and OVOCs. Effective kOH values inferred for C₉–C₁₀ aromatics and small alkenes are lower than literature values by factors of 3 or more, indicating that simple plume-integrated first-order loss does not adequately describe these compounds in aged smoke. A likely reason is that the most reactive species were already preferentially removed during the earliest, OH-rich stage of plume evolution, before the smoke reached Missoula. Additional uncertainty likely arises from dilution, urban mixing, background correction, and species lumping. OVOCs such as methanol, formaldehyde, and acetaldehyde, exhibited weaker and often statistically insignificant ln(EnR)-age trends (p ≈ 0.02–0.21), and in several cases the relationships flattened with age. Such behavior is consistent with secondary production partially offsetting chemical loss, while additional contributions from non-fire sources such as biogenic emissions may further obscure the relationships. Together, these results show that the ln(EnR)-age framework remains informative for relatively long-lived primary VOCs, but becomes progressively less diagnostic for reactive hydrocarbons and especially OVOCs in regionally aged smoke.

To complement the age-framework analysis, we also calculated event-integrated EnRs of PM_2.5 and 27 VOCs relative to CO as bulk descriptors of the aged-smoke event. These values provide an observational reference for VOC-to-CO relationships in regionally aged smoke and are also used for model evaluation in Sect. 8. Species-level values and literature comparisons are provided in Table S3 and Fig. S4.

6 Non-monotonic O₃–PM_2.5 relationship

Figure 6 shows the relationship between maximum daily 8 h average (MDA8) O₃ and daytime mean PM_2.5 at Missoula during September 2020. During the strongest smoke episodes, observed MDA8  O₃ was suppressed by up to ∼ 15 ppb (e.g., 6–7 September; Fig. 6A). These decreases coincide with ∼ 50 % reductions in downwelling shortwave radiation (SW) (Fig. S5A), a first-order proxy for actinic flux and photolysis J(NO₂) and J(O¹D), consistent with reduced radical production and in situ O₃ formation. Consistent with Fig. 4, MDA8  O₃ was only ∼ 3 ppb (7 %) higher during smoke-impacted periods than during the September 2020 background, while NO_x showed no discernible enhancement. In addition, heavy smoke was associated with lower observed daytime maximum temperature and a shallower planetary boundary layer height (PBLH) (Fig. S5B–C), which may further reduce surface ozone by slowing reaction rates and weakening entrainment of O₃-rich air aloft.

Figure 6 Temporal evolution and PM_2.5 dependence of maximum daily 8 h average ozone during September 2020. (A) Time series of daily maximum 8 h average ozone (MDA8  O₃, ppb) from surface observations (black), GEOS-Chem driven by GFAS fire emissions (red), and the AIRPACT forecast system (blue). Shaded bars (right axis) show observed daily mean PM_2.5 (µg m⁻³), providing context for smoke influence. (B) Relationship between MDA8  O₃ and daytime PM_2.5 for observations (circles), GEOS-Chem (triangles), and AIRPACT (squares). Open gray symbols show individual daily values for observations and model simulations. Solid curves denote locally weighted regression (LOWESS) fits applied consistently across datasets (smoothing fraction = 0.5; robust iteration = 3). AIRPACT data were obtained from the Washington State University AIRPACT data portal (https://airpact.wsu.edu/, last access: 30 July 2026). LOWESS curves are shown as descriptive guides.

Across the limited number of BB days, the relationship between MDA8  O₃ and daytime mean PM_2.5 is non-monotonic (Fig. 6B). A locally weighted scatterplot smoothing (LOWESS) curve suggests that MDA8  O₃ increases at low PM_2.5 (< ∼ 20 µg m⁻³), but decreases at moderate PM_2.5 (∼ 20–60 µg m⁻³). The non-monotonic relationship may reflect limited local O₃ production in the absence of a concurrent NO_x enhancement, together with reduced photochemistry at higher aerosol loadings.

The non-monotonic behavior is not unique to Missoula but is also evident across other northwestern sites during the same September 2020 period. At Cheeka Peak (WA) and Eugene (OR), both of which were near BB sources during that month, daily PM_2.5 reached up to 250 and 450 µg m⁻³, respectively. O₃ similarly increased with PM_2.5 under low smoke but leveled off or declined at higher aerosol loadings (Fig. S6). Similar nonlinear O₃–PM_2.5 relationships have been reported in multi-year, multi-site analyses across the western US, with O₃ increasing with PM_2.5 under light-to-moderate smoke and decreasing under heavier smoke (Buysse et al., 2019; McClure and Jaffe, 2018).

7 Chronic and acute health risks

Figure 7 summarizes the upper-limit chronic inhalation risks from wildfire smoke, assuming that a 2020-style wildfire season were to recur annually over a 70-year lifetime. We do not seek to estimate the total public-health burden of wildfire smoke, but instead use several well-defined screening metrics for PM_2.5 and HAPs with available toxicity values as a transparent basis for comparison. Methodological details are described in Sect. S2. Over lifetime exposure (70 years), the total excess cancer risk attributable to wildfire smoke is 100 cases per million people. Roughly 90 % of this smoke-driven risk is attributed to PM_2.5, with the remaining 10 % arising from HAPs combined. Among HAPs, the major contributors are formaldehyde (42 % of HAP cancer risk), benzene (34 %), acetaldehyde (8 %), acrylonitrile (9 %), and naphthalene (6 %). These relative contributions agree with recent aircraft observations of lofted plumes and confirm that PM_2.5 dominates cancer risk both aloft and at ground level (Pye et al., 2024).

Figure 7 Chronic health risk from fine particulate matter (PM_2.5) and hazardous air pollutants (HAPs). (A) Estimated excess cancer risk (cases per million people) attributed to ambient PM_2.5 versus all combined HAPs. (B) Noncancer hazard index (HI) attributed to ambient PM_2.5 versus combined HAPs. In each panel, slices are exploded for clarity, and labels show the percentage of total risk contributed by each component. The bottom annotation reports the overall risk magnitude: total cancer cases per million for panel (A) and HI for panel (B).

For context, the climatological annual-mean of non-smoke background PM_2.5 is ∼ 6 µg m⁻³ in Missoula. If the same BB-impacted days had instead experienced typical non-smoke background concentrations, the corresponding PM_2.5-attributable excess cancer risk would be ∼ 15 cases per million people. Thus, the 2020-style wildfire season increases the PM_2.5-attributable cancer risk by ∼ 7× relative to the PM-related baseline. Even using a 10-fold lower PM_2.5 unit-risk estimate (4.8 × 10⁻⁵ µg⁻¹ m³), PM_2.5 still accounts for ∼ 50 % of the total cancer risk, with the total risk remaining ∼ 20 cases per million people (i.e., ∼ 30 % higher than the baseline). Thus, even when the community's exposure to wildfire smoke is annualized and a lower-bound PM_2.5 unit-risk estimate is used, aged wildfire smoke at the level experienced in September 2020 in Missoula still shows a quantifiable carcinogenic burden.

Non-cancer chronic risk is expressed as a hazard index (HI), with the dominant effects in this study associated with respiratory and cardiopulmonary systems. The full calculation framework is provided in Sect. S2. The calculated HI = 3 indicates appreciable potential for adverse effects. HAPs account for ∼ 90 % of this index, driven by acrolein (60 % of HAP non-cancer risk), formaldehyde (24 %), acetaldehyde (3 %), and all other species (< 1 % each). PM_2.5 contributes to the remaining ∼ 10 %. The result contrasts with aircraft-based assessments, in which PM_2.5 dominated non-cancer risk (O'Dell et al., 2020; Pye et al., 2024). The lower relative contribution of PM_2.5 observed here likely reflects the evaporation of semi-volatile particulate mass as plumes descend from aloft to the surface (Pagonis et al., 2023; Selimovic et al., 2019, 2020). If PM_2.5 were assumed 10 times more toxic, it would instead contribute ∼ 70 % of the HI.

No acute reference exposure levels (RELs) are available for PM_2.5, so the acute HI analysis was restricted to HAPs. During September 2020, 719 of 720 one-hour windows exhibited HI < 1, implying negligible acute non-cancer concern. A single exceedance (19 September, HI = 1.05) was apportioned to acrolein (52 %), formaldehyde (32 %), benzene (12 %), and other HAPs (< 5 %). Expanding to eight-hour windows (Fig. S7), 49 % of periods remained below the HI = 1 threshold, whereas 51 % exceeded it, aligned with smoke incursions. Of the exceedances, 30 % fell within HI = 1–2 and 21 % reached up to HI = 5. The 8 h hazards were driven by formaldehyde (38 %), acrolein (31 %), and benzene (29 %), with all other HAPs contributing < 2 %.

The recurrence of the same three drivers across acute and chronic metrics highlights two practical needs. First, smoke-impacted communities would benefit from continuous, high-time-resolution monitoring of formaldehyde, acrolein, and benzene to track rapidly evolving exposures and guide advisories. Second, because particle-only strategies do not address these gases, indoor interventions should pair particle filtration (HEPA or MERV 13+ for PM_2.5) with sorbent media (e.g., activated carbon) when using indoor air cleaners during smoke events (Maximoff et al., 2022; May et al., 2021).

8 Model evaluation

8.1 Model uncertainties in fire emissions and OH exposure

Figure 8 compares hourly observations of CO, PM_2.5, and four representative VOCs with model simulations; additional VOC comparisons are shown in Fig. S8. The default GEOS-Chem simulation captures the timing of smoke-impacted enhancements during the first two events (5–7 and 11–20 September), with high correlations (R² = 0.8–0.9 for hourly measurements). The agreement indicates that the model generally captures the timing, transport, and location of fire smoke. However, the model fails to reproduce the observed enhancements during Event 3 (22–24 September), likely reflecting overly efficient modeled wet scavenging during the 21 September rainout and/or unresolved smoke transport pathways discussed in Sect. 4.

Figure 8 Time series of observed and modeled concentrations for key pollutants during September 2020 in Missoula, Montana. Observations (black lines) are compared with the base GEOS-Chem simulation (red). Also shown are two model sensitivity tests: one with tripled biomass burning emissions (GEOS-Chem + 3 × BB; blue) and one with biomass-burning emissions turned off (GEOS-Chem + noBB; gray). Panels show hourly averaged concentrations of (A) carbon monoxide (CO), (B) fine particulate matter (PM_2.5), (C) benzene, (D) methanol, (E) formaldehyde, and (F) acetaldehyde. Model outputs are sampled at the observational location and time. Gray shading indicates smoke-impacted days based on the criteria in Sect. 3. The modeled PM_2.5 is calculated offline from the SpeciesConc diagnostic as the sum of inorganic ions (1.10 × $[{\mathrm{NH}}_{\mathrm{4}}^{+}+{\mathrm{NO}}_{\mathrm{3}}^{-}+{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}]$), black carbon (BCPI + BCPO), organic matter ((OCPO + 1.05 × OCPI) × OM:OC), fine dust (DST1 + 0.30 × DST2), sea salt (1.86 × SALA), and secondary organic aerosol (1.05 × (T$\mathrm{SOAS}+\mathrm{SOAP}+\mathrm{SOAIE}+\mathrm{SOAGX})$).

Even during Events 1–2, GEOS-Chem underestimates the magnitude of fire impacts, with low biases of 30 %–40 % for CO and PM_2.5, 60 %–80 % for primary aromatics, and 50 %–90 % for most OVOCs. Model performance during background periods is notably better, suggesting BB-related emissions and/or chemistry drive the majority of these underestimations under smoky conditions. MEK is the only exception, showing an average high bias of ∼ 20 %, reflecting overestimated anthropogenic or biogenic emissions, as pointed out in earlier studies (Chen et al., 2019; Jin et al., 2023).

GEOS-Chem systematically underestimates event-integrated VOC EnRs by 35 % for benzene and 60 %–80 % for other VOCs (Fig. S9), despite ERs being reproduced within ∼ 20 %–40 % in Fig. 5 and other work (Jin et al., 2023). This ER–EnR decoupling points to chemistry-related biases: the model overpredicts plume photochemical age (Fig. 5) and OH exposure (by ∼ 100 %; Fig. S10), which over-oxidizes VOCs relative to CO and depresses EnRs. For OVOCs, the negative bias of EnRs likely also reflects missing secondary production, particularly for methanol, formaldehyde, and acetone, which exhibit multi-day growth in aged smoke in Fig. 5 and previous work (Alvarado et al., 2020; Bates et al., 2021; Holzinger et al., 2005; Jin et al., 2023). By contrast, GEOS-Chem overestimates the event-integrated PM_2.5 EnR by ∼ 33 %, implying excessive PM_2.5 produced (or retained) per unit CO, consistent with the reported biases in the simplified GEOS-Chem SOA scheme used here (Oak et al., 2022; Pai et al., 2019).

Our previous study using airborne observations to constrain total western wildfire emissions showed that GFAS likely underestimates the fuel consumed in these fires by a factor of three (Jin et al., 2023). Thus, we conducted a 3 × BB sensitivity test, which improved overall model agreement with 2020 MSO measurements and reduced modeled OH by ∼ 2×. This lowered the OH-exposure positive bias from ∼ 100 % to ∼ 30 % during Missoula smoke episodes. The remaining modeled OH positive bias likely reflects missing OH reactivity from unrepresented compounds such as furanoids (Coggon et al., 2019; Jin et al., 2026; Permar et al., 2023). Correspondingly, gas-phase NMBs improve on average across smoke-impacted periods (CO: −40 % to +20 %; benzene: −70 % to −10 %; toluene: −80 % to −50 %), although CO appears to be overcorrected during some periods, especially Event 2. Several OVOCs remain substantially underestimated (e.g., methanol −80 % to −50 %; formaldehyde −70 % to −40 %; acetaldehyde −75 % to −35 %; acetic acid −80 % to −45 %; acetone −60 % to −40 %), consistent with missing secondary production as a likely cause of the persistent OVOC low bias. Meanwhile, PM_2.5 shifts from a −30 % bias to a +80 % bias, indicating that a uniform BB scaling cannot reconcile gases and aerosol. It also points to limitations in the default GEOS-Chem “simple SOA” fire treatment, where the SOA precursor tracer (SOAP) is parameterized proportional to BB CO; scaling fire CO therefore scales SOAP and can over-amplify OA and PM_2.5 during smoke events.

8.2 Model uncertainties in ozone formation under high PM_2.5

GEOS-Chem overestimates MDA8  O₃ in September of 2020 by ∼ 15 ppb on average and fails to capture the observed day-to-day variability (Fig. S11), regardless of the BB emission scenario (noBB, base, or 3 × BB). The persistent positive bias in the noBB simulation suggests an overestimate of background MDA8  O₃. Biases further increase on smoke days, reaching up to ∼ 70 ppb at Missoula and scaling with observed daytime PM_2.5 in the base model (r ≈ 0.7) (Fig. S12). Similar smoke-amplified positive ozone biases have been reported in other CTMs (e.g., CMAQ-based simulations; Fig. 6A) and in prior CTM work (Baker et al., 2016, 2018; Zhang et al., 2014), suggesting that regional-to-global CTMs may systematically overestimate background O₃ and/or net O₃ production under BB influence.

However, the main value of this analysis is not simply to document another case of positive O₃ bias in smoke. Rather, Fig. 6B shows that GEOS-Chem does not reproduce the observed non-monotonic O₃–PM_2.5 relationship. At Missoula, the model captures the initial O₃ increase when PM_2.5 is below ∼ 30–40 µg m⁻³, but it predicts increasing O₃ at higher PM_2.5, whereas observations flatten or decline beyond this threshold. A similar failure is evident in AIRPACT, a publicly available regional air quality forecasting system (Chen et al., 2008; Vaughan et al., 2004). Across Missoula and other western US sites (e.g., Eugene, OR; Yreka, CA; Cheeka Peak, WA; Fig. S6), observations show strong O₃ suppression under heavy smoke, while AIRPACT predicts a monotonic increase of O₃ with PM_2.5. Together, these consistent biases across two independent CTM frameworks point to a shared model limitation in representing the nonlinear ozone response under intense BB smoke, which may partly explain the positive O₃ bias commonly seen in CTMs under wildfire influence.

Part of the positive model bias may reflect errors in meteorological processes under smoke conditions. First, GEOS-Chem captures the relative reduction in broadband shortwave radiation on smoke days, suggesting reduced photolysis (i.e., J(O¹D) and J(NO₂)), but it still overestimates the absolute downwelling shortwave flux by ∼ 25 % during smoke-impacted periods (versus ∼ 20 % during non-BB periods) (Fig. S5). These “over-sunny” conditions likely result in excessive radical production, although we lack observed J values for a direct photolysis evaluation. Second, model temperature is biased high by +1.13 °C on smoke days (implying ∼ 20 % faster PAN thermal loss at ∼ 289 K), which could favor enhanced NO_x recycling and contribute to higher modeled O₃. Third, model PBLH tracks the temporal variability of assimilated products (e.g., HRRR, 3 km) but is higher by ∼ 500 m (33 %) during smoke. The elevated PBLH can increase entrainment of O₃-rich air from aloft and thereby contribute to higher modeled background O₃. In addition, heterogeneous ${\mathrm{NO}}_{\mathrm{3}}/{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{5}}$ loss in BB plumes may be too weak (Decker et al., 2019, 2021), which could prolong NO_x lifetime and contribute to overproduction of O₃ under smoke (Shen et al., 2025). We cannot attribute the discrepancy to a single process, but collectively these factors provide a plausible explanation for the spurious modeled MDA8  O₃ increases during heavy smoke and the growth of the O₃ error with PM_2.5.

8.3 Model uncertainties in health-risk estimates

GEOS-Chem captures the order of magnitude of the smoke-attributable excess lifetime cancer risk but still underestimates it by ∼ 40 % (63 vs. 100 per million), primarily due to its low bias in PM_2.5. Because cancer risk is dominated by PM_2.5 whereas HI is driven by a small number of high-potency HAPs, model errors in HAP composition have a much larger impact on HI than on cancer risk. Accordingly, the base simulation substantially underestimates chronic non-cancer risk (HI ≈ 0.3, ∼ 10× lower than observed). For acute risk, GEOS-Chem predicts no 1 h HI exceedances (max 1 h HI = 0.2) and only limited 8 h exceedances (15 h with HI > 1; max 8 h HI = 1.2), far fewer and weaker than observed (476 h; max 8 h HI = 6.3).

In the 3 × BB sensitivity run, the modeled cancer risk increases to 155 per million (exceeding the observation-based estimate), highlighting the uncertainty in modeled PM_2.5 under smoke conditions. Acute non-cancer risk increases in the 3 × BB run (135 h with 8 h HI > 1; max 8 h HI = 3.4), but the model still predicts no 1 h exceedances (max 1 h HI = 0.5). This low bias in HI is driven by incomplete BB VOC representation, including underestimated BB VOC burdens (by ∼ 3×), missing secondary sources, and missing high-potency HAPs (e.g., acrolein, which accounts for > 50 % of the observed HI). Overall, CTMs such as GEOS-Chem can provide a first-order estimate of PM_2.5-driven cancer burden, but the predicted HI is much less reliable without improved HAP speciation, secondary production, and inclusion of high-potency compounds.

9 Conclusion

Missoula, Montana (46.9° N, 114.0° W), experienced persistent and chemically complex wildfire smoke during September 2020, reflecting the widespread influence of wildfires across California and the Pacific Northwest. We report hourly surface observations of PM_2.5, CO, NO_x, O₃, and 75 speciated VOCs in Missoula during the wildfire-smoke-impacted month of September 2020. Of the 75 measured VOCs, 15 were classified as HAPs by the EPA. Leveraging comprehensive measurements, we quantified EnRs of PM_2.5 and VOCs relative to CO in smoke aged for several days, and characterized the temporal evolution of criteria pollutants and VOCs under smoke-impacted and no-/low-smoke conditions. We performed nested GEOS-Chem simulations and compared them with observations to constrain BB emissions and photochemistry in the model. Finally, we assessed public health risks based on regulatory exposure metrics for PM_2.5 and HAPs.

We find that combining gas-phase tracers with conventional particle-based diagnostics improves the identification of surface-level wildfire smoke in urban environments. While traditional aerosol-based smoke diagnostics (HMS, AERONET, and PM_2.5) remain useful for detecting dense, optically thick plumes, they often misclassify events when plumes are aloft, optically thin, or when rainfall removes aerosols. In contrast, VOC tracers, especially furan and maleic anhydride (MA), provide complementary value: furan identifies fresh smoke, while MA captures more aged plumes. We thus develop a balanced multi-indicator rule that requires concurrent enhancements in at least three independent tracers, thereby capturing smoke episodes with higher confidence. Three major smoke events driven by long-range transport resulted in multi-day elevations of CO and total VOC abundances by factors of 2–3, while ozone increased only modestly and NO_x showed little to no change. In addition, PM_2.5 and total HAPs (which comprised approximately half of the measured TVOCs by molar fraction) increased by factors of ∼ 3–8 during smoke periods. Concurrently, the total measured OH reactivity approximately doubled (13.4 vs. 6.5 s⁻¹), shifting the relative OHR contributions from biogenic tracers in background conditions toward CO and BB VOCs during smoke. These results highlight the strong influence of regional wildfire smoke on air quality in the western US.

Photochemical-clock analysis constrains ERs for a wide range of primary VOCs in aged BB smoke, yielding values consistent with those reported for western US wildfires. The observed decay patterns of primary VOCs broadly follow their known OH rate constants (k_OH). Some underestimation of k_OH (by factors of ≈ 3) for reactive VOCs likely arises from rapid depletion, dilution, or background contamination. OVOCs such as methanol, formaldehyde, and acetaldehyde also decrease with age, but their shallow slopes indicate substantial secondary production offsetting primary decay.

Intense smoke episodes in Missoula suppressed surface ozone, with MDA8  O₃ declining by up to ∼ 15 ppb (6–7 September), coincident with 50 % reductions in solar radiation and lower planetary boundary layer heights (PBLH), confirming that reduced photolysis and weaker mixing limit in situ ozone formation. Across smoke-impacted days, the MDA8  O₃–PM_2.5 relationship was non-monotonic: ozone increased under light smoke (PM_2.5 < ∼ 20 µg m⁻³) but tended to flatten or decrease at higher smoke loadings (≳ 20 µg m⁻³).

Chronic and acute health assessments indicate that prolonged exposure to 2020-level wildfire smoke in Missoula poses quantifiable health risks. If a season of similar intensity were to recur annually, the smoke-attributable excess lifetime cancer risk would be ∼ 100 cases per million people (≈ 7× the non-smoke baseline of ∼ 15 per million for the same days), with PM_2.5 accounting for ∼ 90 % of the cancer burden. The chronic non-cancer hazard index (HI) is ∼ 3, with ∼ 90 % attributable to HAPs (driven primarily by acrolein and formaldehyde) and the remaining ∼ 10 % to PM_2.5. For acute non-cancer risk, 1 h HI exceedances are rare, but 8 h averaging reveals frequent smoke-aligned exceedances (∼ 51 % of periods with HI > 1). Across both chronic and acute metrics, the repeated importance of the same few HAP drivers (notably acrolein, formaldehyde, and benzene) indicates that mitigating wildfire-related health impacts should extend beyond particle-only strategies to include gas-phase toxics.

The GEOS-Chem simulation reproduces the timing and transport of major smoke plumes (R² = 0.8–0.9) but underestimates the magnitude of smoke enhancements, with model low biases of ∼ 30 %–40 % for CO and PM_2.5 and ∼ 60 %–90 % for most VOCs. The model overestimates OH exposure by roughly 100 %, leading to excessive chemical loss and artificially low EnRs for reactive VOCs. Persistent low EnRs for oxygenated VOCs further suggest missing secondary production during multi-day aging, whereas the +33 % bias in PM_2.5 EnR points to overly efficient secondary aerosol formation and/or retention in the GEOS-Chem simplified “simple SOA” treatment. Tripling BB emissions (3× BB) largely improves the modeled primary compounds (e.g., benzene from −70 % to −10 %) and reduces the model bias of OH exposure to ∼ 30 % within smoke-impacted periods, but the modeled PM_2.5 is worse (from −30 % to +80 %), demonstrating that uniform emission scaling cannot reconcile errors in both gases and aerosol. Together, these results highlight that reducing model biases in smoky environments requires not only better fire emission magnitudes but also improved representation of missing OH reactivity (e.g., furanoids), secondary OVOC formation, and smoke SOA parameterizations.

GEOS-Chem exhibits a persistent positive bias in September MDA8  O₃ in Missoula (+15 ppb on average) and fails to reproduce observed day-to-day variability, regardless of the BB emissions scenario (base, 3 × BB, or noBB). This is in part due to the modeled O₃ being too high for general urban background, and the model bias further amplifies on smoke days, reaching +70 ppb and increasing with observed daytime PM_2.5 (r ≈ 0.7). More importantly, the model fails to reproduce the observed non-monotonic O₃–PM_2.5 response, where observations flatten or decline beyond ∼ 30–40 µg m⁻³. Instead, GEOS-Chem continues to increase O₃ with PM_2.5. The same monotonic behavior is also present in the independent CMAQ-based AIRPACT forecast system across Missoula and other western US sites, pointing to a shared CTM limitation in representing nonlinear ozone responses under intense BB smoke. These model biases are likely related to overestimated radiation, warm temperature biases, PBLH/entrainment errors, and potentially too-weak heterogeneous ${\mathrm{NO}}_{\mathrm{3}}/{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{5}}$ loss in CTMs, although isolating the dominant drivers will require additional observational constraints.

GEOS-Chem underpredicts the smoke-attributable lifetime excess cancer risk by 40 % (63 vs. 100 cases per million), due to the underestimated PM_2.5. The model's capability to predict non-cancer HI is even worse (∼ 0.3 in GEOS-Chem vs. 3 in observations). The model also underrepresents acute hazards, predicting limited 8 h exceedances (15 h with 8 h HI > 1; max 8 h HI = 1.2), far fewer and weaker than observed (476 h; max 8 h HI = 6). Overall, these results indicate that GEOS-Chem can provide a defensible screening-level estimate of PM_2.5-dominated cancer burden, whereas HI estimates remain unreliable without improved BB VOC speciation, secondary production pathways, and inclusion of high-potency toxics such as acrolein.