Methane and volatile organic compounds and their influence on air quality in Boulder, Colorado

The Northern Colorado Front Range (NCFR), which encompasses the area east of the Rocky Mountains from the wider city of Denver to the Wyoming border to the north, has faced challenges over the past 30 years in achieving air quality standards for several primary and secondary air pollutants. Flocke et al. (2019) provide an overview of historical Denver air quality issues, starting from the 1960s. While in the 1980s and 1990s, carbon monoxide and particulates were the main recognized pollutants (Groblicki et al., 1981; Wolff et al., 1981; Sloane et al., 1991; Neff, 1997), during the most recent decades, secondary pollutants, that is, ozone, photochemically formed fine particulate matter, and wildfire smoke, have become a primary concern. Surface ozone at NCFR monitoring sites has been regularly exceeding the National Ambient Air Quality Standard (NAAQS) for more than 2 decades (Evans and Helmig, 2017; Bien and Helmig, 2018). In 2007, the NCFR was officially designated a nonattainment area for the 1997 standard of 80 parts per billion (ppb) maximum daily 8-h average O₃ mixing ratio (MDA8 O₃). The NAAQS standard was lowered to 75 and 70 ppb MDA8 O₃ in 2008 and 2015, respectively, and after a series of attainment deadlines failed to be reached, the NCFR was designated a severe nonattainment area for the 2008 NAAQS standard in 2022. Furthermore, there have been numerous exceedances of the 24-h 35 μg m⁻³ NAAQS for fine particulate matter, PM_2.5.

The air pollution problems in the NCFR have been associated with a combination of unique topographical and meteorological conditions. Summertime ozone events typically occur under weak synoptic forcing that allows thermally driven mountain-valley circulations to form. Mesoscale circulations such as the infamous “Denver Cyclone” can also lead to stagnation and ozone buildup (Flocke et al., 2019). The relatively high elevation and dry climate of the region foster nighttime radiative cooling that results in a shallow stable nighttime surface layer trapping pollution near the surface (Wilczak and Glendening, 1988; Wilczak and Christian, 1990). Terrain-driven flow regimes promote daytime upslope and nighttime downslope flow patterns that can result in polluted air from densely populated and industry-intensive areas, as well as from oil and natural gas (O&NG) production regions, being transported westward and trapped along the eastern slopes of the Rocky Mountains (Johnson and Toth, 1982). These upslope flows, in essence, capture ozone and its precursors and loft them above the atmospheric boundary layer (ABL) (Brodin et al., 2010; Pfister et al., 2017a). The evening and nighttime downslope flow can recirculate the previous day's ozone and precursors back into the NCFR and Platte River Valley, allowing a further buildup of ozone on subsequent days as long as stagnant synoptic conditions continue (Toth and Johnson, 1985; Reddy and Pfister, 2016). Additionally, the westerlies experienced at summit elevations push back some previous-day polluted air into the free troposphere over the NCFR, which can be entrained into the ABL the following afternoon (Kaser et al., 2017). Frontal passages flush out the air, ending an episode (Flocke et al., 2019). The combination of volatile organic compounds (VOCs) and nitrogen oxides (NO_x) emissions and sunny conditions promotes the photochemical processing of primary pollutants, resulting in high production rates of ozone and secondary particulate matter, which are the two pollutants that have received the most attention because of their NAAQS exceedances (Fehsenfeld et al., 1983; Dingle et al., 2016; Flocke et al., 2019).

During the past decade, there has been an expansion of O&NG well drilling and extraction, with most of this development located in Weld County in the northeastern sector of the NCFR (Figure 1). Weld County, for instance, increased its oil production by 206% from 2014 to 2019, followed by a slight decrease in subsequent years likely due to the COVID-19 pandemic (Colorado ECMC, 2024). Emissions from this industry and their impact on the regional air quality have received notable attention (Flocke et al., 2019; Helmig, 2020). O&NG emissions contribute the largest fraction of VOCs (Gilman et al., 2013; Abeleira et al., 2017), and the atmospheric oxidation of these emissions contributes significantly to the regional production of ozone (McDuffie et al., 2016; Lindaas et al., 2019). Cheadle et al. (2017) estimated a 20-ppb contribution of O&NG emissions to the regional ozone production in and downwind of Weld County on high ozone days. Likewise, modeling results building on data from the Front Range Air Pollution and Photochemistry Éxperiment (FRAPPÉ) field campaign in 2014 reveal an approximately 6–8 ppb contribution to NCFR MDA8 O₃ from O&NG; this contribution can be as high as 40 ppb on select high ozone days (Cheadle et al., 2017; Pfister et al., 2017b; Oltmans et al., 2021). Due to the prevailing east-to-west air flow during daytime, communities located downwind of the O&NG extraction area can be affected as much or even more by the secondary pollutants than the locations where the primary emissions occur (Evans and Helmig, 2017; Bien and Helmig, 2018).

Concerns over regional air quality motivated the implementation of a monitoring program in Boulder County, which lies between Weld County and the eastern slopes of the Rocky Mountains, approximately 40 km northwest of the Denver Metropolitan area. The Colorado Department of Public Health & Environment (CDPHE) relocated the Boulder ozone monitoring site, which previously had been at South Boulder Creek, to the Boulder Reservoir (BRZ) to the north of the city in 2016. Monitoring of methane, VOCs, and NO_x were added by the Boulder County Public Health department in spring 2017 in response to growing concerns of the oil and gas drilling impacting Boulder County’s regional air quality. First results reported on the response of ozone and NO_x to the 2017 solar eclipse (Helmig et al., 2018). Asher et al. (2021) showed that VOC mole fractions were highest nearest the surface in their sampling of vertical gradients near BRZ. Pollack et al. (2021a) provide an in-depth analysis of the causes of elevated ozone: there were 41 days at BRZ in 2017, 2018, and 2019, when the MDA8 O₃ exceeded 70 ppb. Of these 41 high O₃ days, 40 occurred between May and September, 2 were influenced by stratospheric intrusion events, and 13 days were impacted by wildfire smoke. Light alkanes (C₂–C₅), in particular, are dominated by O&NG sources (Pollack et al., 2021a).

Since 2004, there has been a progression of Colorado regulations tailored toward reducing O&NG emissions in light of their increasingly recognized influences on the NCFR air quality (Milford, 2014; CDPHE, 2017; Wells, 2017; CDPHE, 2021). Projections of changes in emissions have mostly been based on bottom-up emission estimates (using inventory data that build on emission factors) prepared by CDPHE. Several field campaigns by university and federal laboratories have derived top-down emission estimates (using atmospheric observations), which have pointed toward the inventory underestimating actual emissions (LaFranchi et al., 2013; Helmig, 2020; Lu et al., 2023).

This study further investigates the sources of methane and VOCs at BRZ for evaluating the State of Colorado’s O&NG VOC emissions inventory and the projected emissions change building on the first 4 years of monitoring data. Two studies (Lyu et al., 2021; Oltmans et al., 2021) previously reported trend analyses of O&NG emissions building on earlier (2013–2016; 2008–2016, respectively), surface and tower datasets. Since the data collections that these studies rely on, there have been at least 6 regulatory measures affecting O&NG emissions (Supplement Table S1). Our research builds on higher time-resolution data compared to these prior publications (10–20 times more measurements per year compared to Oltmans et al., 2021, and ≈100 times more compared to Lyu et al., 2021) from a surface monitoring station, allowing for more targeted data filtering and more precise analyses. Furthermore, our study investigates a subsequent period (2017–2021), adding 5 years of analysis during a period with substantial growth in O&NG production.

The Boulder Reservoir monitoring site (BRZ; 40.0704°N, 105.2198°W, 1,604 m asl) is located on the southeast corner of BRZ (Figure 1), ≈7 km northeast of the center of the City of Boulder. The monitoring station is ≈25 m from the reservoir shore and has ample open fetch on all sides with only a few isolated trees and buildings within a 500 m radius. A sailing and motorboat storage area is located to the west. Approximately 500 m further to the west are the reservoir recreational beach and parking areas. The park hosts several hundred visitors during summer days. Recreational activities peak in the summer and are relatively minor during the off-season months. A major traffic artery into the City of Boulder, Highway 119, runs diagonally from southwest to northeast ≈400 m east of the site. The monitoring at BRZ is ongoing at the time of this publication with real-time data provided to the public (Boulder AIR, 2024). The lower panel of Figure 1 provides detail on the park facilities that are mostly located to the west of the monitoring station. More site descriptions have been provided previously (Helmig et al., 2018; Pollack et al., 2021a). A map illustrating mostly abandoned well sites within a ≈5 km radius of BRZ is provided in Supplement Figure S1.

During this study, 141,944 methane samples were collected with an SRI Model 310 gas chromatograph every 15 min, with an accuracy uncertainty of ≤5 ppb and precision error of ≤8 ppb for single data points. VOCs were measured at 1–2 h temporal resolution, resulting in 19,335 total samples (see Supplement Table S2 for calibration standards). Quantified compounds include ethane, ethene, acetylene, propane, cyclopropane, propene, i-butane, n-butane, i-pentane, n-pentane, isoprene, n-hexane, benzene, n-heptane, toluene, n-octane, ethyl-benzene, o-xylene, m-xylene, and p-xylene. Detection limits varied by species and were in the range of 9–33 pmol mol⁻¹ (ppt) (Supplement Table S3). Nitrogen oxides are measured at 5-min resolution, and winds (speed and direction) are measured at 1-min resolution (Supplement Figure S2). The monitoring and data collection methods have been described previously (Helmig et al., 2018; Pollack et al., 2021b). For convenience, we repeat the most important information in Supplement Text T1.

Methane and VOC data from the Niwot Ridge (NWR) reference site were collected at the T-Van sampling location (40.05°N; 105.58°W; 3,523 m asl) (NOAA GML, 2024). NWR is ≈33 km west of BRZ in the Rocky Mountains and 1,919 m higher in elevation. Figure 1 shows both sites in relation to O&NG wells and county borders. During most times, NWR remains well above the plain’s ABL height (Fahey et al., 1986; Brodin et al., 2010).

Orthogonal distance regression (ODR) analyses were calculated using Python stats and ODR Scipy packages. Specific formulas and residual variance (percentage of variance not explained by the correlation of variables, roughly analagous to 1 − r², where r is the Pearson correlation coefficient) calculations can be found in Boggs et al. (1992). Potential Source Contribution Function (PSCF) analyses determine the probability that an emissions source contributing to elevated mole fractions of atmospheric constituents is located at a grid cell with latitude i and longitude j (Pekney et al., 2006). PSCF modeling was conducted utilizing NOAA’s 3 km High-Resolution Rapid Refresh (HRRR) meteorology output and the Hybrid Single-Particle Lagrangian Integrated Trajectory model (HYSPLIT) to compute back trajectories. Trajectories were modeled at a 5-m starting elevation and calculated in hourly increments over 12 h backward in time. The modeling relies on latitude, longitude, height, and pressure for each hour backward in time provided by the HRRR output. BRZ VOC data were aligned with the closest back trajectory in time for the PSCF analysis. More detail about the modeling application has been provided in Oltmans et al. (2021).

For wind rose plots, wind velocities were paired with VOC measurements by vector averaging the 1-min wind observations during the 10-min sampling period. For methane measurements, the 1-min wind observations closest to each sampling moment (every 10 min) was used.

Results for methane, ethane, and acetylene from the in-situ monitoring at BRZ are compared with the daily sample collection at NWR in Figure 2 for periods when overlapping data were available. Please note that the NWR data are all daytime samples, whereas the BRZ data are sampled 24 h per day every 1–2 h (see Supplement Figure S3 for propane, i-butane, n-butane, and n-pentane comparisons).

Averaged BRZ VOC diurnal cycles for ethane, a primary O&NG emission, and benzene, which has its strongest contribution from transportation sources (Pollack et al., 2021a), are presented in Figure 3, using data from the start of the monitoring (April 2017) through April 2019 to calculate pre-pandemic annual cycles. Ethane mole fractions increase steadily during evening hours, peaking in the early morning, and declining throughout the day. Ethane emissions are likely mostly fugitive emissions that occur day and night, therefore the diurnal ethane behavior is most likely a reflection of atmospheric mixing dynamics. The benzene cycles have two peaks, one in the early morning (2–3 h after the ethane peak) and a second one in the evening. These diurnal mole fraction cycles do not closely follow the diurnal cycle of automobile traffic (Figure 3C) that shows a sharp morning and a wider noon to afternoon peak. This divergence is likely because the atmospheric concentrations are dependent on the changing emissions and the confounding dilution from atmospheric mixing.

Other prior work has shown that diurnal cycles of primary and secondary pollutants that arise from NCFR sources become progressively weaker with altitude (Brodin et al., 2010; Rossabi et al., 2021). At NWR’s elevation, diurnal cycles are expected to be relatively weak. Therefore, there should only be a negligible bias from the difference in the daytime (NWR) versus 24-h (BRZ) sampling. NWR measurements may slightly overestimate species concentrations of free-tropospheric background air on summer afternoons with strong easterly upslope flow (Brodin et al., 2010; Oltmans et al., 2021).

There is greater variability and higher concentrations in observations at BRZ compared to NWR. While the low values in the BRZ measurements were in close agreement with the NWR background data, there was an abundance of occurrences of elevated concentration spikes at BRZ. This behavior was seen for all measured gases, including methane. For methane, enhancements reached up to ≈150% of the ≈1,950 ppb background. The atmospheric lifetime of methane is ≈8 years (Sonnemann and Grygalashvyly, 2014), while C₂–C₆ alkane lifetimes range between 3 and 180 days (Helmig et al., 1996). Because of their shorter lifetime compared to methane, VOCs background mole fractions are about 3 orders of magnitude smaller. With the lower background values, the spikes reached much higher mole fractions relative to their background. For instance, spikes for ethane reached values of up to 100 times the background seen at NWR. The much higher variability and abundance of high mole fraction spikes resulted in an overall higher spread in the statistical distribution, shown in the side-by-side comparison of daytime results in Figure 4. While the lower percentile values were relatively similar, the differences between the 2 datasets become larger for the higher percentile values. The high mole fractions at BRZ were largely caused by frequent but relatively short spikes, lasting between 10 min and a few hours. These elevated mole fraction events occurred throughout the year; however, their abundance and observed maximum mole fractions were higher during the winter. These plumes carried both high methane and high concentration VOCs. Between July 2017 and June 2018, there were 2,187 out of 4,243 samples (51%) in the BRZ record that exceeded the 99th percentile value of ethane at NWR (5.55 ppb). During the colder months of the year, due to weaker solar irradiance and convective mixing, nocturnal boundary layers are stronger and longer lasting, and daytime boundary layer depths are lower. These are most likely the explanations for the more abundant and higher concentration winter spikes, as emissions released at the surface do not mix as vigorously and maintain their elevated VOC signatures for longer distances during transport. In contrast, NWR remains well above the boundary layer on the plains during most times and is subjected primarily to lower free tropospheric air transported into the region (Fahey et al., 1986; Brodin et al., 2010). Seasonal cycles at BRZ are most pronounced for ethane and propane and weaker for heavier VOCs. Seasonal summary plots are shown in time series format in Figure 5, and as monthly box-whisker summer graphs in Supplement Figure S4.

Benzene and toluene deviate from the n-alkanes by showing a secondary summer maximum that is consistent throughout the years of measurements (see more discussion below). A summary of the numerical results is provided in Supplement Table S4. Methane shows its highest seasonal mole fraction during the winter months. The methane seasonal amplitude is on the order of ≈100 ppb. This is well above the seasonal amplitude of ≈30 ppb seen in the background atmosphere. This stark difference underscores the large influence of the reduced mixing in the winter, which exerts a stronger trapping of nearby emissions at the surface and drives the ≈3 times higher concentration enhancement seen in the winter data (BRZ vs. NWR).

BRZ diurnal cycles and weekday/weekend comparisons of methane and selected VOCs are presented in Figure 6. Mean and median values were calculated and binned for each hour of the day in Mountain Standard Time (MST). We chose to present the median results as they are less sensitive to extreme high values and are therefore deemed to provide a more representative result for the general behavior in the data (for a comparison of the median with mean results see Supplement Figure S6). Monday through Friday were considered as weekdays, and Sunday and federal holidays were binned together as holiday-Sundays. Saturdays were excluded. There are approximately 5 times more weekdays feeding into the analysis, which explains the smoother diurnal cycle curves for weekdays.

To our knowledge, there are no plausible reasons for significant diurnal cycles in methane emissions. Therefore, the methane diurnal cycles (Figure 6A) can be used as an indicator of the influence of the diurnally changing boundary layer dynamics on atmospheric mole fractions. The diurnal changes in methane from the suppressed mixing in the nocturnal boundary layer, and increased mixing during the day from convective turbulence, result in a multiyear averaged ≈50 ppb diurnal methane cycle amplitude in surface air. The slight difference seen between weekdays and holiday-Sundays is within the noise of the determination and does not point toward a significant difference. A similar diurnal behavior is seen for ethane and propane. Their data are more noisy, likely due to the less frequent VOC sampling. Using the methane data as reference, the most likely explanation for the higher levels at night and declining values during daytime is again the forcing from mixing/dilution. Both ethane and propane results point toward slightly higher levels on weekdays, particularly during daytime hours, although these differences are not statistically significant. Nonetheless, this possibly points toward higher rates of natural gas releases from venting and other well site operations during weekday work hours. Diurnal cycles for benzene, toluene, and acetylene are distinctly different. These show maxima during early morning and evening hours when vehicular traffic is highest. The second maximum in the evening is around 21:00 MST. This is 3–5 h after the rush hour peak. It may indicate the accumulation of these evening emissions during the time of rapidly weakening boundary layer mixing. Weekday values are mostly higher than weekend/holiday values throughout the day for benzene and toluene, although the difference is less striking for benzene, and both cycles are within 2 standard errors of each other. The behavior seen in the data may be attributed to substantial contributions from traffic and other activities toward total emissions. A Wednesday–Sunday data comparison of nitrogen oxides (NO_x) at BRZ pointed to a similar magnitude of the contribution of traffic on these primary mobile source emissions, as NO_x levels were on average 20%–40% lower on Sundays (Pollack et al., 2021b).

During the summer, the VOC samples showed a prominent presence of the biogenic hydrocarbon isoprene. Isoprene values regularly exceeded 1 ppb, and at its maximum, isoprene occasionally had the highest mole fraction of the monitored VOCs (Figure 7). There is a sudden onset of the isoprene occurrence in spring and a sudden cessation in fall. Outside of this summer growing season, isoprene remained well below 0.1 ppb. During the 4 years of observations, the spring onset of observed isoprene (first value >250 ppt) varied from May 19 to June 10 (2017: June 3; 2018: June 10; 2019: June 3; 2020: May 19), and the date of the last >250 ppt isoprene occurrence varied from September 21 to October 17 (2017: September 21; 2018: September 30; 2019: September 26; 2020: October 17). Taking the isoprene occurrences as an approximate indication of the foliage/growing seasons, this is a remarkable variation in the length of the season, from a minimum of 103 to a maximum of 151 days. During the peak of the summer, monthly mean isoprene can occasionally constitute up to half of the OH reactivity within the spectrum of quantified VOCs (Supplement Figure S7). It is likely that emissions from several large cottonwood trees near the shorelines of the reservoir contribute to the observed isoprene. Cottonwood trees have been labeled as a moderate isoprene emitter, with a standardized isoprene emission rate of 9,050 ng g⁻¹ h⁻¹ (Villanueva-Fierro et al., 2004). The larger area surrounding BRZ is a combination of agricultural land, the Western Great Plains Foothill, and Piedmont Grassland, which are unlikely isoprene sources. A comparison of available concurrent isoprene monitoring data from BRZ with monitoring sites in Longmont and Broomfield for the latter part of the summer of 2020 shows that BRZ daytime maximum values are on average 2–5 times higher than at these comparison sites (Supplement Figure S8). It is therefore unlikely that these isoprene data constitute representative observations for the surrounding area at large. According to this data comparison, atmospheric isoprene shows a high sensitivity to the vegetation species type and distribution at particular monitoring sites. More research on the regional isoprene occurrence is needed and care should be taken extrapolating individual site data for regional atmospheric chemistry observations.

Methane and VOC abundances show a strong dependency on wind speed (Figure 8). All species show a similar pattern with more variability and highest mole fractions seen at lower wind speed, and less variability and lower mole fractions seen with increasing wind speeds. This behavior indicates more dilution of the NCFR air with lower concentration background air as wind speed increases, that is, stronger wind speeds allow pollution plumes to be more quickly dispersed by advection (Pfister et al., 2017a). At the highest observed winds, BRZ mole fractions converge to the background values seen in the NWR data. An interesting and unexpected feature in the BRZ methane graph is that the lowest percentile values at the lowest wind speeds are lower than the methane values at the highest wind speeds (Figure 8A). This is likely a consequence of the annual wind patterns in the region: The strongest winds are mostly westerlies that occur in the winter months, when the seasonal background of methane is at its highest. This is confirmed by NWR data that show methane values that are ≈30 ppb higher during the winter in comparison to the summer minima (Supplement Figure S9).

The dependence of observed mole fractions on winds was further investigated using wind rose and bivariate polar plot analyses. Concentration-wind rose analyses show a high sensitivity to the applied wind speed filter. While calmer winds are relatively evenly distributed, stronger winds mostly originate from the west and north sectors. Data filtered for stronger winds show an increasingly narrower definition of the primary upwind source sectors, with winds from the northeast sector transporting the highest levels of methane and VOCs to BRZ (see Supplement Figure S10 for the dependency of methane and ethane wind roses on the wind speed filter). A clear distinction between different contributing VOC sources can be derived from the analysis of the isomeric pentane ratio (i-pentane/n-pentane). Several prior studies have applied this tool for elucidating VOC source categories. While O&NG VOCs are characterized by a relatively narrow range of ≈0.9–1.1, values of >1.5 are more typical for urban and photochemically processed air (Gilman et al., 2013; Thompson et al., 2014). The distribution of pentane isomeric ratio values in air transported from the northeast sector is distinctly lower and centered in the range of the O&NG signature window compared to values seen in air transported from other sectors (Figure 9), similar to findings at other sites in this region (Swarthout et al., 2013; Halliday et al., 2016; Oltmans et al., 2021). The differences become starker when only data during higher wind speed conditions are considered. These results illustrate that emissions source characteristics are best defined in subsets of the data that have been filtered for higher wind conditions.

Bivariate polar plots (Figures 10, 11, S11–S13) similarly parse the wind and species concentration datasets into wind directional sectors, but instead, color the radial dimension with concentrations corresponding to given wind speeds. The default binning, interpolation, and smoothing settings are used from the OpenAir R package, with direction bins of 10 degree sectors and 30 total speed bins (Carslaw and Ropkins, 2012). The minimum bin (min-bin) parameter specifies the minimum number of values in each bin needed to be considered in the interpolation and shown in the figure. Gray areas represent wind velocities that occurred but with insufficient frequency to meet the min-bin parameter. An evaluation of the wind measurements and the procedure for filtering of data as a function of wind conditions is provided in Supplement Text T2. Results for the bivariate polar plot analyses for methane and five selected VOCs are presented in Figure 10 (see Supplement Figure S13 for a seasonal breakup of these data). These graphical outputs point to a similar conclusion as the wind rose graphs: Elevated mole fractions of methane and light petroleum hydrocarbons originate predominantly from the north to east of BRZ. Elevated pollutant levels were observed even at high wind speeds from this sector despite the stronger dilution. Benzene, toluene, and acetylene were less clearly defined. Relatively high acetylene occurred with northerly winds. This suggests a source from the Reservoir, such as from motor boating, or north of the Reservoir. Similar results for i-pentane and n-pentane suggest that these may similarly have contributions from recreational activities at the BRZ (Supplement Figure S13c).

Figure 11 presents bivariate wind analysis results for five selected ratios of VOC pairs. The most revealing results are for propane/ethane and i/n-pentane where there is a clear signature of higher propane/ethane and lower i/n-pentane ratio values from the O&NG sector. Air transported to the site from the west has higher values for benzene, i-pentane, and toluene, particularly during the summer (Supplement Figure S13). Our best interpretation for these signatures is that these potentially stem from recreational equipment, possibly outboard engines used on boats on the Reservoir.

Observations of the quantified individual VOCs were multiplied by their OH rate constant and added together for a bivariate analysis of the OH reactivity from these VOCs. This analysis was conducted for daytime data in summer to investigate source regions of this partial OH reactivity during maximum ozone formation. Isoprene was excluded from the analysis to focus on anthropogenic contributions. Results shown in Figure 11, panel F, identify two primary source sectors: One reactivity hot spot is the sector extending from the north to the northwest. This sector overlaps with the footprint of the Reservoir and the recreational areas, suggesting that emissions from activities within the park, as hypothesized in the preceding section, constitute a notable contribution to the total VOC reactivity. This interpretation is supported by the bivariate polar wind analysis results for toluene, benzene, i-pentane, and n-butane that all show signs for elevated emissions in the direction of the Reservoir Park footprint during the summer (Supplement Figure S13). Ethene, acetylene, and propene show enhanced mole fractions in air flow from the north during the summer months, which suggests that these alkenes may result from motorboat engine exhaust on the Reservoir (Supplement Figure S13). The other sector that exhibits a signal with increased OH reactivity lies to the northeast in the direction of the center of the O&NG area. It is remarkable that air masses originating from the southeast and the southwest, which is the direction and closest distance to Highway 119 and the City of Boulder, have relatively low OH reactivity. Caution should be exercised as these OH reactivity results exclusively rely on the reactivity from the subset of VOCs that were measured here. The total OH reactivity will have contributions from additional VOCs and other reactive trace gas species (Coggon et al., 2021).

PSCF analysis offers an independent approach for examining wind and transport dependency. PSCF does not rely on the local wind measurements, but instead builds on the paring of observed mole fractions with HYSPLIT back trajectories; in our case, 12-h back trajectories were used. The PSCF output shows the probability that a trajectory that was paired with a measurement result was above the overall median mole fraction when it passed over a given grid cell surrounding the site, assuming a Lagrangian framework, that is, that an emission follows an air parcel without diffusion or chemical loss. The probability scale reflects a normalized count of the trajectory passes over each grid cell. Source region results shown in Figure 12 are similar for methane, ethane, and propane. The likelihood for air parcels to pick up enhanced levels of these gases were highest when trajectories passed over the area from approximately the northeast to the southeast of BRZ.

The source footprint for total VOC carbon (Figure 12D) is similar to the one for the alkanes because alkanes constitute the bulk of the total VOC carbon. Air originating from the direction of Denver to the southeast has a comparatively minor contribution despite Denver being located at the same approximate distance as the O&NG area in Weld County. Results for the i/n-pentane ratio and nitrogen oxides were included in Figure 12 primarily for comparison and as a method sensitivity verification. The likelihood to observe above median values of the pentane ratio, reflecting urban emissions, is highest in the transport from the south, whereas lower values, reflecting O&NG signatures, are observed in transport from the northeast. Results for NO_x are similar: the strongest source region is also to the south and southeast, overlapping with the transportation sector between Boulder and Denver. The breakdown of these PSCF analyses in seasonal segments does not show notable differences in source footprints (Supplement Figure S14). The general conclusions from the bivariate wind and the PSCF analyses are rather consistent as both methods point toward air with elevated VOCs and an O&NG VOC signature transported from an area extending 40–50 km mostly to the north to southeast. These primary source footprints show strong overlap with the O&NG well distribution, mostly in Weld County (Figure 1). Mobile source and urban emissions, as identified by NO_x and a higher i/n-pentane ratio, are more prominent in air transported to BRZ from the south or southeast directions. These footprint analyses and conclusions agree with similar analyses conducted on an earlier (2008–2016) dataset from the Boulder Atmospheric Observatory (BAO) (Oltmans et al., 2021), located approximately 25 km to the east of BRZ within the southwest corner of Weld County.

The highest observed ethane mole fraction over the 4-year reporting window was 117.0 ppb on December 22, 2017. On only one other occasion did ethane exceed 100 ppb (January 5, 2018); 5–10 plume events between 60 and 80 ppb were observed during most winters. The nature of the December 22, 2017, plume is further investigated as a case study to illustrate the chemical composition and possible origin. Three-day time series graphs spanning the event are plotted in Figure 13. This plume was observed at BRZ for approximately 4 h. Three VOC samples were recorded in which VOCs were elevated 30–50 times above the prior runs. For instance, the 12:58 MST sample (recorded time is the middle of the 10-min sample collection period) had an ethane mole fraction of 3.39 ppb. The results for the 3 following runs at 14:52, 16:47, and 18:41 MST were 106.0, 117.0, and 95.8 ppb, respectively. All these recordings were 15–18 times higher than the highest value (6.4 ppb) seen in the available data from NWR. Similar enhancements were seen for other VOCs, that is, peak values were 108.9 ppb for propane, 49.6 ppb for n-butane, and 15 ppb for n-pentane. Methane increased from 1,953 ppb prior to the event to over ≈2,900 ppb, an enhancement of ≈1,000 ppb. The maximum benzene value was 0.83 ppb, which was the second highest benzene recorded in 2017. The mean isomeric pentane ratio in the 3 plume samples was 0.82, characteristic of a natural gas signature, and similar to the Denver-Julesburg Basin (DJB) raw O&NG signature value of 0.86 (Gilman et al., 2013). This plume had a characteristic “wet” (>15 weight % total VOCs) natural gas signature, with methane being the most abundant constituent, followed by ethane, propane, the butanes, and so forth, with relative abundances of these petroleum hydrocarbons decreasing with increasing carbon number. Acetylene, often used as a combustion tracer, reached 1.94 ppb, which is a more moderate enhancement in comparison to the alkanes (i.e., equivalent to a 92.8 percentile value for the full dataset).

December 22 was a cold winter day with a high temperature of −7°C and a low temperature of −15°C. A strong temperature inversion was present at 400–500 m above ground as seen on morning and afternoon Denver radiosonde data, associated with a regional high-pressure system and stable synoptic meteorological conditions. These conditions limit the vertical transport and mixing, enhancing pollutant concentrations near the ground. Winds (Figure 13E) were from the north to northeast prior to and during the onset of the event, but then shifted to the west at 16:00 MST. Wind speed averaged 2–3 m s⁻¹ (equivalent to an approximate 7–11 km surface transport distance per hour), but then decreased to approximately half those values for the remainder of the day. Winds at Longmont were also northeasterly and between 3 and 5 m s⁻¹ (Supplement Figure S15). The back trajectories for the hours surrounding the event (15:00–17:00 MST) show transport from the north to northeast (Supplement Figure S16), in agreement with the wind measurements at BRZ (Figure 13) and the surrounding airports (Supplement Figure S15). These flow conditions transported air from the north and northeast in Weld County for ≈6 h toward BRZ and were preceded and followed by more rapid air transport from the west (Supplement Figure S16). It is therefore likely that the plume originated from a further distant well site beyond the City of Longmont. Further beyond the City of Longmont, the back trajectories passed near the town of Windsor, 54 km to the northeast of BRZ, 6–9 h prior to arriving at BRZ. On that same day, a well site caught on fire within that area after a large explosion at 20:46 MST (Denver Post, 2017). This raises the question if these elevated VOCs possibly originated from this well pad and that the VOC enhancements were from leakage or venting that might have been occurring prior to the time of the explosion. With the ≈4–5 h duration of the plume as observed in the BRZ data, under the December 22 wind conditions, the plume must have had a geographical extent along the trajectory of several tens of kilometers and passed over an area populated by well over 100,000 people, including the City of Longmont.

The graphical comparison in Figure 14 shows the BRZ ethane data in 2017 relative to canister sampling results from 28 other American cities from 1999 to 2004 (Baker et al., 2008). Note that Boulder is by far the smallest city in this comparison (population of 106,000 in 2019, 450 people per square mile). It is also important to note that the BRZ monitoring is outside the Boulder city footprint, whereas the comparison data stem from sampling within those cities. When filtering BRZ data for samples taken during August and the first week of September and between 10 and 19 Mountain Daylight Time (MDT) to compare with Baker et al. (2008), the mean ethane results from Boulder, despite its smaller size, are the fourth highest in this comparison, trailing the Pittsburg, Charleston, and Los Angeles data. Another aspect to consider is the 13–19 years’ time difference in these data. VOCs have steadily decreased in U.S. urban areas throughout the past decades (Warneke et al., 2012; Rossabi and Helmig, 2018). For instance, an average −6.2% yr⁻¹ decline in ambient ethane was reported for greater Los Angeles from 2001 to 2008 (Warneke et al., 2012). If that rate of decline is applied to the Baker et al. (2008) data and the time lag between the datasets is used to correct for the decrease in urban concentrations, then the BRZ ethane results constitute the highest values in this 28-city comparison (Supplement Figure S17). Comparison results for other O&NG tracers were similar, with a tendency of the BRZ VOC results steadily dropping in their ranking with increasing carbon number (not shown).

These analyses underscore that atmospheric VOCs that result from O&NG well operations are remarkably high not just in close proximity to the extraction areas as has been shown in a number of shorter campaigns (Pétron et al., 2012; Gilman et al., 2013; Thompson et al., 2014; Oltmans et al., 2021), but also in downwind populated regions throughout the year.

Atmospheric methane can arise from a multitude of emission sources, that is, livestock, landfills, wetlands, and O&NG, with many of these present and abundant in the NCFR. Previous work concluded that O&NG emissions dominate, contributing ≈75% of total methane emissions in Weld County (Pétron et al., 2014; Peischl et al., 2015). Ethane is a more selective tracer for O&NG emissions because other emission categories make relatively minor contributions to the ethane flux. Building off this relationship, atmospheric ethane observations, combined with methane–ethane O&NG emission ratios, have been used to estimate the O&NG-associated methane flux and its trend (Aydin et al., 2011; Simpson et al., 2012; Helmig et al., 2016; Tzompa-Sosa et al., 2019).

Quantified hydrocarbons show a strong correlation with methane, indicating that they share common sources, footprints, and/or common forcings that influence their ambient air concentrations (Supplement Figure S18). The correlation is the strongest for natural gas hydrocarbons, that is, ethane, with a residual variance of 19%, followed by propane and i-butane. Correlation results become weaker with increasing molecular weight, with the residual variance increasing to 44% for toluene. A similar residual variance value was determined for benzene and acetylene. Positive matrix factorization modeling previously associated the alkanes overwhelmingly to O&NG emission sources, while acetylene and toluene were found to have only a minor association with O&NG emission. The total C₂–C₅ alkanes association with O&NG sources was determined to be 79% ± 1% and 84% ± 20% during the winter and the summer, respectively. In contrast, for acetylene, the main contributing emission factor was traffic (45% ± 6% in winter and 87% ± 32% in summer) (Pollack et al., 2021a).

Figure 15 investigates the dependence of VOCs on the O&NG tracers methane, ethane, and on the i/n-pentane ratio as O&NG indicator. Most of the high mole fraction VOC observations, in particular high ethane mole fractions, are observed at an isomeric pentane ratio of 0.9–1.1. As the pentane ratio increases, reflecting dilution with background or urban air, ethane values decrease. This abundance of high ethane measurements in the pentane isomeric O&NG signature window reaffirms how ethane is dominated by O&NG sources. The predominance of high mole fraction values at an isomeric pentane ratio of 0.9–1.1 is highest for the petroleum hydrocarbons, and highest overall for ethane and propane in comparison to other analyzed VOCs. Using ethane and the pentane ratio as an O&NG tracer, the analysis in Figure 15, panels C–E investigates the dependency of the total identified VOC carbon flux on the isomeric pentane ratio and on methane and ethane. VOCs that were considered for the total atmospheric VOC carbon were ethane, propane, i-butane, n-butane, i-pentane, n-pentane, n-hexane, ethene, propene, benzene, toluene, o-xylene, m- and p-xylene, ethyl-benzene, and acetylene. Total VOC carbon and ethane correlate tightly; there is only a 4% residual variance, compared to 22% for methane. These results add additional evidence suggesting that the sources contributing to the total volatile carbon of the included VOCs are overwhelmingly associated with O&NG emissions.

Correlations between pairs of VOC/VOC (Supplement Figure S19) show narrow correlations between isomeric pairs and VOCs that are closest in carbon number. For instance, the i-butane with n-butane relationship shows a mere 1% residual variability; the result for the pentane isomers is 4%. The ethane–propane correlation shows a 6% residual variability. Correlations are in general narrower among alkanes than between alkanes and aromatic compounds. Residual variability for benzene with alkanes is on the order of 30%–40%. The result for benzene–toluene is 28%, which is relatively weak compared to other results, indicating that both emissions have a significant fraction from non-common emission sources. Acetylene has the highest correlation with benzene, and benzene the second highest with acetylene, confirming their relatively high portion of shared emissions sources. Toluene was found to be mostly attributed to the regional anthropogenic background, while benzene has a relatively even mix of sources from short-lived O&NG, traffic, and regional anthropogenic background sources (Pollack et al., 2021a).

Assuming that most of the VOCs and methane result from a shared O&NG emissions source, and that emissions from a common source are responsible for most of the elevated values, allows using methane and VOC relationships to estimate the approximate composition of the natural gas source. In order to increase the selectivity toward O&NG-influenced samples, the data were filtered for wind direction from the primary oil and gas sector, that is, 0°–120°, and for ethane values over 5 ppb, and an i/n-pentane ratio range of 0.75–1.2. Results in Table 1 show that on average, the O&NG plumes contained 17.8% ethane and 18.1% propane relative to methane. The combined C₂–C₅ alkanes (ethane, propane, butanes, and pentanes) accounted for 59.7% of the methane flux. With this composition, the fractional contribution of these VOCs to the total O&NG hydrocarbon flux is 0.597/(1 + 0.597) ≈37%, with methane contributing the remainder of ≈63%. Inclusion of additional VOCs (>C₅, cyclic compounds) would further increase the VOCs and decrease the methane percentage. Natural gas with <85% methane is defined as “wet” gas. Consequently, the chemical signature in the BRZ data, with a 37% mass contribution of these listed non-methane VOCs alone, clearly points to the wet character of the natural gas source. This characterization agrees qualitatively with the results from airborne measurements where just the ethane/methane emission (mass) ratio for the DJB was estimated to be ≈30% (Peischl et al., 2018).

A top-down DJB ethane flux estimate of 61.3 ± 9.6 Gg yr⁻¹ was estimated from basin-wide aircraft profiling in April 2015 (Peischl et al., 2018). A more recent aircraft study from October 2021 resulted in a lower ethane flux of 27.2 ± 12.3 Gg yr⁻¹ (Fried and Dickerson, 2023). Using correlations between ethane and other VOCs from BRZ data filtered for O&NG emission plumes offers an opportunity to estimate other O&NG VOC fluxes based on these observed dependencies. Since ethane is the more selective O&NG tracer, we developed a flux estimate based on ethane rather than methane, as scaling to the methane flux would add uncertainties in the O&NG emissions allocation because of the stronger contributions from non-O&NG sources to the total methane flux. The principle of using ethane as an anthropogenic O&NG tracer and its relationship with methane and other VOCs for determining O&NG fluxes has been extensively applied in prior research (Katzenstein et al., 2003; Aydin et al., 2011; Simpson et al., 2012; McKain et al., 2015; Roscioli et al., 2015; Ramsden et al., 2022). The VOC/ethane slopes from the O&NG-filtered BRZ data were scaled to both the April 2015 and October 2021 ethane DJB flux estimates. The filtering for isolating O&NG-influenced samples used in Section 3.6 left 2,248 data points out of the overall 19,335 samples. ODR correlation analyses results are provided in Supplement Figure S20, and standard linear regression correlation analysis results are shown in Supplement Figure S21. Supplement Text T3 explains how the nonidentified VOCs were estimated. The slope values of the correlations are presented in Table 2, together with the VOC flux estimates resulting from multiplying the VOC/ethane mass ratio values with a given aircraft ethane flux estimate. 1σ uncertainties are calculated by standard error propagation with the slope uncertainties in the BRZ data and the uncertainties from the top-down ethane flux estimates that were reported in the referenced field campaigns. These calculations result in a 245.0 ± 16.9 Gg yr⁻¹ total VOC flux estimate for April 2015, and 108.5 ± 21.4 Gg yr⁻¹ estimate for October 2021. It is important to note that while we report our estimates as annual fluxes for a comparison to inventory and other literature, these flux estimates reflect the specific time period from each top-down aircraft study, and our estimates rely on the assumption that VOC/ethane ratios have remained constant during the 2015–2021 time window. As mentioned above, these calculations are based on the slope values of ODR analyses which minimizes both errors in the ethane and other VOC data. A sensitivity analysis that was performed for the April 2015 results using a linear regression analysis instead of ODR (Supplement Figure S21, Supplement Table S6) yielded results that were ≈12% lower, that is with a total VOC flux estimate of 229.9 ± 15.9 Gg yr⁻¹.

These flux estimates offer a rare opportunity to evaluate the inventory data that are considered in the recently completed State of Colorado Serious State Implementation Plan for the Denver Metro and North Front Range Ozone Nonattainment Area (Regional Air Quality Council [RAQC], 2020). This report provides inventory estimates for total oil and gas VOC emissions. The VOC inventory estimate excludes ethane as it is considered nonreactive for ozone photochemistry. The VOC inventory flux for the years 2011, 2017, and 2020 projects an overall declining trend of −6.5% yr⁻¹ (relative to 2011), with a slight decrease in this rate for 2023 (RAQC, 2020; 2023). A linear interpolation through these 3 inventory estimates yields April 2015 and October 2021 O&NG VOC (excluding ethane) flux values of 72 and 40 Gg yr⁻¹, respectively (Supplement Figure S22). Our best VOC (excluding ethane) flux estimate from the ODR analysis yields 183.6 ± 12.6 and 81.3 ± 16.1 Gg yr⁻¹ for these respective dates, which is 2.0–2.5 times higher than the respective inventory values. We remind the reader that each ethane flux estimate was obtained during only 2 research flights conducted on 2 days (March 28 and 29, 2015, then October 1 and 5, 2023). Carbon Mapper data suggest that methane emissions from O&NG in the DJB vary by a factor of 2 from season to season (Cusworth et al., 2022). Any flux estimate based on only a few flights may not represent annual data, which is why our high-time resolution and continuous measurements are a useful complement. This flux estimate also hinges on a constant VOC/ethane ratio in O&NG emissions across the DJB. There are several other datasets that allow evaluating the validity of this assumption, using the ethane/n-butane ratio as an indicator. In the BRZ data, this ratio is 3.03, that is, the reverse of the n-butane/ethane regression line slope in the O&NG filtered data (Supplement Figure S20). This value agrees well with results from several other previous studies, that is, with values of 3.01 (ratio of medians; Gilman et al., 2013) and 3.86 (ratio of medians; Swarthout et al., 2013) that were determined in 2011 winter data at the BAO. The range of this ratio was 2–4 in boundary layer aircraft sampling during multiple flights across a wide geographical range of the DJB during the FRAPPÉ (Eric Apel, National Center for Atmospheric Research, personal communication, 2021). These ethane/n-butane data define the relative uniform composition of emissions from O&NG sources across the DJB and the representativeness of the BRZ data. Additionally, an analysis of VOC/ethane ratios for O&NG-filtered data in the BRZ dataset reveals no statistically significant slope over 2017–2021 for any individual VOC (not shown) or total VOCs (Supplement Figure S23), increasing confidence that these ratios are reasonably representative for 2 years earlier in 2015. Thus, the dominant source of error in our VOC flux estimates is likely the uncertainty in the top-down ethane aircraft flux estimates, followed by small temporal fluctuations in VOC/ethane ratios from plumes, which are both captured in the uncertainties reported in Table 2.

Several studies have investigated emission trends in the DJB between 2008 and 2021 yielding inconsistent results. Lyu et al. (2021) noted that ambient O&NG VOCs collected in canister samples in Platteville were lower in 2014 and 2016 than in 2013, with an overall propane and ethane decrease of 49% and 45% for 2013–2016, respectively. Oltmans et al. (2021) cautioned about using such a short time window for determining trends as the high year-to-year variability seen in ambient data is susceptible to misinterpretation. Their study was based on a longer, that is, 8-year NOAA data record from the BAO (2008–2016), filtered by wind direction and time of day to best reflect mixed boundary layer observations and DJB O&NG sources. It did not show any significant trends in methane and propane, which made these authors conclude that O&NG emissions had remained relatively constant.

Two recent investigations determined DJB methane and VOC trends from satellite-based remote sensing. Reddy and Taylor (2022) found a 55% and 73% decrease in O&NG methane and ethane emissions, respectively, between 2013 and 2020, with most of the decrease occurring after the data were collected that were considered by Oltmans et al. (2021). Lu et al. (2023) reported an 11% yr⁻¹ methane flux decrease in the DJB between 2010 and 2019.

In the following, we investigate the BRZ record for more recent changes in O&NG emissions between June 2017 and February 2021. With 4 years of data available for BRZ as of the writing of this article, the record is similar in length to the Lyu et al. (2021) study. However, there is a significant difference in the sampling strategy and data record. While Lyu et al. (2021) had to rely on a single sample collection every 6 days, yielding some 50–55 data points per year, the continuous sampling at BRZ has provided some 5,000–7,000 samples per year. This large body of data allows careful filtering of data for meteorological conditions for a more representative year-to-year comparison. Slope analyses for the BRZ data were explored for methane, ethane, and propane, which are deemed the most suited O&NG tracers. We refer to our investigation as rate of change and slope analyses, rather than trend analysis, as the term “trend” is reserved for instances when longer periods of data are available.

Two initial approaches were pursued to investigate the relative rate of change. In the first approach, to reduce the influence of the seasonal cycle and the start and ending date of the record, observations from particular months were compared over the record. Rates of change were then determined for the median values for each month (Supplement Figure S24), and the median and mean slope and intercept values for the 12 monthly determinations were then calculated. Methane at BRZ is 58.9 ppb higher than in the mostly background air at NWR (difference in the medians; Figure 4). This 58.9 ppb difference reflects contribution from the regionally added atmospheric methane. The overall best estimate of the methane rate of increase is 7.2 (median) and 7.6 (mean) ppb yr⁻¹ (with a standard deviation of 7.7 ppb yr⁻¹). This growth rate is on the same order as the methane growth of 7.0–7.7 ppb yr⁻¹ during 2016–2017 that has been reported for the global background (Nisbet et al., 2019), although the growth accelerated up to 15 ppb yr⁻¹ at the very end of our observation period (Lan et al., 2022). In the second approach, we apply a NOAA GML trend analysis tool (Thoning et al., 1989; Helmig et al., 2016) to the data record. This yields a best estimate of a 14.5 ppb yr⁻¹ rate of change (Supplement Figure S25). As these BRZ methane slope results fall above (result from the NOAA algorithm), and well within (result from the monthly rates of change) the global background range, there is no indication of a decline in regional methane emissions. Furthermore, the difference between BRZ and NWR methane did not show an obvious change over the 4-year time span. In summary, these BRZ versus NWR data comparisons do not point toward a change in the relative enhancement of methane at BRZ versus the (NWR) background, which suggests overall constant methane sources in the NCFR during this time window.

We scrutinized the BRZ methane data further by segregating them into samples collected when winds were from the O&NG sector (NE, 0°–120°) versus “clean” or background samples collected during westerly winds (240°–340°) and compared several percentiles at monthly binning intervals. These results were further broken up by increasing wind speed (Supplement Figure S26). The difference, that is, enhancement in the east sector data, increases with wind speed, which likely stems from more firmly defined and further extending footprint influences at higher winds, reflecting a stronger segregation of the 2 different source sectors. For data collected at winds >3 m s⁻¹, there is a ≈100 ppb enhancement in the air originating to the east, reflecting the contributions from stronger methane emissions within that sector.

Figure 16 shows monthly binned O&NG versus clean sector methane data using the 3 m s⁻¹ wind filter (N = 11,229 for O&NG and N = 9,252 for clean sectors), with box and whiskers showing the 5^th, 25^th, 50^th, 75^th, and 95^th percentile values for each month. Table 3 shows the linear regression best-fit slopes through these percentiles, along with the standard errors of the slopes. Additionally, the monthly averaged time series of the O&NG/clean sector ratios at each percentile are analyzed. We split this analysis into data from the entire year and data from the warm season only (June through September) to remove potential artifacts resulting from the high variability of winter methane. If there were a decrease in O&NG methane emissions, we would expect to see a convergence of clean and O&NG methane values with time in Figure 16, and the O&NG-sector fitted slopes to be lower relative to clean sector slopes. Our analysis shows no evidence of such a convergence, and in fact, most O&NG slopes in Table 3 are greater than clean sector slopes, albeit within one standard error of each other. Additionally, most of the O&NG/clean sector ratios are within one standard error of zero slope. The offset between the two datasets from these sectors shows no indication of change, that is, no evidence of declining methane concentrations that would indicate emission reduction from the O&NG sector between 2017 and 2021 in any of the filtered data.

O&NG emissions are reflected in the time series data primarily by the abundance of and resulting peak concentrations in plumes encountered at the monitoring station (see Section 3.4). Changes in emissions can therefore also be examined by an analysis of the variability in the data. This approach is explored in Figure 17. Data were again segregated into O&NG and clean sectors. Monthly median, standard deviation, and relative standard deviation were determined and plotted over the record. The analysis was conducted for all wind speeds, and for data progressively filtered by increasing wind speeds. The methane variability depicts a seasonal cycle for both groups of data, with higher variability in the winter, reflecting the higher number and higher concentrations in plumes encountered during the winter. While both the O&NG and clean sectors show similar behavior in the full data record, the difference in variability becomes progressively larger in the higher wind speed filtered data. This likely stems from more firmly defined and further extending footprint influences at higher winds, reflecting a stronger segregation of the 2 different source sectors. Most importantly, there is no evidence for a decrease in the methane variability in the NE sector data. This is further confirmed in a plot of the ratio of the NE/W variability (Supplement Figure S27). In conclusion, the wind sector segregated analysis of the BRZ data does not provide any signs of a decrease in methane emissions in the NCFR as a whole, and specifically from the DJB O&NG footprint area over the available record.

The methane trend results from the surface observation are inconsistent with the above-mentioned satellite-based analyses. There are several possible contributing causes to the discrepancies, for instance, the two measurements may represent different footprint areas of the DJB. Additional challenges of satellite studies may include heterogeneous land-use characteristics and complex terrain (Zhang et al., 2020), different smoothing algorithms applied to NWR surface samples versus the determination of monthly means from the Atmospheric InfraRed Sounder satellite data (for the Reddy and Taylor, 2022, study), and transport sources possibly dominating over the local emissions signal in larger grid cells (for the Lu et al., 2023, study). Surface monitoring data are collected during all weather and time of day conditions, including nighttime observations and capture conditions when ambient concentrations on average are highest and are therefore a sensitive tool for evaluating nearby emissions (Figures 3 and 6). Satellite observations largely rely on daytime and clear sky conditions when, due to boundary layer dynamics, absolute concentration values are mostly lower. More work is needed to fully resolve these discrepancies and characterize the biases between these methodologies. Considering the roughly 6-fold increases in O&NG production during the 2010s, the little or lack of decrease in overall methane emissions likely reflects a reduction in emissions per unit of produced oil and gas (Mead et al., 2024). However, for air quality considerations, there would have been little to no improvement because of the relatively small changes in atmospheric concentrations.

All of these trend results should be interpreted with caution. Four years of data, as used in our study, is a relatively short record, and we did not account for dispersion and transport characteristics in a comprehensive manner. The methane monitoring at BRZ has been continuous since 2017, future studies building on the expanded dataset will allow narrowing the uncertainties of the relative trend determinations, and would also benefit from applying Lagrangian transport analysis, and inverse modeling techniques, such as demonstrated by Lin et al. (2021) for the Uintah Basin.

Changes in observed VOCs were investigated using similar tools as described in the preceding section for methane. A slope analysis for ethane and propane using the NOAA algorithm is shown in Supplement Figure S28, regression analyses through monthly binned data in Supplement Figure S29. A summary of these 3 different approaches for determination of the rate of change, along with uncertainties where available, is provided in Table 4. These results show a consistent pattern of overall declining atmospheric VOC mole fractions of these O&NG tracers over the 4-year window, with rates of change, relative to 2017 mole fractions, in the range of 5%–9% yr⁻¹. Decreases in mole fractions were also found for other VOCs, for example, for benzene (−0.007 ± 0.003 ppb yr⁻¹) and toluene (−0.012 ± 0.004 ppb yr⁻¹) (Supplement Figure S29). The investigation of changes in the variability of the data (reflecting the enhancements in plume events) yielded inconclusive results as the month-to-month and year-to-year variability was too large to deduce systematic changes in the data (Supplement Figures S30 and S31).

The decreasing VOCs match trends seen in other U.S. urban areas (Warneke et al., 2012; Rossabi and Helmig, 2018; Liu, 2021). The data for the O&NG tracers may point toward decreasing emissions from the industry. However, care should be excercised in the interpretation of these results. Year-to-year variabilty can have a strong influence on slope calculations from short records. Further, it is plausible that emission changes during the COVID lockdown toward the end of the record may place a bias on these slope determinations. Reduction of other, less O&NG-influenced VOCs may also be explained by reductions in VOC emissions from the mobile sector (transportation, primarily automobiles) that have been observed throughout the United States (Rossabi and Helmig, 2018). It is not obvious why methane was relatively stable while ethane decreased between 2017 and 2021. Future studies should investigate how changes in agricultural or landfill related methane emissions from the O&NG wind direction sector at BRZ may be influencing the methane rate of change.

The VOC changes seen in the BRZ observations offer an opportunity for comparison with the State’s inventory emissions rates for the O&NG sector. In a similar recent effort, Oltmans et al. (2021) provided an evalation of the RAQC’s projected decline in O&NG VOC emissions by comparing BAO monitoring data with the relative rate of change in the inventory emissions. The projected inventory emissions decline is 6.5% yr⁻¹. Relative to the 2011 inventory emissions, during 2011–2020, emissions dropped to 43% of their 2011 levels. The best fit linear regression determined by Oltmans et al. (2021) through ambient measurement data of the O&NG tracer propane resulted in a slope value of −1.5% yr⁻¹ for their 2008–2016 data, which accounts for less than a quarter of the inventory rate of decline. Results from this study align more closely with the projected inventory emissions trend and may indicate a possible increasing rate of decline in the O&NG VOC emissions in the more recent years. However, for the reasons already emphasized above, care should be taken in the interpretation of the results from this realitvely short and COVID lockdown-influenced record.

This study utilized 4 years of continuous, high-time resolution methane and VOC monitoring data from BRZ to investigate ambient concentrations, contributions, source regions, and footprint area, including year-to-year variability and their relative rates of change. These data comprised the most extensive VOC measurements obtained in Colorado, with 4,000–6,000 VOC samples and over 50,000 methane measurements per year, and over 20,000 (200,000) individual air samples analyzed for VOCs (methane) in total to date. This abundance of data allowed the application of stringent filtering to increase the selectivity and to discern specific features, particularly those related to O&NG activities.

VOC mole fractions at BRZ are significantly greater than those seen in 25 of 28 urban areas studied by Baker et al. (2008). The VOC spectrum is dominated by light alkane petroleum hydrocarbons. Mole fractions of O&NG hydrocarbons exceed those of VOCs from other source categories, including traffic. Bivariate polar wind dependency analyses and PSCF analyses gave consistent conclusions: There is a stark difference in the VOC loading, with much higher average total VOC mole fractions present in air transported from the DJB O&NG sector, located to the north to southeast of BRZ, compared to all other directions. Elevated VOCs are mostly determined by the occurrence of concentration spikes, reflecting air plumes that can be enriched up to a factor of ≈100 over the background. Plumes last anywhere from one sample (≈1 h) to half a day. Several hundreds of these plumes were discernable in the BRZ data over the 4 years of observations.

Mole fractions of the most abundant VOCs follow a diurnal cycle due to boundary layer dynamics, with mean levels building up throughout the night, reaching a maximum during the early morning hours, and then gradually decreasing throughout the day. Buildup of emissions during nighttime hours result in ≈2–3 times higher average mole fractions than during afternoon hours. During winter, VOC concentrations are on average ≈2–3 times higher than during the summer. Rulemaking for setbacks of homes from O&NG wells and related operations should consider the full seasonal and diurnal range of VOC concentrations, rather than short-term daytime observations from summertime measurements.

We demonstrated a new approach for estimation of the total regional O&NG VOC emissions. VOC/ethane ratios in O&NG plumes were determined from the monitoring data and then scaled to DJB-wide top-down ethane flux estimates. This resulted in estimated year 2015 and 2021 total O&NG VOC fluxes (excluding ethane) of 183.6 ± 12.6 and 81.3 ± 16.3 Gg yr⁻¹, respectively. These estimates are approximately 2.0–2.5 times higher than the inventory projection for those years. More basin-wide aircraft observations of ethane, which now can be done more easily given that several brands of continuous ethane monitors are available, are highly desirable to reduce the uncertainty and update these flux estimates. More data on the spatial VOC/ethane distribution would help reduce the uncertainty of the total VOC flux determination. While a 4-year record is marginal for trend analyses, stringent filtering of data for wind sectors and wind speed allowed reducing the relative variance imposed by those variables. These analyses show a downward rate of change in the O&NG tracers ethane and propane, at a higher rate than recently reported for an 8-year earlier time span (Oltmans et al., 2021). It should be cautioned that reductions in traffic and industrial emissions from COVID lockdown measures may have contributed to somewhat lower emissions than what would have normally occurred. With 2021 being the last year of the observation window, any reductions would exert a particularly strong influence on the overall 2017–2021 rate of change calculation. These analyses should therefore be revisited in the future as more post-pandemic observations become available.

The discrepancies between the estimates derived from the observations versus inventory projections for total O&NG VOC fluxes reemphasize the need for critical review and improvement of inventory input data and the importance of continuous and high-time resolution observations and top-down data collection for the inventory assessment.

Data from the BRZ monitoring program are available through Boulder County Public Health (https://www.bouldercounty.org/departments/public-health/). VOC data have also been submitted to the EPA Ambient Monitoring Archive for HAPs (https://www.epa.gov/amtic/amtic-ambient-monitoring-archive-haps). NWR methane and VOC data were downloaded from the NOAA Global Monitoring Laboratory data server (https://gml.noaa.gov/dv/data/).

The supplemental files for this article can be found as follows:

Supplemental Material.pdf.

We thank the Colorado Department of Public Health & Environment (CDPHE) for maintaining the ozone measurements and infrastructure at the BRZ monitoring station and for accommodating the monitoring program described in this work. We thank David C. Carslaw and the openair team for producing and openly sharing openair package (Carslaw and Ropkins, 2012; Carslaw, 2019) in RStudio that was used for the bivariate polar plots and PSCF analyses. We also acknowledge NOAA for the historical 3 km HRRR back trajectory data and HYSPLIT model that were used to compute the PSCF outputs. Valuable insights and comments on an earlier manuscript draft were provided by Emily Fischer and Ilana Pollack, Colorado State University. We thank Benjamin Gaubert at University Corporation of Atmospheric Research for discussions during the review process.

The air quality monitoring at the Boulder Reservoir, the data analyses, and the manuscript preparation were funded by Boulder County Public Health.

The authors have no competing interests relating to the content of this publication. DH is one of the Elementa Editors-in-Chief. He was not involved in the review process of this manuscript.

Contributed to design and conception: DH.

Conducted data processing, quality control, and analyses: JH, GG, LSD, HA, JO, JC, BB, DC, SS.

Drafted and/or revised the article: DH, LSD, GG, JO, HA, DC.

Approved the submitted version of the publication: All authors.

How to cite this article: Helmig, D, Greenberg, G, Hueber, J, Blanchard, B, Chopra, J, Simoncic, S, Angot, H, Darby, LS, Ortega, J, Caputi, D. 2025. Methane and volatile organic compounds and their influence on air quality in Boulder, Colorado. Elementa: Science of the Anthropocene 13(1). DOI: https://doi.org/10.1525/elementa.2023.00117

Domain Editor-in-Chief: Alastair Iles, University of California Berkeley, Berkeley, CA, USA

Associate Editor: Gabriele Pfister, Atmospheric Chemistry Division, National Center for Atmospheric Research, Boulder, CO, USA

Knowledge Domain: Atmospheric Science

Part of an Elementa Forum: Oil and Natural Gas Development: Air Quality, Climate Science, and Policy