Creating Measurement-Based Oil and Gas Sector Methane Inventories using Source-Resolved Aerial Surveys

We present a new framework for incorporating aerial measurements into comprehensive oil and gas sector methane inventories that achieves robust, independent quanti�cation of measurement and sample size uncertainties, while providing timely source-level insights beyond what is possible in current o�cial inventories. This “hybrid” inventory combines top-down, multi-pass aerial measurements with bottom-up estimates of unmeasured sources leveraging continuous probability of detection and quanti�cation models for a chosen aerial technology. Notably, the combined Monte Carlo and “mirror-match” bootstrapping technique explicitly considers skewed source distributions and �nite facility populations that have not been previously addressed. The protocol is demonstrated to produce a comprehensive upstream oil and gas sector methane inventory for British Columbia, Canada, which while approximately 1.7 times higher than the most recent o�cial bottom-up inventory, reveals a lower methane intensity of produced natural gas (< 0.5%) than comparable estimates for several other regions. Finally, the developed method and data are used to upper bound the potential in�uence of source variability/intermittency on the overall inventory, directly addressing an open question in the literature. Results demonstrate that even for an extreme case, variability/intermittency effects can be addressed by sample size and survey design and have a minor impact on overall inventory uncertainty.


S1.1 Determination of Active Facilities
The aerial survey was planned using available activity and production data from the British Columbia Oil and Gas Commission (BCOGC) and the PETRINEX reporting system website (Petrinex, 2022a) with the simultaneous objectives of maximizing survey coverage and minimizing sample size uncertainties in the derived emissions inventory.A set of active oil and gas facilities and wells were first identified using Petrinex volumetric monthly production data obtained from BCOGC and supplemental active and suspended facility lists from the Petrinex website.Because publicly available active facility lists are not always accurate or up to date, facilities and wells required to report under the Petrinex system were instead deemed active or inactive based on the existence of reported activity data during the month of the survey.For compressor stations, which do not directly report monthly activity via Petrinex, an initial active count was derived from public activity lists and subsequently checked and updated where necessary based on review of aerial imagery of survey sites (e.g., in cases where images showed all compressors had been removed from the site) and provincial leak detection and repair (LDAR) reporting.As detailed in Table S1, 1,006 facilities within the province of British Columbia (BC) were identified as active during the aerial survey of which 601 (60%) were measured.
The active status of individual wells was also gleaned from Petrinex data where well production volumes are found under one or more unique well identifiers (UWI) linked to facilities reporting to Petrinex.These UWI represent segments of a well and were aggregated to shared surface-holes (wells) using well authorizations (WA) assigned by BCOGC in  As detailed in Table S1, this analysis identified 8,995 active wells within the province, of which 904 (10%) were captured in the aerial survey.Wells and associated production equipment (e.g., separators, line heaters, pump buildings, etc.) maybe co-located with facilities on common pad ("on-site wells") or reside at completely separate well-site location ("off-site wells", OSW).Since equipment associated with on-site wells may not be distinguishable from that associated with the facility, aerial measured equipment sources at these locations were assigned to the facility.Within the inventory a well was deemed to be off-site if the wellhead surface location was located in a different land grid location (i.e., legal subdivision or NTS quarter unit) than the facility where it reported production.The majority of wells in British Columbia, 97% (8,729 of 8,995), were considered to be off-site; these constituted ~78% (705 of 904) of surveyed wells.
For facilities and wells appearing in Petrinex there were 22 associated activity codes (Petrinex, 2022b) in the monthly volumetric data which track production, flaring, venting, transfers (receipts and dispositions), storage (injection), and losses (shrinkage and metering differences) of produced and processed fluid volumes, as well as a "shut-in" status.Shut-in facilities and wells in the monthly data are sites that are capable of processing, producing, or injecting but were inactive during an entire reporting month.Shut-in sites may have been active in months prior to the survey and may becoming active again or may ultimately be moved to a suspended status and abandoned.
For all results presented in the manuscript, shut-in facilities and wells were conservatively assumed to be non-emitting.However, it is possible that some shut-in sites could remain fully or partially pressurized thus having a potential to emit.Reviewing Table S1, at the time of the survey there were an additional 95 "active" but shut-in facilities (9% increase over the 1006 non-shut-in facilities) and 9978 "active" but shut-in wells (11% increase over the 8995 non-shut-in wells).
To bound the potential impact of these shut-in facilities/wells, the entire inventory analysis was repeated including these facilities as active.As shown below in Table S4 and Table S5, including these sites results in a small increase of 4.1 kt (2.8%) in the total inventory.Most of this difference is a 2.7 kt increase in pneumatic equipment emissions, which could be expected if the pneumatic equipment at these shut-in facilities and wells remains pressurized and emitting at expected steady bleed rates.For this reason, regulatory clarification of what qualifies as shut-in, and if possible, differentiating between pressurized and ready to produce versus sealed at the wellhead, is highly recommended.By contrast, including shut-in facilities and wells in the analysis made little difference in the measured source portion of the inventory (112.9 kt including shut-ins, or a 0.6% increase over the 112.2 kt measured source total presented in the manuscript).

S1.2 Sampling Region and Strata
During the initial survey planning phase, candidate survey location were identified by geo-locating all active facilities and wells in ArcGIS Pro using BCOGC permit data (BCOGC, 2022a(BCOGC, , 2022b)), land grid locations (Dominion Land Survey legal subdivisions and National Topographic Systems quarter units), and BCOCG surface-hole locations for wells (BCOGC, 2022c).A review of these candidate sites identified a significantly lower facility density in the northern region centered around Fort Nelson, Figure S1.Given the sparseness of facilities and budgetary constraints, the initial set of candidate sites was constrained to a smaller survey region of approximately 46,000 km 2 south of approximately 58°N.Limiting the area covered by the aerial survey was also essential for the feasibility of conducting parallel on-site follow-up investigations of detected sources (Johnson et al., 2022).An initial set of polygons bounding sites for aerial study were manually specified for facility and well locations in this sub-region.These initial polygons were chosen considering the density of sites (to maximize the economics of the survey and the sample size) and the underlying distribution of sites across the sampling strata (to maximize the relative sample size within each stratum).In the final sample, all verifiably active facilities that could be reliably located and had adequate/current satellite imagery within the survey region were included.Initial polygons were provided to Bridger Photonics for flight planning and cost estimation; the remaining budget was leveraged to measure as many well sites within the survey region as feasible.Aerial measurements were ultimately performed over 508 geographically unique sites (polygons) during September 11 to October 8, 2021.Figure S1 maps the geographic distribution of active and measured sites (facilities and wells) during the measurement survey.
As introduced in the manuscript, the developed inventory protocol uses stratified sampling, in which common facility and well types are aggregated into separate strata.Given the complex diversity of methane sources in the UOG industry, parsing sources into strata has some significant benefits supporting the broad objectives of this survey.First, aggregation of like sources tends to reduce the variance of desired statistics describing each sufficiently large stratum; this corresponds to improved precision in a stratum's calculated mean emission rate and total emissions (i.e., emissions inventory).Second, stratification can be used to combine similar but uniquely different entities with limited sample and/or population sizes; this artificial enhancement of sample/population sizes may come at the cost of increased variance but can enable consideration of these entities using robust analytical methods that need sufficiently large sample sizes.Finally, if strata are defined such that each potential source across the province is contained within one and only one stratum (i.e., the strata do not overlap and provide comprehensive coverage), then the provincial emissions inventory is simply the sum of each stratum's inventory.This allows independent analysis of each stratum, permitting stratum-dependent methodologies that may leverage prior information about the strata.Moreover, this approach provides the relative contribution of each stratum to the whole, which may be informative for regulatory efforts to mitigate emissions.
Within the Petrinex production data, inventory strata for oil and gas facilities (e.g., batteries, gas gathering systems, gas plants, meter stations, etc.) were naturally defined by existing industry assigned facility subtype codes.For British Columbia there are 53 possible subtypes that group production, process, metering, storage, disposal/waste, and treating facilities into specific categories based on characteristics of the underlying site operations (e.g., single or multi-well batteries, oil or gas production, handling sweet or sour gas, etc.).Active facilities gleaned from the above analysis of monthly production data included 33 unique facility subtypes, which are listed alongside population and sample sizes of the survey in Table S1.Six of these subtypes were combined into two larger categories (meter stations and tank farms) to bolster strata size and simplify analysis; thus, the present inventory analysis considered 29 unique facility-related strata.
Strata for wells considered three well typesgas, oil, and water.Gas well and Oil well strata were derived by combining BCOGC defined wellbore fluid types i) Acid Gas (AGAS), Gas (GAS), and Multiple Gas (MGAS) to Gas; and ii) Oil (OIL), Multiple Oil (MOIL), Multiple Oil and Gas (MOG) to Oil.Wellbore fluid for each surface-hole WA was assigned from the publicly available "Well Surface Hole Locations (Permitted)" file.(BCOGC, 2022c)

S2.1 Stratum-Level Measured Source Inventories
Beginning with flight pass-level data provided by Bridger, measured emission inventories (within quantified uncertainties) for each stratum were computed in a statistical framework considering measurement quantification accuracy, detection sensitivity, and finite sample size effects.This is possible via the nested algorithms described in this section and summarized at a high-level in Figure 1 of the manuscript.
For each iteration of the analysis, probabilistic average emission rates for each detected and quantified source within each stratum are first computed from pass-by-pass aerial data via the algorithm detailed in Sections S2.1.1 to S2.1.3.Briefly, for each source that is detected one or more times during the aerial survey, Bridger-quantified emission rates are randomly perturbed according to the recently developed quantification error model for Bridger's GML by Conrad et al. (2022).As discussed in Sections S2.1.1 to S2.1.3,it is possible that a source detected during one or more pass of the aircraft may not be detected during one or more of the other passes; this could be due to variability/intermittency of the source and/or the finite detection sensitivity of GML.Given the prior knowledge of an existing source, these "missed detections" are randomly perturbed from zero via a Bayesian analysis according to GML's probability of detection function (Conrad et al., 2022), estimated 3-m wind speed at the time of the flight pass, aircraft altitude, and the quantified emissions data during other flight passes over the source.This provides a randomized, true emission rate for each flight pass on each measurement day for the source.
Recognizing that variability of the source rate between observations can be expected to increase with the time between observations, a randomized, true, average emission rate for the source is obtained by first averaging over all flight passes on each unique measurement day, then averaging over measurement days.
For each facility or well site, the average emission rate(s) of all sources(s) are then summed to yield a randomized, true, total emission rate for that surveyed facility or well padwhich could be zero if no emissions are detected.For aerial survey sites containing only wells, total emissions from shared equipment (e.g., a common separator fed by multiple wells) are equally distributed among the unique wells.Aggregating these facility-and well-level data yields a set of facility/well-level aerial survey emissions (randomly perturbed to consider GML quantification accuracy and detection sensitivity and specifically including measured zeros) for each stratum.
Emissions for sampled sites in each stratum were then scaled to the stratum population to yield the stratum's measured inventory.Except in a few special cases as noted below, this scaling was done using Sitter's mirror-match bootstrap algorithm (Sitter, 1992), which provides a robust probabilistic estimate of the stratum's measured inventory considering the actual (non-smooth, typically skewed) distribution of emissions from sites in the sample as well as finite population effects.The latter is particularly important given that in most cases, the sample represents a significant fraction of the population.
For six strata where the entire populations were sampled (i.e., subtypes 321 -Crude Oil Multi-Well Group Battery, 393 -Mixed Oil and Gas Battery, 505 -Underground Gas Storage, 611 -Custom Treating Facility, 676 -Natural Gas Liquids Hub Terminal, and wells of undefined fluid; see Table S1), the measured inventory was directly quantified and no further scaling was required.
Similarly, for six strata where no emissions were detected (i.e., subtypes 501 -Enhanced Recovery Scheme, 504 -Acid Gas Disposal, 701 -Surface Waste Facility, 901 -Water Source, 902 -Water Source Battery, and water wells; see Table S1) the measured inventory was conservatively assumed to be zero despite the potential for emitters at entities that were not surveyed.For three subtypes (204 -Gas Transporter, 403 -Gas Plant with Acid Gas Flaring (> 1 t/d Sulphur), and 451 -Liquefied Natural Gas Plant; see Table S1) with small populations (3-5 facilities) and only one facility each with detected emissions, the measured inventory was assumed to probabilistically follow a uniform distribution bounded by two simple cases: a) all unsurveyed entities had zero emissions and b) all unsurveyed entities had emissions equal to the single surveyed entity.Lastly, for a single facility subtype (407 -Gas Plant; Fractionation) that was not surveyed, emissions were estimated by computing the population size-weighted average of facility-level emissions for 47 other gas plants (subtypes 401-405).
Finally, the measured inventory for the province was computed by summing the measured inventories of the 33 unique strata.To fully resolve the effects of GML's quantification error and detection sensitivity, flight pass-level emissions were randomly perturbed as described above  MC = 10 4 times (where the subscript "MC" represents the Monte Carlo considering effects related to the performance of the GML technology).For each of the BMC sets of randomly perturbed flight pass-level emissions, the measured inventory of each stratum and, hence, the province was computed  BS = 10 4 times (where "BS" implies that this analysis considers effects related to sample sizes via the bootstrapping or alternate approaches noted above).Ultimately, this provided  MC ×  BS = 10 8 estimates of the provincial measured inventory for which statistics could be obtained.Final measured provincial and stratum inventory statistics are summarized in Table S4 and Table S5.Reference mean population emission factors for each facility and well stratum would be obtainable by dividing each stratum inventory by the corresponding population from Table S1.

S2.1.1 Successfully Detected Emissions
During the aerial survey, all sites with detected sources were re-flown at least once, 1-10 days after the initial flight, where each flight contained up to 5 passes over each source.Given finite detection sensitivities coupled with the potential for source variability and intermittency, it is possible that any source may be detected and quantified during some flight passes but missed during others.When a source is successfully detected, a randomized true value of the instantaneous emission rate for that flight pass can be computed from the measured emission rate  S3 of Conrad et al. (2022): where the unitless coefficients , , and  are 0.891, 3.82, and 0.918, respectively, and  is a randomly drawn number from the standard uniform distribution.

S2.1.2 Bayesian Analysis of Sources with "Missed" Detections During One or More Passes
For any detected methane source, the pass-by-pass aerial measurement data may include one or more passes in which the source was not detected ("missed").This may be due to the probabilistic nature of detection success by the airplane, finite detection sensitivities for the GML technology, source variability considering these first two factors, or source intermittency.This section describes a robust Bayesian approach to analyze sources that have been both quantified and "missed" during different measurement passes.The method derives a probability distribution for the true emission rate of a known source that is missed during an individual pass, leveraging passspecific POD data (based on airplane altitude and wind speed for each pass) as well as quantified source rate(s) from all other passes.
Let  ̃ and  represent the Bridger-estimated and true source rate and  ̃ and  represent the Bridger-estimated and (unknown) true 3-m wind speed during a flight pass.Similarly let ℎ ̃ represent the aircraft altitude above ground level during the pass.Finally, define a unitless binary variable  to signify a successful detection ( = 1) or a missed detection ( = 0; denoted as ¬, "not ").With these definitions, the objective is to derive the conditional probability distribution: That is, the probability of the source's true instantaneous emission rate conditional on the Bridgerestimated 3-m wind speed and aircraft altitude for the missed detection.
Using the Bayesian perspective, which treats all variables as probabilistic random variates, the conditional probability of Eq. ( S2) is proportional to the joint distribution of all variables: The righthand side of Eq. ( S3) can be re-stated as As is typical in Bayesian analyses, the normalizing constant of proportionality (necessary to satisfy the law of total probability) can be obtained by integrating the righthand side of Eq. (S6) over  such that Two of three elements on the righthand side of Eq. (S7), the probability of detection function ((, , ℎ ̃)) and the error model for the 3-m wind speed ((| ̃)), are available in the literature (Conrad et al., 2022).Thus, the presented analysis only requires a choice of prior distribution for the true source rate.Theoretically an "uninformed"/non-negative prior, which fixes   () to a constant value for all positive values of  could be used.However, this is overly simplistic since it ignores measurement data from other passes where the same source was detected and quantified.Additionally, the "uninformed" prior does not impose an upper-bound on a missed detection rate.To address this latter point, in the present work missed detection rates are upper bounded by the largest measured emission rate for that source during the measurement survey.
Thus, by defining   () as the uniform distribution from 0 to  ̂, where  ̂ is the largest detected and quantified true emission rate of the source (obtained from data calculated using Eq. ( S1)), Eq.
(S7) simplifies to Finally, as in Section S2.1.1,if  is a randomly drawn number from the standard uniform distribution, then a randomized value of the emission rate during a missed detection given an estimated 3-m wind speed and aircraft altitude can be obtained by integrating the conditional distribution of Eq. (S8) and solving for : which, for the published forms of the probability of detection function and wind speed error model, requires numerical integration and root-finding methods.

S2.1.3 Procedure for Averaging Source Measurements During Different Passes and Flights
Flight pass-level emission rates for each source are averaged as described by Tyner and Johnson (2021) after randomization of flight pass level emission rates as described above.Firstly, randomized flight pass-level emission rates are averaged over each measurement date and these are then averaged over the measurement dates in the survey.Let   be the randomized true emission rate during the n th flight pass on the m th measurement date.Furthermore, define  as the total number of measurement dates and   as the total number of flight passes on the m th measurement day for which the source lies within the GML sensor's field of view.With these definitions, a randomized, true, average source rate during the measurement campaign ( ̅ ) for the source is:

S2.2 Unmeasured Sources -Site-Level Emission Factor Development
As detailed in the main text, unmeasured sources are those that are not detected during any flight pass of the measurement survey.Depending on the underlying source distribution, measurement conditions, and sensitivity of the aerial measurement technology, these sources can represent a significant portion of the total inventory and must be considered.The present methodology calculates the unmeasured inventory by combining site counts with site-level emission factors that are estimated via a Monte Carlo (MC) simulation of aerial survey over sources near and below the aerial technology's sensitivity limit.This procedure is summarized in Figure 1b of the main text and detailed in this section.
Given the inherent lack of emission rate data for sources not detected aerially, the first requirement for this analysis is a representative feedstock dataset that provides a source rate distribution near and below the aerial technology's sensitivity limit.In the current work, highsensitivity measurement data of site/stratum-and source-resolved emissions were available from a 2018 ground survey in British Columbia of 267 facilities and wells (Robinson et al., 2018).
These data include quantified emissions from non-pneumatic equipment and abnormally operating pneumatic equipment detected by optical gas imaging and measured where possible using Hi-Flow sampling.Robinson et al. (2018) also counted and identified (manufacturer and model) pneumatic equipment and estimated expected emissions under normal operation based on prior field measurements and manufacturer-specific bleed rates.
To accurately infer the unmeasured inventory, the simulated aerial survey of the feedstock data must be performed similarly to the actual survey.Recalling that unmeasured sources are defined as never being detected during the potentially many flight passes over the source, it is imperative to simulate an appropriate number of flight passes and assume representative measurement conditions, which inform the probability of detecting any one source during a single pass.To this end, for the present study, an empirical probability mass function (PMF) of the number of flight passes over a source was obtained from the aerial survey data (see Figure S2a).Likewise, recognizing that the best-available continuous probability of detection (POD) function for Bridger's GML technology is sensitive to the estimated local wind speed ( ̃ [m/s]; at 3 m elevation) and aircraft altitude above ground level (AGL; ℎ ̃ [m]) (Conrad et al., 2022), empirical probability density functions (PDFs) for these parameters were similarly derived from flight passlevel data during the aerial survey and are shown in Figure S2b and c, respectively.For a single steady methane source, these distributions enable representative simulation of the probability that the source would be detected during a multi-pass aerial survey using Bridger's GML.

S2.2.1 Non-pneumatic Equipment
With the feedstock data, empirical PMF for number of flight passes (see Figure S2a), empirical PDFs for measurement conditions (see Figure S2b and c Non-pneumatic sources in the feedstock data were parsed into four facility/well categories prior to analysis using the presented method.Each stratum in the provincial inventory (refer to Table S1) was captured by one of these categories: off-site wells (OSW), single-well batteries (SWB; facilities with one linked well), multi-well batteries (MWB; facilities with multiple linked wells), and other (facilities with no associated wells).Strata in the "other" category had no feedstock data and were conservatively assumed to have zero unmeasured non-pneumatic emissions.Average site-level emission factors for the OSW, SWB, and MWB categories were computed using the described methodology and are summarized in Table S2.321, 322, 361, 362, 393, 401, 402, 403, 404, 405, 407, 451, 501, 601, 611, 676, and

S2.2.2 Pneumatic Equipment
Pneumatic instruments (e.g., level controllers, pressure controllers, transducers, positioners, etc.) and pumps are ubiquitous in process equipment buildings (e.g., separators, line heaters, compressor buildings, etc.) across production sites in the upstream oil and gas sector.In bottomup component-based inventories, pneumatic emissions are generally estimated by multiplying average device emission factors by the corresponding device counts aggerated over well and facility population strata (Clearstone Engineering Ltd., 2018;Robinson et al., 2018).In the present inventory, average site-level emission factors for pneumatic equipment were derived using the same method described for unmeasured non-pneumatics (see Section S2.2.1).For pneumatic equipment, however, this approach simplifies the method of Figure S3 because normally operating pneumatic devices generally vent at sufficiently low rates to preclude detection by aerial technologies; this is supported by the present aerial survey and follow-up ground survey (Johnson et al., 2022).Firstly, average vent rates for the pneumatic instruments and pumps in BC derived by (Robinson et al., 2018) using data from (AER, 2018;D'Antoni, 2018;Government of Alberta, 2020;Prasino Group, 2013;Western Climate Initiative, 2013) were less than 0.4 and 0.6 kg/h of methane, respectively, which correspond to a negligible POD even when using Bridger's highly sensitive GML technology.Indeed, at the median wind speed and typical aircraft altitude of the survey (4.5 m/s and 175 m AGL, respectively), the predicted single-pass PODs at these rates are less than 10 −5 .This low POD is evidenced by the present survey data itself, where methane was not detected using Bridger's GML at 603 active surveyed sites (177 facilities and 426 off-site wells) where pneumatic equipment would be expected.Moreover, follow-up ground survey of 195 aerially detected sources at 75 locations (Johnson et al., 2022) implicated pneumatic instruments and pumps as a potential contributor to just 24 sources.However, make and model data from pneumatic devices suggest that, if pneumatics were the cause of these 24 aerially detected sources, the identified pneumatics would need to be emitting approximately an order of magnitude greater than published/manufacturer-rated vent rates.However, this is also unlikely given the field data from the 2018 ground survey in BC (Robinson et al., 2018), where pneumatics deemed to be operating abnormally emitted at rates near the published/manufacturer-rated vent rates.These observations permit the reasonable assumption that normally operating pneumatics were not detected during the present aerial survey using Bridger Photonics GML.
Since probabilistic detection of normally operating pneumatic sources can justifiably be ignored, average site-level unmeasured pneumatic emission factors were derived using the same methodology as unmeasured non-pneumatic equipment (Section S2.2.1 and Figure S3).This procedure is thus greatly simplified as it does not require MC analysis of detection sensitivity; all pneumatic equipment contribute to the site-level emission factor, which is simply the average sitetotal vent rate of pneumatic equipment.In this case, the contribution of unmeasured pneumatic equipment to the provincial inventory reduces to the classic bottom-up approach that combines individual pneumatic counts and typical vent rates (emission factors) to yield site-level emissions.
It should be noted that this calculation is likely conservatively low, since the preceding discussion also suggests abnormally operating pneumatics (at least at the magnitudes observed in the field measurements of Robinson et al. (2018) S3; emission factors for well strata are on a perwellhead basis and applied to all surface wells in the province (see Table S1, footnote e) while emission factors for facility strata are on a per-site basis.Originally reported in units of m 3 /h of natural gas in Table 15 of Robinson et al. (2018), these were converted to kg/h of methane assuming a methane content of 88% by volume and a methane density of 0.6785 kg/m 3 at standard conditions (1 atm and 15°C).For gas well and gas facility strata, venting rates reported for conventional and tight gas production types were averaged and weighted by site type.Finally, a venting rate for mixed oil and gas batteriesmulti-well facilities that report both gas and oil productionwas estimated using the weighted average of all oil and gas multi-well batteries.As in Tyner and Johnson (2021) emissions from pneumatics at facility strata not covered by the 2018 OGI survey (e.g., gas plants, custom treaters, LNG Plant, custom treaters, etc.) were conservatively assumed to be zero under the assumption that pneumatics at the majority of these larger facilities are instrument air-driven.Likewise, methane emissions from pneumatic chemical injection pumps were estimated using a seasonal operation factor of 50%, an average vent rate of 0.973 m 3 gas/h (Clearstone Engineering Ltd., 2018)corresponding to 0.5807 kg/h of methane at the standard conditions noted above and again assuming a methane content of 88% by volumeand mean chemical injection pump counts for well and facility strata derived from the 2018 OGI survey, weighted by site type, for natural gas driven pumps (Robinson et al., 2018, Table 16).Consistent with pneumatic instrument emission factor estimates, the mean number of chemical injection pumps at mixed oil and gas batteries were estimated by averaging pump counts across all oil and gas multi-well battery strata weighted by site type.A seasonal operation factor of 50%, based on an assumed six months of operation per year consistent with the national inventory, discounts the contribution of pneumatic pumps to the unmeasured inventory as it inherently ignores pumps that may operate throughout the year.Importantly, this seasonal operation discount factor is likely to underestimate annual-averaged pneumatic pump emissions as some implementations are designed to operate more frequentlye.g., corrosion inhibitor pumps effectively operate constantly.surrogate of manual operations, two null hypothesis tests were deployed.This section describes the implementation of these tests and the results.
The first hypothesis test was performed on the subset of detected sources for which flights were performed on workday(s) and weekend(s).The difference(s) between flight-averaged emissions on the workday(s) and weekend(s) (Δ [kg/h]; positive if workday emissions exceeded weekend emissions and vice versa) were computed for each detected source.Available data for the parameter Δ were considered in aggregatethat is, across all sourcesand for nine unique source categories: compressors and compressor buildings, dehydrators, flares (lit and unlit), piping infrastructure, power generators, separators, tanks, other (e.g., amine sweetening, line heater, fuel gas, wellhead, etc.), and unknown (not identifiable).One-sample t-tests were performed for each set of Δ data with the null hypothesis that the mean workday-weekend difference (Δ ̅ [kg/h]) was zero (i.e.,   :  Δ = 0;   :  Δ ≠ 0).Results of the tests are summarized by source category in Table S10, which provides the size of the data set, the mean workday-weekend difference (Δ ̅ ), the t-statistic, and p-value of the hypothesis test.The present data did not justify rejecting the null hypothesis for any source category at 5% significance, implying no statistically significant difference between workday and weekend emissions.
The one-sample t-test only identifies statistical significance with respect to the mean and therefore does not necessarily capture effects at the tails of the distribution, where emissions from events like manual liquid unloading could be expected to manifest.Thus, an additional hypothesis test was performed to compare the source distributions between detected emissions on workdays and weekends.This was accomplished using the two-sample Kolmogorov-Smirnov (KS) test, which was applied to the same source categories as the t-test.The two-sample KS test compares the empirical cumulative distribution functions (eCDFs) for two datasetshere, quantified emission rates of all sources detected on workdays or weekendswith the null hypothesis that the

Figure S1 :
Figure S1: Geographic locations of 508 aerial survey sites in British Columbia, Canada overlaid on the identified locations of 1006 active facilities and 8995 active wells at the time of the survey.The approximately 46,000 km 2 red bounding box is the convex hull of areas of interest (polygons) measured during the present aerial survey.The inset map shows the aerial survey region within the province of BC and Canada.
using a quantification error model for Bridger's GML.Recently, Conrad et al. (2022) analyzed data from fully and semi-blinded controlled release experiments to derive a quantification error model for Bridger's GML technology.The error model defines the probability of the true source emission rate () given Bridger's estimate ( ̃) for each measurement pass of the airplane; the error model is in the form of a conditional probability distribution denoted by (| ̃).The inverse cumulative distribution function of (| ̃) provides an efficient means to randomly draw a true pass-level emission rate from Bridger's estimated value.From Table true 3-m wind speed to the joint distribution and subsequently marginalizing it out via integration over the positive real numbers.By assuming that the true source rate and wind speed are statistically independent and the uncertainty on aircraft altitude is negligible, the chain rule can be used to expand the joint distribution in Eq. (S4) to give (| ̃, ℎ ̃, ¬) ∝ ∫ (¬|, , ℎ ̃)(| ̃)  (| ̃) is a probabilistic error model for the 3-m wind speed (i.e., the conditional probability of the true 3-m wind speed given the Bridger-estimated 3-m wind speed) and   () represents a prior distribution for the true emission rate of the source.The conditional probability (¬|, , ℎ ̃) is the likelihood of a missed detection given a true emission rate, 3-m wind speed, and aircraft altitude; this is the complement of the probability of detection function () (| ̃, ℎ ̃, ¬) ∝ ∫ (1 − (, , ℎ ̃)) (| ̃)  (

Figure S2 :
Figure S2: Empirical probability mass (a) and density (b, c) functions for the number of flight passes per source (a), estimated 3-m wind speed (b;  ̃), and aircraft altitude (c;  ̃).

Figure S3 :
Figure S3: Monte Carlo analysis procedure for quantifying unmeasured sources (i.e., sources that are not detected during any pass of the aerial survey).
are likely missed by the aerial survey and are excluded in this bottom-up estimate.Measurement-based, average emission factors were derived by Robinson et al. (2018) from pneumatic counts of individual makes and models, collected during the 2018 ground-based OGI survey of production sites in BC, and combined with measured venting data from pneumatic studies in BC and Alberta(AER, 2018;D'Antoni, 2018; Government of Alberta, 2020; Prasino     Group, 2013; Western Climate Initiative, 2013).Strata-specific average emission factors for pneumatic equipment are summarized in Table underlying data come from the same distribution.Letting  ̂ () and  ̂ () represent the eCDFs at source rate , the test statistic () is simply the maximum difference between the eCDFs:  = max  | ̂ () −  ̂ ()| (S11)

Facility Description a Facility Type b Combined Type c Well/ Battery Type d Population Size Sample Size (% of Population) Entities With Sources (% of Sample)
a Facility description as per Petrinex database.b3-digitfacility "subtype" used to identify facility type/description inPetrinex.cCombinedstrata defined as a union of unique facility types.d