Advanced Leak Detection and Quantification of Methane Emissions Using sUAS

Abstract

Detecting and quantifying methane emissions is gaining an increasingly vital role in mitigating emissions for the oil and gas industry through early detection and repair and will aide our understanding of how emissions in natural ecosystems are playing a role in the global carbon cycle and its impact on the climate. Traditional methods of measuring and quantifying emissions utilize chamber methods, bagging individual equipment, or require the release of a tracer gas. Advanced leak detection techniques have been developed over the past few years, utilizing technologies, such as optical gas imaging, mobile surveyors equipped with sensitive cavity ring down spectroscopy (CRDS), and manned aircraft and satellite approaches. More recently, sUAS-based approaches have been developed to provide, in some ways, cheaper alternatives that also offer sensing advantages to traditional methods, including not being constrained to roadways and being able to access class G airspace (0–400 ft) where manned aviation cannot travel. This work looks at reviewing methods of quantifying methane emissions that can be, or are, carried out using small unmanned aircraft systems (sUAS) as well as traditional methods to provide a clear comparison for future practitioners. This includes the current limitations, capabilities, assumptions, and survey details. The suggested technique for LDAQ depends on the desired accuracy and is a function of the survey time and survey distance. Based on the complexity and precision, the most promising sUAS methods are the near-field Gaussian plume inversion (NGI) and the vertical flux plane (VFP), which have comparable accuracy to those found in conventional state-of-the-art methods.

1. Introduction

Why is methane so important? Methane is a greenhouse gas (GHG) that has a global warming potential 86 times that of carbon dioxide in a 20 year time window and is even larger for smaller time-scales. The mitigation of methane and reducing methane emissions can help reduce global warming in the near term. The first step is improving the way we measure emissions in practice, both in accuracy and in frequency. The overall measurement of methane emissions in oil and gas for example, (top-down vs bottom-up) has been shown to have discrepancies and is often underestimated [1,2].

For example, 190 oil and gas production sites were explored in [3,4], and the measurements indicated that well completion emissions were lower than previously estimated. The data also showed how emissions from pneumatic controllers and equipment leaks were higher than the Environmental Protection Agency (EPA) national emission projections. In a report titled, “Lessons from a decade of emissions gap assessments” [5], the authors argued about where we need to be and where we think we are, including the Paris climate agreement and what steps to take in order to keep global warming below 2 °C. One way to combat this is by detecting super emitters through tiered remote sensing strategies, which is outlined in [6].

This approach aims to focus on the detecting and repairing the largest emitters first, which can contribute a significant portion of the overall emissions of oil and gas systems. Furthermore, pressure to reduce contributions of climate change from customers and investors has been seen, insisting on reducing carbon footprints, including from landfills with much needed debate on inventory methods, direct emission measurements, and accountability [7]. The importance of mitigating methane emissions on a wide array of mitigation strategies is needed to stay on track with the Paris agreement [8].

Who is currently looking at methane? From an anthropogenic point of view (e.g., oil and gas), companies, such as Picarro (Santa Clara, CA, USA), Aerodyne Research (Billerica, MA, USA), Bridger Photonics (Bozeman, MT, USA), SeekOps (Austin, TX, USA), Heath Consultants (Houston, TX, USA), Flir (Global), Scientific Aviation (Boulder, CO, USA), Avitas (Houston, TX, USA), Ventus Geospatial (Houston, TX, USA), Aerometrix (Canada), and many more have provided methane detection and quantification solutions in a variety of technologies.

For example, A quantum cascade laser spectrometer is deployed on a small unmanned aircraft system (sUAS) for measuring facility-scale emissions using a mass balance approach with kriging [9]. For biogenic sources in ecosystems, there has been work looking at permafrost bogs [10,11], lakes [12], small ponds, wetlands [13], and vernal pools [14,15,16,17] to name a few. Seasonal dynamics of methane emissions from permafrost landscapes are explored in [18], specifically a lagoon pingo, and emissions estimated using a Thin Boundary Layer approach.

Porewater samples were analyzed using a Quantum Cascade Laser Spectrometer and combined with high resolution images from sUAS as an input into a neural network for creating a prediction map to upscale methane flux [19]. The spatial distribution of methane in the Artic permafrost bluffs was explored in [20] with a backscatter tunable diode laser absorption spectrometer (bs-TDLAS), namely the Pergam Laser Methane mini.

Flux estimates can be made, typically, using methods based on static measurements, on foot, by vehicle, manned aircraft, and by satellite. Static measurements consist of: (1) Eddy Covariance (EC) towers: A footprint modeling technique that looks at the turbulent exchange with the environment and utilizes meteorological conditions with precision concentration measurements to estimate the flux, (2) Chambers (autochambers): An enclosed chamber is placed over a target piece of land and is sampled with a syringe occasionally (to be analyzed at a later time, typically with gas chromatography) or dynamically sampled within a closed loop (such as GASMET’s Fourier Transform Infrared (FTIR) analyzer [21]).

Measurements made on foot consist of handheld sensors that are used with survey equipment and the surface. Surface emission monitoring (SEM) is typically a routine operation for landfills, done quarterly, to maintain compliance with local regulations to account for emissions lost from gas control systems. SEM is a point-based scanning technique that can take on the order of a week or so to complete. The concentrations are measured with devices, such as the flame ionization detector (FID) (regulated by EPA’s guidance and Method 21), and are integrated along sub divided grids looking for elevated methane levels (greater than 500 ppmv) [22].

Landfill emissions are generally calculated using inventory-type estimates. Measurements by vehicle consist of methods, such as the Tracer Correlation Method (TCM) and the environmental protection agency (EPA)’s other test method 33A (OTM33A) [23]. In [24], they explored vehicle-based advanced leak detection (ALD) with a cavity ring-down spectrometer (CRDS) from Picarro and determined that five to eight drives will capture a majority of leaks (>90%) as well as indicate detection limitations (such as wind and soil conditions and variations in methane enhancements making quantifying emissions difficult).

In [25], a vehicle mounted CRDS (Picarro G2301 and G4302) used empirical formulation to turn elevated concentration levels (or leak indications), C, to emission rates, Q (based on the work of [26], $lnC=-0.988+0.817lnQ$) and used Gaussian plume model to quantify emissions from site-level emissions in Utrecht and Hamburg, Germany. In work by [27], 6650 sites were evaluated using inventory and inverse point source Gaussian measurements, and they found that the methane inventory was underestimated by a factor of 1.5.

A series of campaigns were carried out utilizing TCM and downwind mobile measurements to explore the accuracy of different TCM approaches as well as compare CRDS with FTIR instrumentation in multiple source separation [28]. Measurements from manned aircraft have been done using FID, mounted on both fixed wing (Piper (Vero Beach, FL, USA) Seneca or Piper Navajo twin engine) and helicopter (Bell (Fort Worth, TX, USA) 206 Long Ranger) to detect liquid hydrocarbons from pipelines [29]. In [30,31], they utilize the next generation Airborne Visible/Infrared Imaging Spectrometer (AVIRIS-NG) to retrieve methane, carbon dioxide, and water vapor.

In [32], AVIRIS-NG was used to generate the VIsta-CA geospatial dataset to provide a comparison to the attribution of sources with Calfornia Air Resources Board (CARB) Pollution Mapping Tool (CARB PMT) and the U.S. Environmental Protection Agency (EPA) Facility Level Information on Greenhouse gases Tool (EPA FLIGHT). In [33], they investigate the uncertainty for estimating urban fluxes by an aircraft-based mass balance approach. They assess the sensitivity of the estimated city-wide CO${}_2$ and CH${}_4$ fluxes for several flight experiments, including the regional background concentration, depth of the convective boundary layer, magnitude of the wind speed, and type of interpolation technique.

In [34], they utilized a Sky Arrow Environmental Research Aircraft to measure emissions from multiple landfills and combined steady state Gaussian models to distinguish the emission coefficients for each individual site. A Bridger Photonics Gas Mapping LiDAR (GML) system was deployed on a Cessna 172 and blindly evaluated where the detection limits were as low as 1 kg/h depending on the wind conditions [35].

This method was also introduced to the Fugitive Emissions Abatement Simulation Toolkit (FEAST) [36] and shown to be comparable to OGI-based methods at equivalent survey frequencies for the detection and repair of emissions. Optical Gas Imaging (OGI) was explored in [37,38] and the effectiveness was evaluated in [38]. NASA’s Alpha Jet Atmospheric eXperiment (AJAX) and the AutoMObile greenhouse Gas (AMOG) surveyor were used to fuse airborne and ground-based data together (as part of the GOSAT-COMEX Experiment) using an anomaly approach instead of the typical mass balance approach [39]. Measurements by satellite have been explored in [40], where a ResNet-50 was trained from ESA’s Sentinel-2 data and labeled with a U-Net to detect smoke plumes.

Other works in the literature where emissions are detected, quantified, mapped, or localized include: A mid-wave infrared (MWIR) camera was used to compare eight supervised multivariate methods for detecting oil spills along the coastline in [41]. Using an array of stationary laser fetches, a controlled release emission is estimated using a Bayesian Markov chain Monte Carlo (MCMC) approach in [42]. There have been several works devoted to gas distribution mapping (GDM) using the Kernel DM/V methodology [43], including simultaneous localization and mapping (SLAM) [44].

GDM and gas source localization (GSL) with micro-drones have been explored in [45]. GDM has also been used in olfactory simulations in [46]. In [47], different GSL strategies (spiral, surge-cast, spiral-surge, and particle filter) were evaluated using the GADEN gas dispersion simulator. A mobile ground robot system named ARMEx was used to perform gas distribution mapping with a Heath Consultants remote methane leak detector (RMLD) [48].

In recent years, sUAS-based sensing approaches have become increasingly popular amongst practitioners for a variety of reasons, such as the ability to not be restricted to roadways or land locked areas, the ability to operate within the class G airspace at altitudes that traditional manned aircraft cannot operate at (improving resolution), the low cost, and the ability for high frequency deployment for capturing temporal changes.

Here, we provide an overview of some of the recent literature works utilizing sUAS, such as: a fixed-wing SIERRA sUAS with off axis integrated cavity output spectrometer (OA-ICOS) instrument was deployed in Svalbard, Norway prior to the NASA Characterization of Arctic Sea Ice Experiment (CASIE) [49]; single and multi-sUAS systems for source seeking based on the Luenberger observer were explored in [50]; an open path GHG analyzer based on vertical cavity surface emitting laser (VCSEL) was developed and tested in [51] with an aim to provide improved measurements compared to satellites; and volcanic emissions were captured using thermal cameras [52].

In [53], detection and spatial temporal analysis of a thermokarst lake was done with RGB images taken from a plane and sUAS. They looked at the bubble characteristics of the images to determine methane ebulliton. The use of long wave infrared (LWIR), short wave infrared (SWIR), hyperspectral, and visible cameras were used to detect liquid hydrocarbons with machine learning in [54]. Detection of methane gas from a custom open path absorption spectroscopy, mounted on fixed-wing sUAS, was explored in [55]. Emission factors from a combustion source using the EPA-based sensor, Kalibri, were calculated using a sUAS in [56]. Profiling GHG using sUAS-based AirCore system was analyzed with CRDS in [57].

Terra Sana Consultants developed a sUAS system with a path-integrated laser absorption (10 Hz at 30 m with 1 ppm-m) used in the detection of landfill gas. In a field trial, they compared the sUAS results to ground-based walk-over survey, reporting good correlation between the two [58]. A bs-TDLAS equipped drone with laser rangefinder was used to reconstruct 2D plumes under realistic conditions [59]. SEM, drone emission monitoring (DEM), and downwind plume emission monitoring (DWPEM) with CRDS are used with a genetic algorithm (GA) to estimate methane emissions from a landfill [60].

The AlphaSense electro-chemical sensor suite was used on a DJI 100 series sUAS that conducted ziz-zag and spiral localization flights of a stationary source [61]. In [62], red green blue (RGB), near infrared (NIR), and thermal infrared (TIR) cameras were used to map the topography and create digital elevation maps for identifying problematic areas where localized CH4 emissions were present using a static prototype semiconductor sensor. In [63], a sUAS equipped with a Pergam (Renton, WA, USA) backscatter-based tunable diode laser absorption spectrometer (bs-TDLAS) and OGI camera were used to detect and quantify pipeline leaks. The sUAS traveled 4 m from the pipeline during the surveys and had a minimum detection limit of 0.06 g/s.

In [64], atmospheric particulate matter and carbon dioxide were measured using sUAS sampling and a bag collection system. The bags were collected and analyzed in a lab. characterizing termite mounds using ground and sUAS-based laser scanning [65]. In a recent paper, ref. [66] utilized a NDIR instrument to measure CO${}_2$ flux (characterized, corrected, and validated in laboratory experiment at the Integrated Carbon Observation Station (ICOS) in Steinkimmen, Germany) at an ExxonMobil (Irving, TX, USA) natural gas processing facility in Germany. They used an on-board anemometer (FT-205) that was gain and bias corrected prior to the field experiments, where flux measurements were calculated using the mass balance approach.

General questions one can ask are, “What technologies and methods fall under advanced leak detection and quantification (LDAQ)?” Does this include mobile and sUAS-based approaches? These questions, in practice, are unfortunately up to the owners and operators of natural gas facilities, as they have the choice regarding what becomes adopted. However, the potential impact that LDAQ can have on improving methane mitigation is yet to be seen. Is leak grading quantification? What is accurate enough? What size leaks should we (or can we) care about? In the literature, it is often observed that leak quantification estimates and variability are reported in the place of accuracy and uncertainty. How can we determine the necessary and sufficient conditions for application of these methodologies? To the best of our knowledge, these questions remain unsolved in practice.

In the literature, there have been several reviews conducted on topics that deal with emissions, including remote sensing, source term estimation, and fugitive gas emissions. A remote sensing review paper in [67] describes many applications and topics within remote sensing, including the environmental sensing of volcanic eruptions, soil erosion, and geological related areas. A thorough review paper on source term estimation techniques is presented in [68]. In [69], a review was conducted on chemical sensing drones, which includes the sUAS platforms, sensors, and a brief overview of methodologies.

In [70], a review was conducted measuring fugitive gas emissions from landfills using various methodologies, including surface chambers (closed and open), EC towers, stationary mass balance, aerial mass balance, vertical radial plume mapping (VRPM), differential absorption LiDAR (DiAL), tracer gas dispersion (stationary and dynamic), and inverse modeling approaches (stationary, dynamic, and aerial [71]). In [72], several biogas plants in the UK were evaluated using the point source Gaussian plume model.

In this manuscript, we provide a review on the literature regarding the source rate estimation of continuous emission sources, focused around UAV-based methodology. We provide an overview of theoretical methodology as well as establish a quantitative comparison (for papers that have shown validated accuracy) between existing/current approaches and with UAV-based approaches in an attempt to shed light on the current accuracy of these methods. In Section 2, we discuss the problem overview. In Section 3, we overview some common sensors (chemical and wind). In Section 4, we overview the LDAQ methodology. In Section 5, we analyze the reported accuracy. Section 6 summarizes the methodologies, Section 7 discusses possible future directions, and in Section 8, we conclude the paper.

2. Problem Overview

The general problem in this work looks at methane emissions released into the air in gaseous form. The release mechanism or interface can vary depending on the application or system. For instance, in the oil and gas industry, leaks generally appear from tanks, valves, or hatches in the form of a point source, typically an above ground leak. Underground leaks also occur in practice, and the resulting emissions can manifest on the surface in many ways.

This may also be the case for landfills, where many small sources can be present across a very large area. If the distribution of sources is spatially uniform, we refer to this as an area source. The distribution of emissions may also vary, as in the case of natural ecosystems where the amount of methane may be produced at different rates depending on key factors of natural methane and carbon dioxide production (e.g., temperature, soil properties, water properties, etc.). Examples of different types of leaks are shown in Figure 1.

Landfills, which typically have area sources, are required to do quarterly walkover surveys, based on landfill regulations [73], using SEM on gas collection facilities. Landfills also have to consider the production and control of hydrogen sulfide gas, which is reviewed in [74]. Assuming some level of uniformity, chamber measurements have been taken and compared against atmospheric tracer methods (or TCM) [75]. The mass flux for chamber measurements can be calculated as

$$E=\frac{V}{A}p\left(\frac{\Delta C}{\Delta t}\right),$$

(1)

where V is the volume of the chamber, A is the area covered, p is the gas density given the headspace temperature, and $(\Delta C/\Delta t)$ is the change in mixing ratio, which is derived from linear regression of the temporal observations (four to five headspace measurements to achieve an acceptable correlation coefficient). However, due to the large size of these sites they are difficult to accurately measure. For example, four methods (aircraft mass balance, tracer correlation (TCM), vertical radial plume mapping (VRPM), and static chambers) and the California Landfill Methane Inventory Model (CALMIN) were compared in a landfill study in Indiana [76].

A field study comparison of different landfill methods for assessment of fugitive gas emissions was explored in [77]. This included VRPM, TDM, DiAL, MicroMeteoroligcal (Eddie Covariance method), and Flux chamber, and VRPM (close to the source, ≈10 m) and TDM (≈400 m) performed quite well against DiAL. In a paper from the UK, a review quantification of biogas plants was undertaken with inverse dispersion modeling (e.g., bLS), a tracer dispersion model, and OGI for different feedstock cases [72]. Mass balance approaches have also been applied using UAVs developed in [78].

In [79], chamber measurements were used to compare TIR images to quantify emissions in two landfills. The overall site emissions were verified using TCM (which tends to be the gold standard). The methane flux from different types of surface emissions were explored using chamber and FID measurements in [80]. A point-based scanning method, utilizing a portable gas detector (bs-TDLAS based), was correlated in lab testing using chamber methods and deployed on a landfill experimentally. This study showed a positive correlation between ambient methane concentrations and flux as well as directly proportional to flow rates [81]. Then, using this relation, a spatial map of the emissions was derived.

In a landfill study, the TCM was quantified for a 6 day campaign during different wind conditions and found that the methane emitted accounted for 31% of the generated methane [82]. Based on these findings, it is clear that fast and effective methods for estimating emissions from landfills are needed.

In natural ecosystems, which typically manifest as distributed sources (sometimes point sources distributed across a landscape) are generally much lower emissions than anthropogenic sources. It has been thought that these emissions are small; however, recent research has suggested that they are still not well understood. For example, digital elevation maps with SWIR imagery have been used to detect temporal trends in ombrotrophic peatland [83].

A Patagonia peat bog was examined with a UAV with high resolution color infrared (CIR). The images were classified using chamber measurements and different microforms in an attempt to upscale the plot-scale fluxes [84]. Thawing permafrost, peatland bogs, etc. have been typically measured using chambers [85], autochambers, and Eddy covariance (EC) towers. Commercialized chambers include those from LICOR, Picarro, and GASMET [21].

Aside from types of sources, wind is a direct input into the flux calculation and it can introduce a lot of uncertainty into the emission estimation. There are many important weather related measurements that can provide metrics for quantification methods, such as atmospheric stability. These stability classification schemes can depend on mechanical turbulence (roughness length and friction velocity), convective turbulence (mixing depth, Monin–Obukhov length, and heat flux), wind speed, and wind direction fluctuations [86,87]. These meteorological measures are summarized here: the Monin–Obukhov length is given as

$$L=-\frac{u_{*}^3\overline{T}}{\kappa g\overline{w^{\prime }{T}^{\prime }}},$$

(2)

where $\kappa$ is the von Karman constant, g is the acceleration of gravity, $\overline{T}$ is the average temperature, $\overline{w^{\prime }{T}^{\prime }}$ is the mean covariance between the vertical wind speed and sonic temperature [88]; and the friction velocity is given as

$$u_{*}=\sqrt{-\overline{u^{\prime }{w}^{\prime }}},$$

(3)

where $\overline{u^{\prime }{w}^{\prime }}$ is the mean covariance between the horizontal and vertical wind speed components. The effective plume height, $\overline{z}$, can be determined using the following two equations,

$$\begin{array}{c}\hfill L_{x,eff}+x_0=\\ \hfill \left\{\begin{array}{c}(\overline{z}/{\kappa }^2)[ln(c\overline{z}/z_0)-\mathrm{\Psi }(c\overline{z}/L)]{[1-p{a}_1\overline{z}/\left(4L\right)]}^{-1/2},L<0,\hfill \\ (\overline{z}/{\kappa }^2)[(ln(c\overline{z}/z_0)+2b_{2}p\overline{z}/\left(3L\right))(1+b_{1}p\overline{z}/\left(2L\right))+\hfill \\ (b_1/4-b_2/6)]p\overline{z}/L,L>0.\hfill \end{array}\right.\end{array}$$

(4)

This equation is first initialized by setting the effective distance $L_{x,eff}=0$, the effective plume height to the source height $\overline{z}=z_s$, and solving for the integration constant $x_0$. $L_{x,eff}$ is calculated from the longitudinal distance to the source using the angle to the center of the plume ${\theta }_p$ by, $L_{x,eff}=L_{x}cos(\theta -{\theta }_p)$, where $L_x$ is the longitudinal distance from the source. The stability parameter, $\mathrm{\Psi }$, which is dependent on the effective plume height $\overline{z}$ and Monin–Obukhov length L [89], can be calculated (for a given height) as,

$$\mathrm{\Psi }(z/L)=\left\{\begin{array}{c}{(1-a_{2}z/L)}^{1/4}-1,L<0,\hfill \\ -b_{2}z/L,L>0.\hfill \end{array}\right.$$

(5)

The coefficient c is dependent on the shape function parameter, s, described in [89] and given as

$$s=\left\{\begin{array}{c}\frac{1-a_{1}c\overline{z}/\left(2L\right)}{1-a_{1}c\overline{z}/L}+\frac{{(1-a_{2}c\overline{z}/L)}^{-1/4}}{ln(c\overline{z}/z_0)-\mathrm{\Psi }(c\overline{z}/L)},L<0,\hfill \\ \frac{1+b_{1}c\overline{z}/L}{1+b_{1}c\overline{z}/L}+\frac{1+b_{2}c\overline{z}/L}{ln(c\overline{z}/z_0)+\mathrm{\Psi }(c\overline{z}/L)},L>0.\hfill \end{array}\right.$$

(6)

The remaining coefficients (also from [89]) p, $a_1$, $b_1$, and $b_2$ can be set to 1.55, 16, 5, and 5, respectively (as used in [88]). The speed of the plume is given as

$$U\left(z\right)=\frac{u_{*}}{\kappa }[ln(z/z_0)-\mathrm{\Psi }(z/L)].$$

(7)

For the interested reader, the Monin–Obukhov similarity theory overview is given in [90].

3. Sensors and Equipment

There are many types of sensors that can be used on-board sUAS, granted they are light enough for a given platform’s payload capacity. In this section, we overview a few key sensors used and refer interested readers to a thorough review paper for more on chemical sensing drones [69]. Sensors used in many of the works reviewed here are briefly overviewed.

There are generally two types of sensing modalities, passive and active. Passive sensing encompasses any sensor that receives information from the environment. A common example of this includes optical cameras, such as visual spectrum cameras (e.g., RGB), thermal cameras (e.g., thermal infrared (TIR) [79,91], near infrared (NIR), short-wave infrared (SWIR), mid-wave infrared (MWIR), and long-wave infrared (LWIR)).

TIR cameras tend to span a larger bandwidth of wavelengths, whereas hyperspectral cameras can control which wavelengths to focus on. For example, in [92], Telops (Quebec, Canada) used a standoff tripod mounted hyperspectral camera to estimate the flow rate by integrating the mass per unit area and multiplying by the mean velocity of the gas. They utilized a two layer model to calculate the background radiance,

$$L_{tot}=[L_{bkg}{\tau }_{plume}+L_{plume}(1-{\tau }_{plume})]{\tau }_{atm}+L_{atm}(1-{\tau }_{atm}).$$

(8)

When the hyperspectral camera was optimized for methane detection (as in [93], 7.7 $\mathsf{\mu }$m band), two controlled release tests showed flow rates calculated (measured) were found to be $25.3\pm 2.8$ g/h ($23\pm 2.3$ g/h) and $102.9\pm 5.8$ g/h ($100\pm 10$ g/h). The authors claimed that the approach is 40 to 100 times more sensitive and can potentially be mounted to an aerial platform, remotely sensing from several hundred meters, deeming it suitable for both natural and anthropogenic sources. Active sensing encompasses any sensor that actively transmits information into the environment, probing a response.

An example of active sensing includes tunable diode laser absorption spectroscopy (TDLAS), which can come in several forms. The working principle relies on the gas species entering the sensor region or laser path, such that some of the power is absorbed by the gas, and a power drop is detected. One form of TDLAS is the closed path TDLAS (e.g., sensing region enclosed in controlled environment), where the emitter and detector are apart of the same device at a fixed distance apart, optimized for a desired detection species.

Other variations of the TDLAS include: the open path TDLAS (e.g., sensing region is open to the environment (see the open path laser spectrometer (OPLS) [94])), the backscatter TDLAS (bs-TDLAS) where the laser is reflected off the natural environment before being received at the detector, and long path TDLAS (e.g., used with retro-reflectors not connected to the physical instrument). Several examples of bs-TDLAS include: LiDAR based (Continuous wave laser absorption LiDAR, Pulsed Differential Absorption LiDAR (DIAL)), Pergam Laser Methane Mini [95], RMLD, Gasfinder2 [96], and Gasfinder3 [97]) [98].

Other more sensitive laser based instruments are also used in practice, such as cavity ring-down spectroscopy (CRDS) and off-axis integrated cavity output spectroscopy (ICOS). An example of these types of instruments can be seen from the Los Gatos Research Inc. (LGR) micro greenhouse gas analyzer (MGGA) (also referred to in the literature as the ultra-portable greenhouse gas analyzer (UGGA)) or the Picarro G2301 and G4302. These instruments are typically the gold standard for sensing gas, albeit they are also typically the most heavy as well. Other sensors used consist of non-dispersive infrared (NDIR), ceramic metal oxide sensors (CMOS) [99,100,101,102], photo-ionization detectors (PID, such as the Honeywell MiniRAE^® 3000) and electro-chemical sensors (review of applications [103].

A recent survey paper in [104] outlined new electronic nose technologies and applications. As the need for low cost sensing solutions increases and becomes commercially available, we are faced with evaluating the accuracy and characteristics of these sensors for practical use, such as in the work by [105] where low cost commercially available sensors were evaluated for precision and accuracy in a gas mixing chamber, providing promise for applications in continuous monitoring applications. Once these sensors can be evaluated and integrated on a platform with suitable sensor characteristics, they can be applied in practice.

For example, in [102], they demonstrated a proof of concept using a chemical multi-sensor payload for gas monitoring based on the DJI S900 platform. Or in [106], a semiconductor type sensor (Testo Gas Detector, Testo SE & Co. KGaA, Titisee-Neustadt, Germany) was used on a DJI M600 to analyze the spatial distribution of methane at a landfill, as well as compared different spatial interpolation techniques. This required calibration and consideration of the vehicle’s critical velocity.

For sUAS, lightweight and accurate wind sensors are needed to provide in situ measurements that can be used in the quantification methodology. Some examples of lightweight sensors can range from five hole probes (or multi-hole probes) to ultrasonic anemometers that utilize time of flight (e.g., Anemoment (Longmont, CO, USA) Trisonica used in [107] or the Gill (Lymington, Hampshire, UK) WindMaster used on the OPTOKopter [108]) and resonance based (FT Technologies FT742 and FT205 used in [109]) to more custom micro electrical mechanical systems (MEMS)-based solutions (such as in [110]). These wind sensors can also be applied to general wind profiling and mapping applications. For example, in [111], they used in situ wind measurements on sUAS for understanding the atmospheric boundary layer by developing wind profiling measurements using wind-induced perturbations. Mapping wind distributions over complex terrain was explored in [112], where they utilized a Gill WindMaster 3D ultrasonic anemometer mounted on a octocopter (called the WindLocater).

Wind and temperature profiling from fires were explored in [113] based on wind measurements using the Trisonica. Vertical velocity measurements of aerosol cloud interactions were compared with ground-based radar in [114]. If payload limitations prevent integration of an on-board wind sensor, wind estimation techniques can also be explored, such as in [115,116]. For a more thorough understanding of different sensors and estimation techniques, we refer the reader to [117].

The choice of a platform and ancillary equipment depends directly on the choice payload system (e.g., the collection of on-board sensors) that needs to be integrated onto the sUAS. For LDAQ, this typically includes lightweight methane and wind sensors. Due to sUAS payload capacity limitations, this often leads to integration problems as the weight of methane sensors can vary greatly depending on the desired sensitivity and response time (e.g., from a couple hundred grams to a couple kilograms). Lightweight and low cost methane sensors (such as CMOS) are slower in response and are less sensitive. For high accuracy and fast response sensors, which typically weigh more, this ultimately affects the selection process of sUAS, requiring bigger and more costly platforms to maintain safe stability and control.

For this reason, ground-based wind sensors are often used (placed adjacent to operating area) with smaller sUAS platforms gathering the in situ methane measurements. The data from these two systems are collected by a local data acquisition system for live observation and post-processing. This kind of configuration, implicitly assumes that the average wind, or mean wind field, represents the overall sensing region quite well. This is usually only the case in rural areas, where there are little to no obstructions (i.e., trees, hills, buildings, infrastructure, etc.). Most of the scenarios faced in practice, however, deal with obstructions and require on-board, or in situ, wind measurements.

This decision requires platform specific knowledge (e.g., hardware setup, autopilot, co-pilot software, etc.) as well as desired payload knowledge, which will vary depending on the desired application, sensitivity, measurement mode, and so on. For these reasons, and since this manuscript’s focus is on detection and quantification methodologies, we omit these details in this manuscript, and suggest that the interested reader see [118] for a guide on sUAS platform selection.

Payload integration strategies for methane sensors include several configurations, such as the boom-mounted, bottom-mounted, or top-mounted (see Figure 2). Boom-mounted approaches typically consist of TDLAS based sensors, which are subject to disturbances in the measurement from downwash of the propellers. This is avoided by placing the sensor out front of the aircraft along a boom and sampling when the effective wind speed over the sensor is greater than 2 m/s [119].

Bottom-mounted approaches are typically used with bs-TDLAS or OGI. The sensor is mounted on a gimbal system or sometimes hard mounted to the aircraft frame. Top-mounted approaches on typically only suitable for bs-TDLAS or OGI based methods. Point source measurements with TDLAS will, on average, underestimate the concentration (see [69] for more details). On the contrary, top-mounted wind sensors can provide high accuracy if translational and induced wind velocities can be removed [108].

4. Advanced Leak Detection and Quantification Methods

In this section, we overview the conventional and sUAS-based advanced leak detection and quantification (LDAQ) methodology. The LDAQ methods utilize several concepts and approaches within numerics, control, and optimization as well as approaches based on different available sensing modalities (see Figure 3 for effective length-scales). In this manuscript, we divided these approaches into five general categories, namely: Simulation-based (Section 4.1), Optimization-based (Section 4.2), Mass-Balance-based (Section 4.3), Imaging-based (Section 4.4), and Correlation-based (Section 4.5).

In the Simulation-based approaches, the methods depend heavily on simulation and computational tools for solving dynamic partial differential equations, which are used to determine the source rate estimation. Sometimes other source parameters are also estimated in the process and this is typically referred to as source term estimation (STE) or the source determination problem (SDP). Two methods that show up in the literature are backwards Lagrangian stochastic (bLS) and mesoscale recursive Bayesian least squares inverse (RB-LSI).

The optimization-based methods showcased in this manuscript depend on some form of a parameterized system model, which undergoes a model fitting or recursive optimization (statistical or information based). Many of these methods include several variations of the point source Gaussian (PSG) solution of the classical Gaussian plume model. This is seen in the PSG approach based on the EPA’s other test method (OTM) 33A, where the data is gathered from a single sensor downwind and undergoes model fit of the peak concentration measured.

Next is the conditionally sampled PSG (PSG-CS) approach that utilizes meteorological data in the model fitting process using conditionally sampled concentration data based on the incremental changes in wind direction. Another variation to this is the recursive Bayesian PSG (PSG-RB) that utilizes a moving sensor and meteorological data to condition the models likelihood function and prior for updating the posterior distribution that is used to quantify the source estimate. This approach also considers past knowledge about equipment characteristics if this is known.

A different approach to the Bayesian way of thinking is to solve for the parameters of the model conditioned on the observations. This approach also utilizes a particle filter and Markov Chain Monte Carlo (MCMC) to update the posterior and is referred to as the PSG sequential Bayesian MCMC (PSG-SBM). The last optimization approach mentioned in this manuscript is the Near-Field Gaussian Plume Inversion (NGI) approach.

The NGI utilizes fitting the Gaussian plume model based on sampling of a perpendicular plane downwind of the source. The vertical and horizontal dispersion relations are used to find the center of the plume within the perpendicular plane and minimize, by least square fit, the difference between the modeled concentration and the observed concentration (integrated over the lateral dispersion direction).

The next category is the Mass-Balance-based approaches, which includes methods that utilize equations based on mass conservation and continuity. The simplest approach is the vertical flux plane (VFP), which takes a control volume approach to estimating the emission rate by measuring the flux entering and leaving the control volume. Traditionally, the plume is sampled using a raster scanning approach in a perpendicular plane upwind and downwind of the source. The sparse set of observations within the plane undergo a spatial interpolation process and are combined with the wind to estimate the source rate.

A direct variation to this approach is the cylindrical flux plane (CFP), which the sensing system measures concentrations on successive loops around the source at different altitudes. The flux going into and out of this cylindrical plane is used to estimate the flux. Using different sensing modalities (such as imaging or backscatter-based sensors), a path integrated vertical flux plane (PI-VFP) method can be formulated. Both aircraft and sUAS-based approaches to PI-VFP have been implemented, which rely on horizontal scanning of the area of interest.

For sUAS-based PI-VFP, concentric circles are flown to confirm that sources are contained inside the path before estimating the source rate. A flux plane approach has also been explored using a series of TDLAS-based laser fetches at different altitudes and utilizes the time-average of the line-integral of the instantaneous product of the wind speed and concentration. This is advantageous to other VFP approaches as it provides very good performance and does not take time to scan the plane. However, it is in ways impractical as it requires setup of the laser fetches and knowledge of the source geometry.

The next method is the Gauss divergence theorem (GDT) approach. It utilizes the CFP approach with mass flux continuity as well as the expected time rate of change of the mass within the control volume to estimate the source rate. Another VFP approach was included in this review that uses Gaussian plume model optimization with a general linear model (GLM) to help determine the contributions of multiple sources. This approach is referred to as the VFP-GLM.

The last Mass-Balance-based approach is vertical radial plume mapping (VRPM). The VRPM approach uses a ground based laser with retro-reflectors at different altitudes downwind of the source. The path integrated concentrations are measured at different radial angles and used to estimate the flux.

The next category is the Imaging-based approaches that utilize MWIR, hyperspectral cameras, and absorption spectroscopy (such as iterative maximum a posteriori differential optical absorption spectroscopy (IMAP-DOAS)).

The last category covered in this manuscript is the correlation-based approaches, which includes the traditional Eddy covariance (EC) method (in brief) and the tracer correlation method (TCM). The TCM has also been referred to in the literature as the tracer dispersion method (TDM) and atmospheric tracer method (ATM).

4.1. Simulation-Based

4.1.1. Forward Modeling

Forward modeling is typically used for projecting or forecasting dispersion. Forward modeling is not directly used in emission quantification by itself, but rather paired with feedback in the optimization sense. This can include numerically solving a governing set of equations, such as the advection diffusion equation (ADE) or applying a parameterized general model (such as the Gaussian plume). It is also common in practice to utilize existing numerical models, such as the WindTrax 2.0, WRF model, FLEXible PARTicle-Weather Research and Forecasting (FLEXPART-WRF), SCIPUFF, QUIC, and others that can be Lagrangian-base, include turbulence e.g., Large eddie simulation (LES), and Reynolds averaged Navier Stokes (RANS). Interested readers can check the review paper from [122] on dispersion models.

4.1.2. Backward Lagrangian Stochastic (bLS)

The accepted backward modeling approach used in the draft OTM- 33A document [23] and in several applications (e.g., Dairy Farm [123], etc.) is the backwards Lagrangian stochastic (bLS) approach by [124]. The bLS approach aims to answer the general questions: What is the proper form of the LS trajectory model? As well as, how can source estimates be extracted from the particle’s backward LS trajectory? The forward model, formulated as a generalized Langevin equation, is evolved jointly as a Markov process,

$$d{u}_i=a_i(\mathbf{x},\mathbf{u},t)dt+b_{i,j}(\mathbf{x},\mathbf{u},t)d{\xi }_j,d{x}_i=u_{i}dt,$$

(9)

where the particle position is given by $\mathbf{x}=(x_1,x_2,x_3)=(x,y,z)$, and $d{\xi }_j$ is a random increment governed by Gaussian process. The functions $a_i$ and $b_{i,j}$ have to be specified such that the velocity probability density function, $g_a(\mathbf{x},\mathbf{u},t)$, satisfies the Fokker–Planck equation (FPE) [124],

$$\frac{\partial g_a}{\partial t}=\frac{\partial }{\partial x_i}(u_i,g_a)-\frac{\partial }{\partial u_i}\left[a_i(\mathbf{x},\mathbf{u},t)g_a\right]+\frac{\partial }{\partial x_i}\left[B_{i,j}(\mathbf{x},\mathbf{u},t)g_a\right].$$

(10)

This method provides a source estimation for an area source given the source location (with unknown source rate) and assuming horizontally uniform surface source atmosphere in horizontal equilibrium (see Figure 4). To make an emission estimate using bLS, the method utilizes the dispersion model relation,

$$\frac{UC}{Q}=n=f(z_m,z_0,L,h,G),$$

(11)

where L is the Monin–Obukhov length, h is the depth of the mixing layer, G describes the set of parameters characterizing the plume, and $z_m$ represents the measurement height. As the particles from the back trajectories touchdown in the source area, the vertical velocities, $w_0$ are logged and used to estimate n,

$$n\left(z_m\right)=\frac{C\left(z_m\right)U\left(z_m\right)}{Q}=\frac{1}{N}\sum |\frac{2}{w_0/U\left(z_m\right)}|.$$

(12)

Once n is known, an estimate of the source rate can be determined using the measured concentration and wind speed, $Q=n/\left(CU\right)$. In this approach, due to the time-averaged ensemble, the accuracy improves over time (nominal averaging period of 15 min [124]). An alternate expression for the emission estimate is given as

$$Q_{bLS}=\frac{C-C_b}{{(C/Q)}_{sim}},$$

(13)

where $C_b$ is the background concentration and ${(C/Q)}_{sim}$ is calculated using

$${(C/Q)}_{sim}=\frac{1}{N}\sum |\frac{2}{w_0}|.$$

(14)

This Monin–Obukhov similarity theory (MOST)-based bLS emission estimation methodology was validated against the mass balance approach (given the along-wind distance of the source D),

$$Q_m=\frac{1}{D}{\int }_0^{\infty }C\left(z\right)U\left(z\right)dz,$$

(15)

and field tested in [96,125,126]. A lagoon environmental leak was simulated and explored with the bLS approach by constructing a large 45 m by 45 m emission source on a pond. The accuracy was shown to be lower during the summer period due to more frequent unstable atmospheric conditions [127].

4.1.3. Mesoscale Recursive Bayesian Least Squares Inverse (RB-LSI)

Utilizing the NOAA P-3 aircraft and a wavelength-scanned CRDS, ref. [71] used a mesoscale bayesian least squares approach to solve the inverse problem of estimating emissions. They use the FLEXPART-WRF to model the forward problem, which was compared to physical observations and minimized on an iterative cost function that assumes lognormal distributions,

$$\begin{array}{c}\hfill J=\frac{1}{2}{(ln\left(y_0\right)-ln\left(Hx\right))}^T(ln\left(y_0\right)-ln\left(Hx\right))\\ \hfill +\frac{1}{2}\alpha {(ln\left(x\right)-ln\left(x_b\right))}^T(ln\left(x\right)-ln\left(x_b\right)),\end{array}$$

(16)

where the observed concentration enhancements are given as $y_0$, posterior solutions are x, FLEXPART-WRF outputs are H, prior fluxes are $x_b$, error covariance matrix from observations are R, and error covariance matrix from prior fluxes are B in the lognormal space.

4.2. Optimization-Based

In this section, we discuss the emission quantification techniques that utilize some form of optimization in the methodology that fits a model.

4.2.1. Point Source Gaussian (PSG)—OTM33A

In [88], the point source Gaussian (PSG) is discussed. The measurement involves a vehicle with a concentration measurement instrument (CMI) to park downwind of the known source with the vehicle off. The CMI (such as Picarro or LGR UGGA) collects data at roughly 2.5 m above ground at a known distance from the source. The variations in the wind direction are measured using a sonic anemometer (e.g., R.M. Young). The PSG calculations are based off of enhanced emission levels and can be calculated as the fifth percentile of the concentration time series signal [88]. The PSG estimate then becomes a simple 2-D Gaussian integration with no reflection term,

$$S_E=2\pi {\sigma }_y{\sigma }_z{U}_m{C}_p,$$

(17)

where $C_p$ is the peak concentration from the Gaussian fit, $U_m$ is the mean wind speed, ${\sigma }_z$ and ${\sigma }_y$ are the vertical and lateral plume dispersion that can be determined from the meteorological conditions, such as the Pasquill–Gifford stability classification curves [128] (see Figure 5). The accuracy of the OTM33A method is explored in [129,130].

4.2.2. Conditionally Sampled PSG (PSG-CS)

To capture the ensemble mean of the downwind plume behavior, a dispersion model is used in [88]. The model is a function of downwind distance and dispersion factors $D_y(x,y)$ and $D_z(x,z)$, given as

$$C_m(x,y,z)=\frac{S}{\overline{U}}{D}_y(x,y)D_z(x,z).$$

(18)

This method essentially aims to determine the source rate, S, using the conditional mean concentration data, $C_m$, of the downwind plume. The lateral dispersion downwind of a continuous point source can be shown to have a Gaussian distribution such that it can be represented as

$$D_y(x,y)=\frac{1}{\sqrt{2\pi }{\sigma }_y}[-\frac{1}{2}({\frac{y}{{\sigma }_y}}^2)].$$

(19)

However, the vertical dispersion (assuming vertical eddy diffusivity and wind speed that scales vertically to a power law) can be formulated as a parameterized stretched exponential (originally expressed in [131]),

$$D_z=D_z(x,z)=\frac{A}{\overline{z}}exp[-{(\frac{Bz}{\overline{z}})}^s].$$

(20)

The parameters $\overline{z}$, s, A, and B are functions of the atmospheric stability and downwind distance, x. A and B can be described using the usual Gamma function, $\mathrm{\Gamma }(·)$ as

$$A=s\mathrm{\Gamma }(2/s){\left[\mathrm{\Gamma }(1/s)\right]}^2,$$

(21)

$$B=s\mathrm{\Gamma }(2/s)\mathrm{\Gamma }(1/s).$$

(22)

The conditional averaged concentration can be calculated using

$$\langle C|\theta \rangle =\frac{1}{n}\sum _{{\theta }_i\in \mathrm{\Theta }}^{n}C\left({\theta }_i\right),$$

(23)

where the set $\mathrm{\Theta }\left(\theta \right)=\{{\theta }_i:|\theta -{\theta }_i|<\Delta \theta /2,\forall i=1,2,...,n\}$ and $\Delta \theta =2^{°}$. The basic idea is to capture the plume geometry in the crosswind direction, which is further used to derive the least squares source estimate,

$$S=\left[\sum _{i=1}^N\frac{D_y{D}_z}{\overline{U}}\langle C|{\widehat{Y}}_i\rangle \right]/\left[\sum _{i=1}^N{(\frac{D_y{D}_z}{\overline{U}})}^2\right].$$

(24)

As shown in [88], the lateral dispersion can be determined in two ways: classically, using atmospheric stability (for constants $a_y$ and $p_y$) [132],

$${\sigma }_y=a_y{z}_{0}1.9{(L_x/z_0)}^{p_y};$$

(25)

and by reconstructing the lateral dispersion,

$${\sigma }_y=\sqrt{\frac{1}{N}\sum _{i=1}^N{\widehat{\underline{Y}}}_i},$$

(26)

where the N is the number of values in $\langle C|\widehat{\underline{Y}}\rangle$, and $\widehat{\underline{Y}}$ are $\widehat{Y}$ values that are greater than the minimum concentration (i.e., background) and $\pm {40}^{°}$ off the plume center ${\theta }_p$. The distance $\widehat{Y}$ is calculated as

$$\widehat{Y}\left(\theta \right)=L_{x}sin(\theta -{\theta }_p),$$

(27)

with ${\theta }_p=arg{max}_{\theta }\langle C|\theta \rangle$ (see Figure 6).

4.2.3. Recursive Bayesian Point Source Gaussian Method (PSG-RB)

In work from [133,134], a moving sensor measured a point source concentration that can be formulated as

$$C(x,y,z)=\frac{S}{\overline{U}}{D}_y(x,y)D_z(x,z).$$

(28)

The source rate is given as S, the effective wind speed is $\overline{U}$, and the lateral and vertical dispersion is characterized by $D_y(x,y)$ and $D_z(x,z)$, respectfully. The equation is formulated such that the downwind distance, x, is aligned with the predominant wind direction. Since the measurement is taken at closer distances to the source, the lateral dispersion is taken as a random function such that

$${\int }_{-\infty }^{\infty }{D}_y(x,y)dy=1.$$

(29)

This can be advantageous for instantaneous plumes. The integrated lateral concentration can be written as

$$C^y(x,z)=\frac{S}{\overline{U}}{D}_z(x,z).$$

(30)

The choice of the vertical dispersion $D_z$ (originally expressed in [131]) is that of a parameterized stretched exponential function,

$$D_z=D_z(x,z)=\frac{A}{\overline{z}}exp[-{(\frac{Bz}{\overline{z}})}^s],$$

(31)

where $\overline{z}$, s, A, and B are functions of atmospheric stability and downwind distance, x. The lateral dispersion is given as

$$D_y=\frac{1}{\sqrt{2\pi }{\sigma }_y}exp[-\frac{1}{2}{(\frac{y_i}{{\sigma }_y})}^2].$$

(32)

Then, by numerically integrating (28) and incorporating the vehicle movement V,

$$C^y=\sum _{i=0}^{\infty }C(x_i,y_i,z_i)\Delta tV=S\sum _{i=0}^{\infty }\frac{\Delta tV}{{\overline{U}}_i}{D}_z(x_i,z_i)D_y(x_i,y_i).$$

(33)

The recursive Bayesian approach described here is based on well pads and oil and gas production, which are used to help inform the path planning of the mobile sensor. For brevity, we will only cover the formulation of the quantification only. Starting with the definition of the posterior distribution,

$$p(S|M,W,\mathrm{\Lambda })=\frac{p(S|W\left)p\right(M|S,\mathrm{\Lambda })}{p(M|\mathrm{\Lambda })},$$

(34)

where M is the concentration data, W is the ancillary information (e.g., well pad characteristics), $\mathrm{\Lambda }$ is the meteorological conditions, $p(S|W)$ is the prior, $p(M|S,\mathrm{\Lambda })$ is the likelihood, and $p(M|\mathrm{\Lambda })$ is the evidence (which can be thought of as a normalization constant for the likelihood [135]). The prior is given as

$$p(S|W)=\frac{1}{\beta }exp[-{(1+\gamma \frac{S-\mu }{\beta })}^{-\frac{1}{\gamma }}]{(1+\gamma \frac{S-\mu }{\beta })}^{-1-1/\gamma },$$

(35)

where the hyperparameters need to be fit to the application (for well-pad source, $\gamma =1$, $\mu =0.19$, $\beta =0.23$ based on [136]). The likelihood function is chosen to be a Gaussian,

$$p(M|S,\mathrm{\Lambda })=\frac{1}{\sqrt{2\pi }{\sigma }_e}exp[-\frac{1}{2}{(\frac{C^y-C^{y,M}}{{\sigma }_e})}^2],$$

(36)

where $C^{y,M}$ is the modeled concentration for a given source rate, and ${\sigma }_e$ is the combined model and measurement error (outlined in [137]). The recursive approach involves replacing the prior with the previous posterior distribution found using the likelihood function,

$$p{(S|W)}_i=\left\{\begin{array}{c}p(S|W),i=1,\hfill \\ p{(S|M,W,\mathrm{\Lambda })}_{i-1},i>1.\hfill \end{array}\right.$$

(37)

As the number of passes increases, the posterior distribution improves and can be used to estimate the source rate,

$$S={\int }_{S_{min}}^{S_{max}}Sp(S|M,W,\mathrm{\Lambda })dS.$$

(38)

Variations of this method were seen in [134], where the measurement noise was assumed to be Gaussian and also included a UAV with sensor noise and utilized the flux plane mass balance method to estimate the source rate, which was used in the calculation of the posterior distribution. Further field tests of this method were carried out in [138].

4.2.4. Point Source Gaussian Sequential Bayesian Markov Chain Monte Carlo (PSG-SBM)

Utilizing the Gaussian plume model for the likelihood of a sequential Bayesian Markov Chain Monte Carlo (MCMC) method, a UAV scans horizontally to update the estimated posterior distribution in [139]. The parameters are given as ${\mathrm{\Theta }}_k={[{\mathbf{p}}_s^T,q_s,u_s,{\varphi }_s,{\zeta }_s]}^T$, were the position is ${\mathbf{p}}_s$, source rate $q_s$, wind speed and direction $u_s$ and ${\varphi }_s$, and the model diffusion coefficients ${\zeta }_s={[{\zeta }_{s1},{\zeta }_{s2}]}^T$. The point source observations, ${\mathbf{z}}_{1:k}=\{z_1,z_2,...,z_k\}$ are used within Bayes rule to update the posterior,

$$p({\mathrm{\Theta }}_{k+1}|{\mathbf{z}}_{1:k+1})=\frac{p({\mathbf{z}}_{k+1}|{\mathrm{\Theta }}_{k+1})p({\mathrm{\Theta }}_{k+1}|{\mathbf{z}}_{1:k})}{p(z_{k+1}|{\mathbf{z}}_{1:k+1})}.$$

(39)

The likelihood model, $\mathcal{M}({\mathbf{p}}_k,{\zeta }_k)$, in [139], based on observational data, ${\overline{z}}_k=\mathcal{M}({\mathbf{p}}_k,{\zeta }_k)$$+{\overline{v}}_k$, was taken to be detection event, $p({\overline{z}}_k|{\mathrm{\Theta }}_k)$, if $z_k>z_{thr}$,

$$p({\overline{z}}_k|{\mathrm{\Theta }}_k)=\frac{1}{{\sigma }_k\sqrt{2\pi }}exp[-\frac{{({\overline{z}}_k-\mathcal{M}({\mathbf{p}}_k,{\zeta }_k))}^2}{2{\sigma }_k^2}],$$

(40)

and a non-detection event otherwise,

$$p({\underline{z}}_k|{\mathrm{\Theta }}_k)=(\frac{p_b}{2}[1+\mathrm{erf}(\frac{z_{thr}-{\mu }_b}{{\sigma }_b\sqrt{2}})])+p_m+(\frac{p_s}{2}[1+\mathrm{erf}(\frac{z_{thr}-({\mu }_b+\mathcal{M}({\mathbf{p}}_k,{\zeta }_k))}{{\sigma }_b\sqrt{2}})]).$$

(41)

The three terms in the non-detection event account for instrument noise, turbulence, and observing concentrations above the threshold, where $p_b+p_m+p_s=1$, and ${\mu }_b$ and ${\sigma }_b$ are mean background noise and standard deviation, respectively. Using a particle filter, the posterior can be approximated by a set of n weighted random samples ${\{{\mathrm{\Theta }}_k^{\left(i\right)},w_k^{\left(i\right)}\}}_{i=1}^n$,

$$p({\mathrm{\Theta }}_k|{\overline{z}}_{1:K})\approx \sum _{i=1}^n{w}_k^{\left(i\right)}\delta ({\mathrm{\Theta }}_k-{\mathrm{\Theta }}_k^{\left(i\right)}),$$

(42)

where $\delta$ is the Dirac delta function. The un-normalized weights are then updated using

$${\overline{w}}_{k+1}^{\left(i\right)}=w_k^{\left(i\right)}·p(z_{k+1}|{\mathrm{\Theta }}_{k+1}^{\left(i\right)}).$$

(43)

Once the weights are determined they can be normalized by dividing by the summation of all the weights. Additionally, an effective sample size must be considered to avoid the degeneracy problem. The new samples undergo a MCMC step that is accepted with the likelihood probability distribution described earlier (see Figure 7).

4.2.5. Near-Field Gaussian Plume Inversion (NGI)

The near-field Gaussian plume inversion (NGI) method [140,141] is a mass continuity model in principle, where the upwind and downwind concentration measurements, combined with wind measurements, of an emission source are differenced to quantify emission flux. The NGI method is typically sampled around 100 m from the source. The sampling aims to capture the time-invariant behavior of the plume, which, under turbulent conditions, may not map out the characteristic Gaussian plume shape.

This is because it is assumed that spatial variability in the time-averaged plume is Gaussian. This method was initially carried out with a DJI S900 equipped with a ultra portable greenhouse gas analyzer (UGGA) by Los Gatos Research Inc. (LGR). The flux estimate is derived by fitting the experimentally measured flux values, $q_{me}$, to the modeled flux values, $q_{m0}$ given as,

$$q_{me}=(C-C_b)U\left(z\right)\rho ,$$

(44)

where the modeled flux is given by the Gaussian model,

$$q_{m0}=\frac{F_e}{2\pi {\sigma }_y\left(x\right){\sigma }_z\left(x\right)}exp(\frac{-{(y-y_c)}^2}{2{\sigma }_y{(x)}^2})(exp(\frac{-{(z-h)}^2}{2{\sigma }_z{(x)}^2})+exp(\frac{-{(z+h)}^2}{2{\sigma }_z{(x)}^2})).$$

(45)

The lateral and vertical dispersion relations are typically looked up in the PGT stability tables, however, in this method, they are assumed to be linearly proportional to downwind distance,

$${\tau }_y={\sigma }_y\left(x\right)/x,{\tau }_z={\sigma }_z\left(x\right)/x.$$

(46)

Trying to solve (45) is not always well constrained, and thus the method proposes to separate (45) and fit the model along the z-direction,

$$q_{me,y}=q_{me}\frac{{\tau }_{z}x\sqrt{2\pi }}{(exp(\frac{-{(z-h)}^2}{2{\left({\tau }_{z}x\right)}^2})+exp(\frac{-{(z+h)}^2}{2{({\tau }_{z}x)}^2}))}.$$

(47)

The spatial variability in the z-direction has to be sampled to determine ${\tau }_z$. The lateral spatial variability ${\tau }_y$ and plume center $y_c$ are determined simultaneously,

$$y_c=\frac{{\sum }_j\left(q_{me,y}{y}_j\right)}{{\sum }_j\left(q_{me,y}\right)},$$

(48)

$${\tau }_y=\sqrt{\frac{{\sum }_j(q_{me,y_j}{(\frac{y_j-y_c}{x_j})}^2)}{{\sum }_j(q_{me,y_j})}}.$$

(49)

Once the unknowns, ${\tau }_z$, ${\tau }_y$, and $y_c$ are determined, the source emission rate, F, can be estimated by minimizing the least square fit between $q_{me}$ and $q_{mo}$, given $F_e$ and ${\tau }_z$. The uncertainty in F and the impact of limiting ${\tau }_z$ are given in [140].

4.3. Mass Balance Based

The mass balance approach aims to estimate an emission source by balancing the mass flux leaving or entering a control volume. Generally, there are two path planning approaches to the mass balance method: (1) rectangular vertical flux plane (or curtain) downwind of the source and (2) a cylindrical flux plane enclosing the source. For a well behaved plume under stable atmospheric conditions, the downwind plume contains all the flux. The sampling distance from the source may vary based on each submethod. The measured flux plane data can be sparse and is typically subject to spatial interpolation.

4.3.1. Vertical Flux Plane (VFP)

The flux plane method generally involves sampling within a plane, vertically or horizontally, upwind and downwind, of an emission source. It has been applied in several works [9,33,37,66,76,78,120,142,143,144,145,146,147,148,149]. The plane is typically sampled using a raster-scanning approach, capturing the plume within the width and height of the plane. The emission rate (in moles $s^{-1}$) can be estimated as,

$$Q_c={\int }_0^z{\int }_A^B{n}_{ij}(C-C_b)\mathbf{u}·{\mathbf{n}}_{f}dxdz,$$

(50)

where $n_{ij}$ is the mole density of air (given standard temperature and pressure), $(C-C_b)$ is the enhanced mole fraction (referenced to air), $C_b$ is the background mole fraction, $\mathbf{u}$ is the wind speed vector, and ${\mathbf{n}}_f$ is the flux plane normal vector (see Figure 8). Since the measurements are sparse, the integral irregularly spaced. To combat this, the sparsely sampled points are spatially interpolated using techniques, such as inverse distance weighting (IDW) [150] or kriging [151]. This is a common problem in geostatistics to interpret unknown data, $z\left(s_0\right)$, from desired spatial locations $s_0$ in domain $\mathrm{\Omega }\in {\mathbb{R}}^2$, only using N sparse sampling points, $z\left(s_i\right)$, based on some optimal weights, ${\lambda }_i$,

$$\widehat{z}\left(s_0\right)=\sum _{i=1}^N{\lambda }_{i}z\left(s_i\right).$$

(51)

For example, in ordinary kriging [151], a semivariogram is used to model the spatial variability and, given a spatial distance, h, is defined as,

$$\widehat{\gamma }\left(h\right)=\frac{1}{2N\left(h\right)}\sum _{i=1}^{N\left(h\right)}{(z\left(s_i\right)-z(s_i+h))}^2.$$

(52)

This experimental semivariogram can be fitted to the model semivariogram with one of several common functions: circular, spherical, exponential, Gaussian, or linear. The weights are determined by solving

$$\sum _{j=1}^N{\lambda }_{j}C({\mathbf{s}}_i-{\mathbf{s}}_j)+\mu \left({\mathbf{s}}_0\right)=C({\mathbf{s}}_i-{\mathbf{s}}_0),\mathrm{for}i=1,2,..N,$$

(53)

where $C(·)$, in this context, represents the point support covariance matrix. This matrix is related to the semivariogram, $\gamma \left(h\right)=C\left(0\right)-C\left(h\right)$ [151], and the mean square prediction error is ${\sigma }_e^2=Var(z\left(s_0\right)-\widehat{z}\left(s_0\right))$, which, for ordinary kriging, is minimized to make the estimated values $\widehat{z}\left(s_0\right)$ optimal. Furthermore, the estimator should be unbiased (e.g., $E\left[\widehat{z}\left(s_0\right)\right]=E\left[z\left(s_0\right)\right]$), which requires $\sum {\lambda }_i=1$ and the spatial mean to be stationary $E\left[z\right(s\left)\right]=\mu ,\forall s\in \mathrm{\Omega }$.

If the kriging process is not stationary, it is considered, at best, an approximate solution to the spatial interpolation problem and incorrect at worst. A better approach could be to apply a spectral method that takes into consideration non-stationarity and higher frequencies, namely, the high frequency kriging method [152]. Consideration of temporal observations could be included as well, see quantile kriging in [153].

An enhanced version of the IDW was proposed in [154] to include an adaptive distance-decay parameter based on the density characteristics of the sampled points. Available tools, such as Kriging Assistant (KA) [155], Golden Software Surfer, or ESRI Geostatistical Analyst for ArcMap have been used in the literature before. For irregular geographical units with different sizes and shapes, the interested reader should consult [156].

A variation of the VFP technique is illustrated in [157], where a path-averaged long open path duel-comb spectroscopy is operated from a ground vehicle to a sUAS with a retro-reflector. A vertical profile is flown downwind of the source to conduct the VFP. This technique is also vary similar to VRPM.

4.3.2. Cylindrical Flux Plane (CFP)

A variation to the VFP is the Cylindrical Flux Plane (CFP). This method has been used with manned aircraft as it is not as easy to raster-scan a rectangular flux plane. The methodology is essentially vary similar to the VFP and can be found in the work by [158], omitted here for brevity.

4.3.3. Path Integrated Vertical Flux Plane (PI-VFP)

A variation of the VFP is the path integrated vertical flux plane (PI-VFP). This method utilizes a bs-TDLAS approach in that the instrument points straight down and scans or circles the emission source (see Figure 9). In [159], the AVIRIS-NG manned aircraft used IMAP-DOAS technique to retrieve methane concentrations and estimated fluxes using a PI-VFP type calculation. This approach was compared with the GDT and Gaussian inverse approaches during a joint-flight campaign.

The emission rates were estimated by, $Q\approx \mathbf{u}·\mathbf{n}{\sum }_i{V}_i\Delta s_i$, were $V_i$ represents the vertically integrated concentration, and $\Delta s_i$ is a path segment along the boundary. The individual measurements are integrated together (referred to as integrated methane enhancement (IME)) such that IME $=k\sum XC{H}_4\left(i\right)·S\left(i\right)$. The value $XC{H}_4$ is the methane plumes that exceed the minimum threshold of 200 ppm/m and k is a conversion factor.

Using an RMLD sensor fitted to a small quadrotor UAV, a circular scanning approach can be applied to sample horizontally a site of interest. The sensor uses a bsTDLAS to measure integrated methane emissions from a known height. The resulting measurements are then combined with wind measurements to estimate the flux [121,160],

$$q={\int }_0^H{\int }_{-W/2}^{W/2}\mathbf{u}\times (X-X_b)dxdz,$$

(54)

where H and W are the vertical and lateral dimensions, and $X_b$ is the background concentration. This calculation encompasses a single circular loop and if the source is encapsulated, multiple passes can be used to estimate the source,

$$Q=\frac{1}{n}\sum q_i.$$

(55)

In practice, the circular flight path is actually made up of line segments that are box-like. The source location was also identified by course raster scanning over the area of interest followed by a more fine flight pattern free approach combined with triangular natural neighbor interpolation. The maximum observed concentration was used for the source location.

4.3.4. Micrometeorological Mass Difference (MMD)

Utilizing the technique from [161], sampling the plume far enough downwind of the source, the averaged MMD can be calculated as

$$Q=\overline{\int \int U_{(y,z)}({\rho }_{(y,z)}-{\rho }_b)dydz}=\int \chi \left(z\right)dz.$$

(56)

The work in [162] utilized the time-average of the line-integral of the instantaneous product of U and $\rho$ in the y-direction. Alternatively, while using a laser fetch, an instantaneous product of a single wind measurement U and line-averaged laser concentration was used,

$$\chi \approx \Delta y\overline{U_{\left(z\right)}({\rho }_{L\left(z\right)}-{\rho }_b)}.$$

(57)

This method can also be used to calculate the turbulent fluxes,

$$\frac{Q_{tur}}{Q}=\frac{(Q_{\overline{U}\overline{\rho }}-Q_{\overline{U\rho }})}{Q_{\overline{U\rho }}},$$

(58)

where $Q_{\overline{U\rho }}$ is calculated from the flux term in (57) and $Q_{\overline{U}\overline{\rho }}$ in (59),

$$\chi \approx \Delta y\overline{U_{\left(z\right)}}\overline{({\rho }_{L\left(z\right)}-{\rho }_b)}.$$

(59)

This prescription of the flux does not capture the turbulent component of the horizontal flux (albeit wrong), is often necessary due to the short time-scale behavior of the wind (e.g., limitations in wind measurement devices).

4.3.5. Gauss Divergence Theorem (GDT)

In the paper by [163], Conley et al. they focused on the continuity equation,

$$Q_c=〈\frac{\partial m}{\partial t}〉+\int \int \int \nabla ·c\mathbf{u}dV,$$

(60)

where m is the mass of the aerosol, $\langle ·\rangle$ is the expectation or average, $c=C+c^{\prime }$ is the concentration (comprised of an average term and a deviation term), $\mathbf{u}$ is the wind speed, and V is the volume of the area of interest. The flux divergence can be expanded as,

$$\nabla ·c\mathbf{u}=\mathbf{u}·\nabla c+c\nabla ·\mathbf{u}.$$

(61)

The surface integral is taken to be a cylinder, which can be broken into several parts: the floor, the walls of the cylinder, and the top. The height of the cylinder is taken such that the emission is encapsulated with the minimum and maximum height. The resulting emission rate can be calculated as

$$Q_c=〈\frac{\partial m}{\partial t}〉+{\int }_0^{z_{max}}\oint c^{\prime }{\mathbf{u}}_h·\widehat{\mathbf{n}}dldz,$$

(62)

where z represents the altitude, and l the flight path. The temporal trend of the total mass ($\frac{\partial m}{\partial t}$) within the volume can be estimated from the measurements. The cylinder passes can be vertically binned and discretely summed up,

$$Q_c=\frac{\Delta m}{\Delta t}+\sum _{z=0}^{z=Z_t}\left(\sum _0^L\rho ·u_n\right)·\Delta z.$$

(63)

4.3.6. Vertical Flux Planes with GLM (GLM-VFP)

In [34], a 3D grid of airborne measurements are collected across multiple landfill sites. The resulting downwind observational points are then spatially interpolated with IDW and used to calculate the total mass flux. The multiple steady state Gaussian dispersion models,

$$C(x,y,z)=\frac{Q}{2\pi {\sigma }_y{\sigma }_{z}U}exp(\frac{-y^2}{2{\sigma }_y^2})(\frac{1}{\sqrt{2\pi {\sigma }_z}})exp[\frac{-{(z-L)}^2}{2{\sigma }_z^2}],$$

(64)

are applied to a fixed grid (50 by 50 m), where the mixing ratios found over each individual landfill was used to calculate a model mass flux (for each site, integrated along the x, y, and z directions). The experimental measurements are then used with simulation measurements and a general linear model,

$$\underset{\alpha }{min}|MF-\sum _{i=1}^{max}(MM{F}_i·{\alpha }_i)|,$$

(65)

to approximate the emission coefficient, ${\alpha }_i$, from multiple landfill sources. The emission findings are further corroborated with a local Eddie covariance tower measurement.

4.3.7. Vertical Radial Plume Mapping (VRPM)

The vertical radial plume mapping approach (compared with other methods in [77]), utilizes a long path TDLAS instrument from the ground. The laser is aimed at retro-reflectors, situated perpendicular and downwind of the source. The height of the retro-reflector constitutes the different radial angles where the path-integrated concentrations are combined with the normal wind component to estimate the flux (similar to VFP or MMD). An illustration of this is seen in Figure 10.

4.4. Imaging-Based

In this section, we overview the imaging-based methodology for quantifying methane emissions. This typically includes techniques that sample images passively, such as TIR, MWIR, or other OGI-based instrumentation. The methods mentioned here that can quantify methane emissions are considered as quantitative optical gas imaging (QOGI).

4.4.1. Mid-Wave Infrared (MWIR) and Hyperspectral

In the work by [164], the detection limits of MWIR band of a hyperspectral data was explored using the Spatially-Enhanced Broadband Array Spectrograph System (SEBASS) airborne instrument. They also provided a comparison between LWIR and MWIR (see Figure 11) using the radiative transfer model,

$$R_s=(R_T^{\uparrow }+R_S^{\uparrow })+t\{{ϵ}_{s}B\left(T_s\right)+(1-{ϵ}_s)[\frac{R_T^{\downarrow }+R_S^{\downarrow }}{1-S(1-{ϵ}_s)}]\},$$

(66)

where $R_s$ is th etotal radiance at the sensor, $R_T^{\uparrow }$ is the upwelling emitted atmospheric path radiance, $R_T^{\downarrow }$ is the downwelling emitted atmospheric path radiance, $R_S^{\uparrow }$ is the scatter path radiance at the sensor, $R_S^{\downarrow }$ total solar radiance that reaches the surface, t is the atmospheric transmittance, ${ϵ}_s$ is the surface emissivity, and $B\left(T_s\right)$ is the blackbody radiation at the surface temperature.

Other works, such as [98], have used MWIR cameras combined with two Pergam Methane Mini G lasers in pipeline leak detection. In [165], a FLIR GF320 and a RMLD were used together to make volumetric flow rate calculations in the laboratory using a data fusion approach. In [166], they utilized a thermal camera and steady state energy balance approach to estimate methane emissions from thermal anomalies in urban landfills.

4.4.2. Iterative Maximum a Posteriori Differential Optical Absorption Spectroscopy (IMAP-DOAS)

The IMAP-DOAS method was applied to the AVIRIS-NG [30,31] aircraft and measures reflected solar radiation between 0.35 $\mathsf{\mu }$m and 2.5 $\mathsf{\mu }$m with 5 nm spectral resolution and sampling. Using a nonlinear iterative minimization of the differences between modeled and measured radiance. The measured concentrations can be applied to the PI-VFP method to calculate fluxes [159]. Variations in this approach for retrieving methane concentrations has been seen in [167] for albedo correction and [39] anomaly-based mass balance.

4.5. Correlation-Based

4.5.1. Tracer Correlation (TCM)

The tracer correlation method, or isolated source tracer ratio method, initially proposed and implemented in works by Lamb et al. [168] and Czepiel et al. [75], aims to quantify the emission rate of an unknown gas species by releasing a tracer gas at a known flow rate while measuring both the tracer and the unknown signals collocated downwind. This method assumes that the location of the source is known and, at the measurement location, the plume is well mixed. The elevated signal downwind also needs to typically be greater than 50 ppb. The authors report uncertainty estimates of $\pm 15\%$. The general equation is given as

$$Q_m=Q_t\frac{C_m}{C_t},$$

(67)

where $Q_t$ is the tracer release rate, and $C_m$ and $C_t$ are the elevated mixing ratios of the unknown source gas and tracer gas, respectively. A comparison study between TCM and other fugitive emission quantification methods are studied in [77]. The effect of wind on accuracy of the TCM was explored for landfills using WRF model [169]. An in situ method was used to evaluate the collection efficiency of gas extraction wells based on tracer gas [170].

Variations of the quantification of TCM were explored in [28], which quantified emission rates based on the plume integration of a transect, peak height of the transect using a scatter plot to calculate the ratio (best fit line), and comparison with fitted Gaussian plume model. A landfill field comparison of methane emission models were compared to measured emissions using TCM [171]. The TCM method was also applied to quantifying emissions from dairy farms in [123].

A dual tracer method was explored in [172]. The second tracer provides for closer downwind measurements that can be refined by assessment of plume position as well as in the far-field measurements the second tracer becomes an internal standard to the measurement. A mobile version of the TCM approach was proposed in [173].

4.5.2. Eddy Covariance (EC)

The Eddy covariance method aims to estimate the emission flux from a footprint area given the boundary layer meteorology. Historical developments and current implementations of this method are summarized in [174]. This method generally assumes stationarity of the measured data and fully developed turbulent conditions [175]. One way it can be expressed is,

$$Q=\frac{1}{t_f-t_i}{\int }_{t_i}^{t_f}(C\left(t\right)-\overline{C})(w\left(t\right)-\overline{w})dt,$$

(68)

where the time-averaged concentration and vertical wind speed is $\overline{C}$ and $\overline{w}$, respectively. There are several assumptions required to make this flux calculation.

5. Analysis of Methods and Assessment

In an attempt to analyze the methods covered in this paper, we decided to use the following metrics: required assumptions, sample distance, survey time, complexity, average precision, average accuracy, and average cost. The required assumptions are meant to inform the practitioner so that the best method can be applied to a given problem. For example, if the source location is unknown, the PSG method may not be directly applicable unless a source location estimate is supplied. The sample distance is defined as the distance from the source at which the required method needs measurements taken from.

The survey time consists of the time required to make a single flux estimate. Understandably, some methods may require multiple flux estimates in order to approximate the emission source to within an acceptable error. Complexity is the measure of how difficult it is to implement any given method. In order to determine a value for complexity, a scheme was developed using figures of merit (FOM) that assigns factors and weights to the metrics (detailed in

Table 1. Figures of merit for defining complexity of an estimation method.

FOM	(%)	Low (2.5)	Medium (5)	Med-High (7.5)	High (10)
Operator skill	30	Little	Moderate	Professional	Expert
Number of operators	25	1	2	3	3+
Equipment cost	15	<$10,000	<$50,000	<$100,000	>$100,000
Setup Time	20	<1 h	<4 h	<8 h	8+ h
Survey Time	10	<0.5 h	<1 h	<2 h	2+ h

). Determining the values for these factors were based on loose estimates, inferred from papers found in the literature.

Ranges were assigned to the metrics to capture variations in the factors due to either the operators or the equipment being used, and are given in

Table 2. Summary of methods and their assumptions, operational details, complexity, cost, average precision, and average accuracy are generalized over implementations found in the literature. The average precision and accuracy are given as unitless values, normalized by the true source rate (estimated source rate). (measured using ${}^1$ fixed, ${}^2$ foot, ${}^3$ vehicle, ${}^4$ manned aircraft, or ${}^5$ UAV; cost with $ = <$10,000, $$ = <$50,000, $$$ = <$100,000, $$$$ = >$100,000; (·) represents precision normalized on estimated source rate.).

Method	Assumptions	Sample Distance	Survey Time	Complexity (1–10)	Avg Precision	Avg Accuracy	Avg Cost
bLS	x${}_s$, horizontally uniform surface source atmosphere in horizontal equilibrium	20–441 m ${}^1$	15 min–4 days ${}^1$	3.8–8.5 ${}^1$	±0.16–0.36 ${}^1$ (±0.07–0.85) ${}^1$	±0.02–0.30 ${}^1$	$–$$$ ${}^1$
PSG	x${}_s$, steady state source rate, point source, plume evolution via ground-level Gaussian dispersion with no obstructions	441 m ${}^1$ 18–500 m ${}^3$ 50 m ${}^5$	4 days ${}^1$ 37–58 min ${}^3$ 7.22–20 min ${}^5$	4.5–8 ${}^1$ 3.6–6.1 ${}^3$ 3.3–5.5 ${}^5$	(±0.30) ${}^1$, ±0.20–0.67 ${}^3$ (±0.19–0.47) ${}^3$ ±0.31 ${}^5$	±0.0022–0.43 ${}^3$ ±0.50 ${}^5$	$$–$$$ ${}^3$
PSG-RB	PSG assumptions, vertical eddy diffusivity and wind speed approximated by power law scheme	20–200 m ${}^3$	6 min ${}^3$	3.6–6.4 ${}^3$	—	±0 ${}^3$	$$–$$$
PSG-CS	PSG assumptions, continuous source emission, constant wind speed, vertical eddy diffusivity and wind speed approximated by power law	18–106 m ${}^3$	20 min ${}^3$	3.6–5.9 ${}^3$	±0.20–0.67 ${}^3$ (±0.19–0.47) ${}^3$	±0.02–0.26 ${}^3$	$$–$$$ ${}^3$
NGI	constant source rate, ${\sigma }_x$ and ${\sigma }_y$ linear functions of distance to source, plume not capped by atmospheric temp. inversion in z-direction	50–82.25 m ${}^5$	7.35–29.62 min ${}^5$	3.3–5.5 ${}^5$	±0.21–0.58 ${}^5$ (±0.06–0.53) ${}^5$	±0.11–0.13 ${}^5$	$–$$ ${}^5$
MMD	x${}_s$	12–27 m ${}^1$	15 min ${}^1$	4.5–8.1 ${}^1$	±0.06 ${}^1$	±0.10 ${}^1$	$$$–$$$$ ${}^1$
GDT	near to no meandering, steady state source rate	3–8 km ${}^4$	1 hr ${}^4$	5.5–7.8 ${}^4$	±0.07 ${}^4$ (±0.08) ${}^4$	±0.13 ${}^4$	$$$–$$$$ ${}^4$
VFP	GDT assumptions	4.875–10 km ${}^4$ 19.08–510 m ${}^5$	1.5–4.5 h ${}^4$ 20–30 min ${}^5$	5.8–8.3 ${}^4$ 3.3–5.5 ${}^5$	(±0.30–0.53) ${}^4$ ±0.17–0.0.37 ${}^5$ (±0.013–0.62) ${}^5$	±0.10–0.50 ${}^4$ ±0.03–0.50 ${}^5$	$$$–$$$$ ${}^4$ $–$$ ${}^5$
PI-VFP	GDT assumptions	3 km ${}^4$ 0–6.77 m ${}^5$	20–30 min ${}^4$ 15–20 min ${}^5$	5.3–7.5 ${}^4$ 3.3–5.5 ${}^5$	(±0.34–0.58) ${}^4$ ±0.82 ${}^5$	±0.27–52 ${}^5$	$$$–$$$$ ${}^4$ $–$$ ${}^5$
CFP	GDT assumptions	3–17.84 km ${}^4$	1 h ${}^4$	5.5–7.8 ${}^4$	—	—	$$$–$$$$ ${}^4$
GLM-VFP	GDT assumptions	0.4–2.2 km ${}^4$	2.5 h ${}^4$	6–8.3 ${}^4$	(±0.21) ${}^4$	—	$$$–$$$ ${}^4$
VRPM	GDT assumptions	10–100 m${}^1$	1 h ${}^1$	7.1–8.6 ${}^1$	±0.18–0.21 ${}^1$ (±0.21–0.33) ${}^1$	±0.05–0.43 ${}^1$	$–$$$ ${}^1$
QOGI	temperature and pressure of gas at leak location are the same, gas plume length in direction of optical path is small	30 m ${}^1$	1 min ${}^1$	3–6.3 ${}^1$	±0.01–0.02 ${}^1$ (±0.02–0.03) ${}^1$	±0.20–0.24 ${}^1$	$–$$$ ${}^1$
TCM	leak plume and tracer plume are well mixed	100–3546 m ${}^3$	0.5–2 h ${}^3$	4.6–7.6 ${}^3$	±0.06–0.24 ${}^3$ (±0.06–0.74) ${}^3$	±0.0056–0.17 ${}^3$	$$–$$$ ${}^3$
EC	stationarity, fully developed turbulent conditions	25–228 m ${}^1$	—	4.5–7.8 ${}^1$	(±0.08) ${}^1$	—	$–$$$ ${}^1$

. For example, some setups may use more expensive equipment or more people for the same method and, as a result, are reflected in the complexity metric.

Evaluating methane quantification techniques is important, and much work has already gone into this topic through controlled release experiments and evaluation frameworks. Examples from controlled release facilities (CRF) consist of but are not limited to the following:

In the Joint Urban 2003 study [176,177], static sensors were distributed in an urban setting to measure the dispersion of tracer particulates. In [178], area-averaged velocity and turbulent kinetic energy profiles were derived from data collected at the Mock Urban Setting Test (MUST). Mock Urban Setting test (MUST) was also evaluated with photo-ionization detectors (PID) [179,180]. MUST was further simulated using MISKAM 6 [181]. In [182], the WRF model was used to model wind and turbulence inside the Quick Urban and Industrial Complex (QUIC) model for comparing simulated and observed plume transport. A test plan for Jack Rabbit II was developed in [183], which aimed to improve chemical hazard modeling, produce better planning for release incidents, improve emergency response, and improve mitigation measures.

More recently, single-blind tests at the Methane Emission Technology Evaluation Center (METEC) in Fort Collins, Colorado evaluated several types of LDAQ sensing modalities as apart of the Standford/EDF Mobile Monitoring Challenge (MMC) and the Advanced Research Projects Agency-Energy (ARPA-E) MONITOR program (such as by vehicle, plane, and drone—shown in Figure 12 and Figure 13). In the Standford EDF MMC it was observed that the drone based technologies performed quite well (e.g., SeekOps) with an $R^2=0.42$ [144].

While the results shown in Figure 12 seem quite promising, there is still exists some improvements in precision that can be made. In the ARPA-E MONITOR program, 6 of the 11 participants tested their technologies at the METEC facility in [184] against six other industry-based participants. Due to confidentiality agreements at the time of testing, the data gathered from the 12 participants were aggregated to compare the methodologies based on measurement type (handheld, mobile and continuous monitoring). However, to the best of authors knowledge, only four of the MONITOR program participants have published data regarding the METEC tests (shown in Figure 13).

In a white paper by Bridger Photonics, a sUAS-based approach using LiDAR-based sensor, also demonstrated promising results even though the uncertainty is not given. In [121], a RMLD was used on a sUAS with the PI-VFP method. In contrast, ref. [185] utilized a portable TDLAS-based instrument and the PSG method to quantify emissions. Lastly, ref. [186] used a dual frequency comb spectrometer (from over one kilometer away) with the non-zero minimum bootstrap method (see [187]) and the Gaussian plume model to estimate the source rate.

Examples from active operations with comparison to conventional OGI-based methods are conducted in the Alberta Methane Field Challenge (AMFC) [147,188,189], which aimed to answer the questions: Are Leak detection and repair (LDAR) programs effective at reducing methane emissions? As well as, Can new technologies provide more cost-effective leak detection compared to existing approaches?

In order to compare the performances of the each of the methods to one another, their performance metrics were garnered from different studies where the method was utilized in either a field study or a controlled release scenarios and recorded in Table 2. Performance values were gathered from the standard deviations of consecutive flux estimates of a singular source leak scenario. Accuracy pertains to error of the flux estimate to the known source rate.

This information was limited primarily to controlled release scenarios. For each method, performances and details were separated into the broad types of sampling strategies: fixed/ static, on foot, mounted on a vehicle, mounted on an aircraft, and mounted on a sUAS. This prevents convolution of performance values between, for example, long aircraft sampling flights at far distances and short sampling flights near the source via sUAS.

6. Summary of Methods

After analyzing the quantification methods, we can separate the methods based on whether they have used sUAS or not. In this manuscript, we observed that the sUAS-based methods consist of near-field Gaussian plume inverse (NGI), vertical flux plane (VFP), and the path-integrated vertical flux plane (PI-VFP). The non-sUAS-based methods consist of backwards Lagragian Stochastic (bLS), point source Gaussian (PSG), recursive Bayesian point source Gaussian (PSG-RB), conditionally sampled point source Gaussian (PSG-CS), micrometeorological mass difference (MMD), Gauss divergence theorem (GDT), VFP, PI-VFP, cylindrical flux plane (CFP), general linear model verticl flux plane (GLM-VFP), vertical radial plume mapping (VRPM), quantitative optical gas imaging (QOGI), tracer correlation method (TCM), and Eddy covariance (EC).

When comparing their performances in Table 2, it can be seen that, when categorizing by means of mobility (i.e., fixed, on-foot, etc.), methods using static sensors show a trend of having higher complexity values while UAV-based methods display generally lower complexity values. For a subset of the methods, the survey times, sample distances, and average accuracies can be seen in Figure 14. This subset was specifically displayed for these methods had both upper and lower bounds for survey times and sample distances along with accuracy data, which allowed for the plotting of these quantities for each method in the form of ellipses on a log–log plot.

When analyzing this plot, it can be seen that the sUAS-based methods are generally lower in sample distances and survey times as opposed to the aircraft-based method being the one of the highest in both. The bLS and TCM methods are shown to have the best average accuracy with several sUAS and mobile methods close in accuracy. The long sample times of bLS method are due to the values reported in [191], and it is possible that these values don’t reflect typical bLS sample times. The advantages and disadvantages of each of the methods can be seen in

Table 3. Summary of the method advantages and disadvantages along with fields of application.

Method	Application	Advantages	Disadvantages
bLS	General, Oil and Gas, Agriculture	Able to quantify area source and point source emissions	Sensor is fixed, multiple measurements, negatively impacted by obstacles
PSG	Biogas, Oil and Gas	No tracer required	Negatively impacted by obstacles and low wind speeds due to plume advection, mobile sensors limited to roads
PSG-RB	Oil and Gas	No tracer required, accurate source quantification in open environments	PSG limitations
PSG-CS	Oil and Gas	No tracer required	PSG limitations
NGI	General, Oil and Gas	Assumes near-field plume turbulence and wind meandering, No assumption of atmospheric stability class	Near-field ( 100 m), prop-wash interference if unaccounted for
MMD	Agriculture	Able to give instantaneous flux estimates	Fixed sensors
GDT	Regional	Capable of giving emission quantifications of large areas	Sample around closed volume typically large areas, unable to capture instantaneous methane flux
VFP	General, Oil and Gas, Landfill	Does not require exact source location, ease of mobility	Stable atmospheric conditions, unable to capture instantaneous methane flux
PI-VFP	General, Oil and Gas	VFP advantages	VFP limitations
CFP	General, Urban	Capable of giving emission quantifications of large areas	GDT limitations
GLM-VFP	Landfill	VFP advantages	VFP limitations
VRPM	Landfill	Able to give instantaneous flux estimates	Fixed sensors
QOGI	Biogas	Able to give instantaneous flux estimates	Gas velocity determined via gas camera–velocity component parallel to image plane, can be difficult to process images
TCM	General, Oil and Gas, Landfill, Urban	Does not rely on meteorological measurements	Mobile sensors limited to roads, application difficulty due to outside methane source interference
EC	General, Oil and Gas	Can capture emission variations due to long time series	Fixed sensor(s), stable atmospheric conditions, sensitive to time of day, typically requires long sampling times

along with what typical application fields that they were applied in.

The final ranking of the methods depends heavily on the desired application, which also depends on factors, such as sample distance, sample time, and desired accuracy. For that reason, it is difficult to rank the methods in general. Thus, we provide a ranking of the methods in terms of complexity (outlined in Table 1) with highlights from the precision and cost in Figure 15. The results indicate that the simplest methods, in terms of complexity, are the sUAS-based NGI${}^5$ and VFP${}^5$ as well as fixed QOGI${}^1$.

The most complex methods include bLS${}^1$ and manned aircraft-based approaches. In terms of precision, bLS${}^1$, NGI${}^5$, GDT${}^4$, VFP${}^5$, QOGI${}^1$, TCM${}^3$, and EC${}^1$ tend to be the best. Thus, for sUAS-based methods, NGI${}^5$ and VFP${}^5$, are the most promising approaches. Additionally, the GDT${}^4$, TCM${}^3$, and EC${}^1$ approaches can be treated as candidate methods for future implementation using sUAS.

7. Future Directions

What areas should we begin to focus and invest in, and where is the field going? One approach is by looking more into smart sensing using sensor arrays and machine learning [192]. Leveraging gas dispersion modeling in path planning with source estimation approaches. A recent paper in [193] showed a joint estimation method (wind and gas) that performed fairly well compared to existing methods at reconstructing plumes within enclosed structures.

Can these approaches be used with sUAS for smarter path planning to improve LDAQ methodology? As these methods develop for outdoor gas dispersion modeling, sUAS could potentially use these concepts for improved sensor placement, which can yield improved quantification results. For techniques, such as SEM of landfills as well as any 2D plume reconstruction problems, the sUAS-based complex tomography techniques outlined in [59] could be applied. Additionally, the combination of gas tomography with mass balance approaches or model-based optimizations could yield improved emissions estimates as well.

As survey areas become larger and harder to capture with single sUAS systems, the use of swarms can also be leveraged to improve VFP mass balance based approaches, such as in [194], or can be applied to the Bayesian inference framework, which aims to maximize the mutual information, such as in [195] to estimate source parameters. Currently regulatory hurdles and cost may prevent these systems from being applied in practice today. Given that low cost methane sensors are being actively researched, it is possible that swarms of sUAS may be used in the near future.

Considering that some non-sUAS methods (overviewed in this manuscript) may be adopted by sUAS that require longer sampling times (such as PSG and TCM), the use of power-over-tether may become desirable for increasing the survey time of sUAS. This is demonstrated for meteorological applications in [196] and has nearly indefinite flight times. There has also been advances in the digital transformation of technological applications and control, where concepts, such as Digital twins are being applied to perform smart control engineering or industrial artificial intelligence (IAI).

These techniques utilize modeling, machine learning, edge computing, and or internet of things (IOT) approaches to create digital representations of physical assets that evolve system parameters over time. They can be used to understand remaining useful life (RUL) of equipment and perhaps be coupled with LDAQ to understand when to survey equipment that is projected to fail in the near-term. This can allow the limited resources of companies and practitioners to focus on problematic areas in an attempt to detect and mitigate super emitters.

In a previous work, [120], plume modeling was applied to not only improve on methodology but also aide in smarter path planning (as mentioned above). As low cost methane sensors become more sensitive, they can be integrated into existing infrastructures to give system status updates that can be fed back into modeling approaches, such as with PSG-RB method’s a priori well information [133] and PSG-SBM [139]. Furthermore, having access to a priori information and digital twins of the plume (e.g., a model of the system), can allow for improved autonomy of the sUAS. Ultimately, by applying sUAS in this context, early detection, and repair of methane leaks can be better approached.

8. Conclusions

Overall, this manuscript serves to capture the majority of sUAS-related emission quantification strategies as well as provide some accuracy comparisons to more conventional and non-sUAS quantification strategies. LDAQ methods based on sUAS can provide accuracy close to the state-of-the-art conventional methods, while improving the sampling distances and sampling times (see Figure 14). The advantage to using sUAS in some cases allows for better localization of emission sources and provides more flexibility in deployment.

Taking into consideration the operator skill, number of operators, equipment cost, setup time, and survey times, the complexity of the methods were derived. The complexity ranking of the methods indicated that NGI${}^5$, VFP (sUAS-based), and QOGI${}^1$ have the simplest complexity, while bLS${}^1$ and the manned aircraft-based approaches have the highest. Comparing the precision of each method indicated that bLS${}^1$, NGI${}^5$, GDT${}^4$, VFP${}^5$, QOGI${}^1$, TCM${}^3$, and EC${}^1$ have the most precise estimations, while VFP${}^4$ and PI-VFP${}^4$ are the least precise. To conclude, for sUAS-based methods, NGI${}^5$ and VFP${}^5$ are the most promising approaches.

Additionally, the GDT${}^4$, TCM${}^3$, and EC${}^1$ approaches can be treated as candidate methods for future implementation using sUAS. Lastly, sUAS-based quantification approaches, outlined in this manuscript, can be combined with new modeling and control approaches and a priori information (e.g., digital twins, machine learning, or the joint estimation method) to improve autonomy and estimation. For interested readers, the papers and bib file can be made available upon request.