Source characterization of volatile organic compounds at Carlsbad Caverns National Park

ABSTRACT Carlsbad Caverns National Park (CAVE), located in southeastern New Mexico, experiences elevated ground-level ozone (O3) exceeding the National Ambient Air Quality Standard (NAAQS) of 70 ppbv. It is situated adjacent to the Permian Basin, one of the largest oil and gas (O&G) producing regions in the US. In 2019, the Carlsbad Caverns Air Quality Study (CarCavAQS) was conducted to examine impacts of different sources on ozone precursors, including nitrogen oxides (NOx) and volatile organic compounds (VOCs). Here, we use positive matrix factorization (PMF) analysis of speciated VOCs to characterize VOC sources at CAVE during the study. Seven factors were identified. Three factors composed largely of alkanes and aromatics with different lifetimes were attributed to O&G development and production activities. VOCs in these factors were typical of those emitted by O&G operations. Associated residence time analyses (RTA) indicated their contributions increased in the park during periods of transport from the Permian Basin. These O&G factors were the largest contributor to VOC reactivity with hydroxyl radicals (62%). Two PMF factors were rich in photochemically generated secondary VOCs; one factor contained species with shorter atmospheric lifetimes and one with species with longer lifetimes. RTA of the secondary factors suggested impacts of O&G emissions from regions farther upwind, such as Eagle Ford Shale and Barnett Shale formations. The last two factors were attributed to alkenes likely emitted from vehicles or other combustion sources in the Permian Basin and regional background VOCs, respectively. Implications: Carlsbad Caverns National Park experiences ground-level ozone exceeding the National Ambient Air Quality Standard. Volatile organic compounds are critical precursors to ozone formation. Measurements in the Park identify oil and gas production and development activities as the major contributors to volatile organic compounds. Emissions from the adjacent Permian Basin contributed to increases in primary species that enhanced local ozone formation. Observations of photochemically generated compounds indicate that ozone was also transported from shale formations and basins farther upwind. Therefore, emission reductions of volatile organic compounds from oil and gas activities are important for mitigating elevated O3 in the region.


Introduction
Over the last five years, the frequency of elevated groundlevel ozone (O 3 ) mixing ratios has increased at Carlsbad Cavern National Park (CAVE; Figure 1).High levels of O 3 have detrimental impacts on both human health and vegetation, which is of significant concern in US National Parks (Keiser, Lade, and Rudik 2018;Kohut 2007).Since 2018, the maximum daily 8-hour average (MDA8) O 3 at CAVE (Figure 1b) has regularly exceeded the current National Ambient Air Quality Standard (NAAQS) of 70 ppbv (US NPS 2022).For example, in 2022, there were 21 days with MDA8 O 3 � 70 ppbv.To address the emerging air quality issues in southeast New Mexico, a series of special studies was conducted to investigate the sources of air pollutants and their precursors, with an emphasis on O 3 formation.
The formation of tropospheric O 3 occurs when volatile organic compounds (VOCs) are oxidized in the presence of NO x (NO + NO 2 ) and sunlight (Chameides et al. 1992).NO x emissions in the US are dominated by fuel combustion (US EPA 2021).In the western US, however, long-range transport of reservoir species such as peroxyacetyl nitrate (PAN) is also a significant contributor to NO x (Fischer, Jaffe, and Weatherhead 2011;Lin et al. 2017).VOCs can be emitted from both natural and anthropogenic sources (e.g., (Gu, Guenther, and Faiola 2021)).Natural VOC sources include vegetation (e.g., (Geron et al. 2006;Rasmussen and Went 1965;Zimmerman, Greenberg, and Westberg 1988)), oceans (Fischer et al. 2012), biomass burning (Andreae 2019;Jin et al. 2022), and soil and plant litter (Leff and Fierer 2008).Anthropogenic VOC sources in the US include transportation, oil and gas (O&G) development and production activities, other industrial activities, residential wood burning, and volatile chemical products (Cooper et al. 2012;McDonald et al. 2018;McDuffie et al. 2016;US EPA 2021).Contributions of these emissions to VOC mixing ratios and their capabilities to initiate O 3 production vary significantly across regions (Baker et al. 2008;Benedict et al. 2020).Identifying the relative importance of different sources is critical for designing effective strategies to mitigate elevated O 3 levels at CAVE.
Located in southeastern New Mexico, CAVE is within the Permian Basin and adjacent to Delaware shale play within the Permian Basin.The closest O&G well pad is 6 km to the east.CAVE is sometimes also downwind of the Eagle Ford and Barnett Shale formations located several hundred km to the southeast and east, respectively.The Permian Basin is the most productive O&G regions in the US, and O&G production has increased rapidly from 2013 to 2022 (Figure 1c) (US EIA 2023).Satellite observations show that NO 2 levels are increasing in the region (Dix et al. 2020).Benedict et al. (2020) measured VOC concentrations at four national parks in the southwestern US and found strong influence of anthropogenic sources on VOCs at CAVE and in the surrounding region.This screening study lacked sufficient temporal resolution to capture the full variability of VOCs at CAVE for source characterization and motivated a follow-up study with VOC measurements and NO x observations at a higher temporal resolution (i.e., hourly).
Here, we present the source characterization of VOCs at CAVE during the Carlsbad Caverns Air Quality Study (CarCavAQS) conducted between July 25 th and September 5 th , 2019.The campaign was designed to investigate the impacts of regional sources on fine particle formation (PM 2.5 ), VOC composition, O 3 and PAN formation, and nitrogen deposition.The composition, sources, and visibility impacts of PM 2.5 at CAVE have been discussed in Naimie et al. (2022), while O 3 , NO x , and PAN observations are discussed in Pollack et al. (2023).In this work, we characterize sources of VOCs at CAVE using positive matrix factorization (PMF) analysis and consider the contribution of different sources to VOC reactivity.

Monitoring site and VOC sampling
During CarCavAQS, measurements were made at the CAVE Biology Building (32.18°N, 104.44°W; 1355 m above sea level), located within the park and approximately 0.5 km from the visitor center.The major urban center of El Paso-Ciudad Juarez is located 200 km to the west.Smaller cities are located 140 km to the east (Odessa and Midland), 140 km to the north (Roswell) and 30 km to the east (Carlsbad).
Meteorological and O 3 data were collected at the same location by the National Park Service.Mean ambient temperature and relative humidity were 33°C and 36%, respectively, during the campaign.From July to September, the summer monsoon generally transports air masses out from the Permian Basin to the northwest (Figure S1).Therefore, prevailing winds at the site were southeasterly, with occasional gusts from the northwest (Figure S2).When the summer monsoon transport is weak, the site can be more strongly impacted by local and terrain-driven flows, with southeasterly daytime upslope flows and nocturnal northeasterly downslope flows (Crosman 2021).During the summer of 2019, the site was impacted by a subtropical ridge over southern New Mexico, which suppressed the monsoon transport out of the Permian Basin to the west (Figure S1) and led to more frequent northwesterly winds at CAVE (Figure S2).To qualitatively distinguish the origins of air masses, observations are grouped into south-tosoutheast (S-SE) and north-to-northwest (N-NW) wind sectors (Figure S1).The S-SE wind sector, the direction of the Permian Basin, ranges from 90° to 225°, with 0° representing north and 180° representing south.Wind directions larger than 270° or less than 45° are within the N-NW wind sector, and air masses from this sector typically represent background conditions.
To investigate VOC composition near and away from sources, 66 air samples were also collected at various locations using grab canister samples.Stainless steel Silonite®-coated 1.4 L canister samples (Entech Instruments) were used.A comprehensive description of canister cleaning and sampling procedures is available in Benedict et al. (2020).Sampling locations were selected to reflect background conditions (blue squares in Figure 1) and air masses impacted by urban/traffic emissions (blue crosses), O&G development activities (blue triangles), and O&G production activities (blue circles).Grab canister samples collected within 10 km of an O&G site being developed during CarCavAQS, were considered influenced by O&G development (n = 10).Samples were deemed impacted by O&G activities if they were collected within 10 km of an O&G producing well pad and were not in urban areas or heavily impacted by traffic (n = 25).Note that for some locations it is challenging to distinguish impacts from urban and O&G emissions.

Measurements
Measurements discussed in this work include O 3 , methane (CH 4 ), NO x , total reactive nitrogen (NO y ), PAN, and non-methane VOCs (NMVOCs).Detailed descriptions of each measurement can be found in Text S1 of the supporting information (SI).Briefly, O 3 was measured by a Thermo 49i O 3 analyzer.CH 4 was measured via cavity ringdown absorption spectroscopy (Picarro, G2508).NO, NO 2 , and NO y were measured using a single-channel commercial NO analyzer (Eco Physics, model CLD 780 TR) employing NO-O 3 chemiluminescence detection and outfitted with two custombuilt inlets (Pollack, Lerner, and Ryerson 2010).PAN was measured using a custom-built gas chromatograph (GC) instrument using a commercial electron capture detector.
A five-channel, three-GC analytical system was used to analyze 56 individual VOCs in real-time (hourly time resolution), including light alkanes (C 2 -C 7 ), heavy alkanes (C 8 -C 10 ), C 1 -C 2 halocarbons, C 1 -C 5 alkyl nitrates, alkenes, branched alkanes, branched alkenes, cyclic alkanes, aromatic VOCs, and biogenic VOCs (Abeleira et al. 2017).A similar lab-based system was also used to analyze VOCs in the grab canister samples (Benedict et al. 2020).In addition to the GC system, a proton transfer reaction-mass spectrometer (PTR-MS; PTR-MS HS, Ionicon Analytik) was deployed to measure acetonitrile, isoprene, and select oxygenated VOCs (OVOCs), including acetone, methyl ethyl ketone (MEK), and methyl vinyl ketone (MVK).A complete list of VOCs measured and corresponding measurement precisions and method detection limits (MDLs) are provided in Table S1.
The inlets of these instruments were positioned approximately 6-7 m above the ground.All measurements were averaged to one hour to synchronize the data, although instrument sampling frequencies ranged from one minute to one hour.

Positive matrix factorization
PMF has been widely used as a source apportionment technique for VOCs (Abeleira et al. 2017;Bon et al. 2011;Brown, Frankel, and Hafner 2007;Pollack et al. 2021;Ulbrich et al. 2009;Yuan et al. 2012).PMF is a receptor model that decomposes a matrix of speciated sample data (X) into factor profiles (F) and their contributions (G): The sample data matrix (X) has a dimension of m � n, where m is the number of species and n is the number of samples.The factor profile matrix (F) has p � m elements, which are often interpreted as fingerprints of VOC sources or the stages of photochemical processing (Yuan et al. 2012) and are assumed to be constant over the duration of the campaign.Factor contributions vary and reflect the relative importance of different sources.The factor contribution matrix (G) has n � p elements, representing the contributions of the factors.The residual matrix (E) represents data points not fully fit by the model.
The EPA PMF Model (v5.0) is used in this study, which solves Equation (1) using the constraint that no sample can have significantly negative source contributions (Paatero 1997(Paatero , 1999)).The tool uses both sample concentrations and measurement uncertainties to weight individual data points when calculating and minimizing the objective function Q.Samples and uncertainties used in our PMF analysis were prepared following Hopke (2016) and Polissar et al. (1998).VOC measurements with mixing ratios above the corresponding MDL were used directly as inputs and assigned uncertainties listed in Table S1.However, some VOC species had unclear or indiscernible chromatogram peaks for certain periods and were deemed below their MDLs.In this case, the corresponding MDLs divided by three were used as inputs, and 5/6 times the MDLs were assigned as uncertainties.Weighting data points this way allows the level of confidence in the measurement to be accounted for in the PMF model solution.In addition, signal-to-noise ratios (S/N) were considered before running the PMF model.Following Hopke (2016), 11 VOC species with S/N < 0.5 were flagged as "bad" and excluded from the PMF analysis.In total, 47 species and 573 samples from the real-time GC and PTR-MS data were included in the PMF analysis.
Determining the number of factors used in the PMF model is important so that the solution can represent source characterization appropriately.A sufficient number of factors is required to explain the variability in the dataset.However, when too many factors are included, physically relevant factors can be separated into multiple unrealistic factors, known as "factor splitting" (Brown, Frankel, and Hafner 2007;Ulbrich et al. 2009).In this work, PMF solutions with numbers of factor ranging from 3 to 9 were evaluated, and the solution with seven factors was selected because it significantly reduces the model error without unreasonable splitting of factors.Details for selecting the sevenfactor PMF solution are provided in Text S2 and Figure S3-S7.
We estimated the uncertainties of factor profiles and contributions from 100 bootstrapping runs (Ulbrich et al. 2009).Figures and tables throughout this work show factor profiles and factor contributions from the PMF solution, and uncertainties represent the mean and the 25 th and 75 th percentiles from the 100 bootstraps.We note that the mean and median values determined from bootstrapping are nearly identical in the PMF solution (Figure 2).We also found that no additional matrix rotation or constraints were needed to optimize the solution (Figure S8) (Paatero and Hopke 2009).

HYSPLIT backward trajectories and residence time analysis
Backward trajectories were calculated using the Hybrid Single Particle Lagrangian Integrated Trajectory (HYSPLIT) model (Stein et al. 2015).A detailed description of the model setup is provided in Naimie et al. (2022).HYSPLIT trajectories and source-receptor relationships were investigated using residence time analysis (RTA) (Gebhart et al. 2018(Gebhart et al. , 2021)).The residence time is the amount of time an air mass spends in a horizontal grid cell (0.25° × 0.25°), which is represented by the everyday probability (EP).Highconcentration residence time or high probabilities (HP) are calculated using trajectories from the top 10% concentration values at the receptor site.These two residence times can be compared using the incremental probability (IP = HP -EP) to determine if the probability of an endpoint being in a grid cell is more likely during the high concentration periods or throughout the dataset.

VOC-OH reactivity
The potential of a VOC species to initiate the O 3 formation cycle is partly driven by its reactivity with the OH-radical, R OH;VOC (s −1 ) R OH;X is calculated as the inverse of the lifetime of OH reaction with species X as: where k OH;X is the rate constant at 298 K (obtained from Atkinson and Arey (2003) and Atkinson (1986)) and X ½ � is the number density of X in molecule cm −3 .Although this simplified estimate does not consider reactions with other oxidants, nor does it consider chain propagation or termination steps resulting from reaction with OH, R OH;VOC has been widely used to provide valuable insights into the potential for a VOC to contribute to O 3 production, as the initial reaction with OH is often the limiting step in the sequence of reactions to form O 3 (Abeleira et al. 2017;Benedict et al. 2020;Gilman et al. 2013;McDuffie et al. 2016;Rutter et al. 2015).

Photochemical aging
To understand the sources of secondary VOCs and their relationships with O 3 formation, we also estimated photochemical air mass ages using ratios of concentrations of two alkyl nitrates to their corresponding parent hydrocarbons (RONO 2 /RH).Alkyl nitrates are formed when n-alkanes are oxidized (Simpson et al. 2003).Therefore, the ratio of a daughter alkyl nitrate to its parent n-alkane increases as an air mass ages (Benedict et al. 2020;Bertman et al. 1995;Russo et al. 2010;Simpson et al. 2003).We use 2-pentyl nitrate/n-pentane and 2-butyl nitrate/n-butane ratios to construct a photochemical clock based on measured and modeled ratios.The modeled ratios are calculated using laboratory kinetic data (Bertman et al. 1995), assuming OH ½ � ¼ 10 6 molecules cm −3 and that OH-initiated photochemical production is the only source of the given alkyl nitrates.Since the estimated photochemical ages using this method scale with [OH], they are only approximate, allowing for comparing relative photochemical processing among different air masses.In addition, deviations from the modeled relationship are expected because additional sources and sinks are not considered.More details of the model are provided in Text S3.  different wind directions are shown in Figure 3.The contribution of OVOCs to VOC OH reactivity might be underestimated since only a few OVOCs were measured (i.e., acetone, MEK, and MVK). Figure 3 also shows the mean VOC concentrations and R OH;VOC in canister samples collected in different regions.We note that OVOCs and acetonitrile were not analyzed in the canister samples.

VOC observations and VOC OH reactivities
Observations made at CAVE are categorized into S-SE and N-NW wind sectors as defined in Section 2.1 and accounted for 59% and 33% of the time with valid observations.High total VOC episodes (>100 ppbv) occurred mostly with wind from the S-SE wind sector except for one event on 08/19 during which the locally observed wind direction shifted between the two wind sectors.Short-term wind direction shifts occurred mostly as a part of the wind diel pattern.VOC mixing ratios remained low when CAVE experienced consistent N-NW winds (e.g., 08/25-08/27).High O 3 (purpleshaded areas in Figure 2) only occurred when winds were from the S-SE wind sector (pink-shaded areas in Figure 2) but did not always align with high VOC periods.
To better illustrate VOC composition, VOCs are grouped into nine categories: light (C 2 -C 7 ) alkanes, heavy (C 8 -C 10 ) alkanes, aromatics, OVOCs, alkyl nitrates, alkenes, ethyne, biogenics, and halogenated compounds.Tables S1 and S2 list compounds in each group and summary statics for concentrations at CAVE and in grab canister samples.Consistent with Benedict et al. (2020), light alkanes dominated VOCs observed at CAVE and in the grab canister samples collected in the area.Although certain light alkanes can be emitted from combustion sources, light alkanes are generally considered good markers for O&G sources (Abeleira et al. 2017;Benedict et al. 2020;Gilman et al. 2013;Swarthout et al. 2013).VOC mixing ratios during periods with N-NW winds were slightly higher than those in the canister samples from background and urban regions.The background and urban canister samples also had higher biogenic but lower light alkanes.The difference in VOC mixing ratios and compositions indicate that locally measured winds may not fully represent regional transport patterns, especially when the diel wind pattern was dominated by terrain-driven flows.Thus, further analyses of canister samples and backward trajectories are needed to better link VOCs observed at CAVE to sources.
The i-to n-pentane ratio can be used to distinguish emissions between O&G or combustion emissions.An i-to n-pentane ratio much larger than one implies strong influence from combustion sources, including vehicular emissions (Russo et al. 2010b), while a ratio less than one indicates substantial influence from O&G emissions (Evanoski-Cole et al. 2017;Gilman et al. 2013;Swarthout et al. 2013).Measured i-to n-pentane ratios at CAVE ranged between 0.8 and 1.5, and values below one were observed mostly during periods with S-SE winds (Figures 1 and 3).The pattern of i-to n-pentane ratios observed at CAVE is consistent with patterns observed in grab canister samples collected in background areas and urban regions showed higher i-to n-pentane ratios with larger variability compared to those collected in O&G development and production regions.This consistency suggests that wind directions can be used to distinguish air masses impacted by O&G activities and other sources.At CAVE, light alkanes also dominated VOC OH reactivities.Air masses from the S-SE wind sector passing through the Permian Basin had higher total OH reactivity for measured VOCs, indicating strong influence of O&G emissions.Although reactive OVOCs and acetonitrile were not measured in the grab canister samples, grab canister samples collected in O&G development areas showed the highest total OH reactivity, which corroborates the importance of O&G emissions for O 3 formation.Concentrations of biogenic compounds were low, but they were still the second most important contributor to OH reactivity at CAVE because they are highly reactive with OH.In the canister samples collected in background and urban areas, biogenics were the main contributor to total OH reactivity despite being a desert environment.Alkenes were typically the third contributor to OH reactivity at CAVE, followed by OVOCs and heavy alkanes.
The wind direction, tracer ratio, and grab canister sample analyses identify the differences in VOC composition and VOC OH reactivities of different sources.However, they do not explain why O 3 increases were associated with both high total VOC (08/28) and substantially lower total VOC (08/05-09, and 08/11).A PMF analysis could help illustrate the differences between these cases.

PMF factors
Using the NMVOC observations from the on-line fivechannel GC system and the PTR-MS at CAVE, the PMF model identified seven factors, which were labeled as follows: (1) longer-lived O&G, (2) shorter-lived O&G, (3) process-specific O&G, (4) longer-lived secondary, (5) shorter-lived secondary, (6) alkenes, and (7) background.Figure 4 shows factor contributions to each VOC species considered in the PMF analysis, and Figure 5 shows their diel and wind direction profiles.NO x , O 3 , and PAN were left out of the PMF analysis to be used as tracers to facilitate the characterization of the PMF output, and their diel profiles are also presented in Figure 5.Time series of the contributions of the PMF factors are shown in Figure S10, and reconstructed time series of selected compounds from different groups using the PMF factors are shown in Figure S11.
Factors 1, 2 and 3 are consistent with emissions from O&G activities.The longer-lived O&G factor (factor 1) accounts for 20-50% of the observed C 2 -C 7 alkanes and is consistent with contributions of O&G factors to these alkanes reported in previous studies in other O&G-producing regions (Abeleira et al. 2017;Gilman et al. 2013;Pollack et al. 2021;Swarthout et al. 2015).The C 7 alkanes, especially branched ones, contribute significantly to the shorterlived O&G factor (factor 2), along with C 8 alkanes.Emissions of these alkanes have been closely linked to specific activities associated with O&G extractions, such as flashing from oil and condensation tanks (Warneke et al. 2014).The process-specific O&G factor (factor 3) is dominated by C 9 -C 10 alkanes, ethylbenzene, and xylenes.High emissions of these species have been observed during O&G development processes, such as drilling, fracking, and flowback (Hecobian et al. 2019).In studies in the Denver-Julesburg (DJ) Basin, elevated emissions of C 8 -C 10 alkanes have been reported specifically during welldrilling operations (Lachenmayer 2022).All three factors are associated with S-SE winds from the Permian Basin.Thus, the diel patterns of contributions of these factors are similar.Their contributions decreased in the afternoon, indicative of dilution resulting from boundary layer deepening.Contributions of these three factors correlate well with each other (R 2 ranging from 0.37 to 0.76; Figure S9) and with CH 4 (Figure S12).The correlation between the processspecific O&G factor and CH 4 is lower than for the other two O&G factors; this is not surprising given that CH 4 emissions can vary significantly relative to VOC emissions across well development activities.
Secondary VOC species (e.g., OVOCs and alkyl nitrates) produced by photochemical reactions map into two factors.Both factors have an i-to n-pentane ratio lower than one, suggesting their O&G origins.A longer-lived secondary factor (factor 4) is dominated by OVOCs and C 1 -C 4 alkyl nitrates.Yuan et al. (2013) showed that longer-lived species like C 3 -C 5 alkanes can be a principal source of acetone in some environments.Thus, we postulate that the C 1 -C 4 alkyl nitrates and OVOCs associated with the longer-lived secondary factor arise from similar photochemical reactions of light alkanes associated with the longer-lived O&G factor.This hypothesis is supported by the fact that the longerlived secondary factor always consists of both a fraction of light alkanes and C 1 -C 4 alkyl nitrates regardless of the number of factors (Figure 4 and S4-S7).The shorterlived secondary factor is enriched in larger (e.g., C 4 and C 5 ) alkyl nitrates, which have shorter lifetimes than smaller alkyl nitrates (Bertman et al. 1995).Diel patterns of the two secondary factors have a pronounced daytime rise.However, the shorter-lived secondary factor has a similar diel pattern to those of PAN and solar radiation, whereas the longer-lived secondary factor peaks later in the afternoon, suggesting they were likely aged air transported from farther upwind and entrained from above the mixed layer.The shorter-lived secondary factor is largely associated with S-SE winds, while the longer-lived secondary factor is largely correlated with easterly winds.
The sixth factor is dominated by light alkenes and the seventh factor is mostly comprised of isoprene, halogenated compounds, and acetonitrile.Unlike the O&G and secondary factors, the alkene factor has a mean i-to n-pentane ratio larger than one (1.2),suggesting stronger influence on this factor from combustion sources.Light alkenes have been linked to vehicular emissions in urban areas and O&G production regions (Baker et al. 2008;Swarthout et al. 2013) but also can be emitted by industrial activities (Washenfelder et al. 2010), vegetation (Rhew et al. 2017), and biomass burning (Akagi et al. 2012;Hecobian et al. 2011).CAVE was not significantly impacted by wildfire plumes in 2019 during CarCavAQS (Pollack et al. 2023).Since CAVE is more than 140 km away from major urban centers and high contributions of the alkene factor occurred mostly with S-SE winds, we hypothesize that the alkene factor may reflect influence from vehicular or other combustion emissions associated with O&G activities in the Permian Basin.However, vehicular emissions from the Facilities area near the site could also contribute to the alkene factor.The seventh factor is comprised of species with long lifetime or natural sources and shows little correlation with other trace gases (Figure S9) and weak wind direction dependence.This factor is interpreted as representing the regional background.
On average, the long-lived O&G factor had the largest contribution to observed VOCs and their OH reactivities at CAVE during CarCavAQS (Figure 6),  followed by the process-specific O&G factor.During periods with S-SE winds, contributions of the three O&G factors increased.The two secondary factors also contributed significantly to the measured VOCs, especially during periods when O 3 was high, whereas contributions of the three O&G factors decreased during high O 3 events compared to other periods with S-SE winds, especially the process-specific O&G factor 3 (Pollack et al. 2023).

Linking PMF factors to sources
Although the wind direction analyses shown above roughly point to the source regions of different VOC species and PMF factors, in-situ wind measurements near the surface at CAVE may be impacted by local terrain effects and not correctly represent regional transport patterns.In addition, confirmation is needed for the attribution of the PMF factors based mostly on previously reported tracers observed near different sources, especially for the process-specific O&G factor.Emission profiles for these processes in the Permian Basin might be different from those reported for the DJ Basin.Therefore, we further applied the PMF results to grab canister observations and conducted RTAs to better link VOC observations and PMF factors to different source regions.
To check if the PMF analysis correctly captured the characteristics of different sources in the regions, we solved Equation (1) for all grab canister samples using a non-negative least squares solver with the factor matrix from our PMF analysis based on CAVE observations (Lawson and Hanson 1987).More details are provided in Text S4.Note that acetonitrile and OVOCs, including acetone, MEK, and MVK, were not measured in the canister samples and were removed from the fitting.The estimated contributions are normalized using all canister means to show the relative importance of a factor from different sources, and mean normalized contributions are shown in Figure 7.The highest contributions of the shorterlived O&G, process-specific O&G, and shorter-lived secondary factors were observed near active drilling sites.Samples from O&G production regions showed the highest contribution of the longer-lived O&G factor.The alkene factors were observed in urban/traffic, O&G development, and O&G production canister samples.The background factor was relatively consistent near different sources.Thus, PMF contributions derived from source-dominated air samples are generally consistent with our attribution, providing additional support to the PMF analysis.
VOCs sampled in grab canister samples collected near O&G development and production sites show large variability and relatively large residuals when fitted with the PMF factors (see Table S3 and Figure S13).The samples strongly influenced by local sources actually have larger fitting errors (Figure S13), indicating that the PMF factors derived from CAVE observations represent aged air masses.The large variability likely arises from the challenges associated with obtaining representative samples for these sources in close proximity because the wind direction and atmospheric conditions could be non-ideal while emissions and plume interception can both be intermittent.Therefore, instead of using the VOC observations from air samples to conduct source attribution with the PMF method or a chemical mass balance approach directly, we only used the canister results as a check on factor assignments in our PMF analysis.In addition to applying the PMF results to the grab canister observations, we also conducted RTA to investigate the transport pathways associated with different factors.Figure 8 shows the IP distributions using 48-hour backward trajectories for the seven PMF factors.Air masses with a high contribution (top 10%) of a factor were more likely to have passed through the red grid cells than the blue cells.High contribution periods of the three O&G factors (Figure 8a-c) and the alkene factor (Figure 8f) had a similar transport pattern, predominantly passing through O&G development and production areas in the Permian Basin and Eagle Ford Shale formation to the southeast of CAVE.The high contribution of the shorter-lived secondary factor was associated with transport from the southeast and the east sectors (Figure 8e).However, trajectories through central, north, and northwest Texas, Oklahoma, and Kansas contributed more to periods with high contribution of the longer-lived secondary factor (Figure 8d).These air masses might have been impacted by multiple sources at different times, including urban and industrial emissions as well as O&G emissions in the Permian Basin, the Eagle Ford Shale, the Barnett Basin Shale, or the Anadarko Basin in Oklahoma.High contributions of the background factor were associated with transport from the west (Figure 8g).

Photochemical aging
Figure 9 shows 2-pentyl nitrate-to n-pentane ratios versus 2-butyl nitrate-to n-butane ratios derived from the observations at CAVE and grab canister samples as well as the relationship between the two ratios modeled based on laboratory parameters (Bertman et al. 1995).The relationship between the two ratios can be used as a photochemical clock to infer the photochemical aging of air masses (Benedict et al. 2020).The modeled photochemistry matched the values derived from canister samples collected near O&G development and production sources.Canister samples with a shorter photochemical age often had a higher ethane (C 2 H 6 ) mixing ratio, a good indicator for O&G emissions (Benedict et al. 2020, Evanoski-Cole et al. 2017;Swarthout et al. 2013Swarthout et al. , 2015)).The good agreement between the modeled and observed ratios provides support for using the photochemical clock for air masses strongly influenced by O&G emissions, whose alkyl nitrates are mostly generated through OH-initiated reactions.For the canister samples representing background and urban conditions, there were larger discrepancies between the modeled and observed ratios, suggesting factors not included in the model played a more important role, such as other possible sources and removal mechanisms for 2-pentyl and 2-butyl nitrates, and varying branching ratios (Benedict et al. 2020).
The photochemical clock suggests that air masses observed at CAVE might have undergone two different modes of photochemical processing.Fresh air masses generally had a high C 2 H 6 mixing ratio, consistent with results for grab canister samples, and low PAN and O 3 mixing ratios.For aged air masses (>6 hours), there were two clusters.One cluster matched the calculated ratios, and PAN and O 3 mixing ratios increased as they went through photochemical aging.Air masses in the other cluster deviated from the calculated ratios.Despite having similar O 3 mixing ratios, their C 2 H 6 and PAN levels were much lower than those of air masses in the first cluster, suggesting longer photochemical ages than the ones predicted by the photochemical clock.The discrepancies could result from varying OH concentrations and/or photolysis rates during the history of the air mass (Benedict et al. 2020).

Conclusion and implications
During CarCavAQS, C 2 -C 7 alkanes dominated both observed VOC mixing ratios (98% on average) and VOC OH reactivities (79% on average), consistent with the findings from Benedict et al. (2020).These alkanes were generally associated with air masses from the S-SE and with i-to n-pentane ratios < 1.These air masses had a similar VOC composition in the grab canister samples collected near O&G development and production sites.Three of the PMF factors were attributed to O&G activities, comprising VOC profiles previously reported for O&G production and development activities.Contributions of these factors were significantly higher during periods with S-SE winds (69% to VOC mixing ratios) than periods with N-NW winds (45% to VOC mixing ratios).Our RTA further confirmed that the O&G factors were associated mainly with transport from the south and southeast, an area of significant O&G development and production activities.
Our PMF analysis also identified two secondary factors with different atmospheric lifetimes, suggesting that air masses observed at CAVE might have different photochemical ages.Both factors are likely oxidation products of O&G emissions as indicated by their i-to n-pentane ratios less than one.The contribution of the shorter-lived secondary factor was higher from the S-SE based on our wind direction analysis and the RTA.Its contribution typically peaked around noon (10:00-14:00) with the highest solar radiation.Therefore, the air masses with an increased contribution of shorter-lived secondary factor might have been transported during the day directly from source regions near CAVE (e.g., Permian Basin) with freshly produced secondary species.Similar diel patterns have been observed in Northern Colorado in summer 2015 (Abeleira et al. 2018).The diel pattern of the longer-lived factor peaked later in the afternoon (14:00-18:00), which could be related to boundary layer development and entrainment.O 3 increases were observed during periods with high contributions of both shorter-lived and longer-lived factors, and the O 3 diel pattern had a broad peak (10:00-18:00).The RTA showed that air masses with a high contribution of longer-lived secondary factor might have been transported from source regions further upwind from the southeast, the east, and the northeast.The longer-lived secondary factor also consisted of less reactive light alkanes but lacked more reactive heavy alkanes, suggesting the longer-lived secondary species were photochemical products of VOC and NO x emissions from O&G activities.
Our photochemical clock model shows that O 3 enhancements were related to both locally and regionally formed O 3 .A fraction of the air masses with both high O 3 and PAN mixing ratios had a photochemical age between 6 to 12 hours, showing the impacts of locally formed O 3 .However, a fraction of high O 3 plumes had a photochemical age longer than the daytime duration (12 hours).For these plumes, the transport time might be significantly longer than the photochemical age since the [OH] of 10 6 molecules cm −3 used in the photochemical clock model represents daytime conditions with high photolysis rates and nighttime [OH] is usually one or two orders of magnitude lower (Bey, Aumont, and Toupance 1997).Therefore, these plumes may reflect air masses stagnated near CAVE overnight or sources located more than one day upwind of CAVE.
The source apportionment shown here establishes the dominant influence of O&G activities on both primary and secondary VOC at CAVE, which likely affects local and regional O 3 formation.A more comprehensive modeling effort is needed to investigate local and regional photochemistry and their impacts on elevated O 3 at CAVE and in the region, including Carlsbad, NM which has also been experiencing O 3 exceedance.

Figure 1 .
Figure 1.(A) map of the study region, (B) O 3 mixing ratios measured in Carlsbad Cavern National Park (CAVE), and (C) O&G production in the Permian Basin.Grey, pink, and green areas in panel (A) are urban regions, O&G development and production regions (US EIA 2023), and national park units, respectively.Brown line shows the extent of the Permian Basin.Brown triangles show active O&G drilling sites during the 2019 CarCavAQS campaign.Blue squares, crosses, triangles, and circles represent air sampling locations using grab canister samples that reflect air masses in the background, urban, O&G development, and O&G production regions, respectively.See Section 2.1.For how the grab canister samples were categorized.Panel (B) shows O 3 mixing ratios observed at CAVE, and the gray area in panel (B) shows the CarCavAQS period.Red numbers in panel (B) are the number of days with a MDA8 O 3 >70 ppbv.Blue and orange lines in panel (C) show O&G production in the Permian Basin (US EIA 2023).

Figure 2
Figure 2 shows meteorological conditions, O 3 mixing ratios, and NMVOC mixing ratios at CAVE during CarCavAQS.C 2 -C 7 alkanes (orange areas in Figure 2) dominated NMVOCs during the campaign, and time series of other species are provided in the bottom panel of Figure 2. Mean VOC concentrations and R OH;VOC at CAVE for the whole campaign and periods with

Figure 2 .
Figure 2. Meteorological conditions, O 3 mixing ratios, and NMVOC mixing ratios at CAVE during the CarCavAQS.The top panel shows wind speed (black) and wind direction (wind dir; blue) observations.The filled areas and the purple line in the bottom panel show VOC and O 3 mixing ratios measured at CAVE.The bottom panel shows the VOC mixing ratios for species other than the C 3 -C 7 alkanes.The pink-and purple-shaded areas show periods with southerly and southeasterly winds and high O 3 (70 ppbv), respectively.

Figure 3 .
Figure 3.(A) Mean VOC mixing ratios, i-to n-pentane ratios, and (B) mean OH reactivities observed at CAVE under different winds and in the grab canister samples.The boxes and whiskers in panel (A) show the 5 th , 25 th , 50 th , 75 th , and 95 th percentiles of the i-to n-pentane ratios.The triangles in panel (A) show the mean i-to n-pentane ratios.Note that acetonitrile and OVOCs, including acetone, MEK, and MVK, were not measured in the canister samples.

Figure 4 .
Figure 4. PMF factor profiles.The bars represent the contribution of the factor to each species, and error bars represent the 25 th and 75 th percentiles of the contributions from the 100 bootstrapping runs.Dots show the mean values of the bootstrapping runs.

Figure 5 .
Figure 5. Diel and wind direction profiles of factor contributions. Colored triangles and lines in panels (A-C) show 2-hr means of factor contributions. Boxes and whiskers in panels (A-C) are the 5 th , 25 th , 75 th , and 95 th percentiles of the 2-hr composite distributions of the factors.Diel patterns (2-hr means) of NO, NO y , O 3 , PAN (multiplied by 200), wind direction (WD), and solar radiation (SR) are also shown in panels (A-C).Panel (D) shows 30°-wind-direction means of factor contributions.

Figure 6 .
Figure 6.Mean factor contributions to measured VOCs at CAVE during summer 2019 their OH reactivities for the whole duration of CarCavAQS and periods with southerly and southeasterly (S-SE) winds, and northerly and northwesterly (N-NW) winds.

Figure 7 .
Figure 7. Mean normalized contributions of the PMF factors to air sampled in different areas.The factors are the same as the ones shown in Figure 4 but without acetone, MEK, and MVK.The contributions are normalized using all canister means to show the relative importance of a factor in different areas.The black bars show the mean normalized residuals of the fittings.

Figure 8 .
Figure 8. Incremental probability maps of the PMF factors from the RTA (A-G) and a map of potential sources (H).The IPs are based on 48-hour back trajectories and are distance normalized by multiplying each grid cell value by the distance from the origin site (yellow star).The grid cell resolution is 0.25 degrees square.Grids in panel (H) are colored by the number of O&G well pads in the cell (US EIA 2023).Grey lines and areas show major highways and urban areas, respectively.Brown lines show major boundaries of the O&G basins and shale formation.Green areas show national park units.Green crosses in panel (H) show locations of active drilling sites during CarCavAQS.

Figure 9 .
Figure 9. Photochemical clocks using ratios of 2-pentyl nitrate (2-PenONO 2 )/n-pentane and 2-butyl nitrate (2-BuONO 2 )/n-butane for (A) grab canister air samples, (B-D) all observations at CAVE with data points colored by C 2 H 6 , PAN, and O 3 mixing ratios.The lines represent the modeled data.Colored squares, crosses, triangles, and dots in panel (A) represent ratios derived from grab canister samples collected near background, urban/traffic, O&G development, and O&G production areas, respectively.