Changes in summertime ozone in Colorado during 2000 – 2015

In 2016, the Denver Metro Area (DMA)/Northern Colorado Front Range (NCFR) was reclassified from a Marginal to a Moderate O3 Non-Attainment Area due to the prevalence of high summer ozone (O3) occurrences. Hourly surface O3 data collected during 2000–2015 from a total of 80 monitoring sites in the State of Colorado were investigated for geographical features in O3 behavior and O3 changes over time. We particularly focus on summer O3 (June, July, August), which is the time when most exceedances of the O3 National Ambient Air Quality Standard (NAAQS) have been recorded. Variables investigated include the statistical (5th, 50th (median), and 95th percentile) distribution of O3 mixing ratios, diurnal amplitudes, and their trends. Trend analyses were conducted for 20 site records that had at least ten years of data. The majority of Colorado ozone monitoring sites show an increase of the 5th (16 total; 11 of these are statistically significant (p-value ≤ 0.05) trends) and 50th (15 total; 4 statistically significant trends) percentile values. Changes for the 95th percentile values were smaller and less consistent. One site showed a statistically significant declining trend, and one site an increasing trend; the majority of other sites had slightly negative, albeit not statistically significant declining O3. Ozone changes at the two highest elevations sites (>2500 m asl) are all negative, contrasting increasing O3 at U.S. West Coast sites. NCFR urban sites do not show the rate of decreasing higher percentile O3 as seen for the majority of urban areas across the U.S. during the past 1–2 decades. The amplitudes of diurnal O3 cycles were studied as a proxy for nitrogen oxides (NOx) emissions and the diurnal O3 production chemistry. The majority of sites show a decrease in the median summer O3 diurnal amplitude (19 total/10 statistically significant). This is mostly driven by the increase in nighttime O3 minima, which is most likely a sign for a declining rate of nighttime O3 loss from titration with nitric oxide (NO), indicating a change in O3 behavior from declining NOx emissions. Since median and upper percentile surface O3 values in the DMA have not declined at the rates seen in other western U.S. regions, thus far the reduction in NOx has had a more pronounced effect on the lower percentile O3 distribution than on high O3 occurrences that primarily determine air quality. An assessment of the influence of oil and gas emissions on Colorado, and in particular DMA O3, is hampered by the sparsity of monitoring within oil and gas basins. Continuous, long-term, high quality, and co-located O3, NOx, and VOC monitoring are recommended for elucidating the geographical heterogeneity of O3 precursors, their changing emissions, and for evaluation of the effectiveness of O3 air quality regulations.


Introduction
Surface ozone (O 3 ), first identified in the 1940s and 1950s as an air pollutant that adversely impacts vegetation, human health, and crop yield (National Research Council, 1991;Lefohn et al., 1997;The Royal Society, 2008;Jerrett et al., 2009), is a widely recognized air quality problem in many regions around the world (Dentener et al., 2010).The globally averaged O 3 tropospheric lifetime is approximately 23 days (Young et al., 2013), with the O 3 lifetime being shorter near the surface and in urban areas because of deposition and chemical destruction.The resulting global tropospheric O 3 distribution is highly variable by season, location, and altitude (Cooper et al., 2014).Elevated surface O 3 and exceedances of the O 3 National Ambient Air Quality Standard (NAAQS) are primarily determined by transport of O 3 from outside the region, and local and regional production from photochemistry in the lower troposphere.Emissions and resulting concentrations of the O 3 precursor gases nitrogen oxids (NO x ) and volatile organic compounds (VOC), play a key role in local O 3 production, as well as temperature and sunlight (Haagen-Smit, 1952;McKeen et al., 1991;Ryerson et al., 2001).
In the U.S., compliance with the O 3 NAAQS is determined by the 'Design Value', which is calculated as the mean of three consecutive years of the 4 th highest annual value for the maximum daily average 8-hour O 3 value (MDA8) (CFR, 2017).The O 3 NAAQS has been progressively lowered in recent years in consideration of the increasingly recognized 1.The region has experienced remarkable growth, with an estimated 30% population increase for the NCFR from 2000 to 2015 (Supplemental Materials SM_Table_1).This population growth has been asso-ciated with urban sprawl, increase in transportation, power generation, and associated industrial stationary and mobile emission sources.2. The NCFR is subject to naturally elevated O 3 and downfolding transport events that bring elevated O 3 from the mid to upper troposphere to near the surface due to its elevation of 1500-1600 m asl and location on the eastern flanks of the Rocky Mountains.Occasions when such events led to significant enhancements in surface O 3 and high MDA8 values have been reported (Langford et al., 2009).Musselman and Korfmacher (2014) showed that even remote high elevation areas of the southern Rocky Mountains can experience exceedance of the NAAQS value, in particular during springtime, when stratospheric intrusions contribute to elevated O 3 at these elevations.3. Emissions from oil and natural gas (O&NG) industries are another complicating influence.In 2017, there were over 53,000 active O&NG wells in Colorado.
The O&NG development is concentrated in a number of lower elevation basins (Figure 1),with 22,000 wells located in Weld County in the Denver Julesburg Basin (DJB), and thereby in the NAA (COGCC, 2014).
Several recent studies have investigated the influence of O&NG emissions on Colorado O 3 .Building on the speciated measurements of O&NG-associated VOC, both Swarthout et al. (2013) and Gilman et al. (2013) found that emissions from the O&NG industry alone contribute to more than half of the local O 3 production potential at the Boulder Atmospheric Observatory (BAO) in northern Denver.McDuffie et al. (2016), applying a box model, found that O&NG operations contribute approximately 20% to regional photochemical O 3 production in the NCFR.These inferred influences have been confirmed by O 3 monitoring results: A study conducted by the CDPHE in 2008, using O 3 data from four sites along the NCFR, showed that for elevated O 3 events during mid-May to mid-August 2006, air transport from the DJB was associated with the highest O 3 values, whereas transport from surrounding areas, including the DMA, brought in air with lower O 3 levels (CDPHE, 2008b).Using a correlation analysis of ambient O 3 observations and winds, Evans and Helmig (2017) found that for the BAO and a site south of Boulder, during 2009-2012, an average of 65% of elevated O 3 events were associated with transport from O&NG production regions.Based on case study comparisons near Greeley during the 2014 summer, Cheadle et al. (2017) estimated that O&NG emissions contribute up to ~20 ppbv to O 3 production on high O 3 days.
Applying modeling of the same year's data for the wider region, Pfister et al. (2017) concluded that, on average, O&NG emissions contribute to 30-40% of the O 3 produced on high O 3 days.4. During recent years, the region has been subjected to an unusual frequency and extent of nearby wildfires and transported wildfire plumes from outside of the state.These wildfire emissions are contributing to occurrences of elevated O 3 production.Given the irregular and sporadic occurrences, thus far there is no quantitative assessment of the wildfire contribution to Colorado O 3 .1 characterizes sites as rural, suburban, and urban.The amount of data varied widely among sites, from ~2 months to a maximum of 15 years; the fraction of data that were available for each site and year is summarized in SM_Table_2.

Data processing
Data were downloaded as zip files, converted to .xlsxfiles, and then imported into Matlab.All data are reported as a molar ratio in air, i.e. units of parts per billion by volume (ppbv [nmol mol −1 ]).All data were converted to Mountain Standard Time.A few stray data values <−5 ppbv were removed from records.During the review of this manuscript, deviations in the O 3 data for the southwestern Colorado sites Bondad and Ignacio were pointed out; after further investigation, the 2000-2004 data for these two sites were excluded from further analyses (SM_Text_1).Data were plotted as time series, and boxwhisker plots were prepared when at least 80% of data were available for a given year.For box-whisker plots, the center line represents the data median, and the upper and lower horizontal bars of the box the 25 and 75 percentiles of the data, the whiskers are drawn to a length of 1.5 times the interquartile range away from the top and bottom of the boxes, and values outside the whisker range are marked as red crosses.As we particularly focus on summer O 3 (June, July, August), when most exceedances of the O 3 NAAQS occur, hourly O 3 values from 1 st June (12:00 am) until 31 st August (11:00 pm) were extracted and used to compute the summer 5 th , 50 th , and 95 th O 3 percentile values for each year.When at least 80% of summer data were available and a site had at least 10 years of data, trends were calculated for the 5 th , 50 th , and 95 th percentiles with a linear regression fit.The statistical significance of the linear regression result was evaluated from the p-value, tested against the null hypothesis, R 2 = 0, which indicates no linear relationship by using the standard F-statistic.The null hypothesis was rejected with a confidence level ≥95% when the probability P associated with the F statistic p-value was ≤0.05.For the MDA8 and Design Value calculation, the Code of Federal Regulations (2017) protocol was used, with a linear regression line fit through the calculated annual data points for each year.In the maps showing the geographical distribution of trend analyses, slope values are given by the color scale, and the significance of each trend result is represented by the size of the marker: Statistically significant trend results are shown in 4 times larger size than statistically non-significant values.
For calculation of daily amplitudes, the hourly minimum O 3 value occurring in the first 12 hours of a day was subtracted from the following daily maximum.Amplitudes were only calculated when at least 20 hourly data points were available for a given day.The seasonal change of O 3 diurnal amplitude during the year was calculated from the 50 th percentile amplitude for every month when at least 80% of data were available.For each site with at least 10 years of data, box-whisker plots show the year-to-year variability and change of diurnal amplitudes.As for the O 3 mixing ratio, we focused again on the summer months' O 3 amplitudes, and computed the 5 th , 50 th , and 95 th percentile values for each summer when at least 80% of data were available.For sites with at least 10 years of data, a linear regression fit was calculated to evaluate trends in the summer amplitude.
Additionally, all summer data that were available for the 5-year period of 2011-2015 (plus three 2007-2008 Boulder County site data sets) were used to calculate the 50 th and 95 th percentile summer O 3 mixing ratio, and the median summer daily amplitude.In the maps displaying these results, the marker area is proportional to the square of the amount of available data.

O 3 summer mixing ratios
The full time series records for all 80 sites are available in SM_Figure_1.69 data sets have data for 2011-2015.The geographic distribution of the 50 th and 95 th percentile summer O 3 mixing ratios for this time window is visualized in Figure 3.The amount of data considered for each point is reflected by the marker size; the smallest marker, i.e. site No. 52 (Elk Springs) represents 735 hourly data points, and the largest marker, site No. 47 (Palisade), had 10977 hourly data points.Overall, summer median O 3 results span a range from 32-62 ppb across these Colorado sites (Table 2).Within the DMA, relatively low median mixing ratios of ~35 ppbv were found, whereas relatively higher values were calculated for suburban Denver (43-47 ppbv).The highest summer median O 3 mixing ratios were seen at the high elevation mountain sites, with the highest elevation site, Mount Evans (4300 m), having the highest value of 62 ppbv, ~30 ppbv above those for the DMA sites.
A different spatial pattern was found for the 95 th percentile summer O 3 mixing ratios, where the the high elevation locations on the eastern slopes and the greater DMA sites exhibited the highest values.Most of the DMA sites in the blue enlargement box have values greater than 67 ppbv.95 th percentile values were on the order of 20 ppbv or less above the 50 th percentile vales for many of the rural low elevation, and for high elevation sites.In contrast, a more dynamic O 3 range was found for sites located in the NCFR.These sites all show a wider distribution, with 95 th percentile values exceeding the median by ~30 ppbv or more.The relationship between the 50 th and 95 th percentile mixing ratios and the site altitude is shown in Figure 4.There is a statistically significant increase of the summer O 3 50 th percentile value with altitude of 6.5 ± 1.6 ppbv km −1 (p-value < 0.05).In contrast, the 95 th percentile O 3 values are much more variable and do not show a statistically significant dependence on elevation.

O 3 hourly time series data
Three sites with different characteristics were chosen for illustrating the data records and emphasizing main features; their results are displayed in Figure 5: 1. Mesa Verde National Park (MVNP), in the southwestern part of Colorado, representing a rural reference site, 2. Welby (WEL), a suburban site located northeast of Denver, at the southern edge of the DJB, and 3. the Golden -National Renewable Energy Laboratory (NREL) site, close to and west of Denver, in a densely populated area within the NCFR O 3 NAA.
Of the sites shown here, MVNP is the one with the lowest total O 3 variability over the year, indicated by the  Summer O 3 trends were investigated by conducting a linear regression fit to the entire summer data record.Results for the three sample sites are shown in Figure 6; graphs for all other sites are available in SM_Figure_2.MVNP shows signs for decreasing O 3 at all percentile values, with the 95 th percentile value of −0.35 ± 0.34 ppbv yr −1 being a statistically significant negative trend.The O 3 mixing ratios at NREL show a steeply increasing trend of 0.55 ± 0.44 ppbv yr −1 for the 5 th percentile values.Changes in the median values (increasing) and 95 th percentile values (decreasing) are not statistically significant.The WEL results are characterized by statistically significant increasing trends for the 5 th , 50 th , and 95 th percentiles.Details of the trend analyses results for sites that had long enough records to fulfill the selection criteria are listed in Table 3. Maps showing the geographical distribution of trend results for the 5 th , 50 th , and 95 th values are given in Figure 7.In Table 3 we also include trends for the annual 4 th highest MDA8 and Design Value, (time series graphs for the Design Value and MDA8 trend determinations are shown in SM_Figure_3) for comparison of those regulatory metrics with the 95 th percentile results that are presented throughout this paper.Given the longer data coverage needed for calculation of the Design Value (i.e. a three-year running mean through continuous data), this calculation was only possible for a subset of the site records.
For the 5 th percentile values, most sites (16 out of 20) show increasing O 3 values, with 11 of these passing the statistical significance test.Five of these results are >0.63 ppbv yr −1 , equivalent to a remarkable >10 ppbv increase in the 5 th percentile O 3 value over the 16-year record.For the 50 th percentile, 15 sites show increasing O 3 results, 4 being statistically significant.For the 95 th percentile values most (15) sites show a decrease in O 3 ,   .At all sites, diurnal amplitudes are highest during the summer months.Absolute values differ largely between sites, with MVNP having the lowest values, ranging from (median) i.e. ~10 ppbv in winter to ~20 ppbv in summer.These are on the order of 30-50% of those observed at NREL.WEL has an even larger seasonal cycle than NREL, and interestingly, the seasonal cycle shows a wider spread summer maximum.

Daily summer amplitude for 2011-2015
The same sites for which the 50 th and 95 th percentile summer O 3 mixing ratios were determined (Figure 3) were subjected to a statistical analysis of the summer diurnal O 3 amplitude: Tabulated results were included in Table 2, and the geographical distribution is shown in  Results from WEL contrast the other two sites.Here, amplitudes have been increasing at all chosen percentile values; but none of these changes are statistically significant.
The results for all datasets are listed in Table 4.The geographical distribution of the slopes is shown in Figure 12.Notably, the majority of sites exhibit declining amplitudes, with those being reflected at all percentile values.Ten of the negative median slope values passed the statistical significance test.Negative values were determined for most of the sites in the DMA and NCFR.WEL is the only DMA site that shows positive slopes at all percentile values.

Year-to-year variability
Surface O 3 is dependent on numerous influences, causing significant year-to-year variability (Strode et al., 2015).The range (spread), representing the year-to-year variability, of the median O 3 summer values within individual records was between 8-18 ppbv for these Colorado sites.Study of causal relationships determining O 3 have been an active field of research, motivated by the necessity to elucidate the effect of human influences.Processes determining year-to-year changes include prevalence of high pressure weather systems (Reddy and Pfister, 2016), occurrences of stratospheric intrusions (Langford et al., 2009;Langford et al., 2015), wildfires (Jaffe et al., 2008;Jaffe and Wigder, 2012;Jaffe et al., 2013;Baylon et al., 2016;Lu et al., 2016), long-range to trans-continental transport events (Jaffe  , 2007;Zhang et al., 2008;Reidmiller et al., 2009;Zhang et al., 2009), and changes in anthropogenic O 3 precursor emissions.In Colorado, elevation and proximity to urban centers are two factors further influencing absolute levels and variability.The individual contribution of these effects on the O 3 data considered here is beyond the objective of this study.However, recognition of these multiple influencing factors is of importance for evaluating the trend results presented here.While the 10-year minimum cutoff and 95% p-test confidence interval provide robust trend results, it cannot be completely excluded that trend results are influenced by multi-year extremes of any or a combination of the above listed factors.

Dependency of ozone on elevation
Highest O 3 median values were observed at the high elevation sites.The summer elevation gradient of 6.5 ± 1.6 ppbv km −1 found in the data considered here is significantly lower than the 15 ppb km −1 reported by Brodin et al. (2010) for a ~30 km horizontal/2 km elevation eastwest transect west of Boulder.These different results are likely due to the abundance of considered network sites within the NCFR for the study presented here, and the wide spread seen in the low elevation data (Figure 4).While Brodin et al. (2010) referenced their data to a single downtown Boulder site, where O 3 values were relatively low, many of the lower elevation sites that went into the calculation here are in suburban areas, where median O 3 levels are higher overall because these sites have less nighttime O 3 loss from titration by NO.While median O 3 increases with elevation, there is not as much of a dynamic diurnal O 3 change at the high elevation sites, as a larger fraction of O 3 at these sites is due to long-range transport and the increase of O 3 with elevation.The lower elevation sites, in contrast, are influenced more by local O 3 production, causing larger diurnal O 3 amplitudes.Those two dependencies determine the inverse relationship and scatter in the data in Figure 11.

Ozone trends at rural sites
For rural locations in the U.S., O 3 trends have been inconsistent and smaller than in urban areas.The behavior in the western U.S. contrasts that in eastern regions of the country, where more homogenous, downward O 3 trends have occurred (Cooper et al., 2012(Cooper et al., , 2014)).Results from the analyses presented here provide further evidence for this conclusion, as most of the rural and higher elevation sites showed neutral to only slightly downward summer O 3 for the 50 th and 95 th percentile values.

Ozone trends at high mountain sites
Seasonal studies have reported an increase in O 3 at western U.S. mountain sites (Parrish et al., 2012;Cooper et al., 2014;Gratz et al., 2015).This has mostly been associated to increasing O 3 in transcontinental transport from Asia, with this effect being most influential in spring.Two high mountain sites on the U.S. west coast that have been central in demonstrating O 3 increases stemming from outflow from Asia are Lassen NP (1756 m asl), and Mt Bachelor (2763 m) (Jaffe et al., 2003;Jaffe and Ray, 2007;Parrish et al., 2012;Gratz et al., 2015).Incorporating O 3 sonde data in their analyses, Cooper et al. (2012) calculated that springtime free tropospheric O 3 over western North America increased at a rate of 0.41 ± 0.27 ppbv yr −1 between 1995-2011.Available trends from two sites in the Colorado network above 2500 m asl, i.e.Niwot Ridge C1 (3035 m), and Rocky Mountain NP Longs Peak (2742 m) offer a comparison with these west coast high elevation sites.Remarkably, the mean (of these two, ± standard deviation) 5 th /50 th /95 th percentile O 3 trends are all negative, i.e. −0.23 ± 0.40/−0.33± 28/−0.48 ± 0.11 ppbv yr −1 , and inconsistent with the observations from the U.S. west coast.This difference is probably due to the high elevation west coast sites distribution of the summer amplitude data for every year that had at least 80% of data available for a given year.Linear regression fits trough the 5 th , 50 th and 95 th percentile of the data for each year show the trend behaviour.Regression analysis slope results (in ppbv yr −1 ) are given in the legends.DOI: https://doi.org/10.1525/elementa.300.f8being relatively isolated from North American influence, whereas the two Rocky Mountain sites are more responsive to the changes in background ozone in continental North America and the NCFR.

Ozone regulatory considerations
There are two groups of sites with MDA8 and Design Values over the >70 ppbv threshold (SM_Table_3, and SM_Figure_7 showing the geographical distribution): 1. high elevation mountain sites and 2. suburban NCFR sites.Sunlight Mountain, the 7 th highest elevation site at 3225 m, has the highest average number of exceedances per year (32).This site record, however, is relatively short (3 years), which makes these data more susceptible to biases from interannual variability, and places a relatively high uncertainty on this ranking.There are several other high mountain sites where the MDA8 regularly exceeds 70 ppbv, such as Niwot Ridge, and Rocky Mountain NP Longs Peak.These observations are in line with previous elaborations that pointed out the high O 3 at these mountain sites, and the significant contribution from longrange continental and upper troposphere/lower stratosphere transport on elevated O 3 (Langford et al., 2009;Lin et al., 2012;Cooper et al., 2015).The second group are NCFR sites, and here, in particular, the suburban locations on the western periphery of the City of Denver.This chain of stations along the eastern foothills of the Rocky Mountains shows a consistent behavior of a year after year high number of high O 3 days, with an average of 10-27 days per year when the O 3 MDA8 exceeds 70 ppbv.
A contributing factor appears to be the diurnal air flow pattern along the NCFR, where predominantly upslope flow conditions during mid-day to late afternoon bring air that has been enriched in O 3 precursors towards the mountain slopes, causing elevated O 3 in the afternoon at the Plains-Rocky Mountain transition zone (Evans and Helmig, 2017).The analyses presented here expand upon the observational study from two sites near Boulder (Evans and Helmig, 2017), and recent modeling work (Pfister et al., 2017), demonstrating this spatial elevated O 3 distribution along a ~30 km wide belt extending some 150 km north-south (Fort Collins to Chatfield) along the NCFR foothills.

Ozone changes in relation to NO x and VOC
Mandated by the Clean Air Act, emission reductions of O 3 precursors in NAAs have been remarkably successful in many U.S. urban regions, in particular in the mid-western U.S. (Lurmann et al., 2015;Sather and Cavender, 2016).For instance, reductions for the MDA8 ranging from 18-38 ppbv have been reported for Baton Rouge, Houston, El Paso, Houston, and Dallas-Fort Worth over the past 30 years.The most dramatic improvements have been achieved in Southern California, where on average the O 3 MDA8 has decreased by 2.6 ± 0.2% yr −1 between 1973-2010 (Pollack et al., 2013).Nationwide, a 17% decrease has been reported for 2000-2015, with a regional decrease of 10% for the Southwest in 2000-2015, which includes the state of Colorado (EPA, 2016).Averaging the O 3 trend results for DMA sites (#1, 2, 5, 14, 18, 33, 34, 35, 36) yields median 50 th /95 th summer O 3 trends of 0.29 ± 0.41/−0.01± 0.34 ppbv yr −1 .None of these results are indicative of  3, Figure 7).These changes in the distribution of absolute mixing ratios are paralleled by declining diurnal O 3 amplitudes (Table 4, Figure 8).Since nighttime O 3 minima are largely determined by O 3 loss from titration with NO, the changes in nighttime O 3 are indicative of gradually declining NO x concentrations in the urban corridor.Similar, albeit weaker, nighttime O 3 changes are observed at the rural Colorado locations, suggesting a declining NO x signature on the regional scale as well.In these rural areas, increasing nighttime O 3 levels are likely due to 1. a declining nighttime O 3 destruction in air transported to the sites, and 2. a decline in the local nighttime O 3 destruction.These indirectly inferred changes in NO x are confirmed by surface NO x monitoring data in the state.All five of the 12 records overall (SM_Figure_8) that are long enough for a trend analysis show declining NO x (SM_ Figure_9) at all percentiles, with 11 out of these 15 NO x percentile changes being statistically significant downward trends (SM_Table_4).An overall 26.3 ± 5.4% 2005-2013 reduction in tropospheric NO x has been estimated from satellite observations for the Denver area (Lamsal et al., 2015;Lurmann et al., 2015), with indications that the rate of decline has become smaller during the most recent years of the record.These declines in Colorado NO x are consistent with results from other nationwide studies, building on surface monitoring, remote sensing, emission inventories, and modeling (Cooper et al., 2014;Lamsal et , 2015;Strode et al., 2015;Lin et al., 2016;Sather and Cavender, 2016).
Changes in VOC and VOC reactivity are more difficult to assess than changes in NO x .Tailpipe emissions, constituting a major to dominant source of VOC in urban air, have been decreasing throughout U.S. cities, including in the DMA (Bishop and Stedman, 2008).Despite the growth of the automobile fleet and miles driven per vehicle, the signature of these declining emissions can be seen in declining urban air VOC (Warneke et al., 2012;Rossabi and Helmig, 2018).In the NCFR, VOC emissions from the O&NG industries are the dominant VOC source within and downwind of O&NG activities (Gilman et al., 2013;Swarthout et al., 2013;Thompson et al., 2014).Several studies have demonstrated increasing atmospheric VOC trends linked to O&NG fugitive emissions from the rapid growth of this industry across the U.S. (Schade and Roest, 2015;Vinciguerra et al., 2015;Helmig et al., 2016).Within Colorado, CDPHE has been conducting VOC monitoring at four locations.While some of these records indicate potentially declining VOC, at this time these records are too short, and/or the sampling has been inconsistent, and data are too variable for obtaining robust trend results.
In Figure 13 we plot available paired median summer O 3 amplitudes against the summer median NO x values from the records shown in SM_Figure_8 to investigate the O 3 -NO x /VOC relationship further.Data markers for WEL and Denver Camp (DCP) are color-coded by the year of observation.There are several dependencies: 1. Lowest NO x and diurnal amplitudes were observed at rural and suburban locations; highest values were measured at WEL and DCP. 2. At sites with lower NO x (0-5 ppbv summer median), the median O 3 amplitude scales close to linearly with increasing NO x .This behavior then reverses to a slightly negative relationship between 10-25 ppbv of NO x .3. Despite DCP having the overall highest NO x , at equivalent NO x , O 3 amplitudes at DCP are consistently lower (by ~10 ppbv) in comparison to WEL. 4. The colorcoding of the WEL and DCP data shows that while for both sites NO x has been decreasing over the available record, O 3 amplitudes have changed very little, possibly increasing slightly at WEL.

Ozone at Welby
Data in Figure 13 point towards overall increasing diurnal O 3 changes at WEL.Given the decreasing rate of nighttime O 3 destruction by NO (as reflected in the statistically significant positive 5 th percentile O 3 trend), these data argue for an increase in O 3 production chemistry at WEL.This conclusion is confirmed by the O 3 trend analyses: WEL is the only site with (statistically significant) increasing O 3 trends at all (5/50/95) percentiles (Table 3), and the only site where O 3 diurnal amplitudes have been increasing at all percentiles.This behavior clearly shows that, thus far, there have been no improvements in O 3 at WEL despite the achieved statewide reductions in NO x .Conversely, it appears that local O 3 production has increased under declining NO x .Similar conclusions were derived from a recent analysis of trends in weekday/weekend differences in O 3 and NO x at selected NCFR sites by Abeleira and Farmer (2017).Their study found that "most sites in the NFRMA were NO x -saturated, but are transitioning to, and in one case may already have reached, the peak P(O 3 ) cross-over point between NO x -saturated and NO x -limited regimes".Our analyses, relying on further site data and different analytical approaches, confirm the first part of their conclusion.

Ozone and O&NG development
WEL is situated on the northeastern outskirts of the DMA, within 2 miles of I-76, and with other busy roadways a bit further away.The site is also within the southwestern edge of the DJB (Figure 1).A map showing well pad locations surrounding the site is provided in SM_Figure_10; 20 well pads are located within a 10-km distance, and 1445 wells within a 20-km radius.The association of increasing O 3 with the proximity to O&NG operations might suggest that emissions from this growing industry have contributed to increases in O 3 .Atmospheric VOC have been shown to increase along transects from the periphery towards the center of the DJB (Thompson et al., 2014;Helmig et al., 2015;Rossabi et al., 2017), and VOC have been reported at highly elevated concentrations, largely driven by O&NG emission within the DJB and downwind (Gilman et al., 2013;Swarthout et al., 2013;Thompson et al., 2014).
Unfortunately, there are no VOC monitoring data available for WEL.The seasonal cycles of the diurnal amplitudes (Figure 9) show a wider summer maximum than at NREL.This may be linked to differences in VOC/NO x ratios at both sites.There is most likely more NO x , causing lower VOC/NO x ratios at NREL, or higher VOC, causing higher VOC/NO x at WEL. Lower VOC/NO x may inhibit O 3 production more during the spring/fall shoulder season when there are fewer photochemically productive hours than during mid-summer (at NREL) in comparison to possibly higher VOC/NO x conditions at WEL. Greeley -Weld County Tower (site #67) is another site within the DJB to further investigate this dependency.The monitoring location is surrounded by intense O&NG development in all directions.While the site, with a summer 95% percentile value of 71 ppbv, and a median diurnal O 3 amplitude of 51 ppbv, is ranked among the top ten of the Colorado sites for highest O 3 and diurnal production, trend analyses for the O 3 higher percentile values and O 3 diurnal amplitudes show signs for decreasing local O 3 production (Tables 3 and 4).VOC monitoring in Platteville ~25 km to the southwest has shown highly elevated O&NG associated VOC (Thompson et al., 2014;Halliday et al., 2016).There are, however, indications of potentially declining O&NG related VOC concentrations during recent years at Platteville (Pierce, 2017).
Other recent research building on campaign observations has provided observational and modeling evidence that O&NG emissions significantly contribute to summer O 3 production within the DJB and NCFR downwind regions (Gilman et al., 2013;Swarthout et al., 2013;McDuffie et al., 2016;Pfister et al., 2017;Cheadle et al., 2017).This association is difficult to prove with the data and tools that were applied here.Unfortunately, the DJB and its periphery are not well covered by the Colorado surface O 3 monitoring network, and regrettably, a valuable site in the SW sector of the DJB, i.e. the Boulder Atmospheric Observatory, was decommissioned in 2017.Further, given the lack of parallel O 3 , VOC, and NO x monitoring, it is impossible to deconvolute the influence of precursor emissions and their changes over time on the contrasting behavior seen in the data from these DJB sites.

Summary and conclusions
Ozone behavior in Colorado is rather diverse, both along the horizontal scale and the ~2000 m elevation change from the Eastern Rocky Mountain Plains to the peaks of the Continental Divide.Elevated O 3 is prominent at higher elevation sites, and in the NCFR.However, in the NCFR, highest O 3 is not recorded within the central DMA, but instead mostly at suburban sites to the north and west of Denver along the eastern foothills of the Colorado Rocky Mountains.
There is high year-to-year variability in O 3 in the NCFR which adds to the uncertainty of surface O 3 trend results.Ozone trends at considered non-urban sites, as well as at the highest mountain sites, do not show increasing O 3 during the observation period.This contradicts the argument that NCFR elevated O 3 can (in part) be explained by increasing O 3 transport into the state.Analyses presented here underscore that NCFR O 3 trends are primarily determined by local/regional emissions and production.
Our analyses for DMA sites do not show the clear O 3 decreases at the higher percentiles as seen for the vast majority of urban areas across the U.S., including the southwestern states during the past 1-2 decades.Only five locations (collocated with O 3 ) in the State of Colorado have NO x data records suitable for trend analyses.Declining trends seen in those NO x data, and the clear signature of decreasing O 3 amplitudes, are in agreement with remote sensing and modeling results.Together, these records provide solid evidence for a steady decline of NO x in the region.The lack of the higher percentile O 3 response to the NO x decline indicates that O 3 in the region has not been sensitive to these NO x reductions thus far.
While we were able to consider data from an impressive 80 O 3 monitoring sites for this study, there is only one site that provides a long-term surface O 3 record from within the DJB, Colorado's largest O&NG basin.Sites on the southern and western periphery of the DJB are subject to influence by urban development and changes in emission from non-O&NG emission sectors, which complicates assessing the O&NG influence and how its importance is changing over time.Further O 3 monitoring within the DJB and at its northern and eastern boundaries would be beneficial for addressing these questions.With the predominant daytime upslope (east to west) air circulation conditions during summer days, those sites could provide data from upwind of the DJB, which, together with the existing chain of sites along the foothills, would allow assessing daytime O 3 production as air travels over the DJB and becomes enriched with O&NG emissions.This would largely enhance capabilities for studying the role of O&NG operations in the NCFR O 3 production, and the effectiveness of regulations for curbing O&NG O 3 precursor emissions.
The sparseness of direct surface NO x observations poses a severe limitation for understanding O 3 responses to emission changes and for directing O 3 policy decisions.Continuous, high quality, year-round observations of VOC and NO x , co-located with O 3 monitoring, are needed for a more comprehensive evaluation of changing influences and regional differences in the Colorado O 3 production chemistry.Here, we explored the use of summer diurnal O 3 amplitudes as an indirect indicator for site comparisons and for NO x trend analyses, and found this variable to be a sensitive and valuable tool for this research.The lack of response in O 3 amplitudes to declining NO x at two DMA sites indicates that O 3 daily production dynamics

Figure 1 :
Figure 1: Location and abundance of oil and natural gas wells in the State of Colorado in 01/01/2017 (data from the Colorado Oil & Gas Conservation Commission, interactive map).Brown dots depict individual well locations.Also shown, by the red lines, are the borders of the DMA/NCFR ozone NAA.Red circles are used as indicators for Mesa Verde National Park, Golden -National Renewable Energy Labs (NREL), and Welby, which are the sites chosen for the data examples shown in this paper.DOI: https://doi.org/10.1525/elementa.300.f1

Figure 2 :Figure 3 :Figure 4 :
Figure 2: Location of ozone sites within the State of Colorado that were considered in this study.Site numbers refer to the listing in Table 1.The color of the dots represent the data-source: Yellow -AQS, blue -INSTAAR, green -NOAA, and red -NPS.Some site numbers were shifted slightly to avoid marker overlap.Red circles show sites which are used for the case study in this paper.DOI: https://doi.org/10.1525/elementa.300.f2

Figure 5 :
Figure 5: Time series examples for the full record of 1-hour ozone at the rural background site Mesa Verde National Park (top), the urban Golden NREL (middle), and Welby (bottom).Statistical analyses of these data with box-whisker plots (see data processing section for box-whisker format explanation), representing each year of data, are shown to the right.DOI: https://doi.org/10.1525/elementa.300.f5 Seasonal cycle of O 3 diurnal amplitudes Time series results for diurnal O 3 amplitudes for the full data records are presented in SM_Figure_4.The full record of O 3 diurnal amplitudes for the case study sites is shown in Figure 8.The same data, binned by month for illustrating the mean seasonal cycle, are shown in Figure 9 (see SM_Figure_5 for results for other sites)

Figure 10 .
Urban areas consistently show higher daily amplitudes than rural areas.Median DMA amplitude values are on the order of 55 ppbv.Values in suburban Denver are somewhat lower, on the order of 35-45 ppbv.Within the inner city, on high O 3 production summer days, the diurnal O 3 increase is on the order of 60-80 ppbv.Rural and high elevation sites show significantly lower diurnal O 3 changes, with median daily amplitude values of 10-15 ppbv for the most remote locations.There are likely two contributing processes determining the larger urban amplitudes.Some portion is due to the low nighttime values from titration with NO, and mixing of air from aloft with higher O 3 during the breakup of the nocturnal surface layer.The remainder is due to daytime photochemical production.The linear regression analyses between amplitude and site elevation yielded statistically significant negative slope relationships for both the 50 th and 95 th percentiles, as shown in Figure 11a, b.The O 3 diurnal amplitude is inversely correlated with median O 3 (Figure 11c).
3.2.3.Summer daily amplitude trend analysesTrend analyses on the diurnal amplitude data were prepared to investigate changes in the summertime dynamic O 3 behavior.Results for the case study sites are shown in the right side set of graphs in Figure8, results for all other sites are presented in SM_Figure_6.At MVNP, summer diurnal amplitudes have been moving towards slightly lower values over the study period.While slope values are small, the median and 95 th percentile results (−0.16 ± 0.14/−0.37± 0.25 ppbv yr −1 ) are statistical significant changes.At NREL, there has been a downward trend of dirnal amplitudes, with statistical significant decreases of the median and 95 th percentile slope values of −0.67 ± 0.37 and -0.94 ± 0.54 ppbv yr −1 , respectively.

Figure 6 :Figure 7 :
Figure 6: Box-whisker plots of the hourly summer ozone distribution.Linear regression fits trough the 5 th , 50 th and 95 th percentile of the data for each year show the trend behavior, with slope results in the legend in ppbv yr −1 .DOI: https://doi.org/10.1525/elementa.300.f6

Figure 8 :
Figure 8: Left column: Ozone time series of the daily amplitude.Right column: Box-whisker plots of the statistical distribution of the summer amplitude data for every year that had at least 80% of data available for a given year.Linear regression fits trough the 5 th , 50 th and 95 th percentile of the data for each year show the trend behaviour.Regression analysis slope results (in ppbv yr −1 ) are given in the legends.DOI: https://doi.org/10.1525/elementa.300.f8

Figure 10 :Figure 11 :
Figure 10: 95 th (top) and 50 th percentile (bottom) summer O 3 daily amplitude of all data available for 2011-2015 (please note that for Longmont (site 7) and Lyons (site 73) data are summer 2007 values).The marker size is a function of the amount of data which was used for the plot (see SM_Table_2).Please note the different color bar scales in the two graphs.DOI: https://doi.org/10.1525/elementa.300.f10

Figure 12 :
Figure 12: Changes in the summer O 3 daily amplitude for the 95 th percentile (top), 50 th percentile (middle), and 5 th percentile (bottom).Shown are all sites with at least 10 years of summer data.The colors of the marker represent the value of the slope.Statistically significant values (p-value ≤ 0.05) are plotted four times larger than statistically not significant values.DOI: https://doi.org/10.1525/elementa.300.f12

Figure 13 :
Figure 13: Ozone median summer daily amplitude as a function of median summer NO x at 12 sites.Data are colorcoded by year according to the scale on the right for Welby and Denver Camp.DOI: https://doi.org/10.1525/elementa.300.f13 Garfield County, the US National Park Service (NPS), the Southern Ute Indian Tribe Reservation (SUIT), Meteorological Solutions Inc. (MSI), the Colorado Bureau of Land Management (CBLM), the City of Aspen, CH2M Hill, as well as from monitoring by the EPA.Additional data from three sites were obtained from the National Oceanic and Atmospheric Administration (NOAA) Global Monitoring Division (GMD).Data from six sites operated during the summer 2014 Front Range Air Pollution and Photochemistry Experiment (FRAPPE), and from the NSF-funded Long-Term Ecological Research (LTER) Niwot Ridge program, were collected by the Institute of Arctic and Alpine Research (INSTAAR) at the University of Colorado, Boulder.Data for Dinosaur National Monument (DNM) were received from the NPS.NO x data for 12 sites were downloaded from AQS. Site information is summarized in Table 1, and Figure 2 depicts a map with all site locations.AQS sites are represented by yellow markers, NOAA sites by green markers, INSTAAR measurements by blue markers, and data obtained directly from the NPS in red.The NPS Dinosaur National Monument (DNM) site coordinates given in the data record placed the site within Colorado, but we later discovered that the actual measurement location is ~22 km west in the State of Utah.Nonetheless, we retained the DNM record for our analyses.The location setting column in Table 2.1.1.Location of study and data sources We were able to retrieve surface O 3 data from a total of 79 sites in Colorado for 2000-2015.Data for 70 sites were downloaded from the EPA Air Quality System (AQS) archive in September 2016.AQS includes Colorado O 3 data from the Colorado Department of Public Health and Environment (CDPHE), the US Forest Service (USFS), Boul-der County Public Health (BCPH),

Table 1 :
Summary and information on the location, operating agency, and classification of ozone monitoring sites considered in this study.DOI: https://doi.org/10.1525/elementa.300.

Table 2 :
Overview of the median and 95 th percentile summer O 3 , and the daily median and 95 th percentile summer O 3 amplitude of every site with data available for 2011-2015.DOI: https://doi.org/10.1525/elementa.300.t2

Table 3 :
Linear regression results for the trend analyses of summer O 3 for sites with at least 10 years and 80% of data for a given summer, categorized by the 5 th , 50 th , and 95 th percentile levels.Statistically significant trends (p-value ≤ 0.05) are highlighted in bold font.Additionally, the linear regression results for the trend analyses of the 4 th highest MDA8 and the Design Value are listed.DOI: https://doi.org/10.1525/elementa.300.t3

Table 4 :
Linear regression results of the trend analyses of the daily summer O 3 amplitudes derived from the yearly summer data, broken up the 5 th , 50 th , and 95 th percentile levels.Statistically significant trends (p-value ≤ 0.05) are highlighted in bold font.DOI: https://doi.org/10.1525/elementa.300.t4