Spatial patterns in summertime surface ozone in the Southern Front Range of the U.S. Rocky Mountains

Ozone is a criteria air pollutant that threatens human health and biological activity. High levels of surface ozone put vulnerable populations like children, elderly, and individuals with preexisting respiratory conditions at risk of negative health effects including decreased lung function, worsening asthma, and higher rates of respiratory-related premature death (Berman et al., 2012; Fann et al., 2012). For this reason, the U.S. Environmental Protection Agency (U.S. EPA) established an ozone National Ambient Air Quality Standards (NAAQS) in 1997 requiring that the fourth highest maximum daily 8-h average (MDA8) ozone episode per year does not exceed 80 ppb (nmol mol^–1) over 3 consecutive years. This standard was lowered to 75 ppb in 2008 and to 70 ppb in 2015 (EPA, 2015).

While criteria pollutants across the United States, including ozone, have decreased due to the regulation and mitigation of anthropogenic emission sources, surface ozone has had spatially and temporally variable results on local and regional scales (Austin et al., 2015; EPA, 2020a). Most high elevation and rural sites across the country show either decreasing trends or no significant change in spring and summertime ozone since 2000 likely due to precursor emission controls, but the declines in Eastern United States ozone are larger in magnitude (Jaffe et al., 2018). Different trends at eastern and western sites have been in part attributed to contributions from baseline (ozone transported into the United States from outside sources; Parrish et al., 2017) or background (ozone formed from anthropogenic sources outside the United States plus natural sources; Jaffe et al., 2018) ozone, with both of these categories representing sources that are uncontrollable on the local level. The Western United States is especially susceptible to such inputs due to the interception of air masses that undergo trans-Pacific transport, regional wildfires, and propensity for episodic stratospheric intrusions (Parrish et al., 2010, 2012; Cooper et al., 2012; Lin et al., 2012a, 2012b; Jaffe et al., 2018). Synoptic-scale meteorological conditions and associated terrain-driven circulations also contribute to interannual variability in summertime ozone at sites across the region (Reddy and Pfister, 2016). Collectively, these factors can complicate the efficacy of local emission controls and the ability of some locations to maintain compliance with the recently lowered ozone NAAQS. Parrish et al. (2017) show that baseline ozone has decreased since the early 2000s, perhaps in part due to recent reductions in nitrogen oxides (NO_x) emissions from China, and likely aiding the ability for Western United States locations to achieve greater ozone reductions in the future. Nevertheless, deciphering the sources and controlling factors for surface ozone remain challenging in the Western United States.

The Front Range of the Rocky Mountains exemplifies this complex set of circumstances, where the Northern Front Range Metropolitan Area (NFRMA) has been designated as an ozone NAAQS nonattainment area (Flocke et al., 2020). Given the complex mixture of anthropogenic precursor emissions, particularly those from oil and gas wells and agricultural activities, background ozone contributions, and terrain-induced meteorology, a number of studies have explored the sources, chemistry, and distribution of ozone around the NFRMA (Langford et al., 2009; Brodin et al., 2010; Gilman et al., 2013; Reddy and Pfister, 2016; McDuffie et al., 2016; Abeleira, 2017; Abeleira and Farmer, 2017; Cheadle et al., 2017; Bien and Helmig, 2018; Flocke et al., 2020; Helmig, 2020). In contrast, surface ozone in the Southern Front Range, which has notably fewer oil and gas wells and where much of the agricultural activity is confined to the southernmost portion of the region along the Arkansas River Valley, has received less attention despite being a region with a growing population and facing similar meteorological and topographical complications as the region to the north.

Colorado Springs, situated in the Pikes Peak region of the Southern Front Range and approximately 100 km south of Denver, is the second largest metropolitan area in the state of Colorado. In the past decade, El Paso County, home to Colorado Springs, has grown by more than 100,000 people, with another 300,000 expected by 2050 (State Demography Office, 2019). The resulting increase in motor vehicle use, as well as the city’s proximity to the Denver metropolitan area and the broader Denver–Julesburg basin, has the potential to intensify ozone production in this region. Notably, emissions of ozone precursors from inventoried point sources in El Paso County have declined in the past decade, with reported reductions of 51% and 47% in NO_x and volatile organic compounds (VOCs), respectively (Figure 1; EPA, 2013, 2014, 2018, 2020b). Simultaneously, NO_x from moving polluters decreased by 27%, and VOCs from moving polluters decreased by 10% (Figure 1; EPA, 2013, 2014, 2018, 2020b). According to the 2017 National Emissions Inventory (NEI), the highest point source emitters of NO_x in El Paso County are two coal-fired power plants, while local airports are the largest emitters of VOCs. All point sources reported to the NEI for El Paso County are shown in Figure 2. In 2018, the 3-year average of the fourth highest MDA8 ozone in the Pikes Peak region reached 70 ppb for the first time since the ozone NAAQS was lowered in 2015, while in 2019, the region remained slightly below the standard (Colorado Department of Public Health and Environment, Air Pollution Control Division, 2019). Based on observations from 2018 to 2020, the region is now at risk of nonattainment (Colorado Department of Public Health and Environment, Air Pollution Control Division, 2020).

Here, we analyze summertime ozone measurements from Colorado Springs and the surrounding Pikes Peak Region, beginning with an analysis of long-term trends in MDA8 at the two permanent monitoring sites in the area. We then analyze data from a special study during the 2018 ozone season, in which the Colorado Department of Public Health and Environment (CDPHE) added six temporary ozone measurement stations in order to better characterize the temporal behavior and spatial distribution of summertime ozone in this area of the Southern Front Range. We further consider the relationships between surface ozone and observed meteorological conditions through the use of generalized additive modeling (GAM) for the two permanent monitoring sites. We then assess the representativeness of these two permanent monitoring sites during high ozone events in the Pikes Peak region.

The CDPHE operates two permanent ozone monitoring stations in the Pikes Peak region. Ozone measurements began at the Air Force Academy (AFA) site in 2000 and at the Manitou Springs (MAN) site in 2004 (Colorado Department of Public Health and Environment, Air Pollution Control Division, 2019). We use data from June 1 through August 31 for each year between the starting years and the year 2020 for the calculation of long-term trends in summertime MDA8 ozone (Section 3.1). The choice of June to August for this trend analysis offers consistency in comparison with other trend analyses reported in the literature. We also use MDA8 ozone measurements from these two sites from June 1 through September 30, 2014–2020, in GAMs to assess the relationships of summertime ozone with meteorology (Section 3.3). The choice of June to September for this analysis was for consistency with the special ozone study performed in the region during 2018.

From June 1 through September 30, 2018, CDPHE added six temporary ozone measurement sites as part of a special study designed to characterize the spatiotemporal distribution and behavior of summertime ozone in Colorado Springs. The six sites are Black Forest (BLF), Cascade (CAS), Colorado College (CC), Monument (MON), Navigator (NAV), and Rampart (RAM). Together, these six sites plus the two permanent ones represent a variety of elevations and proximities to ozone precursor emitters such as highways (Table 1; Figure 2). We use data from all eight sites and from June 1 through September 30 to examine spatial patterns and diurnal variability in summertime ozone (Sections 3.2), case studies of high ozone (Section 3.3), and representativeness of the two permanent monitoring stations during high ozone events (Section 3.4).

The CDPHE operates Teledyne API (TAPI) model 400 ozone analyzers at the permanent air monitoring sites (AFA and MAN). Hourly data are readily available online through the US EPA Air Quality System (AQS). The TAPI 400 has a detection limit of 0.4 ppb, and a measurement precision that is greater than 0.5 ppb or 0.5% of the reading. Calibration and quality assurance and quality control (QA/QC) procedures follow the CDPHE Quality Assurance Project Plan (QAPP; CDPHE, 2015).

For the 2018 special study, CDPHE operated 2BTech model 205 ozone analyzers at all temporary sites except CC where a TAPI 400 was temporarily installed. Data were collected using a Campbell Scientific data logger that generated hourly mole fractions averaged from ozone measurements detected with 6-s time resolution. The 2BTech ozone analyzers used in this study have a 1-σ measurement precision that is greater than 1.0 ppb or 2% of the reading (2B Technologies, 2018). At all stations, sampling heights were 2 m above ground, and inlet filters were changed either biweekly or monthly. The 2BTech analyzers, power systems, data loggers, modems, and zero-air charcoal containers were contained in weatherproof enclosures. A CPU-style fan located near the top of the structure beneath a protected opening was activated when the internal temperature was warm enough to warrant power consumption. All instruments were placed above the bottom of the enclosure to avoid any potential exposure to water.

The 2BTech analyzers are special purpose monitors by CDPHE’s Air Pollution Control Division’s (APCD) guidance. As such, stations were calibrated following the APCD 2BTech analyzer standard operating procedures, which are approved under the CDPHE QAPP (CDPHE, 2015). Both temporary and permanent monitors have federal equivalent methods (FEMs) designation. All monitors are independently audited following CDPHE QAPP methods and are held to the same passing criteria. In accordance with the CDPHE QAPP guidance (CDPHE, 2015), stations were calibrated at least once per month during the 2018 study period using a combination of 2BTech model 306 and Teledyne model 703 sources. Nightly baseline checks were also performed at all sites with the 2BTech analyzers, resulting in no data available at the 23:00 h at those sites. Tables S1a and S1b provide results of nightly baseline checks and temporary site audits. Due to the verified stability of the 2BTech analyzers over the length of this study period, these instruments were not post-processed for drift, and calibration corrections were not applied. While the data collected using temporary monitors by the APCD are considered comparable in quality to those collected by permanent monitors, the temporary monitors do not fall under the quality assurance umbrella necessary for regulatory decisions. Argent ADS-WS1 weather stations with wind vane and cup anemometer (Argent Data Systems, 2010) were placed at four of the eight sites (BLF, CAS, MON, and NAV) during the summer 2018 study to measure hourly averaged wind speed and direction.

We generated summary statistics, linear regressions, and independent samples t tests using a combination of Excel v15.13.1, IBM SPSS Statistics v25, and RStudio v1.2.1335. We consider statistical significance at p values < .05 in this study. We used ArcMap 10.6 GIS software to spatially visualize and interpret several data sets, including point sources of precursor emissions in Colorado Springs, wind roses, air mass back-trajectories for high ozone events, and spatial interpolation of ozone mole fractions for selected events.

We used the “open-air” package in RStudio to generate wind roses for ozone by direction at the four sites with local wind speed and direction data in the 2018 special study (BLF, CAS, MON, and NAV). We binned hourly ozone data between 7:00 through 18:00 MST to generate daytime wind roses and 19:00 through 6:00 MST for nighttime wind roses at each site. All daytime wind roses have greater than 50% data availability, while nighttime wind roses with less than 50% data availability are not shown. We also generated wind roses for the 0:00, 4:00, 8:00, 12:00, 16:00, and 20:00 h of the day for the same sites (Figure S3). Ozone measurements in Figure S3 are plotted only for each labeled hour (e.g., ozone measurements at hours 1:00–3:00, 5:00–7:00, 9:00–11:00, 17:00–19:00, and 21:00–23:00 and corresponding wind direction are not shown).

Statistical modeling and machine learning have been used previously to examine the complex relationships between ozone and various meteorological factors. For example, Camalier et al. (2007) used a generalized linear model to look at these relationships. GAM (Wood, 2017) is a similar approach but can incorporate linear, nonlinear, and categorical predictors. Gong et al. (2017) and McClure and Jaffe (2018) used GAMs to quantify the contribution to the ozone MDA8 from wildfire smoke, using data from numerous sites in the Western United States and we follow a similar approach here.

We use GAMs to explore relationships between MDA8 ozone and meteorological conditions at the two permanent monitoring sites, AFA and MAN, from June to September 2014–2020. We use hourly meteorological data, including temperature, relative humidity (RH), wind speed, and wind direction, from CDPHE’s Highway 24 permanent monitoring station in Downtown Colorado Springs. The site is 10 m north of Colorado State Highway 24 and 400 m west of Interstate 25 (Figure 2). Meteorological measurements began on August 20, 2014, at this site; therefore, the GAMs do not include data for June 1 to August 19, 2014. We computed daily minimum, maximum, and mean RH as well as daily minimum and maximum temperature using hourly data from 0:00 to 23:00 MST. We computed average wind speed and vector average directions from 6:00 to 14:00 MST. The GAMs reported here use three meteorological predictors (daily maximum temperature, daily mean RH, daily mean 500 hPa geopotential height—obtained from the National Oceanic and Atmospheric Administration [NOAA], Physical Sciences Laboratory NCEP/NCAR Reanalysis data product; Kalnay et al., 1996), and two categorical predictors: month and local wind direction quadrant (Table S3). The numeric variables are fit using penalized cubic regression splines and using an identity link with a Gaussian distribution. The GAMs were computed using the “mgcv” package in R software. To evaluate the GAMs, we randomly subset the data into a training set (90% of data) and a cross-validation set (10% of data). The cross-validation data provide an important confirmation of the model predictions, using data that were not part of the training set. We also subset the data by smoke and nonsmoke days. Classification of smoke days is described in Section 2.4.

We use an inverse distance weighted (IDW) Gaussian model within ArcMap 10.6 to interpolate ozone between sites for selected periods and assess the representativeness of the permanent sites with respect to the rest of the monitored region (Section 3.4). IDW interpolations are computed as a function of the distance between real observations and the site where a prediction is being made (Wong et al., 2004). The tool relies on the inverse distance between points raised to a mathematical power (Environmental Systems Research Institute, 2016). For our IDW interpolations, we use a power of 2 to mimic the behavior of a secondary pollutant such as ozone and keep the influence of a single point spatially broad. The IDW interpolation is appropriate for assessing the spatial distribution of ozone because ozone is less localized than other directly emitted polluters (Wong et al., 2004; Berman et al., 2015). In a study of interpolation techniques for air quality parameters, Jha et al. (2011) found IDW to be a more suitable evaluation method compared to kriging, especially for assessing a network of stations such as that in the Pikes Peak region. Prior knowledge of the behavior of ozone with respect to elevation is helpful as the IDW tool does not take variable topography into account and thus functions best when looking at horizontal spatial distributions. Because of these limitations, IDW is a viable approach when strategizing about monitoring station placement for the protection of human health rather than directly modeling ozone by taking topography and meteorological variables into account. We performed tests of representativeness on the highest ozone events of summer 2018, including those with identified underlying causes such as smoke from regional wildfires, stratospheric intrusions, or outstanding meteorological events because representativeness of a site only concerns ozone values as they relate to one another, not as they relate to the overall cause of the ozone values.

Smoke days from June 1 to September 30 of 2014 through 2020 are identified using a combination of satellite data and surface PM_2.5 (particulate matter < 2.5 μm) observations. We use the NOAA Hazard Mapping System (HMS) Fire and Smoke Product (AirNow-Tech, 2021) as described by Kaulfus et al. (2017) and Buysse et al. (2019), which is a satellite-derived product to identify overhead smoke plumes. Surface PM_2.5 observations are from CDPHE’s permanent monitoring station at CC (available through the EPA AQS; data for 2014–2015 are collected every third day and for 2016–2020 are collected daily). The HMS product is first used to identify days with overhead smoke in Colorado Springs. Using this information, we compute the mean and standard deviation of surface PM_2.5 by month for nonsmoke days; days with overhead smoke according to the HMS and with PM_2.5 concentrations more than one standard deviation above the nonsmoke day mean are considered to have smoke-impacted surface air.

To identify days with elevated ozone during the 2018 ozone study, we looked for instances when seven of the eight sites had an 8-h running mean greater than 70 ppb. When that criterion was met, we considered the total high ozone events to span from the time when the first site exceeded 70 ppb on a running 8-h average to when the last site fell below 70 ppb on a running 8-h average. We only required seven of eight sites in our identification of high ozone events because the CC site never exceeded 70 ppb on a running 8-h average. This metric reveals five high ozone events during the 2018 study (June 12, July 6, July 14, July 18, and July 31 to August 2, 2018). Requiring six of eight sites adds one additional high ozone event, while requiring five of eight sites adds one further episode. Although this short-duration, single-season study cannot address the Pikes Peak region’s compliance (or lack thereof) with the NAAQS, we use the aforementioned 70 ppb metric as a broad representation of the standard and to identify high ozone events. For each high ozone event in 2018, we use the NOAA Hybrid Single-Particle Lagrangian Integrated Trajectory model (Draxler and Rolph, 2014) to calculate 72-h back-trajectories using the Eta Data Assimilation System 40 km gridded meteorological data. We use the coordinates of AFA as the starting location and initiate trajectories from three starting heights: 500 m; 1,000 m; and 1,500 m above ground level (AGL). Starting times for trajectories are the peak ozone hour of each high ozone event.

The aforementioned smoke analysis reveals that three of the five high ozone events during 2018 were impacted by regional smoke (June 12, July 6, and July 31 to August 2, 2018). Meanwhile, in 2020, all of the top five MDA8 days at AFA and MAN were influenced by smoke. The impacts of smoke will be discussed further when we discuss the GAM results. Although smoke events may be considered as exceptional events, we consider them in our study because they demonstrate the behavior of ozone at a critical level that is important for public health in the Pikes Peak region.

Figure 3 displays the 5th, 50th, and 95th percentiles of summertime (June, July, and August) MDA8 ozone at the MAN (2004–2020) and AFA (2000–2020) permanent monitoring sites. Trends are not statistically significant at either site (Figure 3). Median MDA8 ozone at both sites ranges between approximately 50 and 60 ppb from year to year. These values are similar or slightly lower than median daytime (10:00–16:00 MDT) values (50–70 ppb) at sites in the NFRMA during summertime between 2000 and 2015, many of which also do not exhibit significant trends; the exceptions were the CAMP and Welby sites in metropolitan Denver where positive trends were attributed to long-term decreases in NO_x emissions in a NO_x-saturated environment (Abeleira and Farmer, 2017). Bien and Helmig (2018) similarly found significant positive trends in summertime hourly (0:00–23:00) ozone at only four of 20 sites in the NFRMA from 2000 to 2015. The 5th percentile is often considered in long-term trend analyses because at rural background sites, it is shown to be more sensitive to changes in baseline ozone or changes in NO_x emissions from nearby urban areas (Cooper et al., 2012). At AFA and MAN, the 5th percentiles of summertime MDA8 values have increased over the past 20 years, but the trends are not significant. The magnitude of the 5th percentile at each of the two sites (approximately 40 ppb) is similar to the 5th percentile of daytime ozone at urban and semiurban sites during summer in the NFRMA where most sites also did not show significant trends (Abeleira and Farmer, 2017). Meanwhile, Bien and Helmig (2018) found significant increases in the 5th percentile of summertime hourly ozone at 11 of 20 NFRMA sites. Values for the 95th percentile of summertime MDA8 ozone in the Pike Peak region (approximately 60–80 ppb) are again similar or slightly lower than 95th percentiles of daytime ozone at sites in the NFRMA (60–90 ppb) considered by Abeleira and Farmer (2017). Bien and Helmig (2018) showed no significant trends at 18 of 20 sites for the 95th percentile of summertime hourly ozone, results that are similar to the findings from Strode et al. (2015) for the Mountain West.

We acknowledge that interannual variability in regional meteorology can affect ozone levels in any given year and thus influence long-term trends. Normalizing ozone levels against synoptic-scale meteorological conditions can be a particularly important tool for agencies seeking to determine the efficacy of control measures when they are put in place (Camalier et al., 2007; Reddy and Pfister, 2016). We do not apply such an approach here as our work is primarily concerned with trends in surface ozone as compared to EPA standards, the exposure levels of local populations, and determining where monitoring is warranted for the protection of human health. We do, however, minimize the effects of interannual variability in the trends reported in Figure 3 through the use of multidecadal data sets (21 years at AFA and the 17 years at MAN) and seasonal (summer) averages, which are shown to be valuable strategies in other long-term ozone studies (Cooper et al., 2012; Parrish et al., 2012; Bien and Helmig, 2018).

Mean (±1σ) ozone mole fractions at the eight sites during the 2018 Colorado Springs ozone study ranged from 34 ± 19 to 60 ± 9 ppb (Figure 4; Table S2). In the entire data set, the highest 1-h measurement was 96 ppb at RAM. Table S2 provides summary statistics for all eight sites during the 2018 special ozone study. Six of eight sites exhibit normal distributions in hourly ozone observations; the exceptions are AFA and CC, which have higher frequencies of low ozone counts (Figure S2). All hourly data distributions are significantly different from one another based on independent samples t tests, with the exception of MON and NAV (p = .87), and all other site means differ by a range greater than the reported instrument precisions.

For all sites, ozone is, on average, lowest between 4:00 and 6:00 MST and highest between 11:00 and 13:00 MST as expected due to surface deposition and photochemically driven ozone production and loss cycles (Figure 5). On average, the CC site consistently has the lowest ozone across the day, while RAM has the highest. We also detect variability in the amplitudes of the diurnal patterns among sites and their respective elevations. The diurnal amplitude at the highest elevation RAM site (2,840 m ASL) is only 9 ppb and is comparable at BLF (2,233 m ASL; 11 ppb). The CAS (2,445 m ASL), NAV (2,220 m ASL), MON (2,206 m ASL), and MAN (1,960 m ASL) sites display diurnal amplitudes of 15–20 ppb. In contrast, AFA and CC (both 1,971 m ASL) display more pronounced diurnal curves of 35 and 40 ppb, respectively. Although these two sites are closest to the Interstate 25 traffic corridor and city center (Figure 2; Table 1), both follow a typical summer diurnal ozone pattern and do not deviate during rush hour periods in comparison to the lowest elevation site in the Boulder, CO, vertical transect that saw two daily maxima (15:00–16:00 and 1:00–3:00 MST) and two daily minima (7:00–8:00 and 19:00–22:00 MST) driven by the confluence of daily cycles in traffic volume as well as ozone sources and sinks and diurnal boundary layer evolution (Brodin et al., 2010).

The average rates of change in ozone with respect to each hour of the day at each site (Δ[O₃]/Δt; Figure 6) show that AFA and CC, two of the lowest elevation sites that are also closest to the Interstate (Table 1), have faster rates of ozone increase in the morning (9 and 10 ppb h^–1, respectively, from 6:00–7:00 MST). Increasing ozone begins at MON between 4:00 and 5:00 and at NAV, BLF, and CAS between 5:00 and 6:00, but their peak rates of change between 6:00 and 8:00 MST are notably smaller (3–4 ppb h^–1) as compared to AFA and CC. In contrast, MAN and RAM experience larger increases in their rates of change later in the morning than the other sites, with peak ozone increase occurring between 9:00 and 10:00 MST (6 and 3 ppb h^–1, respectively). In the evening, AFA and CC again display larger average rates of change in ozone (e.g., –6 and –7 ppb h^–1, respectively, at 18:00–19:00 MST, while other sites range from –1 to –3 ppb h^–1).

Wind direction also shows pronounced diurnal patterns during the summer 2018 study, with anabatic (upslope) winds during daytime hours and katabatic (downslope) winds during nighttime hours, particularly at sites with mountainous topography sloped toward Colorado Springs’ city center (e.g., CAS and NAV; Figures 7 and S3). Three of four sites (CAS, MON, and NAV) have high frequencies of elevated ozone coming from the city center under daytime upslope flow (Figures 7 and S3). In contrast, BLF shows a high frequency of elevated ozone originating from the southeast, perhaps due to localized topography where a ridge northwest of the site may alter the orientation of daytime upslope winds; however, there is not sufficient nighttime data to fully characterize diurnal wind patterns at BLF. Daytime hourly ozone is only weakly correlated (R² < .2, p < .001) with hourly wind speeds at each of the four sites that have local wind measurements, indicating that wind speed is not a particularly strong predictor of surface ozone levels in these areas.

These observations reflect the controls on diurnal ozone variability discussed in similar studies, including nighttime ozone depletion via surface deposition and NO titration, rapid midmorning ozone increases via boundary layer entrainment followed by in situ photochemical production, and diurnal winds facilitating horizontal advection (Jaffe et al., 2018; Trousdell et al., 2019; Oltmans et al., 2019). These effects can be seen to different degrees across the range of site elevations. For example, the relatively lower overnight minimum ozone and large diurnal amplitudes at AFA and CC point to the increasing effect of surface deposition at lower elevation valley sites that may experience strong temperature inversions in combination with ozone loss via reaction with local NO; in contrast, higher elevation sites like RAM that are further from NO_x emission sources see higher amounts of overnight ozone survival (Cooper and Peterson, 2000; Tong et al., 2006; Brodin et al., 2010; Bien and Helmig, 2018; Jaffe et al., 2018). Stronger inversions prevent nighttime boundary layer mixing at sites in the inversion layer, allowing mechanisms like surface deposition to become a greater surface ozone sink than at higher elevation sites (Aneja et al., 2000). Higher elevation sites can be above the nighttime inversion and so experience significantly less ozone loss at night (Aneja et al., 2000; Ambrose et al., 2011). Although our study design does not permit us to directly determine the magnitude of depositional loss processes or how it may vary across sites, Trousdell et al. (2019) propose that in the San Joaquin Valley, much of the ozone loss at dusk, in the absence of photochemistry and entrainment, is in fact due to surface deposition and not simply titration by rush hour NO emissions. Caputi et al. (2019) found overnight losses of approximately 1.2 ppb h^–1 in California’s San Joaquin Valley.

Meanwhile, rates of ozone increase occur at most sites beginning around 5:00–6:00 MST (Figure 6) due to boundary layer growth and daytime turbulence causing ozone from the residual layer to be mixed with air at the surface across at most sites (Brodin et al., 2010; Fine et al., 2015; Bien and Helmig, 2018; Oltmans et al., 2019; Trousdell et al., 2019). In the NFRMA and the San Joaquin Valley of California, entrainment has been observed to cause morning increases in ozone values with minimal to no input from photochemical production from 6:00 to 9:00 MST (Trousdell et al., 2016; Oltmans et al., 2019), and we assume similar dynamics govern the diurnal profiles in Colorado Springs. Early morning boundary layer mixing followed by daytime photolysis are accompanied by diurnal mountain wind patterns; for example, at NAV and CAS, the nighttime downslope flow shifts to daytime upslope flow between 5:00 and 8:00 MST (Figures 7 and S3). Ozone and its precursors may be transported upslope and continue to build on ozone from entrainment from midmorning until the evening, when surface deposition and residual NO emissions act as a sink for ozone at the lowest elevation sites nearer the city center (Aneja et al., 2000; Brodin et al., 2010; Bien and Helmig, 2018; Caputi et al., 2019). For the highest elevation RAM site, the delay in ozone increase to 8:00–10:00 MST may suggest it is situated above the nighttime boundary layer, while also indicating the time necessary for the transport of precursor emissions and morning photolysis to begin (Mueller, 1994; Brodin et al, 2010; Oltmans et al., 2019). As noted above, mountain observation sites may be above the boundary layer at night depending on the elevation and local topography (Ambrose et al., 2011). At MAN, the delay in morning ozone increase may be due to local topography. Although it is the lowest elevation study site, it is further from the Colorado Springs city center than AFA or CC (Table 1) and sits in a small basin in the foothills of Pikes Peak (Figure 2), which may inhibit the transport of local NO_x emissions; however, without local wind or ozone precursor measurements at MAN, this interpretation is only speculative.

Among all eight sites, the 5th percentile, mean, median, and 95th percentile of hourly ozone increase significantly by (slope ± 95% CI, CI = confidence interval) 3.4 ± 3.2, 2.0 ± 1.8, 1.8 ± 1.6, and 1.1 ± 0.9 ppb, respectively, with every 100 m of elevation (Figure 8). The mean and median vertical ozone gradients are slightly higher than the range of values (1.2–1.7 ppb per 100 m) reported in other studies of vertical transects in primarily rural mountainous areas (Peterson et al., 1999; Cooper and Peterson, 2000; Ribas and Penuelas, 2006; Brodin et al., 2010). More specifically, Brodin et al. (2010) reported an increase of mean ozone with elevation at seasonal rates of 1.1–1.7 ppb per 100 m of elevation at 10 sites near Boulder, CO; the average increase for summertime was 1.5 ppb per 100 m but was not monotonic, with a smaller rate of change found above 2,400 m ASL. Although we do not observe such a pattern in Colorado Springs, our study consists of fewer sites that are more horizontally dispersed compared to those in Brodin et al. (2010) and only two of which are above 2,400 m ASL. The overall observation of higher ozone with increasing elevation is consistent with other studies and reflects the background tropospheric gradient of increasing ozone with elevation, largely due to reduced surface deposition and longer ozone lifetime at higher elevation sites (Cooper and Peterson, 2000; Oltmans et al., 2008; Brodin et al., 2010; Cooper et al., 2011; Jaffe et al., 2018; Caputi et al., 2019). For example, median values from an ozonesonde site in the NFRMA show an increase of approximately 10 ppb from 100 to 2,000 m AGL (approximately 0.5 ppb per 100 m; Cooper et al., 2011). A similar gradient is reported by Jaffe et al. (2018) for ozonesondes at remote sites in California, Oregon, and Washington. Fine et al. (2015) also showed a similar gradient at sites in Nevada and attributed the higher ozone at higher elevations to global and regional pollution. So, while there may be local influences due to increasing distance from the traffic corridor and local anthropogenic precursor emissions, the observed elevation gradients in Colorado Springs are likely driven by the same factors that produce these vertical gradients at other comparable locations.

We use GAMs to investigate the relationship among summertime ozone, meteorology, and smoke influence in Colorado Springs from June to September 2014–2020, using MDA8 ozone from AFA and MAN, the two permanent sites with long-term ozone records. Only three numeric variables gave significant predictive skill: daily maximum temperature, daily mean RH, and daily mean 500 hPa geopotential height. Two categorical variables, month and surface wind direction quadrant, also had a small degree of predictive skill (Table S3). The data indicate a weak positive, but statistically significant (p < .001), correlation with MDA8 at both AFA and MAN for daily maximum temperature (R² = .32 and .35) and daily mean 500 hPa geopotential height (R² = .15 and .12) and an inverse relationship with daily mean RH (R² = .14 and .10). The GAM integrates across these relationships to provide an overall prediction for the MDA8. The GAMs were set up as described in Section 2.3.2, and detailed results are shown in Tables S4–S7.

For both AFA and MAN, model quality control included the examination of the residuals for heteroscedasticity and bias against the input predictors. Further model predictions using the cross-validation data set should show similar predictive skill as the training data set. Finally, the mean residuals (observed MDA8–predicted MDA8) for any GAM analysis must be very close to zero. For both AFA and MAN, all QC tests were favorable, and results using the cross-validation data gave similar or better predictive skill compared to the training data set (see Tables S4–S7). For the remainder of this discussion, we reran the GAMs using all data (training and cross-validation).

The R² comparing the observed and predicted MDA8s was .45 and .42 for MAN and AFA, respectively. These R² values are lower than values reported for more urbanized environments, where R² values of between .52 and .81 were reported by Gong et al. (2017) for a number of cities in the Western United States. We interpret the lower R² to reflect a greater influence from regional sources in the Pikes Peak region compared to local photochemical production. Additionally, we find that despite the higher mean ozone at altitude, all sites, with the exception of CC, show strong correlations with one another using the 2018 MDA8 values (Table 2). These correlations also likely reflect the importance of regional control on ozone rather than local control. This relationship points to the fact that the MDA8, by definition, looks at the peak ozone period each day (generally the midday hours), and as shown in Figure 5, these are more similar than the daily mean across all sites. In the case of local production, variables such as daily maximum temperature and wind speed are expected to show stronger predictive power than we find for the Pikes Peak region. Nonetheless, the GAMs can still provide important insight into the mechanism of high ozone days.

The standard deviation of the GAM residual, 6.4 ppb for both MAN and AFA, provides a measure of the overall model performance. While the overall mean GAM residual is essentially zero (<1 × 10^–12), the smoke days have a significant positive mean residual of 4.3 and 4.7 ppb for MAN and AFA, respectively, suggesting that wildfire emissions of ozone precursors add 4–5 ppb to the MDA8. Given the large number of smoke days identified in 2018 and 2020, this has a significant implication for attaining the NAAQS. For both AFA and MAN, three of the five highest MDA8 days in 2018 had smoke influence, and for 2020, all of the top five highest MDA8 days had smoke influence.

For nonsmoke days included in the 2014–2020 analysis, we find that the 95th percentile of MDA8 ozone days at AFA occurred with significantly (p < .001) higher daily maximum temperatures (mean ± 1σ = 30 ± 2.9 °C), lower daily mean RH (43 ± 12%), and higher daily mean 500 hPa geopotential height (5,911 ± 50 hPa) as compared to the 5th percentile of nonsmoke MDA8 (20 ± 7.5 °C; 66 ± 19%; 5,845 ± 60 hPa). Similar values are determined for the 95th and 5th percentile nonsmoke MDA8 days at MAN (95th percentile: 30 ± 2.9 °C; 43 ± 12%; 5,917 ± 52 hPa; 5th percentile: 17 ± 6.7 °C; 71 ± 18%; 5,828 ± 53 hPa). Additionally, we calculated Δ[O₃]/ΔT slopes using reduced major axis regression on the MDA8 and daily maximum temperature values. The slopes for AFA are 1.8 and 2.1 ppb °C^–1 for nonsmoke and smoke days, respectively. The values for MAN are very similar at 1.8 and 2.2 ppb °C^–1 for nonsmoke and smoke days, respectively. The higher slope on smoke days is indicative of stronger ozone production on these days. These results are consistent with other studies and confirm that in the absence of regional smoke, elevated summertime MDA8 ozone in the Pikes Peak region tends to occur on warmer, drier days associated with upper atmosphere ridging (Jaffe, 2021).

Using the metric described in Section 2.4, which considers hourly ozone at AFA and MAN as well as the six temporary sites, we identify five high ozone events during the 2018 ozone study: one in June (June 12, 2018) and four in July (July 6, July 14, July 18, and July 31 to August 2, 2018). Events range from 9 to 67 h in length, and peak ozone mole fractions at the eight sites range from 85 to 96 ppb (Table S8). Three episodes (June 12, July 6, and July 31 to August 2, 2018) have evidence of smoke at the surface due to regional wildfires (Figures S4, S5, and S8). On July 6, 2018, the HMS smoke product shows five active fire complexes in Colorado alone and others across the Mountain West, with smoke covering the entire northeast half of the state and elevated ozone observed along the Front Range (Figure S5). Back-trajectories depict south–southeasterly winds related to a large surface high-pressure system over the Midwest, northeast of the study area, which by the following day had moved further east and a surface low-pressure system moved into the area from the northwest. The July 31 to August 2, 2018, event (Figure S8) was a similar region-wide smoke event, with stagnant conditions on August 1 and August 2 and elevated ozone again observed along the Front Range, particularly on August 2, 2018. In contrast, the June 12, 2018, smoke event (Figure S4) occurred under rapid westerly flow with less stagnation, and sites in Denver reported slightly lower ozone levels than in Colorado Springs. According to CDPHE meteorologists, there was a stratospheric intrusion during the June 12, 2018, event in addition to surface smoke.

On July 14, 2018, and July 18, 2018, when the HMS products do not show smoke over the study area (Figures S6–S7), surface meteorology is consistent with the aforementioned pattern of high ozone occurring on warmer, drier days. On these two days, for example, daily maximum temperature was above 30 °C, daily mean RH was around 50%, and daily mean 500 hPa geopotential height was approximately 5,950 hPa. Upper air soundings on the afternoons of July 14, 2018, and July 18, 2018, show a very dry lower atmosphere with a subsidence inversion around 500 hPa, while the 500 hPa geopotential height map shows upper air ridging over the U.S. Southwest. On July 14, 2018, a high-pressure system was situated over Colorado, which is depicted by the clockwise flow near the Pikes Peak region; the sharp turn in the trajectory around 12:00 UTC on July 13, 2018, is related to a cold front moving over the Great Plains, while another cold front advancing southward on July 15, 2018, coincides with the decline in regional ozone. Northwesterly flow on July 18, 2018, is similarly related to a surface high-pressure system southwest of Colorado Springs. During both the July 14, 2018, and July 18, 2018, events, 8-h average ozone was elevated throughout the Front Range (Figures S6–S7). These conditions appear to exemplify conclusions drawn from the larger data set, that in the absence of smoke, high ozone days tend to occur on days with higher temperatures, lower RH, and upper air ridging—conditions that affect surface ozone beyond the local environment of the Pikes Peak region.

Here, we consider the representativeness of the AFA and MAN permanent sites during the highest ozone events of the 2018 study. We use an IDW Gaussian model interpolation based on-site location and peak hourly ozone during the five 2018 high ozone events to visualize the spatial distribution of ozone relative to monitors in the region (Section 2.3.3). We then independently remove the AFA and MAN data points and rerun the IDW analysis. The difference between the actual permanent site ozone mole fraction data and the predicted ozone mole fraction data by the IDW model is used to assess whether or not the data recorded at permanent sites are representative of the surrounding region during episodes of high ozone across the Pikes Peak area. It is important to emphasize that the results from the IDW tool are intended to aid in assessing the location of monitoring sites relative to one another to be more protective of human health and not to directly predict ozone as a model tailored to local topography and meteorology might do.

Overall, for the five high ozone events of the 2018 study, we find that the IDW overpredicts ozone mole fractions at permanent sites by 1–5 ppb. This range is within the established CDPHE quality control criteria (Colorado Department of Public Health and Environment, Air Pollution Control Division, 2015). We also use linear regressions of peak hourly ozone during high ozone events with respect to elevation to test the vertical representativeness of permanent sites, again removing the permanent sites from the regression one at a time and observed the change in slope. We find that in all five instances, removing each of the permanent sites results in a steeper slope, with an average change of 0.003 or 0.3 ± 0.1 ppb per 100 m. This indicates that the linear regression slightly underpredicts ozone values for the elevations of AFA and MAN by 0.3 ± 0.1 ppb per 100 m of elevation. This is also within established CDPHE quality control criteria (Colorado Department of Public Health and Environment, Air Pollution Control Division, 2015). Figures 9 and 10 show the July 14, 2018, to July 15, 2018, event as an example of this exercise.

These models suggest that within the scope of the study, the permanent sites are necessary to provide accurate spatial representations of ozone in the region during high ozone events and are accurately representative of ozone values at their respective elevations. Additionally, continued monitoring at these permanent sites is essential for assessing long-term trends. However, beyond the study, when there are no temporary sites to monitor ozone at other locations and elevations, it becomes problematic that both permanent sites report on the lower end of what is observed across the region. When comparing MDA8 values at each of the six temporary sites from the 2018 study to the two permanent monitors for the same time period, we find that the MDA8 values are highly correlated (Table 2) with slopes ranging from 0.8 to 1.0, consistent with the conclusion of the GAMs analysis that the Pikes Peak region’s ozone is strongly impacted by regional background air. However, most sites display a detectable intercept when regressed against either of the two permanent sites. While the CC site reports lower MDA8 values than AFA or MAN (regression intercepts of −7 and –2 ppb, respectively), likely due to the impact of motor vehicle NO_x emissions given the site’s proximity to the interstate (Section 3.2), all of the other sites display positive intercepts (+3 to +15 ppb) when regressed against AFA or MAN (with the exception of CAS vs. AFA, where the intercept is near zero), further demonstrating that the two permanent sites generally report lower MDA8 ozone than most temporary sites in the 2018 ozone study.

One important implication for Colorado Springs, and potentially for other similar urban environments, is the need for additional long-term ozone monitoring for public health at higher elevation populated sites. As an example from this study, BLF is at an elevation of 2,233 m ASL, and the area just south of it has a high population density (Figure 11). Although Figure 11 relies on population data from the 2010 census, housing development has increased in areas such as the BLF as population increases in the Pikes Peak region as a whole (State Demography Office, 2019). Considering the relationships among ozone, elevation, and proximity to the highway corridor (Section 3.2), an additional higher elevation permanent site may provide currently missing data about the upper end of ozone levels in populated areas of the Pikes Peak region for the protection of human health. Although the FEMs designation deems the 2BTech temporary monitors to be of equal accuracy and precision to other FEM monitors including the permanent monitors in this study region, it should be noted that this study was run for only 4 months, and therefore, any changes made based on observations from the temporary monitor must take this abbreviated time frame into consideration.

We have explored the spatial and temporal behavior of summertime surface ozone in the Pikes Peak region of Colorado, United States. Our study is unique in that it focuses on an urban environment in the Southern Front Range of the Rocky Mountains, where ozone has been less extensively studied in comparison to the adjacent NFRMA. More broadly, our work contributes to the body of knowledge about ozone in the Western United States. We identify the effects of diurnal wind patterns, proximity to the highway corridor, and elevation on the mean diurnal variability of ozone across the Pikes Peak region. These factors can complicate the ability of lower elevation permanent sites to be representative of the broader region and may emphasize the need to place additional long-term monitors in higher elevation areas that are also relatively densely populated as well as further from high-traffic areas. Meanwhile, we find that MDA8 ozone in this location is strongly influenced by regional background air masses, particularly those containing wildfire smoke. Despite the interannual variability of fires, projections of increased fire frequency and severity in the Western United States due to climate change (Barbero et al., 2015; Westerling, 2016) highlight the need for continued consideration toward the impact that wildfires may have on air quality and surface ozone in the Pikes Peak region and other similar areas of the Western United States.

Data from Colorado Department of Public Health and Environment’s permanent monitoring stations can be accessed via the U.S. Environmental Protection Agency Air Quality System (AQS) database (https://aqs.epa.gov/api). Data from the summer 2018 special study are available in the Supplemental Material.

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

Figure S1. Time series of continuous hourly ozone observations and boxplot distributions of summertime hourly ozone at permanent Colorado Springs ozone monitoring stations. (A) Time series of continuous hourly ozone observations at AFA (left). Boxplot distributions of summertime (June–July–August) hourly ozone by year at AFA (right). (B) Time series of continuous hourly ozone observations at MAN (left). Boxplot distributions of summertime hourly ozone by year at MAN (right). Boxplots were generated using R Studio. The center line represents the median ozone mole fraction, the box represents the interquartile range, and upper/lower whiskers are either maximum/minimum value or the upper/lower quartile value plus/minus 1.5 times the interquartile range. Data after October 31, 2020, is considered preliminary and may be flagged or removed before being certified. AFA = Air Force Academy; MAN = Manitou Springs.

Figure S2. Distributions for hourly ozone frequency by site from June to September 2018. Histogram distributions for hourly ozone mole fractions at each site with permanent sites marked in red. All sites except MON and AFA have normal distributions. Table S2 displays descriptive statistics for the hourly ozone measurements for each site. AFA = Air Force Academy; MON = Monument.

Figure S3. Wind roses for ozone frequencies. Ozone frequencies by direction for the 0:00, 4:00, 8:00, 12:00, 16:00, and 20:00 h of the day at the MON (NNW of city center), BLF (NNE of city center), NAV (E of city center), and CAS (W of city center) sites during the summer 2018 ozone study (see Figure 2 for map). Ozone measurements are plotted only for each labeled hour (e.g., ozone measurements at hours 1:00–3:00, 5:00–7:00, 9:00–11:00, 17:00–19:00, and 21:00–23:00 and corresponding wind direction are not shown). MON = Monument; NAV = Navigator; BLF = Black Forest; CAS = Cascade.

Figure S4. NOAA HMS active fires (red triangles), smoke product (gray shading), and surface 8-h ozone observations on June 12, 2018, with 72-h HYSPLIT back-trajectory and ozone time series. (A) Map of fires burning and NOAA HMS smoke product for near peak ozone hour on June 12, 2018, and (B) the 72-h HYSPLIT back-trajectory for June 12, 2018, at peak ozone hour (18:00 UTC) for 500; 1,000; and 1,500 m AGL. (C) Time series plot of hourly ozone at eight sites in Colorado Springs during the high ozone event. NOAA = National Oceanic and Atmospheric Administration; HMS = Hazard Mapping System.

Figure S5. NOAA HMS active fires (red triangles), smoke product (gray shading), and surface 8-h ozone observations on July 6, 2018, with 72-h HYSPLIT back-trajectory and ozone time series. (A) Map of fires burning and NOAA HMS smoke product for near peak ozone hour on July 6, 2018, and (B) the 72-h HYSPLIT back-trajectory for July 6, 2018, at peak ozone hour (21:00 UTC) for 500; 1,000; and 1,500 m AGL. (C) Time series plot of hourly ozone at eight sites in Colorado Springs during the high ozone event. NOAA = National Oceanic and Atmospheric Administration; HMS = Hazard Mapping System.

Figure S6. NOAA HMS active fires (red triangles), smoke product (gray shading), and surface 8-h ozone observations on July 14, 2018, with 72-h HYSPLIT back-trajectory and ozone time series. (A) Map of fires burning and NOAA HMS smoke product for near peak ozone hour on July 14, 2018, and (B) the 72-h HYSPLIT back-trajectory for July 14, 2018, at peak ozone hour (21:00 UTC) for 500; 1,000; and 1,500 m AGL. (C) Time series plot of hourly ozone at eight sites in Colorado Springs during the high ozone event. NOAA = National Oceanic and Atmospheric Administration; HMS = Hazard Mapping System.

Figure S7. NOAA HMS active fires (red triangles), smoke product (gray shading), and surface 8-h ozone observations on July 18, 2018, with 72-h HYSPLIT back-trajectory and ozone time series. (A) Map of fires burning and NOAA HMS smoke product for near peak ozone hour on July 18, 2018, and (B) the 72-h HYSPLIT back-trajectory for July 18, 2018, at peak ozone hour (22:00 UTC) for 500; 1,000; and 1,500 m AGL. (C) Time series plot of hourly ozone at eight sites in Colorado Springs during the high ozone event. NOAA = National Oceanic and Atmospheric Administration; HMS = Hazard Mapping System.

Figure S8. NOAA HMS active fires (red triangles), smoke product (gray shading), and surface 8-h ozone observations on July 31, 2018, to August 2, 2018, with 72-h HYSPLIT back-trajectory and ozone time series. (A) Maps of fires burning and NOAA HMS smoke product for near peak ozone hour on July 31, 2018; August 1, 2018; and August 2, 2018, and (B) the 72-h HYSPLIT back-trajectories for July 31, 2018, to August 2, 2018, at peak ozone hour (18:00–22:00 UTC) for 500; 1,000; and 1,500 m AGL. (C) Time series plot of hourly ozone at eight sites in Colorado Springs during the high ozone event. NOAA = National Oceanic and Atmospheric Administration; HMS = Hazard Mapping System.

Table S1. A summary of results from (A) nightly baseline checks at temporary sites that used the 2BTech model 205 and (B) audits of temporary sites by a CDPHE independent QA team at three ozone levels during the 2018 Colorado Springs ozone study.^a,b,c

Table S2. Summary statistics for ozone mole fractions (ppb) during the June to September 2018 special ozone study.

Table S3. Variables used in the GAMs for both MAN and AFA sites.^a

Table S4. Training and cross-validation results for MDA8 values at MAN.

Table S5. GAM count and residuals (ppb) for nonsmoke and smoke days for MAN.^a

Table S6. Training and cross-validation results for MDA8 values at AFA.

Table S7. GAM count and residuals (ppb) for nonsmoke and smoke days for AFA.^a

Table S8. High ozone event dates, durations, and peak ozone during the 2018 ozone study.

Data S1. Hourly ozone data collected for the 2018 CDPHE special study of the Pikes Peak region.

The authors thank Phillip Stauffer of Colorado Department of Public Health and Environment for contributions to the deployment of instrumentation, data collection, quality assurance and quality control, and initial data assembly and assessment for the summer 2018 special study. The authors also thank Matt Cooney and Francis Russell of the Colorado College GIS Laboratory for their assistance with geospatial analysis.

Funding was provided by the Colorado Department of Public Health and Environment (E. J. Mattson) and by a Colorado College Student-Faculty Collaborative Research (SCoRe) grant (M. T. Flynn and L. E. Gratz).

The authors have declared that no competing interests exist.

Contributed to the conception and design: EJM.

Contributed to the acquisition of data: EJM.

Contributed to the analysis and interpretation of data: MTF, LEG, DAJ, EJM.

Drafted and/or revised the article: MTF, LEG, DAJ, EJM.

Approved the submitted version for publication: MTF, LEG, DAJ, EJM.

How to cite this article: Flynn, MT, Mattson, EJ, Jaffe, DA, Gratz, LE. 2021. Spatial patterns in summertime surface ozone in the Southern Front Range of the U.S. Rocky Mountains. Elementa: Science of Anthropocene 9(1). DOI: https://doi.org/10.1525/elementa.2020.00104

Domain Editor-in-Chief: Detlev Helmig, Boulder AIR LLC, Boulder, CO, USA

Associate Editor: Samuel J. Oltmans, NOAA/ESRL Global Monitoring Division, Boulder, CO, USA

Knowledge Domain: Atmospheric Science