Testing the near-field Gaussian plume inversion flux quantification technique using unmanned aerial vehicle sampling

Abstract. Methane emission fluxes from facility-scale sources may be poorly quantified, leading to uncertainties in the global methane budget. Accurate atmospheric measurement based flux quantification is urgently required to address this. This paper describes the test of a new near-field Gaussian plume inversion (NGI) technique, suitable for facility-scale flux quantification, using a controlled release of methane gas. Two unmanned aerial vehicle (UAV) platforms were used to perform 22 flight surveys downwind of a point-source release of methane gas from a regulated and flow-metered cylinder. One UAV was tethered to an instrument on the ground, while the other UAV carried an on-board high-precision prototype instrument, both of which used the same near-infrared laser technology. The performance of these instruments from UAV sampling is described. Both instruments were calibrated using certified standards, to account for variability in the instrumental gain factor. Furthermore, a modified approach to correcting for the effect of water vapour applied and is described here in detail. The NGI technique was used to derive emission fluxes for each UAV flight survey. We found good agreement of most NGI fluxes with the known controlled emission flux, within uncertainty, verifying the flux quantification methodology. The lower NGI flux uncertainty bound was, on average, 17 % ± 10(1σ) % of the controlled emission flux and the upper NGI flux uncertainty bound was, on average, 218 % ± 100(1σ) % of the controlled emission flux. These highly conservative uncertainty ranges incorporate factors including the variability in the position of the plume and the potential for under-sampling. While these average uncertainties are large compared to methods such as tracer dispersion, we suggest that UAV sampling can be highly complementary to a toolkit of flux approaches and may perform well in situations where site access for tracer release is problematic. We see tracer release applied to UAV sampling as an effective combination in future flux quantification studies. Successful flux quantification using this UAV sampling methodology demonstrates its future utility in identifying and quantifying emissions from methane sources such as oil and gas infrastructure facilities, livestock agriculture and landfill sites, where site access may be difficult.



Introduction
The global methane budget is subject to significant uncertainties (Kirschke et al., 2013;Saunois et al., 2016b;Nisbet et al., 2019), particularly from inventory uncertainty in facility scale sources such as landfill sites (Scheutz et al., 2009), herds of cattle (Blaxter and Clapperton, 1965) and oil and gas infrastructure (Brantley et al., 2014), which collectively contribute 35 significantly to global methane emissions (Dlugokencky et al., 2011;Saunois et al., 2016a). These uncertainties can be reduced through the accurate source identification and subsequent quantification of methane emission fluxes using top-down (atmospheric measurements based) methods, in order to validate bottom-up (inventory based) emission flux estimates (Lowry et al., 2001;Nisbet and Weiss, 2010;Allen, 2016;Desjardins et al., 2018).
Accurate top-down flux quantification from facility scale sources requires a combination of wind vector measurements along with in situ measurements of atmospheric methane mole fraction (Dlugokencky et al., 1994;Rigby et al., 2017). Facilityscale emission fluxes can be derived from near-field sampling (less than 500 m from the source), which may be acquired from an unmanned aerial vehicle (UAV) platform (Gottwald and Tedders, 1985). UAVs are cheap, versatile and relatively easy to use (Villa et al., 2016), compared to large manned aircraft (Illingworth et al., 2014;Lehmann et al., 2016). They can 45 fly near to source and can be directed automatically using waypoints, to enable even and unbiased spatial sampling (Greatwood et al., 2017;Feitz et al., 2018). There are three principal approaches for measuring methane mole fraction from a UAV in situ: on-board air samples can be collected for subsequent analysis (Chang et al., 2016;Greatwood et al., 2017;Andersen et al., 2018), air can be pumped through a long tube to a sensor on the ground for analysis (Brosy et al., 2017;Wolf et al., 2017;Shah et al., 2019) or air can be analysed live using a sensor mounted on-board the UAV (Berman et al., 50 2012;Khan et al., 2012;Nathan et al., 2015;Golston et al., 2017). Yet, a key limitation to accurate source identification and flux quantification is the precision of methane mole fraction measurements (Hodgkinson and Tatam, 2013). Miniaturised sensors suitable for UAV sampling are emerging (Villa et al., 2016), but high precision lightweight in situ sensors, featuring superior techniques, such as off-axis integrated cavity output spectroscopy, have previously failed to materialise.

55
Some studies have also used UAV remote sensing measurements to derive emission fluxes Yang et al., 2018). However, to our knowledge, only Nathan et al. (2015) have derived methane emission fluxes using UAV in situ measurements. In that study, a UAV with an on-board in-situ low precision sensor (±0.1 ppm at 1 Hz) flew in orbits around a gas compressor station, using mass balance box modelling, with geospatial kriging for interpolation, to derive the emission flux. However this method was not tested for UAV sampling using a (known) controlled release of methane gas, beforehand.

60
It is crucial that novel flux quantification techniques are tested by sampling a controlled flux release, prior to investigating unknown emission sources (Desjardins et al., 2018;Feitz et al., 2018). Shah et al. (2019) were the first to test an in situ flux quantification technique using UAV sampling downwind of a controlled methane release. In that study, a UAV was connected to a high precision methane analyser on the ground using 150 m of tubing. Two-dimensional downwind sampling on a vertical flux plane was used to develop the near-field Gaussian plume inversion (NGI) technique for flux quantification 65 (Shah et al., 2019). Fully manual UAV piloting was employed in this previous study to actively pursue the position of the emission plume on the sampling plane, using mid-flight knowledge of its position. This resulted in calculated emission fluxes that were significantly positively biased compared to known emission fluxes; this represents a source of vulnerability in fully manual UAV sampling, which we address in this work.

70
Here we test the application of the NGI method with unbiased UAV sampling of controlled methane emission sources, by flying two UAVs downwind of the release. In this work, the causes of positive flux bias reported in Shah et al. (2019) were addressed in our sampling strategy, by flying a UAV without prior knowledge of the position of the emission plume. One UAV was connected to a commercially available instrument on the ground and the other carried a lighter prototype on-board instrument (sect. 3). Both instruments were characterised and calibrated (sect. 2). Our modified approach to water vapour 75 correction is also outlined in sect. 2. Sampled data was then used to derive NGI flux uncertainty ranges (sect. 4) for each of 22 flight surveys. In sect. 5 the success of the NGI method is assessed overall and its sampling constraints are summarised. given in units of parts-per-million (ppm) throughout this paper, which are defined here as the number of moles of methane per million moles of air (10 -6 • mol methane mol -1 ), with parts-per-million (ppb) defined as the number of moles of methane per billion moles of air (10 -9 • mol methane mol -1 ). In this section, the ABB Los Gatos Research, Inc. Micro-portable Greenhouse Gas Analyzer (MGGA) and a lighter prototype MGGA (pMGGA), designed for UAV use, are compared and characterised to assess their performance. The technical specifications of both instruments are compared in Table 1. Both instruments use 85 off-axis integrated cavity output spectroscopy (ICOS) to derive simultaneous measurements of methane, carbon dioxide and water mole fraction, from the absorption of a near-IR (1650 nm) laser. The pMGGA uses an additional laser to measure carbon dioxide mole fraction more accurately. Off-axis ICOS techniques reflect a tuneable laser between two mirrors in a high-finesse optical cavity, to obtain high-precision mole fraction measurements (see Paul et al. (2001) and Baer et al. (2002) for further details on off-axis ICOS).

90
The e-folding time of the high-finesse cavity in both sensors was measured here by fitting an exponential decay function to the transition from a high to low mole fraction standard gas (see Table 1 for results). This represents the time taken for 63.2% of the contents of the high-finesse cavity to be replaced. The Allan variance precision of each sensor was also derived at various integration times, by measuring [X] from a dry cylinder of compressed air for at least 17 hours of continuous 95 sampling (see Fig. S1 and Fig. S2 for Allan variance plots). The Allan deviation uncertainty factor (σ AV ), defined here, represents the Allan deviation at the maximum sampling frequency. σ AV for the MGGA and pMGGA are 2.71 ppm (at 10 Hz) and 5.44 ppm (at 5 Hz), respectively, with the Allan deviation at 1 Hz and 0.1 Hz given in Table 1. The optimum Allan variance integration time was also assessed for each sensor ((20±3) s for the MGGA and (70±10) s for the pMGGA); this represents the maximum sampling time before instrumental drift begins to dominate over instrumental noise. 100

Water vapour correction
Raw wet methane mole fraction measurements ([X] 0 ) recorded by each instrument were corrected for the influence of atmospheric water vapour, on mole fraction retrievals. Water vapour influences measurements of [X] for three main reasons O'Shea et al., 2013;Rella et al., 2013). First and most significantly, dilution effects occur, where the bulk presence of water reduces the quantity of methane in the cavity at a given pressure. Second, strong, broad infrared 105 absorption bands of water can lead to interference with the absorption spectrum of methane. Third, pressure broadening of absorption peaks can occur, where collisional interaction between water and methane molecules changes the shape of the methane spectral absorption band. These three effects have a net effect of decreasing [X] 0 in both instruments.
To account for these effects, both the MGGA and pMGGA use an internal retrieval algorithm to derive methane mole 110 fractions, which includes empirically-derived estimates of the pressure broadening coefficient in the presence of water vapour. The instruments output both raw dry mole fraction measurements, which have additionally been corrected for the effect of mole fraction dilution by water vapour, and raw wet methane mole fraction measurements ([X] 0 ) which have not been corrected for this dilution effect (but are still calculated using the same empirically derived pressure broadening is some variability from unit to unit. Therefore, to obtain a more accurate correction for the influence of water vapour on the individual instruments used here, we apply a further post-processing correction factor to the [X] 0 measurements (without the dilution correction) reported by the instruments, using reported measurements of water mole fraction ([H 2 O] (1), where a is the water baseline offset, b is the water baseline coefficient and w is the water baseline decay factor. The plotted data used to fit [H 2 O] 0 is given in Fig. S3 and Fig. S4.
130 a, b and w for both instruments are given in Table 2. w for the pMGGA was assumed to be the same as w for the MGGA.
The same w value was used for the pMGGA as gas was sampled by the MGGA at more [X] 0 dry points, resulting in an improved fit.
Having established a well characterised water baseline, a post-processing water correction factor (ν) was derived by 135 sampling gas from a single cylinder, which was humidified to 9 fixed dew points (from 0 °C to 18 °C), using a dew point  Fig. S6 and Fig. S7). A quadratic fit was applied to both curves, with the intercept forced to unity. The first order coefficient (α) and second order coefficient (β) of the quadratic fit, given in Table 2, were then be used to derive ν 145 using Eq. (2), as a function of [H 2 O]. (2) (1) can be substituted into Eq.
The uncertainty in the water correction was quantified using the water correction residual (R) from Eq. (2), to derive a water uncertainty factor (σ ν ) for each instrument (see Table 2).
N is the total number of residuals. Using the above to correct for the effects of water vapour in the measurement cavity, we can increase the accuracy of [X] 0 measurements. For example our correction was used to increase the MGGA measurement 160 accuracy of [X] 0 (at 2 ppm) by +0.27%, at a humidity of 0.001 mol water mol -1 , and by +1.8%, at a humidity of 0.01 mol water mol -1 .

Calibration
In order to convert [X] 0 into [X], both instruments were calibrated by sampling gas from two cylinders: one contained a low The average offset (C) and offset uncertainty (σ C ) was calculated by taking the average and standard deviation, respectively, of individual offsets, calculated using Eq. (7) and Eq. (8) (Pitt et al., 2016).
A key advantage of this calibration procedure is that the uncertainty in G is well quantified up to [X] high (for stable 185 temperature and pressure) and can be incorporated in the measurement uncertainty. Separate in-field calibrations would be preferable to enhance measurement accuracy, by accounting for variability in temperature and pressure. However there are logistical challenges with in-field calibrations, such as the need for calibration gases and the time required to perform calibrations in dynamic atmospheric temperature and pressure conditions. The laboratory calibrations described here required at least three hours of sampling in order to characterise the variability in G: this may be impractical in field 190 conditions. Therefore the calibration coefficients presented here are useful as they account for variability in G under laboratory conditions.

Methane enhancement and uncertainty
The calibration and water vapour correction procedures described above show that G is almost equal to 1, C is almost equal to 0 (see Table 3) and ν is almost equal to 1 (at 5 ppm methane and a very high 0.01 mol water mol -1 humidity) for both 195 instruments. This means that both instruments record raw [X] 0 measurements with very little systematic error, even when uncalibrated. Thus for most methane measurement purposes, [X] 0 may not need to be corrected. However in this work, G and ν were applied to [X] 0 for optimal instrumental accuracy.
The uncertainty in E (σ E ) can be calculated by combining σ b with the precision and accuracy uncertainty components of [X], using Eq. (11). Precision is characterised by σ AV and accuracy is characterised by σ G and σ ν terms.  Fig. 1. The two adjacent grass fields, in which all UAV sampling took place, belong to a fully operational dairy farm. Methane was released from within the operating site at one of two controlled flux rates (F 0 ), from 0.25 m above ground level (see SI for controlled release details). F 0 was undisclosed during flux analysis, prior to the comparisons shown later in this paper, allowing for blind method testing.
Two adapted DJI Spreading Wings S1000+ octocopter UAVs (labelled UAV1 and UAV2) were used to sample the methane plume on a downwind vertical plane, roughly perpendicular to mean wind direction (see Table 4 for UAV details).

225
Measurements of [X] from both platforms are given in Fig. 2. UAV1 was operated using pre-programmed waypoints and  Table S1 and Table S2.

245
A lightweight wind sensor (FT205EV, FT Technologies Limited) was mounted on-board UAV1, on a carbon fibre pole 305 mm above the plane of the propellers. It recorded wind speed and direction at 4 Hz. This data was used to model the change in wind speed with height above ground level (z). A two-dimensional stationary sonic anemometer (WS500-UMB Smart Weather Sensor, G. Lufft Mess-und Regeltechnik GmbH) was also situated on the southern boundary of the operating site (see Fig. 1), (3.30±0.03) m above ground level. This provided wind speed, wind direction, temperature and pressure 250 measurements every minute. Wind measurements from both sensors were combined to derive the average absolute wind speed as a function of z (WS(z)), for the duration of each flight survey. This is described in detail in the SI.

265
Satellite-derived altitude was corrected to obtain the height of the air inlet above ground level, by taking into account takeoff altitude and the height of the air inlet when on the ground. This step ensures that the data represent the true point of sampling. After converting longitude and latitude from degrees into meters, metric longitude and latitude were projected onto a plane perpendicular to and a plane parallel to mean wind direction. Mean wind direction was derived from the stationary anemometer for the duration of each flight survey. The coordinate projection procedure is described in further 270 detail by Shah et al. (2019).
In order to calculate flux, flux density, q, (in kg s -1 m -2 ) was derived. To achieve this, each geospatially mapped E measurement was combined with WS(z), using Eq. (12).
Geospatially mapped q, on a plane perpendicular to mean wind direction, for each flight survey, is given in Fig. 3 for UAV1 and in Fig. 4 for UAV2. Measurements of [X] (see Fig. 2) were not used in the flux analysis, but are nevertheless of interest, as they show [X] to generally reduce with z, as expected, to support observations of enhancements in q shown in Fig. 3 and  Fig. 4 can also be seen at a much greater lateral distance from the source, projected onto the sampling plane. This is likely a consequence of many UAV2 flight surveys sampling at a greater distance from the source than UAV1 flight surveys, which gave the plume 290 more time to spread out.

Flux quantification
Calculated flux density (q) from each flight survey was used to derive an emission flux (in units of kg s -1 ) using the nearfield Gaussian plume inversion (NGI) flux quantification technique (see Shah et al. (2019)). In principle, the NGI method accounts for turbulent variations in wind using Gaussian statistics. The method also takes into account sampling on a slightly 295 offset sampling plane (compared to the plane perpendicular to mean wind direction) by introducing a third dimension to the traditional two-dimensional Gaussian plume model. The NGI method uses a least-squares approach to compare measured and modelled values of q. Residuals in q are minimised to output model parameters, which include an initial flux estimate (F e ).

300
Full details of the NGI method can be found in Shah et al. (2019). We provide a brief overview here. The size of the plume is assumed to increase linearly with distance from the source, by assuming q to decrease according to the inverse square law with distance (an assumption which is valid over short distances). Therefore instead of using constant crosswind and vertical dispersion terms, these terms are allowed to increase with distance from the source, with both terms being fixed at a one metre distance. The crosswind dispersion term (at 1 m) is characterised using measurements of q, rather than assumptions of 305 atmospheric stability, as these assumptions are valid for plumes characterised by dispersion, rather than turbulent advection.
In addition, the centre of the plume in the crosswind direction is derived from measurements of q, as the precise position of the source may be unknown. The vertical dispersion term (at 1 m) and F e can then be acquired by inverting modelled values of q, derived by minimising residuals, as described above. and limited spatial sampling extent. This simulation step therefore gives an important indication of the systematic error due to potential under-sampling. All F e , σ -, σ + and σ F values for each flight survey are given in Table S5.

4 Flux results and discussion
Calculated NGI emission fluxes were compared to the known (controlled) emission fluxes, using the ratio between the NGI flux uncertainty range and F 0 (see Fig. 5). As this was a blind flux analysis, F 0 was not revealed to the analysis team researchers prior to calculating the flux uncertainty range. Fig. 5 shows that the NGI flux uncertainty range agrees well with

325
Although no flux uncertainty range exceeded F 0 , T2.3 spanned a large flux range, much of which fell above F 0 . Flux underestimation may be explained using the plots shown in Fig. 3 and Fig. 4, which demonstrate the following: a limited sampling duration made it possible to almost entirely avoid the emission plume, thus resulting in low flux results; similarly, some flights intersected the emission plume multiple times resulting in flux overestimation in cases, although large NGI uncertainty ranges can conservatively account for this effect. Therefore it is clear that the F e value obtained using the NGI 330 method must not be taken at face value and the full NGI flux uncertainty range must be considered. Furthermore, the flux ranges in Fig. 5 represent uncertainty bounds of one standard deviation; it is statistically realistic to expect some discrepancy between F 0 and the NGI flux uncertainty range.
The flux uncertainty ranges given in Fig. 5 are asymmetric, although the magnitude of this asymmetry was different for 335 flight experiments conducted by the different UAVs. σ + was (0.33±0.13(1σ)) times larger than σfor UAV2 but was only (0.06±0.03(1σ)) times larger for UAV1. This is because UAV2 sampled further from the source, on average, on a similar sized sampling plane to UAV1. As UAV2 was further from the emission source, the instantaneous plume had a greater likelihood of extending beyond the sampling plane and being missed (beyond the horizontal edges of the sampling plane), due to spatially limited sampling. This potential loss of in-plume sampling may have otherwise contributed towards the 340 overall flux, thus enhancing σ + . Therefore σ + is comparatively larger than σfor flights conducted by UAV2.
It is important to recognise the magnitude of the NGI uncertainty ranges in Fig. 5, relative to F 0 , which are due to the difficulties in inverting sparse spatial sampling to derive an emission flux, following the NGI method. These uncertainties uncertainty ranges in Fig. 5 are large, we believe that the true value of the NGI method with UAV sampling is to derive snap-shot rapid flux estimates at low cost, with an order-of-magnitude level precision, for subsequent flux investigation using more precise flux quantification techniques. Although longer periods of sampling in each flight survey may reduce the uncertainties in Fig. 5, this is practically difficult with limited UAV battery life, with little additional benefit. Tethered power or multiple UAV flights may alternatively be used, as was the case with UAV1, but wind conditions can quickly change 350 when sampling for prolonged periods with too many lengthy intervals between flights.
In order to assess whether multiple flight surveys could be used effectively to capture the known controlled emission flux, within uncertainty, the upper and lower NGI uncertainty bounds were averaged for all surveys (see penultimate row of Fig.   5). The average lower NGI flux uncertainty bound as a fraction of F 0 (F -̅ ) was 0.2±0.1(1σ) and the average upper NGI flux 355 uncertainty bound as a fraction of F 0 (F + ̅̅̅ ) was 2±1(1σ), for all surveys. Thus F 0 (i.e. 1 in Fig. 5)  are large, we emphasise that our methodology has been adapted for rapid flux analysis, rather than precise flux estimates for inventory publication.

365
The ability of the NGI method to calculate a target emission flux was further assessed by calculating the central flux estimate as a fraction of F 0 (F c ) for each flight survey, using Eq. (13). F c is distinct from F e (as a fraction of F 0 ), in that F c finds the centre of an asymmetric flux uncertainty, whereas F e is an initial flux estimate calculated using the NGI method, which does not take into account the potential effects of under-sampling, which may result in a potential negative flux bias. To conclude, UAV sampling can be used to practically derive unbiased snap-shot emission fluxes with the NGI method, with an order-of-magnitude precision, by sampling on a plane perpendicular to wind direction from at least approximately 50 m 390 away from the source. Although typical flux uncertainties were high, NGI UAV fluxes serve as an important tool for snapshot source identification and flux quantification. Our UAV methodology fills an important gap between cheap leak detection techniques (such as infrared cameras), which do not provide fluxes, and reliable flux quantification techniques (such as the tracer dispersion method), which require expensive instrumentation and may be more difficult to organise. For example, tracer methods can be problematic in cases where site access for tracer release is impossible or in cases where the 395 plume may be lofted. The UAV methodology we describe is highly suitable for leak detection and source isolation, for regulatory leak detection, with the added capability to gauge the severity of flux leaks, for subsequent investigation using other approaches. We anticipate a combination of UAV sampling with a tracer release, where both a target gas (in this case methane) and a proxy tracer can be measured simultaneously downwind, taking advantage of vertical sampling enabled by UAVs, as a powerful future toolkit for precise flux quantification.

5 Conclusions
Two UAVs were used to test the near-field Gaussian plume inversion technique for flux quantification. One UAV was connected to the MGGA on the ground using a tether, while the other carried a new ABB pMGGA prototype instrument onboard. Both instruments measured atmospheric methane mole fraction, which was calibrated and corrected for the influence of water vapour, following laboratory testing.

405
The flux approach was tested for 22 UAV flight surveys, by deriving fluxes from a controlled release of methane gas. This  Table S3