Significance Our results demonstrate that access to high-resolution spatiotemporal activity data and multiscale, contemporaneous measurements is critical to understanding oil- and gas-related methane emissions. Careful consideration of all factors influencing methane emissions—including temporal variation—is necessary in scientific and policy discussions to develop effective strategies for mitigating greenhouse gas emissions from natural gas infrastructure.

Abstract This study spatially and temporally aligns top-down and bottom-up methane emission estimates for a natural gas production basin, using multiscale emission measurements and detailed activity data reporting. We show that episodic venting from manual liquid unloadings, which occur at a small fraction of natural gas well pads, drives a factor-of-two temporal variation in the basin-scale emission rate of a US dry shale gas play. The midafternoon peak emission rate aligns with the sampling time of all regional aircraft emission studies, which target well-mixed boundary layer conditions present in the afternoon. A mechanistic understanding of emission estimates derived from various methods is critical for unbiased emission verification and effective greenhouse gas emission mitigation. Our results demonstrate that direct comparison of emission estimates from methods covering widely different timescales can be misleading.

The US power and industrial sectors now have the lowest historic carbon intensity, primarily due to increased use of natural gas and renewables (1); however, the natural gas supply chain is a significant emitter of methane—the principal component of natural gas and a potent greenhouse gas (2, 3). Substituting natural gas for more carbon-intensive fuels, such as coal, for power generation has climate benefits provided that methane emissions are less than ∼3% of production (4). Many recent studies have sought to quantify anthropogenic methane emissions at global, national, and regional scales (5, 6), with particular emphasis on emissions from natural gas infrastructure due to increased use of horizontal drilling and hydraulic fracturing (7).

Historically, “bottom-up” (BU) emission estimates, such as the US Greenhouse Gas Inventory (2, 3), have been developed by multiplying average emission factors for each known source category by an activity factor for that source category (i.e., per unit emission rate times the number of units) to estimate the annual emissions from a facility or emission source. These emission factors are developed from direct emission measurements at the component or facility level. BU estimates may be made for specific regions, industry segments, facilities, or sets of infrastructure components. Since emission inventories often guide policy decisions (8), it is critical that BU inventories accurately reflect both total emissions and the relative contribution of different source categories.

More recently, methane emissions have also been estimated using “top-down” (TD) approaches at the regional scale [e.g., an individual oil and gas (O&G) production basin]. In TD approaches, including aircraft mass balance (AMB) and other methods (9⇓–11), measured atmospheric methane concentrations are used in models to infer emission rates. In the AMB method, transects are flown upwind and downwind of the study region (see black arrows in SI Appendix, Fig. S9). A flux is then calculated by taking the difference in methane dry mole fraction measured at downwind and upwind transects while accounting for planetary boundary layer height, horizontal wind speed and direction, and other factors (12). Recent basin-scale comparisons of TD and BU estimates have shown persistent differences between results obtained with the two methods (9⇓–11, 13⇓⇓–16), with TD estimates typically 1.5 times larger than BU estimates (16) for the same region. A better understanding of the TD–BU discrepancy is needed to improve confidence in emission estimates from both approaches.

Several studies (17⇓–19) have suggested that BU estimates may be low due to inaccurate emission factors or underrepresented sources (e.g., “superemitters”) whose overall emissions are not adequately characterized by measurement campaigns employing systematic sampling. However, it is well-known that emissions from O&G infrastructure vary both temporally and spatially (20) and that aggregated and annualized activity and emission factors often deviate substantially from local emissions occurring during shorter time periods. For example, maintenance activities, such as manual liquid unloadings (MLUs) or depressurization of equipment (“blowdowns”), are often triggered by human operators during daytime work-week hours and may produce high emission rates for short durations (21). Prior TD–BU comparisons (9⇓–11, 13⇓–15) had little to no site access or activity data to model specific O&G operations occurring contemporaneously with short duration, regional TD measurements. In lieu of such information, those studies utilized BU activity estimates drawn from annualized national or regional activity inventories. Furthermore, component and facility measurements are rarely made contemporaneously with TD aircraft estimates. As a result, BU estimates often use emission factors based on measurements aggregated at the national or subnational level, typically years before the TD measurement. To eliminate confounding effects resulting from spatiotemporal misalignment, multiple studies (8, 11, 16, 18, 22) have suggested that TD estimates be compared with BU estimates assembled using contemporaneous activity data and emission measurements.

This study addressed the issue of spatiotemporal misalignment, which is present in all prior TD–BU comparisons in O&G basins. Six measurement teams (12, 23⇓⇓⇓⇓–28) made a combination of coordinated basin-, facility-, and component-level measurements of key sources within the “study area” during a 4-wk field campaign in September–October 2015. The study area spans ∼150 km east–west and ∼65 km north–south and covers the eastern portion of the Fayetteville Shale play in Arkansas, as shown in SI Appendix, Fig. S9. In addition to site access for measurement teams, study area operators provided an extensive set of activity data detailing operations during the “study period” [October 1 and 2, 2015, the days corresponding to TD AMB flights by Schwietzke et al. (12)]. See Methods for a description of the role of study area operators. In a significant advance over prior TD–BU comparison studies, measurement and activity data for key sources within the study area were combined with literature data to develop a spatially and temporally resolved BU Monte Carlo (29) model. The BU model developed herein combines and extends those employed in Bell et al. (28), Vaughn et al. (26), and Zimmerle et al. (27). The BU model also allows a direct comparison with TD AMB estimates of Schwietzke et al. (12) during TD measurement windows on both days of the study period. Study period meteorological measurements were used to propagate facility-scale emissions downwind to ensure consistent comparisons with TD aircraft measurements. Unlike earlier work, this experimental design enabled kilometer-scale, hourly resolution BU estimates and provides specific methodological advances for future comparative studies.

Results Temporal Variations in Emissions. Hourly results from the gridded BU model for the total study area are shown in Fig. 1. Study area total emissions exhibited significant variability throughout the day. On both days of the study period, modeled emissions peaked during midafternoon hours due to MLUs performed and recorded by production facility operators. Study area operators recorded 107 MLU events over the 2-d study period, which can be compared with the daily average (54 MLUs per d) derived from the cumulative number of MLUs reported by the same operators to the US Environmental Protection Agency Greenhouse Gas Reporting Program (EPA GHGRP) (30) in 2015. Fig. 1. Hourly averaged methane emissions for the study area estimated by the BU model. Sources in the natural gas production and gathering segments together account for 78% of the study period average emission rate. Emissions exhibit a strong diurnal variation, primarily driven by MLUs in the production sector, which are triggered by operators during daytime hours. Error bars shown represent a 95% CI, as given by the 2.5th and 97.5th percentile of Monte Carlo results, on each hourly BU estimate. The spatial distribution of emissions within the study area for the October 1 TD measurement window is shown in Fig. 2. An animated summary of hourly emission rates is provided in Movie S1 and described in SI Appendix, section S1.3. For this study area, MLUs were the largest emission source in the production sector, and gas quantities released during MLUs were large enough to impact basin-level emissions (other basins may exhibit similar or different behavior depending on the population of MLUs in the specific basin). MLUs serve to clear accumulated liquids from a wellbore to maintain a desired gas production rate. MLUs are initiated and terminated by operators during work hours and typically begin in the morning. During the study period, the average duration of an MLU was 4.5 h, with a range from 15 min to 22 h. Consequently, the combined emissions from multiple MLUs created a basin-scale, midday peak in emissions. The peak BU estimate, which was roughly double that of emissions occurring during early morning hours, occurred ∼1–3 h before the TD measurement windows. TD measurements were made when atmospheric conditions were most suitable for the AMB approach. Late-afternoon flights, between ∼1 and 4:30 PM local time, allowed maximum vertical mixing of near-surface emissions throughout the planetary boundary layer (see ref. 12). Natural gas gathering sector emissions did not exhibit the significant diurnal pattern seen in the production sector (Fig. 1). Gathering compressor stations were operated at high utilization throughout the study period and therefore exhibited little temporal variation in emissions. The largest methane emission source for gathering stations in the study area (other basins may be different) was unburned methane entrained in compressor engine exhaust (“combustion slip”). Basin-scale emissions from combustion slip are affected mainly by the number and type of compressor engines operating, which varied little during the 2-d study period (26). Other gathering sector sources with variable emissions included engine starts, blowdowns of pressurized equipment (which occurred when compressors were taken off-line), and abnormal process conditions. These variable sources proved to be small relative to MLUs during the study period. Emissions from remaining source categories shown in Fig. 1 exhibited negligible diurnal variation, either because no variation existed (e.g., fully automated systems or time-invariant natural processes) or the data necessary to model diurnal variations were unavailable (e.g., temperature-dependent changes in methane from livestock manure handling). These time-invariant sources (transmission, livestock, geologic seeps, and wetlands) contributed less than 25% to the study period average BU emission rate. Source categories contributing less than 1% each to the study period average emission rate are omitted from Fig. 1 for clarity (e.g., natural gas distribution, rice cultivation, landfills, wastewater treatment, etc.; SI Appendix, S4). Spatial Variation in Emissions. Operations (and emissions) within the study area vary spatially as well as temporally. Emissions from each source category in the BU model were assigned to their known location on a 0.04° (∼3.8 km) grid as shown in SI Appendix, Fig. S9. Total BU methane emissions during the TD aircraft measurement window on October 1 (Fig. 1) are shown in Figs. 2 and 3 and in greater detail in SI Appendix, Fig. S1. Spatial assignment of emissions varied by source category based on available data; for example, livestock populations were available at the county level, whereas exact coordinates were available for most natural gas infrastructure. Fig. 2. Comparison of TD–BU longitudinal emission profiles. Simulated BU longitudinal emission rate profiles compared with profiles of Schwietzke et al. (12). CIs overlap for 73% of the east–west distance for October 1 and 66% for October 2. Data from ref. 12. Fig. 3. Results from the spatially and temporally resolved BU model agree (i.e., 95% CIs overlap, shown by error bars) on two consecutive days for the total study area as well as the eastern and western subportions; A compares aggregate BU to TD emissions during the AMB flight on October 1, and B compares aggregate BU to TD emissions during the flight on October 2. Potential explanations for differences in mean estimates between TD and BU are explored in Discussion. In a significant advance over prior studies, improvements in the AMB method and favorable meteorological conditions enabled temporally and spatially resolved TD estimates (12) on two consecutive days. To enable a direct comparison with aircraft results, a Gaussian dispersion model (31, 32) was used to simulate the transport of emissions to the downwind transect of the AMB flights. Emission contributions from each grid cell were summed into 0.02° longitudinal bins [consistent with the TD estimates (12)] at the southern edge of the BU model grid, which corresponds to the aircraft’s flight path during measurements. The simulated BU longitudinal emission profiles are strikingly similar to TD estimates during the midafternoon period of TD AMB flights (Fig. 2) but dissimilar at other times of day when MLU emissions are lower (Movie S1). CIs overlap for 73% of the east–west distance for October 1 and 66% for October 2. The similarity in features observed in this TD–BU comparison of longitudinal emission profiles shows that spatially and temporally resolved BU models can reproduce aircraft observations and provides confidence that both (independent) methods captured the same emissions phenomena. Aggregate BU results also compare well to TD estimates provided in the companion article by Schwietzke et al. (12). Modeled BU emissions from each grid cell were aggregated for the eastern subregion, the western subregion, and the entire study area. Results indicate agreement between TD and BU emission estimates (as defined by overlapping 95% CIs) on both days (Fig. 3) for the eastern and western subregions as well as the total study area. Variability present in all input emissions and activity data are propagated through the BU Monte Carlo model, resulting in a distribution of possible emissions for each modeled category. BU 95% CIs are given by the 2.5th and 97.5th percentiles of emissions predicted by the model. TD 95% CIs are derived from uncertainties in measured wind speed, wind direction, background dry methane mole fraction, and planetary boundary layer height, as described in supporting information section 6 of ref. 12. Missing aircraft vertical wind profiles on October 2 result in greater TD relative uncertainty compared with October 1 (41% and 29% for study area total emissions, respectively).

Discussion Although TD and BU results agreed based on overlapping 95% CIs, the BU model predicted lower mean emissions than TD estimates for most comparisons shown in Fig. 3. The degree of underprediction depended on study area subregion (eastern or western) and day (October 1 or October 2). Several factors could explain the remaining difference in means between TD and BU estimates, as well as local differences in longitudinal emission rate profiles (Fig. 2). Activity Data Accuracy and Uncertainty. Access to on-site operational information enabled a mechanistic understanding of the temporal nature of emission sources, which proved critical in reconciling TD and BU estimates. However, despite access to an unprecedented level of operator information, some uncertainty remains in the exact timing of MLUs, engine starts, and other intermittent, worker-stimulated events that are not typically logged electronically. The degree to which this uncertainty affects the accuracy of TD–BU comparisons depends upon (i) the hourly rate of change of the aggregate BU estimate and (ii) the instantaneous emission rate of individual emission sources that contribute to that estimate. For example, an aggregated shift of the flight window by ±1 h (analogous to shifting a mix of simulated intermittent emission events into, or out of, the flight window) changes BU emissions during flight windows by only ±2% on October 1 but by +17%/−9% on October 2, based exclusively on the mix and timing of intermittent emission events. In addition, the mass emission rates of MLUs are assumed constant within modeled time periods despite evidence that MLU emission rates vary at subhourly timescales (21). A specific example can be seen in Fig. 2. The TD longitudinal emission rate profile on October 2 (bottom panel) shows a large source at approximately −91.75° longitude that is absent in the BU profile. This source (or a very similar source) is included in the BU model (Movie S1) but appears before the TD flight window, based on provided activity data. Temporal uncertainty may also influence the TD estimate because TD estimates were based on multiple transects of the study area, which took ∼2 h to complete (1 h, 54 s on October 1 and 1 h, 42 s on October 2); a different mix of intermittent emission sources may have been encountered on subsequent transects. The averaging of multiple individual transects may blur the contribution of distinct intermittent emission sources. For example, the BU model predicts different longitudinal emission rate profiles during each aircraft transect on October 1, as shown in SI Appendix, Fig. S17. The difference is apparent in the western portion of the study area, where a greater number of MLUs were reported by study area operators. Emission Factor Accuracy and Applicability. To investigate potential reasons for the remaining disagreement between the mean TD and BU emission estimates depicted in Fig. 3 for October 1, alternative scenarios of BU estimates were made using adjusted emission factors. Emission factors for individual source categories were adjusted in each scenario, as shown in Table 1, to assess the sensitivity of the TD–BU comparison with the uncertain emission factors. A TD–BU match (difference in total emissions of 1% or less) was achieved by increasing the mean MLU emission rate by 37% or increasing the mean gathering station emission rate by 89%; other source categories required fourfold or greater increases over emission rates used in the base BU model. Adjustments would be lessened if multiple categories were considered simultaneously. Further discussion focuses on the first two scenarios involving emissions from MLUs and gathering stations. Table 1. Emission rates from individual source categories in the BU model were increased to explore potential explanations for differences in the means of TD and BU estimates BU MLU emissions were modeled using actual activity data for the study period from study-partner operators and emission rates from Allen et al.’s (21) measurements of eight horizontal wells located in the same production region (“Mid-continent”) as Fayetteville (21). The mean time-averaged methane emission rate from Allen et al.’s (21) measurements was 513 kg/h. The 37% increase in mean MLU emission rate (i.e., a mean emission rate of 703 kg/h) that results in TD–BU match falls within the range of the underlying measurements (247–1,253 kg/h) and thus appears plausible. Furthermore, the mean emission rate of the one MLU measured (using the downwind tracer method) during this field campaign (810 kg/h, 95% CI 603–1,018 kg/h) was above the mean from Allen et al. (21) (513 kg/h) and the adjusted scenario (703 kg/h). Finally, adjusting the MLU emission rate in the BU model resulted in nearly perfect spatial TD–BU agreement for the eastern and western subportions of the study area. Thus, a plausible adjustment in one critical emission factor in the BU model substantially improves the agreement between mean TD and BU estimates. Matching the TD study area total by increasing gathering station emission rates required an 89% increase from the base case and did not improve east–west agreement in mean emission rate. Further, in contrast to MLUs, none of the emission sources within gathering compressor stations could easily account for an 89% increase in emission rate. For example, the largest component of gathering emissions—methane entrained in the exhaust of compressor engines, commonly called combustion slip—accounted for ∼78% of emissions from normally operating gathering compressor stations in the study area (26). Combustion slip emissions would need to increase by 150% to increase total gathering emissions by 89%. Such an increase does not agree with data from 111 recent (within year, but not during the field campaign) measurements of combustion slip from engines in the study area (26). Underrepresented, Unobserved, and High-Emitting Sources. Another potential source of underestimated emissions in the BU model stems from “abnormal process conditions,” which here refer to undesired equipment or process failures that release natural gas to the atmosphere. Abnormal process conditions that result in large emissions relative to similar processes under normal conditions, termed superemitters, have received extensive attention (17, 19, 33, 34) as a possible cause for disagreement between TD and BU estimates. One common example is a dump valve (a valve used to empty liquid from drop-out tanks on compressors and gas lines) that fails to close, allowing gas to escape through liquid storage tanks to the atmosphere. Abnormal process conditions leading to tank emissions beyond the measurement capabilities of on-site teams were observed at two gathering stations during the field campaign. At the first, a large (∼150–600 kg/h) emission source was identified on multiple days by the aircraft during basin survey raster flights (distinct from the mass balance flights discussed earlier). At the second, on-site measurement teams using optical gas imaging observed significant tank venting, which was estimated using tracer measurements at ∼140 kg/h. An 89% increase in gathering station emissions based on similar abnormal process conditions alone would require a 10-fold increase in their frequency, which is unlikely given field observations in this study and prior studies (35, 36). Abnormal process conditions often cannot be measured during field campaigns owing to their size and infrequency. Prior studies like Zavala-Araiza et al. (17) or Zimmerle et al. (37) used statistical estimators of potentially unobserved emissions to reconcile BU estimates with TD measurements. In contrast, we used available measurements to model emissions from abnormal process conditions and additionally assumed that (i) the frequency of occurrence in each hour matches that observed across the whole 4-wk field campaign and (ii) the site of these emission sources was randomly located within the study area during any given hour (SI Appendix). A combination of the factors discussed above could plausibly account for remaining discrepancies between the TD and BU estimates. A hypothetical scenario shown in SI Appendix, Fig. S4 depicts nine “estimated potential sources” added to the BU model during the first TD transect on October 1. Added sources totaling 6.6 Mg/h improved the match of aggregate TD–BU estimates (west, east, and total) and produced a longitudinal emission rate profile whose 95% CI overlapped with TD for 89% of the East–West distance modeled (SI Appendix, S2). Protocol and Analysis Recommendations. These coordinated basin-, facility-, and component-level measurements, and the subsequent modeling effort, highlight that TD and BU methods can be reconciled solely with contemporary, multiscale measurements and high-resolution spatiotemporal activity data. Short-duration basin-scale TD measurements made in locations with large temporal variability in emissions do not warrant comparison with BU emission estimates based on annual averages. A more robust protocol that utilizes detailed activity data, unfettered and unbiased site access, and time-resolved operations data for the study period (as employed herein) is necessary to ensure a definitive and meaningful comparison. Activity data should represent, at as fine a time resolution as possible, emissions that would propagate downwind and be captured by the aircraft during measurement. Unprecedented access to high-resolution activity data substantially closed the gap observed in many prior TD-to-BU comparisons; however, even hourly resolution was insufficient to resolve all source behaviors captured by aircraft measurements, as emissions from large intermittent sources, like MLUs and blowdowns, can fluctuate greatly at subhourly timescales. Direct measurements of key sources should be made at the time and location of the study (or as close as possible)—as were made for most sources in Fayetteville during our field campaign—to capture study-relevant emission characteristics. Due to resource limitations, MLUs could not be measured directly in this study (Methods), and analysis herein illustrates that differences between regional and study-area emission rates for MLUs alone could explain the difference in TD and BU emission estimates. Further, episodic MLUs are responsible for roughly doubling emissions at the midday peak over other steady-state sources in the study area. Understanding the relative contribution from large episodic sources (like MLUs), which occur as a part of normal operations in active O&G basins, has implications in developing sound policy for emission reductions. Our study area represents a specific test case for a basin-level comparison of TD and BU emission estimates. The mix of emission sources and their temporal behavior will likely vary substantially among basins, due to differences in facility age, gas composition, production volume, infrastructure, operator practices, and other factors. However, the methods presented here [and in the many supporting papers from this overall study (12, 23⇓⇓⇓⇓–28)] are broadly applicable and, when used together, can lead to improved understanding of emission sources from a production region. Understanding the relative contribution of specific sources, and their time dependency, forms a sound basis for prioritizing emission reduction at the regional level. This understanding would be difficult to achieve using conventional TD or BU estimates alone.

Methods Study Area and Field Campaign. The study area (SI Appendix, Fig. S9) was selected due to the compact and isolated nature of the basin and availability of natural gas industry participants to provide site access and activity data for the BU model. The study area covers 10,100 km2 in the eastern Fayetteville Shale play, an active exploration and production subbasin of the Arkoma basin that produced ∼75,000 MMcf (million cubic feet) of natural gas in September 2015. Approximately 80% of 5,652 currently active and producing wells in the study area were drilled after 2008 (38). The gas produced is dry (90–98% methane), with essentially no natural gas liquids or heavier hydrocarbons, although most wells coproduce water. There are no gas storage facilities or processing plants in the study area. Only the western edge of the study area abuts a neighboring production region, providing a clear delineation of the study area. Isolation of the study area was enhanced by persistent north-to-south winds on October 1 and 2, 2015. The study team was supported by three companies that collectively operate 99.1% of active wells (accounting for 99.8% of September 2015 production) and 110 of 125 gathering compressor stations in the study area. Two production partners, representing 88% of production and 84% of producing wells, and 99 gathering compressor stations provided both activity data and site access for facility-level measurements. A third study partner provided activity data only. Two additional study partners provided measurement access to their four transmission compressor stations, out of six known transmission facilities in the study area. The local natural gas distribution company provided access and activity data for all distribution operations (which were not extensive owing to the rural character of the study area). Measured facilities were randomly selected, and in all cases strict procedures were in place to ensure that (i) facility operations were not altered by partners before, or during, facility measurement and (ii) all activity data were objectively and accurately collected and interpreted. Measurement details are provided in companion papers, which document measurement protocols and results for production facilities (28), gathering compressor stations (26), distribution systems and gathering pipelines (27), and downwind (24, 25) and aircraft (23) facility-level measurement methods and results. TD AMB. The AMB method implemented in this study is described in detail in a companion paper (12) and is briefly summarized here. A Scientific Aviation light aircraft (39) was used to measure CH 4 , C 2 H 6 , CO 2 , H 2 O (at ∼0.5 Hz), and horizontal winds (averaged at ∼0.1 Hz) along upwind and downwind transects north and south of the study area. Species measurements were also made during vertical profiles through the planetary boundary layer (PBL). Throughout the duration of the field study, a 915-MHz boundary-layer radar wind profiler (BLRWP) was deployed in the study area to measure horizontal winds in the PBL. The aircraft vertical profiles and BLRWP measurements were used to (i) determine PBL depth, (ii) corroborate that the PBL was well-mixed vertically (such that lateral downwind measurements at distinct altitude levels were representative of the PBL), and (iii) to verify that no airmass recirculation occurred in the study area immediately before each flight (to satisfy the steady-state mass conservation principle). Total emissions of CH 4 were then calculated as E C H 4 = ∫ − b b Δ X C H 4 v c o s θ d x ∫ z G r o u n d z P B L n A i r , d r y d z , where ΔX CH4 and v are the measured CH 4 mole fractions relative to background and horizontal wind speed, respectively. Background CH 4 mole fractions were estimated for each downwind transect based on the CH 4 mole fractions at both edges of the downwind transect plume and along the transect upwind of the study area. ΔX CH4 and v are integrated over each plume width increment across the study area (−b to b), corrected for the mean wind direction normal to the flight track ( cos ⁡ θ d x ). The term n Air,dry accounts for measured dry air molar density (to compute CH 4 emissions in mass units based on measured CH 4 mole fractions) across the vertical from ground elevation (z Ground ) to the PBL top (z PBL ). BU Model. TD–BU reconciliation was achieved by constructing a comprehensive spatiotemporal Monte Carlo model of methane emission sources. The model used emissions measurement data, high-resolution spatiotemporal activity data, engineering calculations, and an understanding of operational characteristics gleaned from study-partner operators. The BU model calculates emissions based upon empirical emission and activity distributions and uses Monte Carlo methods to propagate uncertainty. Methods follow those of Ross (29) and recent national models developed for the transmission and storage (37) and gathering and processing (36) sectors. Spatial modeling of emissions was performed by assigning emissions to the BU model grid by spatial intersection with source activity data. For example, latitude and longitude were known for O&G facilities, livestock activity data (40) were available at the county level, and wetland activity data (41) were available directly in geospatial format. Grid cells for the BU model grid are based on 0.04° longitude increments. Each grid cell on the southern edge of the BU model grid contains two receptor points. These receptor points are used in a Gaussian dispersion model to develop simulated BU longitudinal emission rate profiles. Temporal variation was modeled for AMB flight windows and at a 1-h resolution for the 48-h study period (October 1 and 2, 2015). Results from each modeled period represent the average of emissions that occurred during that period. On each iteration of the Monte Carlo model, calculated methane emissions from each grid cell were propagated downwind, according to Gaussian dispersion theory, based on measured prevailing wind speed, wind direction, and atmospheric stability class during the AMB flights. Detailed descriptions of emission calculations for each category in the BU model are provided in SI Appendix. MLUs. Facility location and basic equipment inventories for most categories of O&G facilities (e.g., well count on well pads) were available from public sources, such as the Arkansas Oil and Gas Commission (38). These public data were extended and augmented by partner data to provide a better description of facility equipment and to capture temporal variation in operations which impact emission estimates from key source categories. For example, the annual number of manual unloading events is reported in the EPA GHGRP, but the timing of the events is not provided. These variations in activity data, and the corresponding impact on emission rates, are integral to the TD–BU comparison and were uniquely captured for this study. The average duration of an MLU was 4.5 h, with a range from 15 min to 22 h. With one exception, emission rates from MLUs were not measured during the field campaign. Emissions from MLUs were modeled in the BU estimate using eight measurements of horizontally drilled wells taken by Allen et al. (21) in the midcontinent area, which is the same National Energy Modeling System region (42) as the Fayetteville Shale. MLU events during this study were also from horizontally drilled wells. In the Allen et al. (21) data, the time average emission rate for MLUs varied by a factor of five (513 kg/h, range 247–1,253 kg/h). In comparison, the average emission rates measured during this field campaign at gathering compressor stations (including emissions from abnormal process conditions) and well pads without vented liquid unloadings were 74.5 kg/h and 0.4 kg/h, respectively. Activity data used to model MLUs were provided by study partners after the field campaign and included the location, estimated start time, and estimated duration of each MLU during the 2-d study period. BU Model Input Data. Input data to the BU model are described in detail in SI Appendix. Additional input data for MLUs at production facilities and compressor engines at gathering stations are included in BU_Input_Data_Study_Period.zip (https://hdl.handle.net/10217/190251), which contains two files. Start times and durations are provided for MLUs in ProductionManualLiquidUnloadingTiming.txt along with counts of production facilities (well pads) and counts of individual wells within BU model grid cells. Compressor engine counts, with horsepower by engine type, are provided in GatheringCompressorEngineCountByType.txt, along with counts of gathering facilities within BU model grid cells. BU Model Output Data. BU model results are summarized graphically for flight windows in SI Appendix, section S1. BU model results corresponding to TD AMB flights on October 1 and 2, 2015 are shown in SI Appendix, Figs. S1 and S2, respectively. Hourly BU model results are summarized in Movie S1 for the 48-h study period spanning October 1 and 2, 2015 as described in SI Appendix, section S1.3. Additionally, results from modeled source categories and subcategories are provided in tabular format in BU_Output_Data_Flight_Windows.zip and BU_Output_Data_Study_Period.zip (https://hdl.handle.net/10217/190251).

Acknowledgments We thank all field measurement personnel and those who aided in data compilation for their significant, coordinated efforts. This work was supported by Research Partnership to Secure Energy for America/National Energy Technology Laboratory Contract 12122-95/DE-AC26-07NT42677 to the Colorado School of Mines. Cost share for this project was provided by Colorado Energy Research Collaboratory, the National Oceanic and Atmospheric Administration Climate Program Office, Southwestern Energy, XTO Energy, a subsidiary of ExxonMobil, Chevron, Statoil and the American Gas Association, many of whom also provided operational data and/or site access. Additional data and/or site access was also provided by CenterPoint, Enable Midstream Partners, Kinder Morgan, and BHP Billiton.

Footnotes Author contributions: T.L.V., C.S.B., C.K.P., G.A.H., G.P., D.J.Z., R.C.S., and D.N. designed research; T.L.V., C.S.B., C.K.P., S.S., G.P., D.J.Z., and R.C.S. performed research; T.L.V., C.S.B., C.K.P., G.A.H., G.P., and D.J.Z. analyzed data; and T.L.V., C.S.B., S.S., G.A.H., G.P., D.J.Z., and D.N. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

Data deposition: Datasets are available at https://hdl.handle.net/10217/190251.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1805687115/-/DCSupplemental.