Decadal Trends and Variability in Intermountain West Surface Ozone near Oil and Gas Extraction Fields

Abstract. Decadal trends in the annual fourth-highest daily maximum 8-hour average (A4DM8HA) ozone (O3) were studied over 2005–2015 for 13 rural/remote sites in the U.S. Intermountain West. No trends were observed in A4DM8HA O3 at two reference sites, which are located upwind of and thus minimally influenced by emissions from oil and natural gas (O&NG) basins. Trends, or a lack thereof, varied widely at other 11 sites in/near O&NG basins resulting from different controlling factors rather than a simplistic, uniform one. The decreasing trends at Mesa Verde (−0.76 ppbv/yr) and Canyonlands National Park (−0.54 ppbv/yr) were attributed to a 37 % decrease in natural gas production in the San Juan Basin and 35 % emission reductions in coal-fired electricity generation, respectively. The decreasing trend (−1.21 ppbv/yr) at Wind Cave National Park resulted from reduced solar radiation due to increasingly frequent precipitation weather. The lack of trends at remaining sites was likely caused by the increasing O&NG emissions and decreasing emissions from other activities. Wintertime O3 stagnant events were associated with the Arctic Oscillation (AO). Box model simulations suggested that both volatile organic compounds (VOCs) and nitrogen oxides emission reductions during negative AO years while VOC emission reductions alone in positive AO years could effectively mitigate high wintertime O3 within the O&NG basins. Our findings suggest that emissions from O&NG extraction likely played a significant role in shaping long-term trends in surface O3 near/within O&NG basins and hence warrant consideration in the design of efficient O3 mitigation strategies for the Intermountain West.



Introduction
Tropospheric ozone (O3) is a short-lived trace gas that either originates naturally from the stratosphere (Stohl et al., 2003) or is produced in situ by photochemical oxidation of nitrogen oxides (NOx) and volatile organic compounds (VOCs) or carbon monoxide (CO) (e.g., Monks et al., 2009).On 1 October 2015, the U.S. Environmental Protection Agency (EPA) lowered the National Ambient Air Quality Standard (NAAQS), namely the O3 design value, to 70 ppbv to improve protection of public health and welfare (EPA, 2015).Unlike in the eastern U.S., where NOx emission reductions have led to O3 declines, in the Intermountain West no decreasing trends were observed over 1988 -2014 and surface O3 could exceed the NAAQS in both winter and summer (Schnell et al., 2009;Cooper et al., 2012;Edwards et al., 2014;Lin et al., 2017).Most studies attributed the increasing background O3 trends over the western U.S. to increasing anthropogenic Asian emissions (e.g., Cooper et al., 2012;Parrish et al., 2012;Lin et al., 2017).Lefohn et al. (2012) found that stratospheric intrusion affected the interannual variability of surface O3 at both high-and low-elevation sites.However, a multi-model study suggested that surface O3 in the U.S. was far more sensitive to domestic emission changes than to global emission changes in spring and summer (Reidmiller et al., 2009).Thus, long range transport from Asia and stratospheric intrusion may not be the sole contributors to the observed increasing trends in western U.S.
Recent expanded use of horizontal drilling and hydraulic fracturing technologies enabled access to more natural gas resources in shale deposits (EIA, 2015).Extraction of oil and natural gas (O&NG) can result in the emission of O3 precursor gases and episodes of elevated O3 mixing ratios in O&NG production basins frequently exceeding the NAAQS.Most notable were the episodes that occurred during the cold season in the Intermountain West (Carter and Seinfeld,  2012; Edwards et al., 2014;Rappenglück et al., 2014;Schnell et al., 2009).Carter and Seinfeld (2012) studied four O3 episodes in the Upper Green River Basin in Wyoming and found one episode highly NOx-sensitive and the others VOCs-sensitive.Edwards et al. (2014) used a box model to simulate a stagnant 6-day high O3 event in the Uintah Basin and found that carbonyl photolysis was the dominant oxidant source.McDuffie et al. (2017) studied O3 in the Colorado Northern Front Range (NFR) in summer 2012 and found that O&NG emissions contributed to 17% of the modeled maximum photochemically produced O3.To the best of our knowledge, no studies have investigated the long-term impact of expanded O&NG extraction activities on O3 over the intermountain U.S. using more than ten years of measurement data.
The effect of emissions from O&NG production on surface O3 is difficult to quantify because a wide range of factors work together to determine the mixing ratio of O3 at a given location, such as other anthropogenic emissions, natural emissions, transport processes, stratospheric intrusion, O3 photochemistry, and changing global climate (Parrish et al., 2013); however, the significance of these factors varies by season.For example, the number of high O3 days in summer was strongly correlated with wildfire burned area over the western U.S. as wildfires are important sources of NOx, CO, and VOCs (Jaffe et al., 2008;Jaffe and Wigder, 2012).The frequency and intensity of wildfires in the western U.S. may be increasing, driven by increasing temperatures, earlier snowmelt, drier conditions, and accumulation of fuels (Westerling et al., 2006;Jaffe et al., 2008).Moreover, all components that can impact the distribution of O3 are associated with varying atmospheric circulation patterns on interannual to decadal time scales (e.g., Lin et al., 2014;Zhou et al., 2017).For instance, more frequent stratospheric intrusion events occurred in late spring, when the polar jet meandered towards the western U.S. following La Nina events (Lin et al., 2014).Therefore, it is fundamentally important, yet also challenging, to extract the impact of increasing emissions from O&NG production on surface O3 levels from all of these various factors.However, these findings could ultimately lead to improved modeling and to better assess the efficacy of emission controls.
A most common approach to assess the influence of a certain factor or source on O3 concentrations is through model simulations (e.g., Rodriguez et al., 2009;Lin et al., 2017).
However, the application of current chemical transport models used to investigate the impact of emissions from O&NG extraction on tropospheric O3 is particularly challenging because of the large uncertainties existing in the emission inventories of O&NG extraction (Pé tron et al., 2014;Helmig et al., 2014;Ahmadov et al., 2014;Thompson et al., 2017;Allen, 2014;Allen, 2016).
Past studies have found that the emissions of methane and VOCs from oil and gas extraction was largely underestimated (Karion et al., 2013;Ahmadov et al., 2014).Therefore, for this study the approach of data analysis was chosen in order to provide measurement-based contribution estimates of various factors.Specifically we investigated the decadal impact of emissions from O&NG extraction on variability and long term trends of surface O3 in the Intermountain West since 2005, when O&NG extraction activities started to expand rapidly (EIA, 2015).In this study, we used long term O3 measurement data obtained from the National Park Service (NPS), Clean Air Status and Trends Network (CASTNET), and Wyoming Department of Environmental Quality (WDEQ), many of which are from Class I areas as designated by the Clear Air Act.However, some of these sites have been reportedly influenced by O&NG emissions since 2005 (EIA, 2015).This study was focused on the summer and winter, when exceedances of NAAQS tended to occur in the Intermountain West (Edwards et al., 2014;Lin et al., 2017).

Site Selection and Data description
Long term surface O3 observations were available at 18 sites in remote and rural areas of the U.S. Intermountain West (Figure 1 and Table S1).Among the 18 sites, 11 sites are located within 100 km of a shale play and are more likely affected by emissions from O&NG extraction activities (Figure 1 Therefore, GRBA, ZION, GRCA, PEFO and BADL were not included in the analysis of this study. Continuous hourly measurements of O3 were available at each site.We calculated the annual fourth-highest daily maximum 8-hour average (A4DM8HA), upon which the O3 design value was based, for the 13 monitoring sites and their trends were examined through ordinary linear least-square regression.The trends were also calculated separately for 5 th , 50 th , 95 th percentiles of daily maximum 8-hour average (DM8HA) O3 in winter and summer.

Box Model
Box model simulations were performed using BOXMOX (Knote et al., 2015) and the NCAR Tropospheric Ultraviolet and Visible Radiation Model (TUVv5.3.1)(Madronich et al., 1998).BOXMOX is an extension to the Kinetic PreProcessor which allows an easy set up for zero dimension box model simulations (Knote et al., 2015).We selected the Master Chemical Mechanism (MCM v3.3) for the model chemistry scheme.The MCM is a near-explicit chemical mechanism including gas-phase tropospheric degradations of 74 VOCs and OVOCs with a total of 5259 species of and 15176 reactions (Jenkin et al., 2015).
All simulations were initialized with observations from field campaigns in the Intermountain West.These campaigns were Uintah Basin Winter Ozone Study (UBWOS) over 2012 -2014 (Edwards et al., 2013), Nitrogen, Aerosol Composition, and Halogens on a Tall Tower (NACHTT) in 2011 (Swarthout et al., 2013), and Front Range Air Pollution and Photochemistry Éxperiment (FRAPPÉ) in 2014 (Pfister et al., 2017).Section S8 provides detailed information about field campaign measurements.The modeled mixing ratios of CO, CH4, NO, NO2, and non-methane VOCs were forced to match the measured diurnal profiles by introducing turbulent mixing to the model (Knote et al., 2015;Section S8.3).Therefore, the "box" can dynamically exchange with its surroundings.The model was integrated forward with a time step of 10 minutes until the concentrations reach a diurnal steady state, when the cycles of simulated species exhibit less than 0.01% variation from the previous day (Edwards et al., 2013).
The last 24 h were used to represent the simulated diel cycles.Section S8 provides further information on the model treatment of all chemical observations.

General Characteristics in Long-term Variations of O3
First, changes in the A4DM8HA were examined for the 13 sites before and after 2005, when O&NG extraction started expanding rapidly.Before 2005, statistically significant increasing trends were observed at the reference site CRMO as well as three other sites, CANY, CNTL, and PNDE, near O&NG basins with no significant trends at the other sites (Figure 2).
Trends in seasonal 5 th , 50 th , 95 th percentiles of DM8HA O3 were further examined (Table S3).At the two reference sites, no trend was observed in any seasonal percentiles of DM8HA O3 at YELL, while statistically significant decreasing trends (-0.14 to -1.18 ppbv yr -1 ) were observed in wintertime 50 th /95 th and summertime 5 th /50 th /95 th percentile values at CRMO.Significant decreasing trends were also observed at 4 sites located in/near O&NG basins.Specifically, at MEVE significant decreasing trends of -0.49 to -0.63 ppbv yr -1 were found in seasonal 50 th /95 th percentile values in both winter and summer (Table S3).In contrast, significant decreasing trends were only found in summertime 5 th /50 th at CANY, summertime 95 th at WICA, and wintertime 50 th /95 th percentile values at GTHC.Overall seasonal median values near the O&NG extraction basins showed relatively more consistent patterns of interannual variation at relatively low elevation sites (i.e., CANY, DINO, MEEK, RANG, CAMP, MEVE, and WICA) in winter, while in summer, large differences were found between sites (Figure 3).Previous studies have examined long term trends in tropospheric O3 for varying time periods over the continental U.S. (e.g., Cooper et al., 2012;Cooper et al., 2014;Lin et al., 2017;Gaudel et al., 2018).Gaudel et al. (2018) assessed surface O3 trends at 2702 non-urban sites from the Tropospheric Ozone Assessment Report (TOAR) database.They found significant increasing trends (0.5 to 1.0 ppbv yr -1 ) in daytime average O3 at 6 sites in the Intermountain West and no trends at 8 other sites including MEVE, PNDE, YELL in winters of 2000 -2014 (Gaudel et al. 2018).However, in this study a significant decreasing trend (-0.49ppbv yr -1 ) in wintertime 95 th percentile values was observed at MEVE after 2005, possible causes for which are discussed further in Section 4. In summers of 2000 -2014, significant decreasing trends (~ -0.5 ppbv yr -1 ) in daytime average O3 were observed by Gaudel et al. (2018) at most sites including CRMO, PNDE, CANY, MEVE, which was consistent with trends we observed (Table S3).
The following questions arose from examining the above characteristics in surface O3 at the 11 sites near/in O&NG basins: 1) Were the decreasing trends in A4DM8HA at MEVE, CANY, and WICA potentially linked in some ways to changes in O&NG production activities nearby?
2) Was there seasonal variation in the influence of emissions from O&NG extraction on seasonal DM8HA O3?
3) Were there changes of O3 production regime as a result of emissions from O&NG basins over the years?These questions are addressed in the following sections with the goal to delineate the role of emissions from O&NG production in long-term trends in surface O3 at the 11 sites.

Contribution of emissions from O&NG production vs other anthropogenic sources
The two reference sites, YELL and CRMO, are located >100 km upwind of O&NG  S4).This is supported by the strong declines in tropospheric NO2 column densities in summer at a rate of ~ -2.8 × 10 13 molecular cm -2 yr -1 (p ≤ 0.05) surrounding MEVE around the San Juan Basin over 2005 -2015 (Figure 4b).Decadal changes were not computed for winter due to scarce wintertime column NO2 retrievals over the Intermountain West.Because O3 production in U.S. rural areas was typically sensitive to changes in NOx emissions (Cooper et al., 2012), the decreasing NOx emissions associated with reduced natural gas production and emission controls in other sectors were most likely the primary cause for the significant decrease  S4).
As a result of these changes, significant decreasing trends in the seasonal 95 th percentile would be anticipated, but were not observed in both winter and summer.This indicates that there are additional factors producing significant influences on O3 in the area, which are further investigated in Sections 5 -6.
No trend was observed in the A4DM8HA O3 at CAMP, which was likely associated with changes in natural gas vs. oil production activities in the Powder River Basin.CAMP is located within the Powder River Basin, which is traditionally known for its coal production and has been one of the fastest growing oil producing regions in recent years.While natural gas production in Campbell County decreased by 62% (Figure 4a), oil production in Campbell County increased from 9×10 8 barrels in 2005 to 2.3×10 9 barrels in 2015 (http://wogcc.wyo.gov/).WICA is located downwind of the Powder River Basin, experiencing a significant decreasing trend in the A4DM8HA at a rate of 1.21 ppbv yr -1 (p = 0.05) over 2005 -2014.Curiously, total NOx emission increased by 29% from 1,327 tons in 2005 to 1,712 tons in 2014 in Custer County, where WICA is located (Table S4).The decreasing A4DM8HA at WICA is attributed to other factors, as discussed in Section 6. DINO, RANG, and MEEK are located within the Uintah-Piceance Basin O&NG fields (Figure 1).Increased emissions, most likely from a 68% rise in natural gas production over 2005 -2015 in the Uintah Basin (Figure 4a), together with snow cover (Edwards et al. 2014;Section xx), could have contributed to very high A4DM8HA O3 at DINO and RANG.Meanwhile, NOx emissions from highway vehicles decreased by 86% in Uintah County (Table S4).MEEK, in Rio Blanco County, is located at the eastern edge of the Piceance Basin and relative lower 4DM8HA The four high-elevation sites, ROMO, PNDE, CNTL, and GTHC, experienced no trends in their A4DM8HA due likely to opposite changes in emissions from O&NG production, urban sources, and stratospheric intrusion (or a lack thereof).ROMO is located to the west of the Denver-Julesburg Basin and has experienced northwesterly and southeasterly wind most frequently over 2005 -2015 (Figure 5a) with higher O3 (>60 ppbv) from the east and southeast (Figure 5e).Natural gas production in Weld County increased by nearly a factor of 3 from 2009 to 2015, which coincided with use of horizontal drilling starting in 2009 (Figure 4a).Meanwhile the NOx emission reductions of ~37% from the urban area of Denver offset the effect of increased NOx emissions from O&NG extraction, as evidenced by the significant declines in tropospheric column NO2 (Figure 4b).PNDE is located to the north of the Green River Basin, which has two U.S. largest natural gas fields, i.e., Pinedale Anticline and Jonah Field.Natural gas production increased steadily over 2005 -2010 followed by a decline in 2011 and onward.At PNDE, while the dominant winds were from the northwest, southeast, and northeast (Figure 5b), O3 mixing ratios corresponding to southwesterly wind, from where expansive O&NG extraction occurred, reached higher than 60 ppbv (Figure 5f).Overall, no trend was observed in natural gas production in Sublette County in the basin (Figure 4a), which was consistent with the relative constant A4DM8HA at PNDE.CNTL was located to the north of the North Park Basin and southeast of the Green River Basin (Figures 1 and 5g).In the past decade, natural gas production in Jackson County in the North Park Basin decreased by 72% (Figure 4a) and oil production increased comparatively by 73%, which likely contributed to constant A4DM8HA O3 at CNTL.Similar to all high-elevation sites (i.e., ROMO, PNDE, CNTL), GTHC is another site located at the mountain ridge.The high O3 concentrations from the north (> 65 ppbv) and southwest (> 52 ppbv) indicate that the site was likely under the influence of aforementioned opposite changes in emissions from the North Park Basin, Denver-Julesburg Basin, and San Juan Basin (Figures 5d and h) ultimately resulting in the site exhibiting no trend in the A4DM8HA.

Contributions of transport from the Arctic and the West Coast and stratospheric intrusion in winter
Winter seasonal median DM8HA O3 showed consistent interannual variation with higher O3 levels (38 -59 ppbv) in 2011 and 2013 whereas lower values (32 -44 ppbv) in 2012 at the two references sites and most sites near/within O&NG basins, except the three high-elevation sites (i.e., GTHC, ROMO, CNTL) (Figures 3c-d).Seasonal 5 th and 95 th percentile values exhibited similar patterns (not shown).Site average O3 mixing ratios exceeded the decadal average 41 ppbv in 2006, 2008, 2010, 2011, 2013, and were below average in 2007, 2009, 2012, 2014, 2015.The difference of the 850 hPa geopotential height between the higher and lower O3 years is shown in Figure 6a.A pronounced positive difference of ~60 geopotential meters (gpm) was observed over the Arctic polar region and a negative difference of ~ 40 gpm over the midlatitudes.This pattern strongly indicates a negative (positive) phase of Arctic Oscillation (AO) associated with a higher (lower) O3 winter in the Intermountain West.Indeed significant negative correlations (-0.91 to -0.58) were found between the AO index and winter seasonal 50 th /95 th percentile DM8HA O3 at most sites, including the reference site YELL over 2005 -2015 (Table 1 and Figure 6b).Significant correlation was not found at the other reference site CRMO, which is located to the south of the Pioneer Mountains and received few air masses (3%) from the north (Figure S1b).Weaker correlations were also observed in the seasonal median DM8HA at CNTL and PNDE.During a longer time period dating back to years prior to 2005, significant negative correlation was observed between seasonal median DM8HA O3 and the AO index at the reference site YELL (-0.45 over 1997 -2005), as well as at CANY (-0.47 over 1993 -2015), CNTL (-0.46 over 1990 -2015), and ROMO (-0.67 over 1988 -2015).This indicates that the impact of AO on wintertime surface O3 has been consistently significant during different time periods over the past decades in the Intermountain West.
The negative correlation between AO and wintertime O3 could be a result of multiple factors, including transport from the Arctic and California, in situ O3 photochemical production, and stratospheric intrusion.While long-term O3 trends in lower percentiles could be influenced by long-range transport, O3 in higher percentiles was more sensitive to local emissions and extreme events (Lin et al., 2017).First, the influence of transport from the Arctic was examined.
Artic influence would be foremost demonstrated in surface temperature.During negative AO years, frigid winter air masses could extend far into the middle of North America, contributing to relatively cold temperatures over the U.S. Intermountain West (Hess and Lamarque, 2007).
Significant positive correlation was found between the AO index and surface temperature at CANY, DINO, MEEK, RANG, and GTHC.Except for GTHC, all sites are located within the basins and weaker positive correlations were observed at other sites, except for CNTL, the highest elevation site in this study (3175 m amsl) (Table 1 and S1).Significant negative correlation was also found between surface temperature and wintertime O3 at sites within the basin (i.e., CANY, DINO, MEEK, and RANG), with the potential effect of stratospheric intrusion removed (Partial correlation in Table 1).Higher correlation coefficients were found in seasonal 50 th /95 th DM8HA at DINO, MEEK, and RANG, indicating that the sites were more likely under local influence (Table 1).During a colder winter with constant snow cover, high concentrations of O3 precursors emitted from O&NG extraction were trapped in a shallow boundary layer within the basins followed by strong O3 photochemical production during daytime leading to the highest O3 mixing ratios of the season (Schnell et al., 2009;Edwards et al., 2014).For sites outside the basins that are exposed to distant sources (Figure S1), the significant correlation between wintertime median DM8HA O3 and the AO index is a strong indication of impacts of long-range transport from source regions afar.
Backward trajectory simulations suggest that long distance transport from the West Coast may obscure the correlation between the AO index and surface O3 at MEVE, ROMO, YELL, WICA, CNTL, and PNDE.A total of 902 five-day backward trajectories were performed at each site during winters of 2006 -2015 (Figure 7).The difference in numbers of trajectories between high and low O3 years suggests that more air masses reached the six sites coming from low altitudes (<900 hPa) over the area of southern California and Arizona, centered around 32.5°N/120°W during high O3 years (Figure 7).At YELL, WICA, PNDE, ROMO, and CNTL, more than 80 backward trajectories passed through the surface of California (Figure 7), where large amounts of O3 could be produced from emissions of facilities in urban areas (Ryerson et al., 2013).Huang et al. (2013) also found that anthropogenic emissions from southern California could contribute 1 -15 ppbv to surface O3 in air masses transported to the Intermountain West.
Additionally, more than 20 trajectories originating at YELL, CNTL, ROMO, and MEVE came from higher altitudes (>700 hPa) in the north (>45 °N) during high O3 years (Figure 7).During negative AO (high O3) years, cold Arctic mid-tropospheric air rich in O3 would cross Canada and reach the U.S. Intermountain West more frequently (Hess and Lamarque, 2007).This is in part evidenced in the significant negative correlation (r = -0.49,p = 0.01) between the AO index and Stratospheric influence on surface O3 was also examined.One of the physical characteristics of stratospheric air is high values of potential vorticity (PV).The differences in PV between high and low O3 years exhibited positive anomalies of ~0.5 PV (10 -6 m 2 s -1 kg -1 K) in the mid-troposphere (315 K) over the Intermountain West (Figure 8).Significant negative correlations of -0.64 (p < 0.01) over 1988 -2015 and -0.76 (p = 0.01) over 2005 -2015 were found between PV and the AO index over the study region (37.5 °N -45 °N, -102.5 °W --110°W).This indicated that stratospheric intrusion over the Intermountain West was associated with the negative AO with a strong impact at sites outside the basins over the region.

Effects of precipitation weather on summertime O3
We found that summertime O3 in the Intermountain West was strongly correlated with relative humidity and the wildfire index total fire index (TFI) (Table S6).Studies have suggested that in the Intermountain West, summertime O3, especially the seasonal high percentile levels, was impacted by photochemical production from wildfire emissions (Table S6; Jaffe and Wigder, 2012;Lu et al., 2016).To avoid repeating the extensive body of literature on this topic, our own analysis of the effect of wildfires on O3 can be found in Section S7.One key point coming out of our analysis was, with fire influence removed as indicated in partial correlation, significant negative correlations between summertime O3 and relative humidity were found at the two reference sites (CRMO and YELL), as well as at the 7 O&NG emission influenced sites (CANY CAMP, MEEK, RANG, ROMO, PNDE, and WICA) (Table S6).Further investigation, detailed as follows, revealed that this correlation essentially illuminated the effect of reduced solar radiation on O3 production caused by cloudiness associated with precipitation weather.
In the Intermountain West, precipitation varies by elevation and latitude among other factors (Williams et al., 1962).The decadal summertime average precipitation abundance was 0.17 kg m -2 for the reference site YELL and 0.15 -0.27 kg m -2 for GTHC, CNTL, ROMO, and WICA, which were higher than that at CRMO (0.08 kg m -2 ) and the remaining sites (Figure 9a).YELL, GTHC, CNTL, ROMO, and WICA are located in the mountains, where localized convective storms could develop quickly in summer leading to higher precipitation amounts (Williams et al., 1962) and naturally, more cloudiness at these five sites compared to the other sites located within the basins.This is consistent with the significant negative correlations between summertime O3 and relative humidity at CNTL, ROMO, WICA, and the reference site YELL (Table S6).Note that no significant partial correlation was found between summertime O3 and relative humidity at GTHC, indicating that summertime O3 at this site was mostly influenced by wildfire emissions.Significant correlations were also found between solar radiation and summertime median DM8HA O3 at CNTL (r = 0.67, p = 0.05), ROMO (r = 0.43, p = 0.09), WICA (r = 0.62, p = 0.05), and YELL (r = 0.72, p = 0.01).These results indicate that reduced solar radiation flux near the surface resulting from increased cloudiness accompanying increased precipitation resulted in less O3 production over the mountain ranges.Situated on the southern section of the Black Hills, WICA received the highest summertime average precipitation of 0.27 kg m -2 (Figure 9a).Total Precipitation showed a significant negative correlation with summertime median DM8HA O3 at WICA (r = -0.65,p = 0.04, Figure 9b).As noted in Section S7, the decadal highest summertime median DM8HA O3 was observed in 2012 at 7 sites because of the decadal maximum wild fire emissions.WICA was the exception with its 2012 summertime median DM8HA O3 being the second largest (58 ppbv) of the decade, compared to the largest value of 61 ppbv in 2006 (Figure 3b), despite the fact that the decadal maximum TFI occurred in 2012 (0.022 g NOx m -3 in Figure 9b).O3 levels were expected to be lower in the summer of 2006 because it had the least influence of wildfire emissions evidenced by a very low TFI value (0.010 g NOx m -3 ) (Figure 9b).However, the decadal minimum precipitation (0.18 kg m -2 ) in the summer of 2006 counteracted the effect of less wildfire emissions (Figure 9b).By comparison, increased cloudiness accompanied more precipitation (0.22 kg m -2 ) at WICA during the summer of 2012, which dominated over the influence of the decadal maximum wildfire emissions (TFI = 0.021 g NOx m -3 ).An increasing trend of 6.7 g m -2 yr -1 (p = 0.17) was found in total precipitation at WICA for the time period of 2005 -2015, indicative of increasing cloudiness and ultimately contributing to the decreasing trend in the A4DM8HA over the decade at this site.
Compared to YELL, CNTL, GTHC, ROMO, and WICA, sites located at the lower parts of the basins (i.e., CANY, DINO, MEEK, RANG) had lower annual precipitation levels ranging from 0.08 -0.12 kg m -2 over 2005 -2015 (Figure 9a).Significant negative correlations were only found between summertime 5 th percentile DM8HA O3 and relative humidity at CANY, MEEK, RANG (Table S6).This indicates a stronger impact of cloudiness associated with precipitation weather on summertime baseline O3 levels at sites within the basins.

In situ photochemical O3 production in winter and summer
To determine the influence of in situ photochemical O3 production on ambient O3 in/near O&NG basins, box model simulations were conducted for the Uintah-Piceance Basin, Denver-Julesburg Basin, and its downwind region.In particular, field measurements of VOCs and NOx in the Uintah Basin during UBWOS 2012 -2014 were used to understand the role of in situ During NACHTT, photochemically produced O3 at BAO ranged from 6 ppbv at 6:00 am MST to 22 ppbv at 15:00 pm MST, contributing 20 -52% of the observed O3 (Figure 10a).
Photochemically produced O3 contributed 30 -60% of observed O3 during UBWOS 2012 (Figure 10b).As seen in Figures 10c and 10d, it is possible for the contribution of photochemically produced O3 during UBWOS 2013 and UBWOS 2014 to exceed the observations.For example, a modeled maximum value of 111 ppbv photochemically produced O3 was simulated for 16:00 pm MST during UBWOS 2013, larger than the 96 ppbv observed at the site, which indicates that instantaneous mixing with surrounding air and fast transport lowered the photochemically produced O3 to the observed level.
The decadal highest winter seasonal median DM8HA O3 was found at DINO (59 ppbv), MEEK (44 ppbv), and RANG (44 ppbv) in winter 2013 (Figure 3c).A positive AO phase was encountered during the winter of 2012 (Figure 5b), where the ambient conditions were relatively warm without snow cover (Edwards et al., 2013), and the resulting conditions exhibited only moderate O3 production (~31 ppbv in Figure 10b).In contrast to 2012, winter 2013 was in a negative AO phase (Figure 6b), placing the Intermountain West under a strong influence of frigid Arctic air (Section 5) and subsequently under meteorologically stagnant conditions and increased snow cover within the basins.As a result, the high levels of VOC emissions from O&NG fields that accumulated in the shallow boundary layer coupled with the increased photolysis rates from the snow albedo contributed to rapid O3 production (~78 ppbv in Figure 10c) within the basin, consistent with results from Edward et al. (2014).This result supports the point that surface O3 at the sites within the basins was attributed mostly to in situ photochemical reactions (Section 5).
In contrast to the four wintertime campaigns described above, FRAPPÉ was conducted in the summer.The O3 measurements from ROMO, located in the center of the FRAPPÉ campaign area, showed a diel cycle with ~40 ppbv nighttime and ~49 ppbv daytime O3 levels.In contrast, photochemically produced O3 at ROMO showed a distinct diel cycle with 0 ppbv at night and ~41 ppbv during the day.The strong diel cycle of photochemical production of O3 at ROMO probably resulted from the combined effect of its high elevation and close proximity to the NFR.
During the day, large amounts of O3 precursors from the NFR transported upward from the southeast (Figure 5a) and contributed to photochemical O3 production at ROMO (Benedict et al., 2018) followed by dramatically decreased photochemical production accompanied by the collapse of the planetary boundary layer during sunset at ~20:00 MST.

VOCs-NOx-O3 sensitivity
The O3 sensitivity to NOx and VOC emissions, which is relevant to the attainment of regional air quality standards, was evaluated in/near O&NG basins.Fifteen simulations for each field campaign were performed with scaled NOx mixing ratios and the base case VOC scenario to test the sensitivity of daily maximum photochemical O3 production.The O3 formation regime was found be associated to the AO.Photochemically produced O3 during UBWOS 2012 and UBWOS 2014 was sensitive to VOCs (Figure 10f), indicating that O3 photochemistry at DINO, RANG, and MEEK was also likely to be radical limited in winter 2012 and 2014.In contrast to these two winters, O3 photochemistry had mixed sensitivity during UBWOS 2013 (Figure 10f), when in situ photochemically produced O3 could be as high as 111 ppbv.Winter 2013 was in a negative AO phase, while winters of 2012 and 2014 were in positive phases.As stated in Section 5, surface O3 at the sites within the basin was attributed mostly to in situ photochemical reactions.Therefore, for areas within the basins, emission reductions in VOCs alone would lead to O3 mitigation in positive AO years, while emission reductions in both VOCs and NOx would be effective in negative AO years.In the NFR region, photochemical O3 production in winter 2011 surrounding BAO was sensitive to VOCs during NACHTT (Figure 10f).Photochemical O3 production at ROMO in summer 2014 had mixed sensitivity (Figure 10f).Edwards et al. (2013) also simulated an average UBWOS 2012 day (15 January to 1 March 2012) and found the same O3 formation regime as our study.They found that the radicallimited O3 photochemistry in 2012 was driven by the very low radical production rates (~2.3 ppbv day -1 ) in comparison to the emission rate of NOx.While an average UBWOS 2013 day ( 23January to 24 March 2013) was simulated in our study, Edwards et al. (2014) simulated a single stagnant six-day (31 January to 5 February 2013) strong O3 event, when DM8HA O3 increased from 67 to 107 ppbv.It was found that O3 in winter 2013 was sensitive to changes of NOx, and the total net radical production rate could be as high as ~19 ppbv day -1 , which was sufficient to prevent NOx saturation (Edward et al., 2014).In winter, seasonal 50 th /95 th DM8HA O3 was associated with the AO at 8 sites, including the reference site YELL.For sites within the O&NG basins, O&NG emissions, colder surface air temperature, and enhanced solar radiation due to snow cover contributed to higher O3, especially in the seasonal 95 th levels during negative AO years, and emission reductions in both VOCs and NOx could lead to effective O3 mitigation.For sites outside the O&NG basins, the high O3 in seasonal 50 th levels were a result of transport from the Arctic or California and stratospheric intrusion.In summer, the interannual variation of O3 was predominantly affected by precipitation weather at 9 sites including the two reference sites.At high-elevation sites in the mountains, more abundant precipitation (0.18 -0.27 kg m -2 ) induced more cloudiness and consequently reduced solar flux leading less O3 production in seasonal 5 th /50 th /95 th levels.However, at sites within the basin, the cloudiness associated with precipitation weather had a stronger impact on summertime baseline O3 levels (i.e., the 5 th percentile).

Summary
This study is the first one to investigate the long term impact of O&NG extraction activities on the distribution and trend of surface O3 over the intermountain U.S. O&NG activities have varied greatly over the past decade.In Wyoming, annual O&NG production levels have fallen significantly in the Jonah Field since 2009 and the Pinedale Field since 2012, because of lower natural gas prices relative to crude oil (http://wogcc.wyo.gov/).In Colorado, the number of active wells in Weld County increased by ~2000 between 2012 and 2014, followed by a decline since early 2015 (http://cogcc.state.co.us/data.html).However, current emission inventories underestimated VOC emissions from O&NG productions by a factor of 2 or more (Pé tron et al., 2014).While data analysis studies provided measurement-based estimates of contributions from various processes to decadal variability of surface O3, we understand their limitation in clearly separating such contributions.Detailed multiyear studies that integrate measurements and three-dimensional chemical transport model simulations with accurate emission inventories are needed to further quantify the contribution of each process identified in this work.2006, 2008, 2010, 2011, 2013) and low O3 years (2007, 2009, 2012, 2014, 2015).2006,2008,2010,2011,2013) and low O3 years (2007,2009,2012,2014,2015).Red stars indicate study sites.
).Specifically, Mesa Verde National Park (MEVE), Pinedale, WY (PNDE), Centennial, WY (CNTL), Gothic, CO (GTHC), Rocky Mountain National Park (ROMO), and Wind Cave National Park (WICA) are located near O&NG Basins, while Canyonlands National Park (CANY), Dinosaur National Monument (DINO), Rangely, CO (RANG), Meeker, CO (MEEK), and Campbell, WY (CAMP) are located within O&NG basins.Of the remaining 7 sites, Yellowstone National Park (YELL) and Craters of the Moon National Monument and Preserve (CRMO) were the only two sites located upwind O&NG fields and were therefore used as reference sites to investigate the decadal O3 change with minimal influence of O&NG emissions.Four sites, Great Basin National Park (GRBA), Zion National Park (ZION), Grand Canyon National Park (GRCA), and Petrified Forest National Park (PEFO), were found, based on backward trajectory cluster analysis, to often receive air masses influenced by nearby anthropogenic emissions from Las Vegas (Section S2).Badlands National Park (BADL), located more than 100 km away from the shale play, was impacted one third of the time by air masses from O&NG activity in the Powder River Basin during the 2005 -2015 decade (Section S2).
Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2019-164Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 19 March 2019 c Author(s) 2019.CC BY 4.0 License.basins (Figure 1) and exhibited no trends in A4DM8HA after 2005.In comparison, decadal trends of O3 at the 11 sites situated in or near the O&NG extraction basins ranged widely, indicating complex influence of varying changes in emissions from O&NG production activities combined with other anthropogenic sources.The strong decreasing trends in A4DM8HA at MEVE and CANY appeared to be driven by significant reductions of emissions from local industry during 2005 -2015.MEVE is located near the San Juan Basin (Figure 1), the second largest natural gas field in the U.S. A 37% decrease in natural gas production was reported in the San Juan Basin over 2005 -2015 (Figure 4a), subsequently leading to reduced emissions of O3 precursors.Meanwhile NOx emissions from other sectors in San Juan County decreased by 65% from 80,734 tons in 2005 to 27,996 tons in 2014 (Table in the A4DM8HA at MEVE (-0.76 ppbv yr -1 ) during 2005 -2015.The significant decreasing trend in the A4DM8HA observed at CANY was possibly a result of the declining extraction of coalbed methane (CBM) in the Paradox Basin and emission reductions in coal-fired electricity generation in Emery County.Historically, Emery County in the Paradox Basin was the main CBM producer.Permitting activities for CBM extraction have declined by 48% over 2005 -2015 resulting from lower natural gas prices (NGI, 2018) (Figure 4a).Meanwhile, NOx Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2019-164Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 19 March 2019 c Author(s) 2019.CC BY 4.0 License.emissions from coal-fired electricity generation, contributing ~88% of the total NOx emissions in Emery County, decreased by 35% from 28,407 tons in 2005 to 18,336 tons in 2014 (Table Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2019-164Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 19 March 2019 c Author(s) 2019.CC BY 4.0 License.O3 (64.1 ppbv) was observed, and trends were not examined for the site's short data record (6 yrs).
photochemistry in determining observed O3 at RANG, DINO, and MEEK during winters of 2012 -2014 (Figure 1).Measurements during FRAPPÉ at ROMO were used to determine photochemically produced O3 in summer 2014, while the NACHTT campaign at the Boulder Atmospheric Observatory (BAO) Tower was used for winter 2011.Constraining the model concentrations of VOCs, NOx, and the radical precursors HCHO, HONO, and ClNO2 using the observed diel profiles reduced the impact of turbulent mixing and transport processes on the in situ photochemical O3 production.Figure 10 compares the observed average diel cycle of O3 mixing ratios with those simulated for each campaign.Differences between observed and simulated diel O3 profiles represent the combined impact of transport processes, stratospheric intrusion, and model uncertainty on surface O3.
Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2019-164Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 19 March 2019 c Author(s) 2019.CC BY 4.0 License.Our study suggested that decadal trends in the A4DM8HA over the Intermountain West during 2005 -2015 were shaped by precipitation weather, wildfire emissions, stratospheric intrusion, and transport from the Arctic and the West Coast facilitated by regional to hemispheric circulation.Trends over the areas near/within O&NG basins were predominantly impacted by changes in emissions from the O&NG basins among other anthropogenic sources.Two reference sites, YELL and CRMO, exhibited no significant trends in the A4DM8HA O3 over 2005 -2015, indicating, without the immediate influence of O&NG emissions, relatively constant O3 levels over the Intermountain West.In contrast, decadal trends in the A4DM8HA varied at the 11 sites located either in or near the O&NG extraction basins.A significant decreasing trend of -0.76ppbv yr -1 at MEVE was associated with the decreasing natural gas production (37%) in the San Juan Basin, while a significant decreasing trend (-0.54 ppbv yr -1 ) at CANY was predominately influenced by declining extraction of CBM in the Paradox Basin and emission reductions of 35% in electricity generation from coal.Increasing precipitation (6.7 g m -2 yr -1 ) resulting from increasing cloudiness in the southern portion of the Black Hills contributed to the decreasing trend in the A4DM8HA (-1.21 ppbv yr -1 ) at WICA.No trends were found in the 4DM8HA at the remaining sites, resulting likely from the combined effect of increasing emissions from O&NG extraction and decreasing emissions from other sectors.

Figure 1 .
Figure 1.Topographic map of the mountain states with tight O&NG plays in the basins highlighted in yellow (https://www.eia.gov/).Eleven sites are located within 100km of a shale play (red dots), while seven sites are located at the periphery of the shale play (black dots).Blue dots show two campaign locations, Horse Pool and Boulder Atmospheric Observatory (BAO).

Figure 2 .
Figure 2. Time series of the A4DM8HA at each site.Dashed lines indicate the current ground-level O 3 standard (70 ppbv).

Figure 3 .Figure 5 .Figure 6 .
Figure 3.Time series of seasonal median values of DM8HA O3 at each site in summer (a-b) and winter (c-d).

Figure 10
Figure 9. (a) Decadal average summertime total precipitation (kg m -2 ) over 2005 -2015, and the sites with low precipitation levels highlighted in red.(b) Time series of summer seasonal median DM8HA O 3 , convective precipitation, total precipitation, and total fire index at WICA.