A new aerial approach for quantifying and attributing methane emissions: implementation and validation

Methane (CH₄) is a powerful greenhouse gas that is produced by a diverse set of natural and anthropogenic emission sources. Biogenic methane sources generally involve anaerobic decay processes such as those occurring in wetlands, melting permafrost, or the digestion of organic matter in the guts of ruminant animals. Thermogenic CH₄ sources originate from the breakdown of organic material at high temperatures and pressure within the Earth's crust, a process which also produces more complex trace hydrocarbons such as ethane (C₂H₆). Here, we present the development and deployment of an uncrewed aerial system (UAS) that employs a fast (1 Hz) and sensitive (1–0.5 ppb s⁻¹) CH₄ and C₂H₆ sensor and ultrasonic anemometer. The UAS platform is a vertical-takeoff, hexarotor drone (DJI Matrice 600 Pro, M600P) capable of vertical profiling to 120 m altitude and plume sampling across scales up to 1 km. Simultaneous measurements of CH₄ and C₂H₆ concentrations, vector winds, and positional data allow for source classification (biogenic versus thermogenic), differentiation, and emission rates without the need for modeling or a priori assumptions about winds, vertical mixing, or other environmental conditions. The system has been used for direct quantification of methane point sources, such as orphan wells, and distributed emitters, such as landfills and wastewater treatment facilities. With detectable source rates as low as 0.04 and up to ∼ 1500 kg h⁻¹, this UAS offers a direct and repeatable method of horizontal and vertical profiling of emission plumes at scales that are complementary to regional aerial surveys and localized ground-based monitoring.

1 Introduction

Methane (CH₄) is the second-largest contributor to anthropogenic radiative forcing, with a global warming potential (GWP-20) greater than 80 × that of carbon dioxide (CO₂) (Solomon et al., 2007; Schneising et al., 2020; Holmes et al., 2013). Methane is produced from many natural and anthropogenic sources which can be further classified as biogenic or thermogenic in origin. The largest biogenic sources result from anaerobic decay such as wetlands, melting permafrost, or the breakdown of organic matter in the guts of ruminant animals. Thermogenic CH₄, “natural gas”, is generated during the breakdown of organic matter at high temperatures and pressure within the Earth's crust. The later process also produces more complex trace hydrocarbons such as ethane and propane, collectively known as natural gas liquids (NGLs). Coal mining, natural seeps, oil and natural gas (O&NG) activities, and supply chain leakage each produce varying amounts of NGLs. While sometimes other gases predominate (e.g., CO₂), all thermogenic sources contain some fraction ethane (C₂H₆) (Kutcherov and Krayushkin, 2010; Glasby, 2006; Etiope and Sherwood Lollar, 2013), and it is the second-most prevalent NGL in natural gas processed for energy generation (Hodnebrog et al., 2018; Solomon et al., 2007; Simpson et al., 2012; Peischl et al., 2018; Karion et al., 2015; Johnson et al., 2019). Biogenic methane emissions do not contain ethane (Masson-Delmotte et al., 2023; Cheewaphongphan et al., 2019), so C₂H₆ is therefore a critical marker for source attribution, particularly when distinguishing between biogenic and thermogenic methane emissions (Solomon et al., 2007; Karion et al., 2015; Johnson et al., 2019). Accurately quantifying contributions from various sources, both spatially and temporally, is an important step towards building local, regional, and global CH₄ emission estimates as well as informing policy decisions.

The accuracy of regional and global methane emission estimates is limited due to the sheer number of unknown variables. There are two primary methods to estimate methane emissions: statistical analysis of source activity inventories, known as “bottom-up” (BU) estimations, and “top-down” (TD) estimates based on observations from instruments deployed on satellite, aircraft, and ground-based platforms (Heath et al., 2015; Schneising et al., 2020; Vaughn et al., 2018; Cheewaphongphan et al., 2019; Frankenberg et al., 2016). A variety of BU and TD methods have been used to detect methane leaks (Frankenberg et al., 2005) and inform policy decisions (Heath et al., 2015), but persistent discrepancies remain between the BU and TD methodologies (Masson-Delmotte et al., 2023; Solomon et al., 2007; Vaughn et al., 2018; Peischl et al., 2018; Cheewaphongphan et al., 2019). In general, BU and TD estimation methods only represent a snapshot in time and do not effectively account for temporal variations due to diurnal, seasonal, and activity cycles. Recent studies have shown that O&NG methane emissions vary significantly in time (Lavoie et al., 2015; Johnson et al., 2019; Vaughn et al., 2018), showing that more frequent and repeated measurements are important in helping to reduce uncertainties in emission monitoring (Frankenberg et al., 2005; Space Studies Board et al., 2019).

High-altitude and space-based systems offer robust methods of monitoring CH₄ and other greenhouse gases (GHGs) at regional and global scales (Sherwin et al., 2024; Schneising et al., 2020). However, most large-area satellites (Sentinel, Landsat) drop off in detection at ∼ 1000 kg h⁻¹ with only the very sophisticated or targeted systems able to quantify sources < 500 kg h⁻¹ (Sherwin et al., 2024). These instruments for TD emission estimates are ideally supplemented with contemporaneous ground-based measurements to constrain the monitoring capabilities – e.g., Total Carbon Column Observing Network, TCCON (Vaughn et al., 2018; Pétron et al., 2020; Kort et al., 2014; Parker et al., 2011; Turner et al., 2015; Gisi et al., 2012).

A common method for estimating GHG emissions is a mass balance approach, which uses GHG concentration measurements upwind and downwind of target sources to isolate the source from background concentrations (Frankenberg et al., 2016; Schwietzke et al., 2017; Johnson et al., 2019). These measurements can be collected from aircraft instruments along upwind and downwind transects (Frankenberg et al., 2016; Schwietzke et al., 2017) or by simultaneous measurements using similar ground-based instruments installed upwind and downwind of the sources (Gisi et al., 2012; Saad et al., 2014; Heerah et al., 2021). The accuracy of the latter method depends on relatively constant wind speed and direction during data collection. Therefore, accurate knowledge of principal wind conditions is required for proper instrument installation, and only those data obtained under favorable conditions can be used to accurately estimate the flux (Gisi et al., 2012; Lavoie et al., 2017). Aircraft quantification methods generally drop off around 10 kg h⁻¹ for the highly sensitive instruments, so it is difficult or impossible to quantify low-emitting sources using combined ground, aircraft, and satellite data to constrain and validate the estimates. Ground-based detection systems are the most straightforward and accessible methods, but effective site monitoring is highly dependent on wind directions. Column-averaged measurement techniques using Fourier transform spectrometers (FTSs) also require boundary layer height estimations and are unable to provide vertical profiles of source plumes.

There is a need for direct, repeatable, and cost-effective methods for detecting and quantifying CH₄ emissions from relatively small sources (< 1 kg h⁻¹), which do not require a priori assumptions. Small uncrewed aerial vehicles (UAVs) offer new approaches to airborne air pollution and emission monitoring over scales and locations which are difficult to detect or access with other regional monitoring systems (Chen et al., 2024; Villa et al., 2016; Burgués and Marco, 2020). While designs can vary dramatically between models, UAVs are either fixed- or rotary-wing platforms. Fixed-wing UAVs are typically able to cover larger areas and generally allow more options for sensor mounting configurations, but they are unable to hover and quickly adapt to environmental conditions while tracking emission plumes. Rotary wing platforms – also known as vertical takeoff and landing (VTOL) or multirotor UAVs – generally have lower required operating velocities and have the ability to hover and can therefore be used for more complex, discontinuous missions at higher spatial resolution (McKinney et al., 2019; Villa et al., 2016; Burgués and Marco, 2020). Multirotor UAVs usually have between four and eight individual propellers to generate thrust and the power consumption results in shorter flight times than their fixed-wing counterparts, but they do not require specialized equipment or runways for takeoffs and landings. Recent advances in control technology have made multirotor systems easier to reliably operate.

A major issue with fixed-wing systems is that the higher operating velocities and minimum height requirements result in low detection probabilities during site surveys (Barchyn et al., 2019). Rotary UAV systems do not require a minimum velocity to stay aloft and can therefore be outfitted with a wider range of equipment for physical and chemical sensing (Hollenbeck et al., 2018; Shah et al., 2019; McKinney et al., 2019; Villa et al., 2016; Burgués and Marco, 2020). For instance, onboard anemometers for in situ wind speed and direction have been shown to be more accurate onboard VTOL platforms, while higher relative winds and aerodynamic flows around fixed-wing platforms often result in less accurate in situ measurements (Hollenbeck et al., 2018). While multirotor UAVs have a significant propeller wash effect below the body, induced winds are negligible when anemometers are mounted above the shallow inflow layer (Barbieri et al., 2019; Hollenbeck et al., 2018; Villa et al., 2016; Barchyn et al., 2019).

The fast development of commercial UAVs and low-weight sensors has driven a multitude of scientific studies which would have been difficult or impossible to conduct over a decade ago. Both inverse modeling and mass-balance approaches have been used to calculate total emissions across a wide range of spatial scales. McKinney et al. (2019) deployed a hexarotor UAV outfitted with adsorbent cartridges to collect biogenic volatile organic compound (VOC) emissions at various locations in central Amazonia. Shah et al. (2019) calculated 3D flux densities by fitting direct measurements of downwind CH₄ plumes to near-field Gaussian plume models. Olaguer et al. (2022) used contemporaneous UAV and mobile ground-based measurements to estimate biogenic CH₄ emissions from a landfill. Bel Hadj Ali et al. (2020) compared multiple ground-based emission monitoring techniques with downwind plume measurements fitted to Gaussian plume models. Gålfalk et al. (2021) piloted a small quadcopter to fully surround a known biogenic methane hotspot for mass balance emission estimates < 200 kg h⁻¹.

This study describes the development and implementation of a new uncrewed aerial system (UAS) to address critical areas of research by (a) accurately and directly quantifying emissions on small spatial and temporal scales, especially in hard to reach places; (b) constraining BU inventories driven by small-scale estimations; and (c) separating and attributing thermogenic and biogenic sources. The New Mexico Tech (NMT) UAS design combines rapid sampling of chemical and meteorological data with a mobile platform capable of vertical and horizontal profiling relative to target sources for direct and high-resolution sampling of emission plumes. A commercially available UAV is outfitted with a multi-sensor onboard payload including a compact mid-IR spectrum analyzer and lightweight anemometer for efficient, repeatable quantification and characterization of various localized anthropogenic sources.

2 Methods

2.1 System design

The full uncrewed aerial system (UAS) is depicted in Fig. 1. It includes four main components: (1) a mobile mid-IR methane and ethane sensor (Aeris MIRA Pico), (2) a lightweight 3D vector wind and environmental sensor (Anemoment TriSonica Mini), (3) an uncrewed aerial vehicle (DJI Matrice 600 Pro, M600P), and (4) an onboard computer for system monitoring and data collection.

Figure 1 Uncrewed aerial system (UAS) design overview. (a) System block diagram: (A) TriSonica Mini weather sensor (TWS), (B) gas sampler inlet, (C) Matrice 600 Pro uncrewed aerial vehicle (M600P), (D) Raspberry Pi 4 onboard computer, and (E) Aeris MIRA Pico methane and ethane sensor. (b) Flight-ready UAS collecting ground data prior to takeoff. TWS and sampler inlet (A and B) mounted on a carbon fiber tube mast above the M600P's propeller inflow layer (“rotor wash”).

2.1.1 Methane and ethane sensor

The MIRA (Mid-InfraRed Analyzer) Pico leak detection system (LDS) developed by Aeris Technologies employs a solid-state laser and multi-pass absorption cell with a spectral band pass between 2.5 and 4.7 µm. The wide spectral range allows for simultaneous measurements of both CH₄ and C₂H₆ mixing ratios with a precision of 1 ppb s⁻¹ and 0.5 ppb s⁻¹, respectively (Aeris Technologies Inc, 2019; Scherer, 2017). The MIRA Pico is outfitted with a flexible length of tubing (see Sect. 2.1.4), and the constant pump flow rate results in a phase lag of ∼ 2 s based on laboratory tests of the response delay from pulsed gas releases (Sect. 2.1.4). The MIRA has been used extensively on the ground, including controlled-release tests at Colorado State University's (CSU) METEC facility to find and quantify leaks (Travis et al., 2020).

2.1.2 Onboard anemometer and weather sensor

In addition to hydrocarbon concentrations, direct flux quantification requires measurements of vector winds, temperature, and pressure. The system includes an ultra-light 3D sonic anemometer: the TriSonica Mini weather sensor (TWS) by Anemoment has a mass of 50 g with a volume of less than 450 cm³ (Anemoment LLC, 2021). Additionally, environmental measurements can also be used as a way to detect and adapt to unsafe conditions during flights (Hollenbeck et al., 2018). The TWS senses vector winds ($\left|\mathit{u}\right|$), temperature (T), and pressure (P) at 5 Hz with an accuracy of δu_m=δv_m = 0.2 m s⁻¹, δT = 2 C, and δP = 10 hPa, respectively (see Table 1).

Table 1 Sensor uncertainties for error propagation (Anemoment LLC, 2021; Dà Jiāng Innovations, 2018; Aeris Technologies Inc, 2019).

^∗ Average 1σ standard deviation of residuals from background estimation (ε); actual uncertainty is dependent on the quality of the background estimate (χ₀) for each of the individual datasets.

2.1.3 Uncrewed aerial vehicle

The DJI Matrice 600 Pro (M600P) is a hexarotor vertical takeoff and landing aircraft that is capable of flying with relatively large payload masses of up to 5.5 kg. This powerful and mobile platform became commercially available in 2014 and has been used in a variety of scientific and commercial applications (McKinney et al., 2019; Villa et al., 2016; Hollenbeck et al., 2018). The M600P can be remotely piloted at a distance of up to 5 km at a maximum altitude of 125 m above ground level. Flight times depend on factors such as payload mass, winds, and flight pattern. Under typical winds of 2–4 m s⁻¹ and the standard total payload mass of ∼ 3 kg, total flight times range between 18 and 25 min. We employ dual battery packs that are cycled between flying and charging in the field; due to battery life and charging limitations, the UAS typically can accommodate about one flight per hour. Additionally, the flux quantification method described in Sect. 2.6 relies on relatively stable wind fields around the source, optimally between 2 and 6 m s⁻¹.

The M600P is controlled via proprietary software and firmware with an expansive application programming interface (API) for telemetry logging and flight control (Dà Jiāng Innovations, 2018). However, the reliance on proprietary software limits system customizability and sensor integration (see Sect. 2.1.5). While it is possible to automate flight plans for the M600P, it is more important for the operator to maintain control of the UAS throughout the entire flight to account for changing or unexpected flight conditions.

2.1.4 Payload sensor mounting

Both the Raspberry Pi (RPi) onboard computer and the MIRA Pico mount to the underside of the M600P as shown in Fig. 1. The M600P achieves flight by funneling air downwards to create thrust, constantly displacing air around the vehicle's body in the process (Dà Jiāng Innovations, 2018). Computational fluid dynamics simulations carried out by McKinney et al. (2019) determined that the Matrice 600 Pro causes disturbances up to 5 m below the M600P's center of mass; however, the magnitude of the disturbances drops off quickly at a height of 60 cm above the M600P body (McKinney et al., 2019).

The TWS mounts to a mast made of lightweight, carbon fiber tubing (0.5 in diameter) anchored near the M600P's center of mass. This places the TWS ∼ 85 cm above the M600P and above the shallow inflow layer generated by propellers during normal flight. At this height above the propellers, air disturbances are below the noise floor of the TWS (∼ 0.2 m s⁻¹) or easily calibrated out during data processing. A ∼ 1 m section of flexible tubing connects the payload's gas cell to a sampling port approximately 80 cm above the body of the M600P. This physical offset between the MIRA and sampling port introduces a phase lag of ∼ 2 s between the measured MIRA concentrations and the TWS, which is accounted for prior to the calculation of instantaneous fluxes. The sampler port is outfitted with a plastic Y-shaped inlet for redundant openings to decrease the chance of damaging the pump due to debris clogging the inlet.

While the TWS has an integrated magnetometer to measure heading, this device is intended for stationary deployment in an open location where averaging of compass readings can result in a relatively accurate direction. However, any magnetic disturbances created from adjacent devices or operating motors perturb the heading measurements. In contrast, the M600P has a triple-redundant positioning system with an accuracy of < 3° at 200 Hz. The TWS must be physically aligned with the more accurate heading data from the M600P's inertial measurement unit (IMU) to accurately calculate static-frame vector winds (more in Sect. 2.2).

In order to ensure proper alignment, the M600P is set up facing north with the TWS mounting point tightened so that the anemometer's north arm is facing the same direction. The MIRA Pico and onboard computer are both mounted under the body of the M600P as shown in Fig. 1. Initial M600P test flights revealed that payload mounting at this location improves flight performance, especially during takeoff and landing, due to the lower center of mass.

2.1.5 Onboard data logging and transmission

Each of the three main devices (MIRA Pico, TWS, and M600P) has its own output data stream that must be separately parsed and then temporally synchronized with the other streams. The M600P does not support the use of DJI's proprietary “Payload SDK” – firmware enabling sensor data to be transmitted along with the radio frequency control signal. Therefore, remote monitoring and control of the fully integrated UAS are handled independently of the M600P via an onboard device for both acquisition and transmission. A fourth-generation Raspberry Pi (RPi) single-board computer is the ideal choice due to its size, weight, and programmability via the Linux Kernel. The RPi is powered independently via a 10 Ah payload battery (see Fig. 1), allowing for ground-based data collection between flights, typically while the M600P batteries are charging.

2.2 Relative wind adjustments

The TriSonica Mini weather sensor (TWS) is designed for static installation with one arm facing north (Anemoment LLC, 2021) and will therefore produce incorrect results when the device's heading (“yaw”, ϕ) is not a multiple of 2π (Anemoment LLC, 2021; Hollenbeck et al., 2018). In addition, the motion of the M600P during data collection induces an apparent wind that is folded into the TWS wind measurement. Transforming the raw TWS vector wind measurements (${\mathit{u}}_{\mathrm{m}}=[u_{\mathrm{m}},v_{\mathrm{m}},w_{\mathrm{m}}]$) to Earth-fixed coordinates therefore requires accurate, real-time measurements of ϕ and UAS velocity (V_s). In order to transform the measured TWS data to a static, Earth-fixed coordinate system, we apply a standard Galilean transformation as shown in Eq. (1).

$$\begin{array}{}\text{(1)}& \begin{array}{rl}\mathit{u}& =\mathbf{R}\left(\mathit{\varphi }\right){\mathit{u}}_{\mathrm{m}}+{\mathit{V}}_{\mathrm{s}}\\ & =\left[\begin{array}{ccc}\cos \left(\mathit{\varphi }\right)& \sin \left(\mathit{\varphi }\right)& \mathrm{0}\\ -\sin \left(\mathit{\varphi }\right)& \cos \left(\mathit{\varphi }\right)& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]\left[\begin{array}{c}{u}_{\mathrm{m}}\\ v_{\mathrm{m}}\\ w_{\mathrm{m}}\end{array}\right]+\left[\begin{array}{c}{V}_x\\ V_y\\ V_z\end{array}\right]\\ & =\left[\begin{array}{c}u\\ v\\ w\end{array}\right]\end{array}\end{array}$$

Here, ${\mathit{u}}_{\mathrm{m}}=[u_{\mathrm{m}},v_{\mathrm{m}},w_{\mathrm{m}}]$ are the TWS-measured vector winds, $\mathit{u}=[u,v,w]$ is the adjusted wind speed in Earth-fixed coordinates, and ${\mathit{V}}_{\mathrm{s}}=[V_x,V_y,V_z]$ is the instantaneous UAS velocity. R(ϕ) is the counterclockwise (CCW) Euler rotation matrix around the vertical axis ($\widehat{z}$), where ϕ is the heading (yaw) of the M600P.

Rotations around the $\widehat{y}$ and $\widehat{x}$ axes – caused by variations in M600P pitch (θ) and roll (φ), respectively – are both generally less than 3° during steady, level flights through the source plume and are therefore neglected in Eq. (1). The instantaneous pitch and roll angle can be significant when the M600P is being maneuvered to different altitudes between transects, but measurements collected during these spurious adjustments are filtered out during the flux quantification process (see Sect. 2.6).

Figure 2 Raw measured vector winds (a) and heading- and velocity-corrected vector winds (b) before, during, and after flight. Pre- and post-flight ground wind speeds at the takeoff location highlighted in gray.

A time series of the TWS-measured raw winds during a typical flight is shown in Fig. 2, along with the derived static-frame vector winds based on the M600P headings and velocities throughout the flight. The raw-wind measurements show clear signatures of back-and-forth, quasi-horizontal motions of the M600P during sampling of a methane plume. Winds adjusted for changes in heading and horizontal translation show a more realistic, continuous structure, but some residual effects can be discerned at turning points in the flight, when the platform's horizontal acceleration is large and the corresponding pitch and roll angles become appreciable. In general, pitch and roll angles during level, steady flight sections are less than 3° (see Fig. A1), but these angles can increase to 10° or more during rapid horizontal accelerations. We estimate that a 15° pitch or roll angle may introduce an error of up to 4 % in measured horizontal winds; therefore data obtained during turning points and deliberate changes in flight altitude are screened from the analysis of methane fluxes (see Sect. 2.6). The adjusted and filtered winds have been compared to tower anemometer measurements during controlled-release validation experiments (Sect. 3), and these show agreement to within ± 15 % over ∼ 1 min averages and with spatial separations of about 100 m. Individual wind vector components have also been compared in the field (see Fig. A2) when the UAS was flown between 10 and 30 m downwind of a static, single-component wind gauge. Comparisons of 5 s averages show agreement to within 10 % over a range of winds between 2 and 7 m s⁻¹.

2.3 Background concentration estimation

While the MIRA Pico is highly sensitive with a large dynamic range (∼ 20 ppb < ${\mathit{\chi }}_{{\mathrm{CH}}_{\mathrm{4}}}$ < 4000 ppm) (Aeris Technologies Inc, 2019; Meyer et al., 2022; Follansbee et al., 2024), there are observable levels of semi-periodic drift in the measured mole fractions for both hydrocarbons at the parts per billion level. This effect is especially noticeable for ${\mathit{\chi }}_{{\mathrm{C}}_{\mathrm{2}}{\mathrm{H}}_{\mathrm{6}}}$, partly because ambient C₂H₆ concentrations are generally ∼ 100 × smaller than CH₄. In addition to these instrumental drifts, the raw concentration measurements (χ) may reflect actual changes in ambient background. Nearby point sources do not directly influence these background variations; instead they may be related to regional-scale emissions and meteorological influences such as winds and stability at the mesoscale level. Both of these effects contribute to structure in the measured background (χ₀), which must be accurately quantified in order to isolate enhancements from the target source. This background variability generally occurs over periods of minutes and is independent between ${\mathit{\chi }}_{{\mathrm{CH}}_{\mathrm{4}}}$ and ${\mathit{\chi }}_{{\mathrm{C}}_{\mathrm{2}}{\mathrm{H}}_{\mathrm{6}}}$, so direct comparison prior to background removal can be misleading, especially over longer period datasets.

Figure 3 Time series of methane measurements during a typical flight, illustrating the procedure for background estimation and removal. The flight involved multiple cross-sections through the downwind plume of a known point emitter, and the plume enhancements in methane are visible as large spikes in the methane time series. (a) Three-point-smoothed raw data. (b) Gradient filter to detect the majority of in-plume samples. (c) Split remaining samples into segments and remove statistical outliers from each segment. (d) Third-order polynomial fit using remaining samples; resultant regression coefficients used with initial (non-smoothed) raw data samples used to estimate background. Estimated CH₄ emission rate for this dataset is 0.378 ± 0.147 kg h⁻¹ (see Sect. 2.6).

The UAS measurement strategy involves repeated crosswind transects through the dispersing plume downwind of a target source (as discussed in more detail in Sect. 2.5). The resulting time series concentration data contain a number of finite-width enhancements superimposed on a slowly varying background, as shown in Fig. 3a. Over the course of these measurements, the solar insolation was increasing and the wind speed and direction were changing, which led to a monotonic increase in the background methane concentration as seen with the increasing background. It should be noted that each plume “spike” in Fig. 3a should not be expected to have the same peak height or width, since these plume intersects all occur at different times, altitudes, and downwind distances.

The background, χ₀, is estimated by fitting a polynomial to the measurements outside the plume. After extensive testing, the procedure described below and outlined in Fig. 3 proved to be the most effective method of filtering inter-plume samples. The raw χ time series is initially three-point smoothed ($\stackrel{\mathrm{̃}}{\mathit{\chi }}$) to decrease the effect of noise before further filtering (Fig. 3a). The gradient of $\stackrel{\mathrm{̃}}{\mathit{\chi }}$ is then used to detect abrupt changes in $\stackrel{\mathrm{̃}}{\mathit{\chi }}$ which are indicative of the UAS entering or exiting the higher-concentration plume, as shown in Fig. 3. Samples where $\mathrm{\nabla }\stackrel{\mathrm{̃}}{\mathit{\chi }}$ is greater than a specified threshold get filtered out (i.e., more than ± γσ, where γ is a predefined constant; ± 0.2σ in Fig. 3b and in the subsequent analysis of measurements presented here). This gradient filter removes the majority of the inter-plume samples, leaving behind the “background” samples which vary at a much lower frequency. However, the gradient filter is imperfect, and occasionally some inter-plume samples are not removed during this step. Therefore, a follow-up outlier filter is used to remove samples which are significantly far from the mean of the samples (more than 1 standard deviation from the mean in Fig. 3c). The remaining samples are taken to be the background with respect to the target plume (background measurements, ${\stackrel{\mathrm{̃}}{\mathit{\chi }}}_{\text{bg}}$), and a variable-order polynomial is fit to these samples. The polynomial order is empirically selected based on the individual dataset, usually between second and fourth order (cubic in Fig. 3d). Higher-order polynomials can affect the edges of the fitted time series, but generally, higher-order polynomials (> 4) are not required after appropriately tuning the gradient and outlier filter threshold parameters. The polynomial regression coefficients fitting ${\stackrel{\mathrm{̃}}{\mathit{\chi }}}_{\text{bg}}$ are then applied to the original χ sampled timestamps to estimate χ₀ for each independent hydrocarbon and dataset. This estimated background is then subtracted from the original time series to get the isolated source plume enhancements, (χ−χ₀), for both CH₄ and C₂H₆.

Figure 4 presents the distribution of CH₄ and C₂H₆ background estimate residuals (${\mathit{\epsilon }}_{\mathrm{0}}={\mathit{\chi }}_{\mathrm{0}}-{\stackrel{\mathrm{̃}}{\mathit{\chi }}}_{\text{bg}}$) from three independent flights around different source types. These three datasets were specifically chosen to highlight the variability in ε₀ caused in part by the source strength and environmental conditions. The average 3σ confidence interval (99.7 %) calculated from more than two dozen flights is approximately 16 ppb for CH₄ and 2.5 ppb for C₂H₆. This minimum detection limit governs the lower bounds on source strengths and fluxes that can be quantified with this UAS (further discussed in Sect. 2.7)

Figure 4 Residuals from the CH₄ and C₂H₆ background polynomial fit (outlined in Fig. 3) for three independent flights. The 3σ confidence interval (99.7 %) for each dataset is shown with a dashed line, highlighting the variability in the background estimates due to different sources and environmental conditions. The average 3σ confidence interval, calculated from 28 flights during 2022 and 2023, is overlaid with a solid ellipse.

2.4 Plume simulations

This section describes the modeling analysis that was conducted to develop optimal flight patterns and sampling strategies for the UAS. In these idealized simulations, emissions from the target source are held constant and surrounding environmental winds are steady during the M600P's 15 to 25 min flight times. Gaussian models offer a reasonable approximation for the structure and evolution downwind of the source's emission plume under these conditions (Shah et al., 2019; Stockie, 2012; Jacob et al., 2022; Seinfeld and Pandis, 2006; Meyer et al., 2022; Follansbee et al., 2024; Woodward, 1998), but it should be emphasized that source fluxes reported in this work do not rely on the results of Gaussian plume simulations or inverse modeling (see Sect. 2.6).

The Gaussian plume equation for the mass density of a gas (C, units of kg m⁻³) downwind of a point-source emitter is shown in Eq. (2).

$$\begin{array}{}\text{(2)}& \begin{array}{rl}C(x,y,z)=& {\displaystyle \frac{Q}{\left|\mathit{u}\right|}}\cdot {\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }{\mathit{\sigma }}_y\left(x\right){\mathit{\sigma }}_z\left(x\right)}}\\ & \times \exp \left[-y^{\mathrm{2}}/\left(\mathrm{2}{\mathit{\sigma }}_y(x{)}^{\mathrm{2}}\right)\right]\\ & \times \exp \left[-(z-H{)}^{\mathrm{2}}/\left(\mathrm{2}{\mathit{\sigma }}_z(x{)}^{\mathrm{2}}\right)\right]\end{array}\end{array}$$

Q is the emission rate (kg s⁻¹); u is the constant horizontal wind speed (m s⁻¹), and σ_z and σ_y are standard deviations for Gaussian distributions in units of meters, which are generally derived using stability classes (Seinfeld and Pandis, 2006; Woodward, 1998). H is the height of the source (m). C is the estimated increase in gas density at position x (directly downwind, centerline), y (horizontal from centerline), and z (vertical from centerline) given in meters. It is important to make the distinction between the mass density, C (kg m⁻³), and the mole fraction, χ (mol mol⁻¹), measured by the MIRA Pico. The conversion between the two is C=ηρχ, where η is the ratio of molar masses of the gas to that of air, and ρ is the air density (kg m⁻³).

Figure 5 Cross-sections of modeled Gaussian plume under conditions of low wind (2 m s⁻¹) and high solar insulation (stability class “A” from Woodward, 1998). Each panel depicts the Gaussian plume cross-section at distance X downwind of the source. The horizontal distance y is with respect to plume centerline, and the vertical distance z is above ground level. The plume model assumes a steady wind field and constant emission rate from the source.

Figure 5 shows cross-sections of a Gaussian plume simulated at three downwind distances from a 25 kg h⁻¹ point source. Mean winds in this case are ∼ 2 m s⁻¹, and the stability class corresponds to high solar insolation. As expected, the simulated plume becomes more dispersed and broadens both horizontally and vertically with increasing distance from the source. Plume enhancements in methane range between 1 and 5 ppm, which are easily detectable with this system, and the horizontal and vertical scales (∼ 50 m, ∼ 30 m) are readily accessible to the UAS.

2.5 Flight pattern

The full UAS setup and pre-flight checks can be completed in under 15 min, and any certified remote pilot, or person under the supervision of a certified pilot (Federal Aviation Administration, 2021), can operate and maneuver the UAS around a target source's emission plume. Part of the setup process involves the verification of MIRA measurements for both CH₄ and C₂H₆. Prior to each flight, a small canister of natural gas with a known composition (∼ 7 % C₂H₆:CH₄ ) is rapidly opened and closed about 1 m upwind of the UAS gas inlet. This pre-flight release “pulse” is measured on the MIRA to test for any lag or unexpected gain offset on either of the channels. As discussed in Sect. 2.3, the raw concentration measurements can be biased due to a quasi-periodic sensor drift which must be estimated and removed before comparing the channels. This step is completed during data processing after each deployment, and it should be noted that the results from each of these controlled pulses has been consistent to within a fraction of the known C₂H₆:CH₄ ratio (< 1 % disagreement) across all deployments between 2021 and the present.

The flight time of ∼ 20 min constrains the types of flight patterns that can be used to quantify source emission rates. In order to maximize the number of in-plume samples within a relatively short flight period, our strategy is to remain downwind of the source and fly horizontal “transects” perpendicular to the mean wind direction. This boustrophedonic curtain flight pattern is perpendicular to the average wind direction and involves multiple cuts through the plume in order to measure both the in-plume concentrations, χ, and the ambient or background concentrations used to determine χ₀. Each individual transect through the downwind plume is at a relatively constant altitude and horizontal velocity between 2 and 5 m s⁻¹, depending on the wind conditions and proximity to the source. Transects are typically 50 m to 1 km in length, depending on terrain, wind variability, and source distribution and downwind distance. Extended sources and measurements collected at larger downwind distances from the source require longer transects to ensure that the dispersing plume gets fully traversed during each transect. The ability to quickly adjust altitude between each of the individual transects allows for direct measurements of the plume's vertical structure in addition to the horizontal dispersion.

Direct measurements of a source plume with the MIRA Pico sensor require the UAS to be physically maneuvered downwind of the source. While a steady-state Gaussian plume model can be used for flux inversions (e.g., Jacob et al., 2022; Shah et al., 2019; Stockie, 2012; Seinfeld and Pandis, 2006; Woodward, 1998; Bhattacharya, 2013), these models represent plume dispersion probabilities which are not generally observed in the superposition of horizontal transects through the plume. During direct measurements, large-scale turbulence in the wind flow results in changes to the plume's size and relative location throughout the measurement period. Based on simple scaling arguments, we expect that plume variability is largest in the horizontal plane due primarily to wind directional variability. Figure 6 shows multiple simulated transects through a Gaussian plume's downwind cross-section in the case of a steady-state flow. We simulate the stochastic variability in the plume by introducing random shifts in the centerline between each of the temporally separated transects, which can be seen at the bottom of Fig. 6. This pattern is generally more consistent with our observations in comparison to the steady-state model. However, the stochastic model shows that total plume flux is unaffected by horizontal plume displacements, provided that the integral is taken along horizontal transects.

Figure 6 Simulated flight through Gaussian plume: (a) cross-section of plume X meters directly downwind of source, with simulated transects shown as dashed lines. (b) Resample of modeled cross-section at constant horizontal transect velocity V_s. (c) Random shifts in the centerline horizontal position, simulating variable wind conditions and temporally separated plume transects.

2.6 Direct flux quantification

The data collected from a single flight are broken into K individual transects through the source's downwind plume, each transect representing data collected at a nearly constant altitude, z, and non-zero horizontal velocity, V_s, while intersecting the plume. The altitude, z, is only adjusted between transects with the UAS outside the plume (measuring ambient or background concentrations).

Mass balance estimation techniques require knowledge of the wind speed perpendicular to the direction of travel (transect crosswind). Optimally, each transect would therefore be perpendicular to the mean horizontal wind direction ($\stackrel{\mathrm{‾}}{\mathit{u}}$), but local topography, changing wind fields, and flight safety considerations often result in an angle other than 90° between V_s and u. Therefore, it is important to first define a unit vector, $\widehat{\mathit{n}}$, perpendicular to the direction of travel.

$$\begin{array}{}\text{(3)}& \left\{\begin{array}{l}\Vert \widehat{\mathit{n}}\Vert =n_x^{\mathrm{2}}+n_y^{\mathrm{2}}=\mathrm{1}\\ {\mathit{V}}_{\mathrm{s}}\cdot \widehat{\mathit{n}}=V_x{n}_x+V_y{n}_y=\mathrm{0}\end{array}\right.\end{array}$$

Equation (3) shows the system of equations used to calculate this unit vector, $\widehat{\mathit{n}}$, using the instantaneous horizontal UAS velocities V_x and V_y defined previously. The transect crosswind, a measure of wind speed perpendicular to the UAS transect, is $(\mathit{u}\cdot \widehat{\mathit{n}})$. Maintaining a constant heading and horizontal velocity during flights proved difficult during periods of irregular and shifting winds. Direct measurements of u and V_s are variable between samples, so $\widehat{\mathit{n}}$ is required to constrain the flux calculation with respect to each transect.

Each of the K horizontal transects gets processed individually to calculate the intermediate transect-integrated flux f_k (Eq. 4), the horizontal integral of the samples along each transect (units of mass flux rate per vertical distance, $\mathrm{kg}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{1}}$).

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{f}_k& =\mathit{\eta }\sum _{i=\mathrm{0}}^{n-\mathrm{1}}{\mathit{\rho }}_i(\mathit{\chi }-{\mathit{\chi }}_{\mathrm{0}}{)}_i(\mathit{u}\cdot \widehat{\mathit{n}}{)}_i\mathrm{\Delta }{s}_i\\ & =\sum _{i=\mathrm{0}}^{n-\mathrm{1}}(C-C_{\mathrm{0}}{)}_i(\mathit{u}\cdot \widehat{\mathit{n}}{)}_i\mathrm{\Delta }{s}_i\end{array}\end{array}$$

$\mathit{\eta }\mathit{\rho }(\mathit{\chi }-{\mathit{\chi }}_{\mathrm{0}})=(C-C_{\mathrm{0}})$ is the background-adjusted mole fraction, $\mathrm{\Delta }s=\left|{\mathit{V}}_{\mathrm{s}}\right|\mathrm{\Delta }t=\sqrt{V_x^{\mathrm{2}}+V_y^{\mathrm{2}}}\mathrm{\Delta }t$ is the distance between samples (Δt≈1 s), and ρ and η have been defined previously.

The total flux, F_tot, is then calculated through integrating the sum of f_k and the vertical distance between physically adjacent transects, Δz_k (Eqs. 5 and 6).

$$\begin{array}{}\text{(5)}& {\displaystyle }\begin{array}{rl}& \mathrm{\Delta }{z}_k=\\ & \left\{\begin{array}{ll}\left|{\stackrel{\mathrm{‾}}{z}}_{k+\mathrm{1}}-{\stackrel{\mathrm{‾}}{z}}_k\right|/\mathrm{2}+\mathrm{3}\left|{\stackrel{\mathrm{‾}}{z}}_k\right|/\mathrm{4}& {\stackrel{\mathrm{‾}}{z}}_k=min({\stackrel{\mathrm{‾}}{z}}_{\mathrm{0}},\mathrm{\dots },{\stackrel{\mathrm{‾}}{z}}_{K-\mathrm{1}})\\ \left|{\stackrel{\mathrm{‾}}{z}}_{k-\mathrm{1}}-{\stackrel{\mathrm{‾}}{z}}_k\right|& {\stackrel{\mathrm{‾}}{z}}_k=max({\stackrel{\mathrm{‾}}{z}}_{\mathrm{0}},\mathrm{\dots },{\stackrel{\mathrm{‾}}{z}}_{K-\mathrm{1}})\\ \left(\left|{\stackrel{\mathrm{‾}}{z}}_{k-\mathrm{1}}-{\stackrel{\mathrm{‾}}{z}}_k\right|+\left|{\stackrel{\mathrm{‾}}{z}}_{k+\mathrm{1}}-{\stackrel{\mathrm{‾}}{z}}_k\right|\right)/\mathrm{2}& \text{otherwise}\end{array}\right.\end{array}\text{(6)}& {\displaystyle }{F}_{\text{tot}}=\sum _{k=\mathrm{0}}^{K-\mathrm{1}}\left(f_k\mathrm{\Delta }{z}_k\right)\end{array}$$

Figure 7 Plume transect flight pattern viewed from east (c), south (d), and above (b). Approximate source location marked with black square. The wind rose plot (a) shows the Earth-fixed wind in the horizontal direction. Estimated CH₄ emission rate from the UAS is 4.1 ± 1.11 kg h⁻¹ (compare to metered emission rate of 3.52 ± 0.01 kg h⁻¹).

For most transects, Δz_k ranges from the midpoint distance between the previous, lower-altitude transect up to the midpoint between the next, higher-altitude transect. Adjustments to the vertical integral step Δz_k are taken at the bottom and top transects to account for extrapolation to the ground and above the flight pattern. The UAS cannot profile below a height of about 1.5 m, but the plume flux is expected to decrease nonlinearly to zero at the ground. Therefore, we extrapolate below the lowest flight transect by assuming a constant mixing ratio and wind speed between the lowest transect and 75 % of the distance to ground. This choice of integral quadrature corresponds to a logarithmic vertical profile for the transect-integrated flux f below level z, $f=f_{\mathrm{0}}\ln (z/z_{\mathrm{0}})$, with z₀=0.018z. In the absence of tall trees, power lines, or other obstructions, the lowest flight transect typically lies between 2–4 m above ground level (including the UAS mast), and the assumed z_o therefore ranges between about 4 and 7 cm. Extrapolation at the top edge of the plume can be more complicated, but ideally, the uppermost flight transect will lie above the plume so that a linear interpolation can adequately account for the top of the plume profile. However, some flight patterns (e.g., Fig. 9) may not completely span the possible vertical extent of the plume, and in these isolated cases we adopt a conservative approach by assuming a constant mixing ratio and wind above the transect altitude, extending only to a height defined by half the vertical distance between the top two transects. Typically, this extrapolation amounts to 2–4 m above the uppermost transect. A typical flight downwind of a natural gas point source is shown in Fig. 7. This flight is composed of 11 horizontal transects over a period of approximately 11 min. The transects were flown approximately 130 m downwind (roughly south-southeast) of the methane point source with a known emission rate of 3.42 ± 0.01 kg h⁻¹ (El Abbadi et al., 2024).

2.7 Precision and uncertainties

Errors are calculated using the standard Eq. (7) described in Taylor (1996) for an arbitrary measured variable, q. Table 1 gives the measurement uncertainties for the major components of the UAS.

$$\begin{array}{}\text{(7)}& \mathit{\delta }q(x,\mathrm{\dots },z)=\sqrt{{\left({\displaystyle \frac{\partial q}{\partial x}}\right)}^{\mathrm{2}}\mathit{\delta }{x}^{\mathrm{2}}+\mathrm{\dots }+{\left({\displaystyle \frac{\partial q}{\partial z}}\right)}^{\mathrm{2}}\mathit{\delta }{z}^{\mathrm{2}}}\end{array}$$

The transect-integrated flux, f_k, given in Eq. (4) is a function of the background-adjusted concentration (C−C₀), wind speed perpendicular to UAS curtain, and the distance between transect samples (Δs). Each of these terms introduces an error in f_k as defined in Eqs. (8) through (10).

The uncertainty on the background-adjusted concentration, C, is dependent on the calculated air density ($\mathit{\rho }=M_{\text{air}}P/\left(RT\right)$ kg m⁻³), measured hydrocarbon concentration (χ), estimated background (χ₀; see Sect. 2.3), and the uncertainty associated with each value (see Table 1).

$$\begin{array}{}\text{(8)}& \begin{array}{rl}\mathit{\delta }(C-C_{\mathrm{0}})& =\mathit{\eta }\sqrt{{\mathit{\rho }}^{\mathrm{2}}\mathit{\delta }(\mathit{\chi }-{\mathit{\chi }}_{\mathrm{0}}{)}^{\mathrm{2}}+(\mathit{\chi }-{\mathit{\chi }}_{\mathrm{0}}{)}^{\mathrm{2}}\mathit{\delta }{\mathit{\rho }}^{\mathrm{2}}}\\ & =\mathit{\eta }{M}_{\text{air}}/R\sqrt{\begin{array}{c}(P/T{)}^{\mathrm{2}}[\mathit{\delta }\mathit{\chi }+\mathit{\delta }{\mathit{\chi }}_{\mathrm{0}}]+(\mathit{\chi }-{\mathit{\chi }}_{\mathrm{0}}{)}^{\mathrm{2}}\\ \times \left[(\mathit{\delta }P/T{)}^{\mathrm{2}}+(P\mathit{\delta }T/T^{\mathrm{2}}{)}^{\mathrm{2}}\right]\end{array}}\end{array}\end{array}$$

Note that the uncertainty on the estimated χ₀ is calculated independently for each flight.

The crosswind uncertainty is dependent on the static-frame horizontal wind speed (u) calculated as described in Sect. 2.2 with the full uncertainty propagation given in Appendix B1. While the crosswind is dependent on the unit vector ($\widehat{\mathit{n}}$), the conservative upper-limit approximation for crosswind is given in Eq. (9).

$$\begin{array}{}\text{(9)}& \begin{array}{rl}\mathit{\delta }(\mathit{u}\cdot \widehat{\mathit{n}})& =\sqrt{n_y^{\mathrm{2}}\mathit{\delta }{u}^{\mathrm{2}}+n_x^{\mathrm{2}}\mathit{\delta }{v}^{\mathrm{2}}}\\ & \approx \sqrt{\mathit{\delta }u(u_{\mathrm{m}},\mathit{\varphi },V_x{)}^{\mathrm{2}}+\mathit{\delta }v(v_{\mathrm{m}},\mathit{\varphi },V_y{)}^{\mathrm{2}}}\end{array}\end{array}$$

The final parameter in δf_k is the separation between the sample along the transect, Δs, which is dependent on the UAS horizontal velocity ($V_{\mathrm{s}}=\sqrt{V_x^{\mathrm{2}}+V_y^{\mathrm{2}}}$); sampling rate (Δt = 1 Hz); and the associated uncertainties as shown in Eq. (10).

$$\begin{array}{}\text{(10)}& \begin{array}{rl}\mathit{\delta }\left(\mathrm{\Delta }s\right)& =\sqrt{\mathit{\delta }{V}_{\mathrm{s}}^{\mathrm{2}}\mathrm{\Delta }{t}^{\mathrm{2}}+V_{\mathrm{s}}^{\mathrm{2}}\mathit{\delta }(\mathrm{\Delta }t{)}^{\mathrm{2}}}\\ & =\sqrt{\begin{array}{c}{\left[(V_x\mathit{\delta }{V}_x{)}^{\mathrm{2}}+(V_y\mathit{\delta }{V}_y{)}^{\mathrm{2}}\right]}^{\mathrm{2}}(\mathrm{\Delta }t/V_{\mathrm{s}}{)}^{\mathrm{2}}\\ +[V_x^{\mathrm{2}}+V_y^{\mathrm{2}}{]}^{\mathrm{2}}\mathit{\delta }(\mathrm{\Delta }t{)}^{\mathrm{2}}\end{array}}\end{array}\end{array}$$

The calculated uncertainty on f_k is therefore calculated as shown in Eq. (11).

$$\begin{array}{}\text{(11)}& \begin{array}{rl}\mathit{\delta }{f_k}^{\mathrm{2}}=& {\left[\mathit{\delta }(C-C_{\mathrm{0}})(\mathit{u}\cdot \widehat{\mathit{n}})\mathrm{\Delta }s\right]}^{\mathrm{2}}\\ & +{\left[(C-C_{\mathrm{0}})\mathit{\delta }(\mathit{u}\cdot \widehat{\mathit{n}})\mathrm{\Delta }s\right]}^{\mathrm{2}}\\ & +{\left[(C-C_{\mathrm{0}})(\mathit{u}\cdot \widehat{\mathit{n}})\mathit{\delta }\left(\mathrm{\Delta }s\right)\right]}^{\mathrm{2}}\end{array}\end{array}$$

The total source flux is the summation of each f_k multiplied by the corresponding Δz_k, as shown in Eq. (6). The vertical range of each independent transect is given in Eq. (5). $\stackrel{\mathrm{‾}}{z_k}$ is the average height of all samples within transect k, and $\mathit{\delta }\left(\mathrm{\Delta }z\right)={\mathit{\sigma }}_{z_k}$ is the uncertainty in altitude.

$$\begin{array}{}\text{(12)}& \mathit{\delta }{F}_{\text{tot}}=\sqrt{\mathit{\delta }{f}_k^{\mathrm{2}}\mathrm{\Delta }{z}^{\mathrm{2}}+f_k^{\mathrm{2}}\mathit{\delta }(\mathrm{\Delta }z{)}^{\mathrm{2}}}\end{array}$$

The measurement uncertainties detailed in Eqs. (8) through (12) can be used to identify the major sources of uncertainty in source emission rate estimates.

While the minimum flux threshold will depend on the specific meteorological conditions and flight patterns during sampling, under favorable flight conditions with high solar insulation, 2 m s⁻¹ mean wind speed, and a downwind plume cross-sectional area of 100 m² (see Fig. 5), the minimum source strength that can be quantified is on the order of 0.0062 kg (CH₄) h⁻¹ and 0.0018 kg (C₂H₆) h⁻¹.

3 Results

This section gives an overview of a number of quantification flights between spring 2022 and fall 2023. Section 3.1 details the results of multiple controlled releases using the system described in El Abbadi et al. (2024). Section 3.2 discusses measurements from various, smaller biogenic and thermogenic emission sources collected between winter 2021 and fall 2023.

3.1 Controlled-release experiments

In fall of 2022, the NMT team participated in a single-blind controlled-release validation campaign in Casa Grande, Arizona, USA (Sherwin et al., 2024; El Abbadi et al., 2024). While this validation campaign focused on larger-scale aircraft and satellite system validation, it was a unique opportunity to quantify the accuracy for this method of direct methane flux quantification. The UAS was deployed multiple times around the controlled-release location during two separate 3 d validation trips: 10–12 October and 14–16 November 2022. Multiple flights at various times throughout the daily testing window (between 10:00 and 14:00 MST most days) measured concentrations downwind of the release stack.

Suboptimal environmental conditions and unforeseen instrumentation issues resulted in the rejection of the majority of October 2022 flights from the analysis process. The November 2022 3 d campaign produced a total of 12 flights, 8 of which were during times of low-variability wind fields and good flight conditions to allow for reliable flux quantification.

It is important to note that, at the time of this campaign, the UAS flight pattern and quantification method were still in active development. While the NMT team quantified and submitted flux estimates during the single-blind and partially blinded phases of the unblinding process (described in El Abbadi et al., 2024), the final unblinded metered emissions were invaluable to debug and update the initial quantification process, ultimately leading to the numerical integration outlined in Sect. 2.6. Here, we present the analysis of the quantification results informed by the unblinded, ground-truth methane flow rates.

Figure 8 Comparison of the UAS CH₄ flux calculation versus metered emissions reported by El Abbadi et al. (2024) and Sherwin et al. (2024). Both plots show eight individual flights during the November 2022 field campaign with 2σ error bars. Linear regression fit shown as solid line for comparison with the idealized one-to-one fit corresponding to perfect agreement between UAS and metered data (dotted line). (a) Logarithmic plot highlighting the UAS's large dynamic range. (b) Linear plot focused on emission rates below detection limits of most aircraft quantification methods; linear regression fit does not include emission rates above 100 kg h⁻¹.

Figure 8 shows the comparison between UAS-calculated flux and the corresponding unblinded metered emissions from eight independent November 2022 flights. The metered rates, taken to be the ground-truth emission rates, range from 1.7 to 1500 kg h⁻¹. This broad range of emission rates highlights the system's dynamic range, and the fitted linear regression shows the quality of emission estimates compared to the one-to-one agreement with metered rates (dotted line). Our method of flux quantification shows reasonable agreement with the expected emission rates above 1 kg h⁻¹. However, the results suggest an underestimation of the emissions which is not correlated to number of transects or wind speeds during flights. This may be due to the limited flight time and under-sampling the downwind plume with the UAS. One interpretation of these results is that our field measurements could represent a reasonable lower bound to the true emission rate of the source.

3.2 Targeted point and distributed sources

The deployment strategy for this UAS is designed to quantify emissions from targeted local sources such as O&NG wells, manure storage, and biogenic lagoons. With a limited operational distance and flight time, larger distributed sources (such as large-scale dairies and agricultural centers) are difficult to properly quantify using the numerical integration technique described in Sect. 2.6. However, smaller-scale municipal waste facilities with heterogeneous emission profiles can be quantified as long as the UAS is able to fully, and repeatedly, transect the complex and irregular plumes downwind of the facility.

3.2.1 Municipal waste facility

During the course of system development and testing, the UAS was deployed around a local municipal waste facility (MWF) in Socorro, New Mexico (NM). This location serves a county population of approximately 16 300 with multiple cells over an area of 45 ha. Larger facilities studied by Olaguer et al. (2022) and Lan et al. (2015) reported emission rates between 85 and > 2000 kg h⁻¹; this much smaller local facility, however, was expected to emit at rates less than 10 kg h⁻¹ (Olaguer et al., 2022; Bogner and Matthews, 2003). The low emission rate and temporal heterogeneity, likely due to changes in cell activity, are evident across the multiple visits between spring 2022 and summer 2023. The flight path from one such deployment can be seen in Fig. 9.

Figure 9 Example flight, split into individual transects, from a landfill in New Mexico. The landfill itself is composed of active and inactive cells distributed over an area of approximately 45 ha, resulting in a non-uniform downwind plume profile. In comparison to the localized “point source” presented in Fig. 7, longer horizontal transects across a wider vertical range are required to measure the multiple distributed and highly variable emitters. Estimated flux rate for this flight is 1.33 ± 0.58 kg h⁻¹.

In comparison to flights downwind of point sources such as O&NG wells or controlled natural gas (NG) releases (Fig. 7), the flight pattern in Fig. 9 was positioned about 0.5 km downwind of the source and involved plume transects up to ∼ 500 m long, and over an altitude range of up to 50 m. In addition, the methane enhancements are about a factor of 10 lower than in Fig. 7, and the plume is distributed over a larger area and centered at a higher altitude. Nevertheless, the overall plume structure allows for the same flux analysis as discussed for the localized NG point sources above. The measured methane emission rate for this flight was 1.33 ± 0.58 kg h⁻¹.

3.2.2 Orphan well

During April 2023, the UAS was also deployed around an “orphaned well” located in Hobbs, New Mexico, USA. This site has been out of use for more than 2 decades and was in the initial phase of being plugged when visited in April 2023. A detailed analysis of data from this field campaign is presented in Follansbee et al. (2024). The observed plume structure was similar to that shown in Fig. 7, but mean wind speeds were much larger and methane plume enhancements were ∼ 400 ppb, roughly between the range of plume enhancements shown in Figs. 7 and 9. Three downwind flights yielded consistent fluxes in the range of 0.3 to 0.4 kg h⁻¹ (Follansbee et al., 2024).

3.2.3 Wastewater treatment plant

Wastewater treatment plants (WWTPs) are a known source of biogenic methane emissions via biodegradation of pollutants by anaerobic bacteria (Song et al., 2023). The small town of Socorro, NM, has a local wastewater treatment plant that processes less than 1 ×10⁶ gal d⁻¹ (44 L s⁻¹), so its contribution to anthropogenic methane is relatively quite low (when compared to O&NG and agricultural operations). However, the low emission rate of this location was useful for the testing of the UAS's lower detection and quantification limits.

A time series of measured CH₄ and C₂H₆ concentrations from a short flight in August 2023 is shown in Fig. 10. The CH₄ and C₂H₆ detection limit (calculated from the background fit residual 3σ confidence interval; see Sect. 2.3) is overlaid on the time series shown in Fig. 9. While the in-plume CH₄ concentration levels peak at around 40 ppb, no corresponding C₂H₆ enhancements are detected while passing through the plume. The lack of a correlation between C₂H₆ and CH₄ is strong evidence that the measured plume is from anaerobic digestion at the WWTP and is not contaminated by another source. However, it is important to note that weaker source strengths lead to a decrease in the signal-to-noise ratio, so C₂H₆ from thermogenic sources may not be detectable. This is shown visually at the bottom of Fig. 10 where the expected C₂H₆ time series for 10 % and 5 % C₂H₆:CH₄ thermogenic plumes are modeled with respect to the measured CH₄ time series. While the C₂H₆ signal for the 10 % thermogenic mixture peaks above the detection threshold (1.5 ppb in this dataset), the 5 % mixture is barely visible above the noise floor. However, it is important to note that weaker source strengths lead to a decrease in the signal-to-noise ratio, so C₂H₆ from thermogenic sources may not be detectable. This is shown visually at the bottom of Fig. 10 where the expected C₂H₆ time series for 10 % and 5 % C₂H₆:CH₄ thermogenic plumes are modeled with respect to the measured CH₄ time series. While the C₂H₆ signal for the 10 % thermogenic mixture peaks above the detection threshold (1.5 ppb in this dataset), the 5 % mixture is barely visible above the noise floor. The emission rate for the flight shown in Fig. 10 is estimated to be $\mathrm{68.7}(+\mathrm{171.5},-\mathrm{61.7})\mathrm{g}{\mathrm{h}}^{-\mathrm{1}}$.

Figure 10 (a, b) Time series of measured and background subtracted CH₄ and C₂H₆ concentrations (Sect. 2.3) from a flight downwind of the city of Socorro wastewater treatment plant in August 2023. (c) Expected C₂H₆ time series if measured CH₄ plume (top) was thermogenic with 10 % and 5 % C₂H₆:CH₄ ratios. The mean winds were ∼ 5 m s⁻¹, and the solar insolation corresponded to stability class “B” due to partial cloud coverage (Seinfeld and Pandis, 2006; Woodward, 1998). The plume was intercepted on two transects separated by about 40 s, as seen in the CH₄ time series (top) about 4–5 min into the flight. No corresponding C₂H₆ enhancements we measured indicated that the plume is biogenic (compare with expected thermogenic C₂H₆ time series). Estimated CH₄ emission rate of 0.061 ± 0.032 kg h⁻¹.

3.3 Source attribution

As discussed in Sect. 1, simultaneous measurements of both CH₄ and C₂H₆ are useful for source classification as either biogenic or thermogenic. In addition, there are varying ratios of trace hydrocarbons found in thermogenic natural gas sources, and C₂H₆ is the second-most dominant compound in refined natural gas (Peischl et al., 2018; Hodnebrog et al., 2018; Kutcherov and Krayushkin, 2010; Glasby, 2006; Space Studies Board et al., 2019; Solomon et al., 2007; Hansen et al., 2000; Meyer et al., 2022). Therefore, the ratio of C₂H₆ to CH₄ can be used to estimate the percentage of non-methane compounds in an O&NG plume and distinguish between sources.

Figure 11 shows the ratio of measured C₂H₆ and CH₄ mole fraction from five different field measurements around various sources, including the municipal waste facility discussed in Sect. 3.2.1. Multiple downwind plume measurements showed a negligible C₂H₆ content (<|0.2 %| C₂H₆:CH₄) from this biogenic source. Consistent with measured concentrations in Fig. 9, CH₄ concentrations from the MWF are much lower in comparison to the sampled thermogenic methane sources.

Figure 11 Ethane and methane enhancement ratios from multiple independent field tests at different sites. “October and November 2022 controlled releases” are two independent multiday field campaigns in Casa Grande, AZ, during a controlled NG release. “San Juan Basin well”: multiple flights near an active wellhead in Cuba, NM. “Permian Basin well”: three independent flights downwind of an orphan well in Hobbs, NM. “Municipal waste facility” is a relatively small waste disposal location near Socorro, NM.

The two different controlled-release campaigns in fall 2022 were about 1 month apart, and these data are segregated in Fig. 11. The measured percentages are in agreement with the contemporaneous CH₄ concentration measurements reported by the controlled-release team (El Abbadi et al., 2024). The same ratio of C₂H₆ to CH₄ was measured across all flights during the first release and similarly during the second campaign, though the ethane content was noticeably different during the second campaign.

Emissions of unrefined NG from the orphan well discussed in Sect. 3.2.2 contained a smaller fraction of CH₄ due to the presence of other compounds such as H₂S, so it is incorrect to approximate the C₂H₆ ratio in Fig. 11 as the non-CH₄ percentage. However, this unique C₂H₆:CH₄ ratio was consistent across all orphan well flights and in good agreement with ground-based systems deployed at the same site (Follansbee et al., 2024); the larger ethane content is due to the unrefined natural gas seeping from the unused well (Stolper et al., 2018). Furthermore, measurements of unrefined NG from a leaking O&NG well in the San Juan Basin, obtained in October 2023, showed a distinct ratio from that of the orphan well from the Permian Basin. The San Juan Basin is primarily a coal-producing region and therefore has a noticeably lower C₂H₆ content compared to the Permian Basin, the latter region being primarily composed of oil. Thus, each of the sources shown in Fig. 11 had a unique and consistent C₂H₆ ratio which can be used to characterize and differentiate multiple sources based on this percentage. This is in agreement with the findings of Meyer et al. (2022) and Germain-Piaulenne et al. (2024).

Figure 12 Comparison of estimated methane flux from various anthropogenic sources. The UAS has a large dynamic range capable of quantifying emission rates from smaller sources (less than 1 kg h⁻¹). Horizontal dashed lines show the absolute lowest detection limits for satellite and aircraft quantification methods (Jacob et al., 2022; Duren et al., 2019) with respect to the estimated lower limit of this system (∼ 0.007 kg h⁻¹; see Sect. 2.7).

Figure 12 highlights the dynamic range of the UAS, which is able to detect CH₄ enhancements of as little at 20 ppb s⁻¹ (after background removal; Sect. 2.3). Additionally, the UAS has demonstrated the ability to measure much larger concentrations of more than 40 000 ppm (empirically determined based on MIRA Pico saturation levels during Arizona (AZ) validation flights).

4 Conclusions

Our results demonstrate the capabilities of this integrated UAS package, along with associated flight strategies and data analysis methodologies, to quantify and characterize methane of point and distributed emission sources. Through direct, in situ measurements of CH₄ concentrations and vector winds downwind of target sources, the UAS can be used to quantify emission rates spanning more than 3 orders of magnitude. The lower-bound limit for CH₄ flux measurements is estimated to be ∼ 0.007 kg h⁻¹, and it is determined primarily by the precision of the concentration measurements and uncertainties associated with determining the background levels with respect to variable plume enhancements.

Environmental conditions largely dictate the plume dispersion and transport with respect to a target source; optimal downwind distances and transect velocities during flight are highly dependent on insolation and wind speed and direction. The results from the controlled releases shown in Fig. 8 and Table B1 were obtained during daytime periods of relatively steady winds and low cloud coverage (high solar insulation leading to large vertical mixing). In all cases, the measured emission rates showed reasonable agreement with metered flow rates, generally within 1σ error. However, the results suggest the potential for a small systematic underestimation of emissions across all scales, depending on environmental conditions. This underestimation may arise from limitations on the horizontal extent or vertical distribution of transects during each flight, resulting in an under-sampling of the continuous emission plume. Therefore, UAS flux quantification may represent an effective lower-bound value for the true emission rate. It is important to note that the results in Sect. 3.1 were informed by the fully unblinded metered rates, whereas the aforementioned aircraft systems used fully blinded or partially unblinded data for the reported estimates (El Abbadi et al., 2024; Sherwin et al., 2024).

The primary limitations of our UAS measurement approach are related to meteorological conditions. Optimum mean wind speeds are in the range of 2–6 m s⁻¹. Wind speeds below 2 m s⁻¹ have been shown to produce less reliable fluxes due to higher variability in the plume position and shape, while winds above 6 m s⁻¹ are too strong for safe flights with this UAS. The steadiness of wind directionality can also be a factor, although it is a much less severe constraint than speed provided that the flight pattern is sufficiently wide to intersect the plume on every transect. Proper deployment also requires access for takeoff and landing within 1 km of the source, and the target source's location must also be known to within a few kilometers. While successful deployment and flux quantification is constrained to a specific range of environmental conditions, the UAS offers a relatively simple method of deployment for quantifying emitters that would otherwise be difficult or impossible for other monitoring techniques to access or detect. Thus, although this system is not optimized for wide-area surveys, it is well suited for site quantification of known sources such as O&NG wells or processing facilities, small dairies, municipal waste facilities, and wastewater treatment plants.

Appendix A: Supplemental information

Figure A1 Instantaneous UAS pitch and roll versus horizontal velocity from all flights presented in Fig. 12 and Table B1 (> 13 500 samples). Plot contains all samples prior to filtering and transect processing as described in Sect. 2.6.

Figure A2 (a) Comparison of vector wind (u) from UAS adjustment (Sect. 2.2) with directional ground anemometer aligned to the west. (b) Height of UAS above ground level (blue) with respect to anemometer (orange). Highlighted regions depict steady and level transects within 20 m horizontal distance from the fixed anemometer.

Appendix B: Supplemental equations

Static-vector wind uncertainty

$$\begin{array}{}\text{(B1)}& {\displaystyle }\mathit{\delta }u=\sqrt{{\left({\displaystyle \frac{\partial u}{\partial u_{\mathrm{m}}}}\right)}^{\mathrm{2}}\mathit{\delta }{u_{\mathrm{m}}}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial u}{\partial \mathit{\varphi }}}\right)}^{\mathrm{2}}\mathit{\delta }{\mathit{\varphi }}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial u}{\partial V_x}}\right)}^{\mathrm{2}}\mathit{\delta }{V_x}^{\mathrm{2}}}\text{(B2)}& {\displaystyle }\mathit{\delta }v=\sqrt{{\left({\displaystyle \frac{\partial v}{\partial v_{\mathrm{m}}}}\right)}^{\mathrm{2}}\mathit{\delta }{v_{\mathrm{m}}}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial v}{\partial \mathit{\varphi }}}\right)}^{\mathrm{2}}\mathit{\delta }{\mathit{\varphi }}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial v}{\partial V_y}}\right)}^{\mathrm{2}}\mathit{\delta }{V_y}^{\mathrm{2}}}\text{(B3)}& {\displaystyle \frac{\partial u}{\partial u_{\mathrm{m}}}}=\cos \left(\mathit{\varphi }\right),{\displaystyle \frac{\partial u}{\partial \mathit{\varphi }}}=-u_{\mathrm{m}}\sin \left(\mathit{\varphi }\right)+v_{\mathrm{m}}\cos \left(\mathit{\varphi }\right),{\displaystyle \frac{\partial u}{\partial V_x}}=\mathrm{1}\text{(B4)}& {\displaystyle }\partial v\partial v_{\mathrm{m}}=\cos \left(\mathit{\varphi }\right),{\displaystyle \frac{\partial v}{\partial \mathit{\varphi }}}=-v_{\mathrm{m}}\sin \left(\mathit{\varphi }\right)-u_{\mathrm{m}}\cos \left(\mathit{\varphi }\right),{\displaystyle \frac{\partial v}{\partial V_y}}=\mathrm{1}\end{array}$$

Table B1 Measured emission rates from various source types shown in Fig. 12. Presented error ranges are the 3σ bounds calculated from sensor uncertainties and background estimation as described in Sect. 2.6.