A near-field Gaussian plume inversion flux quantification method, applied to unmanned aerial vehicle sampling

: The accurate quantiﬁcation of methane emissions from point sources is required to better quantify emissions for sector-speciﬁc reporting and inventory validation. An unmanned aerial vehicle (UAV) serves as a platform to sample plumes near to source. This paper describes a near-ﬁeld Gaussian plume inversion (NGI) ﬂux technique, adapted for downwind sampling of turbulent plumes, by ﬁtting a plume model to measured ﬂux density in three spatial dimensions. The method was reﬁned and tested using sample data acquired from eight UAV ﬂights, which measured a controlled release of methane gas. Sampling was conducted to a maximum height of 31 m (i.e. above the maximum height of the emission plumes). The method applies a ﬂux inversion to plumes sampled near point sources. To test the method, a series of random walk sampling simulations were used to derive an NGI upper uncertainty bound by quantifying systematic ﬂux bias due to a limited spatial sampling extent typical for short-duration small UAV ﬂights (less than 30 min). The development of the NGI method enables its future use to quantify methane emissions for point sources, facilitating future assessments of emissions from speciﬁc source-types and source areas. This allows for atmospheric measurement-based ﬂuxes to be derived using downwind UAV sampling for relatively rapid ﬂux analysis, without the need for access to di ﬃ cult-to-reach areas.


Introduction
Methane is the second most important greenhouse gas in terms of radiative forcing [1], with a direct radiative forcing of +0.61 W m −2 [2]. Though a global mean methane concentration up to 800,000 years prior to 1850 was constrained between 300 ppb and 800 ppb [3], recent global mean measurements (2018) have exceeded 1850 ppb [4]. Following a hiatus at the turn of this millennium, a was used. However, as near-field in-situ sampling was used in the Nathan et al. [73] study, the use of kriging as a form of interpolation could be subject to error due to effects of near-field turbulence. Yang et al. [74] used a downward facing remote sensing instrument attached to a small UAV platform to detect gas leaks, using a new interpretation of mass balance box modelling tailored for remote sensing sampling. Their UAV performed concentric circles around the presumed emission source at a fixed height [74]. Mass balance box modelling can be a highly suitable method for UAV sampling, as shown in a study by Yang et al. [74], whose method integrated the statistical average of remote sensing measurements from above an emission plume.
In this paper, an NGI mass balance method (see Section 2) is developed and tested using a preliminary sample dataset of geospatially mapped in-situ [CH 4 ] measurements obtained from eight UAV flights. The NGI method used here was specially adapted to allow fluxes to be quantified using near-field UAV sampling downwind of (assumed) point sources by tethering the UAV to a high-precision in-situ methane sensor on the ground. An advantage of UAV sampling is that it does not require access to difficult-to-reach areas (for example, access to industrial or restricted sites). In this study, controlled releases of methane were conducted by the UK National Physical Laboratory (NPL) for evaluation. The refinement of flux quantification techniques using known controlled emission fluxes, such as those provided by the NPL in this study, is crucial in order to legitimise the integrity of further work [21,25,34,35,74]. The NPL fluxes were also used to test simultaneously-derived TDM fluxes (described in detail in Fredenslund et al. [75]) by two TDM teams operating independently of the UAV team and each other. In Section 3, we then test the NGI method and account for inherent flux biases associated with the method, using simulated measurements acquired from a random walk simulation sampling a static plume. These simulations also serve as further guidance on improvements to the sampling approach in future use of the method.

Method
The NGI method is a mass continuity model in principle, applied to near-field sampling, where in-situ [CH 4 ] measurements upwind and downwind of an emission source are differenced and combined with a wind vector, resulting in a measured flux density (q me ). This makes the NGI method suitable for near-field UAV sampling. Provided that upwind air is relatively homogenous in comparison to downwind enhancements and well mixed (regardless of whether it is enhanced in methane overall), the NGI method can be used to difference between upwind and downwind methane sampling. The NGI method requires near-to-source sampling (defined here to be~100 m from the source) on a vertical two-dimensional plane roughly perpendicular to mean wind direction (θ), to define a net source flux into a volume of air advecting with the wind. The [CH 4 ] value upwind of the source (background) is assumed to be roughly constant compared with enhancements expected from the source of interest. This background variability can be measured and accounted for in error budgeting (see later). This assumption requires there to be no nearby interfering emission sources that are detectable from the point of sampling. Additionally, loss terms between the source and sampling plane are assumed to be negligible, as is the case for methane as its atmospheric lifetime (~9 years [9]) significantly exceeds the timescale of near-field advection downwind of a facility-scale source. For highly reactive trace gases, such considerations may become important if this method was adapted.
Under turbulent conditions, near-field sampling may not map out a characteristic conceptual time-invariant Gaussian plume for constant sources. Instead, a non-Gaussian distribution of q me may be observed across the flux plane due to sampling different positions of an instantaneous plume, sampled at different times. However, provided that sampling on the vertical flux plane captures the centre of a time-averaged plume and extends sufficiently far outwards, a flux can be derived by assuming that the time-averaged plume also exhibits a Gaussian plume morphology. Sampling duration should also be sufficient to ensure unbiased measurement of the flux plane. This is the basis of the NGI approach, applicable to near-field UAV measurements of point source emissions. In the later discussion, we attempt to assess and account for the limitations of such sampling and examine sources of bias and error in the method using simulations.

Preparing the Sample Data
A DJI S900 hexacopter UAV (flying up to (24 ± 3) m above ground level) was used to obtain [CH 4 ] measurements using a Los Gatos Research, Inc. Ultra-portable Greenhouse Gas Analyzer (UGGA) on the ground, sampling air at a rate of 1 Hz, which was connected using a 150 m length of PFA tube to an inlet above the plane of the UAV rotors (see Appendix A for further experimental details including the sampling location, instrumental set-up, geolocation measurements and time corrections). These data served as a sample dataset of geospatially-mapped in-situ [CH 4 ] measurements, with corresponding measurements of altitude, latitude, and longitude, downwind of the NPL Controlled Release Facility (CRF). Eight flights were conducted for between 13 and 28 min each (see Table S2 for individual flight details). Each flight was conducted with the intention that it would be centred around the time-averaged plume as projected onto the flux plane and was adapted in-flight to capture the maximum extent of the plume. The UAV was directed towards the time-averaged plume mid-flight to maximize in-plume sampling based on mid-flight knowledge of its position estimated from smoke flares released near to the emission source. Such targeted sampling can result in a positive flux bias that (with hindsight) invalidates the assumptions of mapping a time-averaged plume.
Altitude, longitude, and latitude corresponding to each [CH 4 ] measurement were converted into height above ground level (z), the horizontal distance perpendicular to mean wind direction (y) and the distance from the source along the mean wind vector (x). To calculate z, the lowest in-flight altitude measurement was set to the height of the air inlet above the base of the UAV. Thus, the systematic accuracy of geospatial height measurements is not important, provided that they are of a sufficient precision (random uncertainty). To calculate x and y, longitude and latitude were converted into distances in metres (assuming a spherical Earth). These distances were then projected onto a sampling plane parallel to mean wind direction (θ) to calculate x and onto a plane perpendicular to θ to calculate y for the duration of sampling. The assumed y position of the emission source projected onto the sampling plane was set to 0 m as a reference, and the assumed x position of the source was set to 0 m. Figure 1a gives an example of a UAV flight track, with geospatially mapped measured [CH 4 ] as a function of x, y and z ([CH 4 ] me ) on the flux plane (see Figure A3 for all UAV flight tracks). The upwind background methane concentration ([CH 4 ] 0 ) and corresponding background uncertainty (σ b ) was assumed here to be constant for each flight. [CH 4 ] 0 was found by taking the peak of a log-normal distribution applied to a histogram of the lowest [CH 4 ] measurements (see supplement for full details). All sampling uncertainties in this paper are expressed as single standard deviations.
Anemometers positioned at various heights were used to infer continuous wind speed (U(z)) and wind direction profiles, along with corresponding wind speed variability (σ U (z)) profiles (see Appendix A for further details of the anemometers used). A single mean value of θ was derived for each flight: This was calculated from the average of two-dimensional wind vectors at the height of each [CH 4 ] me measurement. The atmospheric methane mass density (ρ), assuming all of air is methane (in units of kg m −3 ), was calculated using Equation (1), where M is the molar mass of methane (0.0164 kg mol −1 ), R is the gas constant (8.314 m 3 Pa mol −1 K −1 ), and T and P are the average temperature and pressure, respectively, for the duration of each UAV flight. Provided that changes in T and P are small, ρ can be assumed to remain constant for the duration of sampling.

Flux Estimation Using the NGI Method
Here, we describe the procedure applied to each UAV flight from the sample data in order to estimate an initial flux using the NGI method. The qme, as a function of y, z, and x, was calculated using Equation (2), assuming constant ρ. The qme is the emission flux per unit area on the flux plane, derived from each individual [CH4]me measurement (see Figure 1b for flight 2, for example). As Equation (2) requires the concentration enhancement over the background, instrumental resolution  4 ] me measurements obtained during flight 2 from the sample data (117 m from the emission source). z is plotted on the vertical axis and y is plotted on the horizontal axis; (b) q me values calculated using [CH 4 ] me values in panel "a"; (c) q mo values modelled on q me values in panel "b." Panel "b" and panel "c" share a mutual colour scale; (d) Residuals between q me values in panel "b" and q mo values in panel "c."

Flux Estimation Using the NGI Method
Here, we describe the procedure applied to each UAV flight from the sample data in order to estimate an initial flux using the NGI method. The q me , as a function of y, z, and x, was calculated using Equation (2), assuming constant ρ. The q me is the emission flux per unit area on the flux plane, derived from each individual [CH 4 ] me measurement (see Figure 1b for flight 2, for example). As Equation (2) Atmosphere 2019, 10, 396 7 of 25 requires the concentration enhancement over the background, instrumental resolution (1 ppb for the sample data) can be a limiting factor in the magnitude of fluxes detectable: If fluxes are too low, concentration enhancements may be undetectable.
Equation (3) was used to calculate the uncertainty in q me (∆q me ), using the instrumental Allan deviation precision (σ AV , specifically 0.00127 ppm for the sample data) and an instrumental uncertainty factor (σ i ([CH 4 ])) associated with each [CH 4 ] me measurement (see Appendix A for details on σ i ([CH 4 ])).
The NGI method assumes that spial variability in the time-averaged plume is Gaussian. As a consequence, the NGI method also assumes that the plume is not capped in the z direction by an atmospheric temperature inversion. Conceptually, over an infinite period of randomised sampling, and assuming that the prevailing meteorology does not vary, a Gaussian morphology of q me would be observed as an average of instantaneous sampling. This essentially smooths out the manifestation of turbulence in the time-averaged plume morphology. However, in practice, limited sampling time constraints (particularly for UAV sampling) may potentially lead to the observation of non-Gaussian q me across the sampling plane due to turbulence (see Figure A3 for the example of the sample data). The NGI method is therefore a trade-off between the need to sample for a long enough time to characterise the time-averaged plume as much as possible and the need to minimise the systematic effect of any potential dynamic changes in atmospheric stability, source strength and prevailing meteorology over the course of measurement (parameters which may be derived from complementary specialised UAV sampling [57]). Such competing needs would need to be considered on a case-by-case basis by considering weather forecasts, the local topography and the nature of the expected sources being studied.
An initial flux estimate (F e ) can be derived by fitting q me values to corresponding modelled flux density (q mo ) values using an adaptation of the two dimensional Gaussian plume model [76] given by Equation (4). Equation (4) does not include wind speed, as this has been incorporated into q me in Equation (2). An advantage of incorporating wind speed into q me in this way is that variations in wind speed with altitude are incorporated into the measured plume morphology in flux space. An advantage of this method is that it does not require us to assume any stability class, as plume spatial variance is implicitly measured.
σ y (x) is the crosswind mixing factor, σ z (x) is the vertical mixing factor, y c is centre of the plume in the y direction (based on the actual sampled q me measurements) and h is the emission height.
On an idealised downwind vertical flux plane, perpendicular to θ, σ y (x) and σ z (x) would be constant for all q me values. However, as σ y (x) and σ z (x) are both functions of x, they may vary depending on the orientation of the sampling plane, if not perfectly perpendicular to θ. In reality, it is difficult to sample on a perfect perpendicular plane, especially with UAV sampling (see Figure S2 for the example of the sample data), as θ may vary throughout the sampling period. To account for this variability in sampling distance from the source, we assume σ y (x) and σ z (x) to increase linearly with x (i.e., q me decreases with increasing distance from the source). Though σ y (x) and σ z (x) may not increase linearly with x over large distances, they can be approximated to be linear functions of x under near-field conditions and where the variation in x remains small (less than 40 m). σ y (x) and σ z (x) for each stability class were plotted up to 500 m, using the equations found in Turner [76], to test the validity of this assumption (see Figures S4 and S5). This assumption of linearity allowed us to define a crosswind mixing factor evaluated at a distance of 1 m (τ y ) using Equation (5) and a vertical mixing factor evaluated at a distance of 1 m (τ z ) using Equation (6). τ y and τ z are both then scalable factors that can be used to derive σ y and σ z as a function of x and are therefore applicable to the variable sampling typical of UAV measurement.
As τ y and τ z are both independent of x, Equation (4) can now be expressed as a three dimensional Gaussian plume, using Equation (7).
One could attempt to solve Equation (7) iteratively for the four unknown constant parameters (F, y c , τ y and τ z ). However, such a solution is not always well constrained by the sampled data. Instead, by utilising the separability of Equation (7), we can place constraints on y c and τ y by equating the mean and standard deviation (in the y dimension) of q mo to the equivalent values for q me . This requires us to remove the vertical (z) component of q me using Equation (8) to define q me with the z Gaussian component removed (q me,y ). As Equation (8) depends on τ z , good spatial sampling in the z direction is required in order to derive τ z .
y c and τ y were then derived analytically from measured data using Equations (9) and (10), respectively, which require each "j" value of z, x and q me,y . Foster-Wittig et al. [39] used an equation similar to Equation (10) to derive σ y (x), rather than τ y , using downwind concentration measurements acquired at a single position on the ground. As τ z in Equation (8) was unknown, Equation (9) and (10) were solved simultaneously along with Equations (7) and (8) as part of the model. A key advantage of calculating y c in this way is that, provided geospatial measurements of a sufficient precision are recorded, the overall systematic accuracy of geospatial sampling is not crucial, as the centre of the plume is derived from spatial measurements relative to each other.
Equations (9) and (10) require good spatial sampling in the y direction, as y c weights the flux contributions to the position of the centre of the measured plume (as sampled), independent of z. This means that if q me at any point in the y direction was equal to 0, it would have no influence in changing the centre of the calculated plume morphology in the y direction. By calculating y c and τ y in this way we constrain the model and optimise the flux calculation procedure, such that the only remaining unknown variables in Equation (7) are τ z and F e .
The initialised values of τ z and F e were then evaluated computationally to minimise the sum of squared residuals between all q me and q mo values by solving Equations (7), (9), and (10) simultaneously (as y c and τ y are both derived from values of q me,y , which ultimately depend on τ z ). This was achieved here using the lsqcurvefit function of Matrix Laboratory (MATLAB) 2016a, which applies a least-squares solution to the input data. During the optimisation procedure τ z and F e were both constrained to ensure good parameterisation, as described below, without which the model may produce unrealistic flux results. It was clear that the maximum constraining value of τ z (τ z,max ) had an impact on F e : If τ z,max was too small, τ z may have been equal (or almost equal) to τ z,max (as illustrated in Figure A4 for example). Furthermore, F e sometimes failed to stabilise as τ z,max increased, regardless of how high τ z,max was set to (see Appendix C for further discussion). Therefore, before a final constant F e result could be produced, the flux evaluation procedure had to be repeated at various values of τ z,max . A constant value of F e was assumed to be produced once incremental changes in τ z,max no longer resulted in changes in F e . Figure 1c gives an example of final constant modelled q mo values for flight 2 from the sample data. Final τ z and F e values were then substituted into Equations (9) and (10) to establish definite values of y c and τ y . All resulting y c , τ y and τ z values for the sample data are given in Table S4.
Once constant values of τ z and F e were acquired for each flight, a measurement uncertainty in F e (∆F m ) could be calculated through Equation (11), by using each "j" value of q me and ∆q me from Equation (3), assuming each ∆q me term is uncorrelated.
A total minimum uncertainty (∆F − ) was then calculated using Equation (12) from residuals between each "j" calculation of q me and q mo (see Figure 1d for flight 2, for example). Equation (12) finds the average q me residual, weighted by q me .
The residuals between q me and q mo in Equation (12) incorporate uncertainty in wind direction. Under stable wind conditions, at a constant wind direction, and with infinite spatial sampling, the residual (q me j −q mo j ) term in Equation (12) would converge to zero. However, under turbulent (or variable) wind conditions, regardless of sampling spatial extent and sampling duration on the flux plane, the residual (q me j −q mo j ) term may be non-zero due to the potential for overlap of various q me measurements at the same point in space at different times. Furthermore, the limited extent of spatial sampling on the flux plane can lead to a flux under-estimate, which must be incorporated into the total maximum uncertainty (∆F + ). This is described in Section 3.

Flux Sensitivity Random Walk Simulations
The capability of the NGI method to derive a prescribed modelled point source target flux (F t ) was tested here by simulating randomised near-field sampling of a time-averaged plume. As well as evaluating ∆F + , these simulations can be used for guidance to better plan future sampling, to yield fluxes within some threshold of acceptable uncertainty in operational settings. The simulations were designed to replicate real sampling in the y, z, and x directions (see Figure S6 for an example). Each simulation was based on parameters retrieved when calculating F e from the NGI model (y c , h, τ y and τ z ), the orientation of the q me sampling plane in the x direction and the spatial distribution of q me sampling. Each simulation was constrained by the maximum bounds of q me sampling.
Each simulation was initialised from a random y and z position on the sampling plane (point 1 in Figure 2, where Figure 2 is a schematic diagram of a random walk). The distance of each step in the random walk was constant (on a plane perpendicular to x), with a length (L) equal to the average distance between successive real measurement points. In order to replicate changes in sampling direction on the plane perpendicular to x, an exponential probability distribution was fitted to the absolute change in angle between each successive measurement step (see insert in Figure 2). This probability distribution was used to model the change in angle (in magnitude but not in sign) between each successive step in the random walk (δθ), from an initial random direction (θ R ). At the sampling plane boundary, the random walk direction was deflected by an angle of either 135 • or 225 • , provided that the random walk remained within the spatial confines of sampling. The orientation of the simulated sampling plane in the x and y direction, relative to mean wind direction, was chosen to match the orientation of the sampling plane from real sampling (see supplement for details). Simulated q me values were then generated for each step in the random walk using Equation (7), where F e was replaced by F t and q mo was replaced by q me .
Atmosphere 2019, 10, x FOR PEER REVIEW 10 of 25 was chosen to match the orientation of the sampling plane from real sampling (see supplement for details). Simulated qme values were then generated for each step in the random walk using Equation (7), where Fe was replaced by Ft and qmo was replaced by qme. Figure 2. A schematic of the random walk simulation with an example δθ exponential probability distribution function (pdf) in the top right hand corner. Each step in the random walk is numbered from 1 to 4. Upon reaching the plane boundary at point 4, a deflection of 225° was applied, since a deflection of 135° (red square) failed to keep the random walk within the confines of spatial sampling. The sampling plane spans from the minimum real y measurement (ymin) to the maximum real y measurement (ymax) and from the minimum real z measurement (zmin) to the maximum real z measurement (zmax).
For each test described in the following sections, multiple independent random walks (and simulated qme values) were generated from scratch, for either a single fixed simulated sampling time (t) or for various values of t. For each random walk in each test, t was adjusted by altering the total number of steps in the random walk. Each individual random walk sampled the flux plane following a different route, making each random walk unique. A virtual flux estimate (Fe, v) was calculated from each random walk (which is a function of t), following the NGI method. As modelled fluxes can fail to stabilise when τz, max is equal to model output value of τz, five conditions were put in place to ensure that the model had stabilised. These conditions are described in detail in Appendix B. Figure 2. A schematic of the random walk simulation with an example δθ exponential probability distribution function (pdf ) in the top right hand corner. Each step in the random walk is numbered from 1 to 4. Upon reaching the plane boundary at point 4, a deflection of 225 • was applied, since a deflection of 135 • (red square) failed to keep the random walk within the confines of spatial sampling. The sampling plane spans from the minimum real y measurement (y min ) to the maximum real y measurement (y max ) and from the minimum real z measurement (z min ) to the maximum real z measurement (z max ).

Upper Flux Uncertainty Bounds
For each test described in the following sections, multiple independent random walks (and simulated q me values) were generated from scratch, for either a single fixed simulated sampling time (t) or for various values of t. For each random walk in each test, t was adjusted by altering the total number of steps in the random walk. Each individual random walk sampled the flux plane following a different route, making each random walk unique. A virtual flux estimate (F e,v ) was calculated from each random walk (which is a function of t), following the NGI method. As modelled fluxes can fail to stabilise when τ z,max is equal to model output value of τ z , five conditions were put in place to ensure that the model had stabilised. These conditions are described in detail in Appendix B.

Upper Flux Uncertainty Bounds
The difference between F t and simulated flux (F e,v ) was assessed at a fixed simulated sampling time corresponding to the duration of actual sampling for each flight in the sample data. Fluxes for this analysis were derived by continuously repeating a random walk simulation to accumulate a large number (180 for the sample data) of virtual fluxes (see green lines in Figure 3a for example). The average value of F e,v (F e, v , solid green line in Figure 3a) was negatively biased from F t due to three independent effects.

1.
The simulated sampling extent was restricted as a consequence of the managed sampling strategy, resulting in an inadequate characterisation of the entire emission plume.

2.
Sampling was physically restricted in the z direction due to the non-zero height of the air inlet (as was the case for our UAV platform), resulting in an under-sampled area very close to the ground. 3.
The random walk was time-limited, resulting in an incomplete exploration of the available flux plane and hence, sampling gaps, leading to a residual negative flux bias.
These effects collectively result in a lack of spatial sampling across the flux plane and thus insufficient sampling with which to calculate τ y and τ z , which are derived here from measurements of q me across the sampling plane. In theory, these effects can be alleviated with infinite spatial sampling on an infinitely sized sampling plane, but this cannot realistically be achieved from a UAV (or any) sampling platform. These three effects were incorporated into ∆F + by adding the deviation of F e,v from F t to ∆F − , using Equation (13).
Thus, Figure 3a shows that as a consequence of the three effects outlined above, the NGI method can result in a large variation in flux estimation with the potential for a large negative flux bias.

Uncertainty Sampling Thresholds
By gradually increasing t in each simulation, the convergence of F e,v towards F t was assessed (blue crosses in Figure 3b, for example). This makes it possible to analyse the limiting uncertainty that exists in F e,v as a function of t, due to the random error associated with the NGI method itself rather than due to the quality of sampling. This test is only valid for relatively non-variable atmospheric conditions, as it assumes the time-averaged plume to be static. For each flight from the sample data, independent random walks of up to four hours were generated. A simulated scale factor as a function of t was applied to these fluxes to correct for the inherent negative bias that occurs as a function of t (see previous section). This scale factor was derived by repeatedly performing random walks at a range of fixed sampling times to quantify the average negative flux bias as a function of t (see supplement for details). Without applying this scale factor, F e,v gradually increased towards F t , as a function of t (see Figure S7 for example). Following the application of this scale factor, the absolute difference between bias-corrected F e,v and bias-corrected F e, v was used to find the flux anomaly (A) as a function of t (red diamonds in Figure 3b, for example). The decrease of A with increasing t represents a reduction in the fundamental methodological uncertainty due to a greater understanding of the time-averaged plume morphology. The flux anomalies in Figure 3b show that with a sufficient sampling duration, the fundamental methodological uncertainty can gradually be reduced until it converges towards zero. of Fe, v from Ft to ΔF−, using Equation (13).
Thus, Figure 3a shows that as a consequence of the three effects outlined above, the NGI method can result in a large variation in flux estimation with the potential for a large negative flux bias.  The flux anomaly was assumed to decrease exponentially with t, from which a fit could be produced using Equation (14) to define an anomaly coefficient, α, and an anomaly decay constant, β (see Table S5 for values for the sample data). This fit was weighted to take into account fewer simulations at longer t, such that the fit was evenly applied as a function of t.
By setting t in Equation (14) to the duration of each UAV flight in the sample data and by setting F t to F e , the limiting uncertainty (A f ) associated with each flight could be derived. A f represents the inherent limiting uncertainty, determined by the sampling strategy (assuming non-variable atmospheric conditions), which decreases as a function of t. Alternatively, Equation (14) could be rearranged to make t the subject. Then, by setting A to a specified target fraction of F t , the minimum time required to surpass a target limiting uncertainty threshold could then be derived. This may be useful to decide future UAV sampling durations within a desired minimum confidence level. For example, setting A to 1% of F t yielded the 1% fundamental uncertainty sampling threshold (t 1% ). The use of limiting confidence timescales is only valid if non-variable atmospheric conditions are assumed. Therefore, the threshold time must be comparably smaller than the timescale of prevailing meteorological change.

Testing the NGI Method
Finally, virtual flux simulations were used to test the NGI method in principle. First, the flux plane was extended to mitigate biases due to the limited size of the sampling plane. The maximum sampling extent in the z direction was doubled, and the left and right x bounds were also doubled (resulting in an increase in the size of the original sampling plane by a factor of approximately 6). Furthermore, the flux plane was extended down to ground level. The total simulated sampling time (t) was then set to an arbitrary (approximately 72 h) sampling time across the extended flux plane to mitigate the effects of limited spatial sampling coverage, as well as to minimise the limiting uncertainty discussed in the previous section. The simulation was then repeated 60 times for each flux. Average F e,v results deviated from F t by (−0.7 ± 0.5)%, (see Table S6 for results using the sample data). This test validates the utility of the NGI method by confirming that it accurately yields F t , in principle.

Results and Future Guidance
The sample data used to develop the NGI flux quantification method described here was acquired using the controlled release of natural gas (supervised by the NPL) at two flux rates ((1.5 ± 0.1) g s −1 and (3.0 ± 0.2) g s −1 ). The NGI fluxes for each of eight UAV flights are displayed in Figure 4, along with flux uncertainty bounds, expressed as a percentage of the NPL controlled release facility (CRF) reported emission flux (see Table S1 for individual flux values for each flight). All values of F e , ∆F + , and ∆F − are given in Table S7. As this was a preliminary study, regardless of uncertainty, we note that most NGI fluxes calculated for the sample data exceeded CRF fluxes due to the biased nature of our sampling. Though potential sources of negative bias due to sampling limitations are incorporated into the NGI flux uncertainty range, the effects of positive bias are clearly dominant here. One such source of positive bias is the targeted sampling strategy described in Section 2.1, resulting in plume chasing. The NGI method inherently assumes that the measurements constitute an unbiased sample of the time-averaged Gaussian plume. However, during the acquisition of the sample data, the UAV was directed towards the perceived position of the instantaneous plume to maximise in-plume sampling.
Future studies adopting this technique should avoid chasing the plume in this way, as it results in the over-representation of in-plume measurements, which biases NGI fluxes. Another possible source of bias may be wind measurements; if the winds at the location of the anemometers were systematically different from those at the UAV sampling locations, this would have a direct impact on the calculated fluxes.

Results and Future Guidance
The sample data used to develop the NGI flux quantification method described here was acquired using the controlled release of natural gas (supervised by the NPL) at two flux rates ((1.5 ± 0.1) g s −1 and (3.0 ± 0.2) g s −1 ). The NGI fluxes for each of eight UAV flights are displayed in Figure 4, along with flux uncertainty bounds, expressed as a percentage of the NPL controlled release facility (CRF) reported emission flux (see Table S1 for individual flux values for each flight). All values of Fe, ΔF+, and ΔF− are given in Table S7. As this was a preliminary study, regardless of uncertainty, we note that most NGI fluxes calculated for the sample data exceeded CRF fluxes due to the biased nature of our sampling. Though potential sources of negative bias due to sampling limitations are incorporated into the NGI flux uncertainty range, the effects of positive bias are clearly dominant here. One such source of positive bias is the targeted sampling strategy described in Section 2.1, resulting in plume chasing. The NGI method inherently assumes that the measurements constitute an unbiased sample of the time-averaged Gaussian plume. However, during the acquisition of the sample data, the UAV was directed towards the perceived position of the instantaneous plume to maximise in-plume sampling. Future studies adopting this technique should avoid chasing the plume in this way, as it results in the over-representation of in-plume measurements, which biases NGI fluxes. Another possible source of bias may be wind measurements; if the winds at the location of the anemometers were systematically different from those at the UAV sampling locations, this would have a direct impact on the calculated fluxes. The flux biases in this work can be compared to a previous study by Brantley et al. [35], which quantified the controlled release of methane using Gaussian plume modelling of stationary measurements; their method both overestimated and underestimated fluxes by between −87% and +184%, respectively. The wide variation in fluxes acquired by Brantley et al. [35] highlights issues  The flux biases in this work can be compared to a previous study by Brantley et al. [35], which quantified the controlled release of methane using Gaussian plume modelling of stationary measurements; their method both overestimated and underestimated fluxes by between −87% and +184%, respectively. The wide variation in fluxes acquired by Brantley et al. [35] highlights issues similar to those experienced in other studies using near-field concentration measurements in turbulent (non-ideal) conditions with time-limited sampling, i.e. that sampling a non-disperse emission plume can result in large flux uncertainty.
Though seven out of eight F e values were above the CRF flux due to the biased sampling, all reported CRF fluxes fell within the full UAV flux uncertainty range, which spanned between (−76 ± 5)% and (+132 ± 15)% of F e . Flight 6 had a particularly high NGI flux uncertainty range due to a high quantity of sampling at the edge of the plume (see Figure A3), rather than the centre. Failure to adequately sample the centre of the plume can result in a much larger false plume being extrapolated, using the NGI method, in certain cases. It is clear from Figure 4 that the ∆F m uncertainty range was dwarfed by the extended uncertainty range, which was higher in all cases due to the dominating effect of the uncertainty in wind direction. To assess the various sources of UAV sampling uncertainty in ∆F m , a sensitivity test was performed (see Table S8). The principal source of uncertainty in ∆F m was from uncertainty in wind speed (see Table S8). An improved understanding of the temporal variability of the two-dimensional wind field on the flux plane, either by on-board UAV wind measurement or by other means, would serve to reduce this dominant source of uncertainty. However, the combined effects of instrumental uncertainty and uncertainty in the background were small compared with effects associated with natural variability in the wind field, which suggests that instruments such as the UGGA are suitable for measurement of fluxes similar to those in our study at similar downwind distances.
We compared the UAV flux uncertainties from the sample data to other methods. Seven of the eight TDM flux estimates which were acquired simultaneously during this study (see Table S1) agreed with the CRF, within uncertainty, with an average bias of (−2 ± 10)%. The average uncertainty in individual TDM fluxes (through comprehensive uncertainty budgeting) was ±21% [75]. In a different study using mass balance box modelling, Krautwurst et al. [42] derived the emission flux from a single landfill site based on both remote sensing column-concentration aircraft measurements and in-situ aircraft concentration measurements, with average flux uncertainties of ±34.5% and ±19.7%, respectively. Elsewhere, the uncertainty in fluxes obtained by Nathan et al. [73], using mass balance box modelling of polar coordinates and kriging to quantify the emission flux from a natural gas compressor station, was ±55%, although we dispute the use of kriging as an appropriate tool in that case. Yacovitch et al. [37] used a Gaussian plume inversion of ground based measurements to quantify emission fluxes from oil and gas infrastructure facilities, which used σ y (x) and σ z (x) values based on stability classes (rather than derived from q me ), with asymmetric 95% confidence uncertainties of between −33.4% and +334%. Therefore, we conclude that the NGI method using UAV sampling can be a competitive tool for rapid flux assessments of emitter facilities, assuming sampling from point source emitters. While the uncertainties we derived are comparatively large, the NGI method is able to derive rapid flux estimates from time-limited UAV sampling.
Though we report large NGI flux uncertainties for the sample data, alternative techniques such as mass balance box modelling cannot be used reliably to quantify fluxes using off-site downwind UAV sampling. To explain these conclusions, we have assessed the suitability of two mass balance box modelling approaches here. Mass balancing using geospatial kriging cannot be used, as kriging requires the interpolation of time-invariant downwind sampling. In this work, the q me plume sampled by a UAV near-to-source was not time invariant due to turbulent fluctuations in the position of the plume on the sampling plane. Though kriging may be used in principle in stable (non-turbulent) wind conditions, stable winds are unlikely when sampling with a UAV near-to-source. Alternatively, mass balancing may be applied by averaging q me into equally sized grid squares on a plane perpendicular to mean wind direction. However, this requires the entire extent of the time-averaged plume to be captured by the UAV, down to ground level. It also requires dense spatial sampling without any large sampling voids. While the method of grid square averaging may again be possible in principle, in order to produce sufficiently dense spatial sampling, a longer UAV sampling duration would be required compared to the sampling presented here. Furthermore, sampling down to ground level can be problematic, as it relies on the position of the air inlet on a UAV being close to the base of the UAV, and effects or prop-wash may be expected. Thus, we suggest that neither mass balance box modelling approach can be used for the type of near-field UAV sampling presented here-rather, the NGI method is a suitable flux quantification technique for near-field UAV sampling.
Future UAV flight sampling should be planned to minimise flux uncertainties, including potential sources of bias. The random walk simulations were used to generate values of A f for each flight from the sample data, which are given in Table 1. They show that regardless of measurement uncertainty, wind variations and the spatial extent of sampling, an average (58 ± 30)% limiting threshold uncertainty remains due to limited spatial coverage on the sampling plane. The random walk simulation could also be applied to obtain t 0.01 required to surpass a 1% uncertainty threshold, assuming good randomised sampling of a static plume. It is clear that the theoretical t 0.01 values in Table 1 (of at least 5 h) are unrealistic for UAV sampling, which cannot practically be achieved by a single light UAV platform over the timescale of a single flight (typically < 30 min using current battery technology constraints). However, using multiple UAV flights or tethered power (for uninterrupted sampling) may increase sampling time. The t 0.01 values are only valid when sampling in stable atmospheric conditions with a relatively non-variable wind direction. In reality, over the t 0.01 times given in Table 1, non-variable atmospheric conditions are unlikely, leading to the potential for large NGI uncertainties from UAV sampling. Therefore, if additional UAV sampling is used, a compromise sampling time must be reached over which atmospheric conditions change very little while the sampling period is maximised to minimise uncertainty.

Conclusions
A near-field Gaussian plume inversion was used to quantify methane emission fluxes from eight UAV flights downwind of a controlled emission source. The new NGI method takes into account the imperfect orientation of a two-dimensional vertical sampling plane perpendicular to mean wind direction. The number of unknown parameters in the inversion was reduced by deriving constraining parameters in the y direction from flux density measurements. Uncertainties were derived from measurement uncertainty as well as the difference between measured and modelled flux density results. The upper flux uncertainty bound was further enhanced by taking into account the lack of spatial sampling extent, using a random walk simulation of sampling. The random walk sensitivity tests also highlight areas of future improvements to the sampling strategy, which will be a focus of future experimental work. They also highlight the importance of evenly-distributed sampling of an adequate duration and the need for a sufficiently large sampling plane.
We suggest that the NGI method is a highly suitable technique that can be used to accurately quantify emission fluxes from static point sources of both methane and other trace gas sources, with measurement potential by UAVs. This allows for the full potential of near-field in-situ spatial sampling that UAVs provide to be realised. Though we acknowledge the magnitude of uncertainty may be high over a single UAV flight, the NGI method allows for rapid flux estimation with a defined maximal uncertainty. Mass balance box modelling may not be suitable for UAV in-situ sampling as it requires either the sampling of a time invariant plume or dense spatial sampling of a time-averaged plume. Though the NGI method yielded highly successful results when using a random walk simulation, the preliminary fluxes from the sample data resulted in high uncertainties. It is suggested that future studies that utilise the NGI method should adopt a sampling strategy that addresses the potential sources of bias highlighted here by ensuring that plume-chasing is avoided and by ensuring that sampling is sufficient to capture the entire extent of the time-averaged plume.
With the possibility of in-situ UAV wind measurements in future, there is further scope to improve flux accuracy. Future progress may be made by mounting an in-situ methane sensor directly onto a UAV, to eliminate issues associated with the use of a tether. Such improvements and further near-field sampling studies using the NGI method can play an important role in efforts to improve emission inventories. Additionally, UAVs may be used to add a vertical sampling dimension to the TDM technique using in-situ tracer gas detection. To conclude, a Gaussian plume inversion flux quantification method, using in-situ [CH 4 ] measurements of a sub-ppm precision, has been characterised here and future improvements are viable. This permits future studies of unknown well-defined sources (such as landfill sites, oil and gas infrastructure facilities, or herds of cattle) to be conducted using the method, depending on the precision of instrumentation.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2073-4433/10/7/396/s1, Table S1: Reported CRF methane flux releases with the time duration of continuous methane emission given in seconds past midnight (spm), local time (GMT); Figure S1: A photograph of the UAV platform at rest (without battery), with the PFA tubing pointing upwards; Table S2: UAV flight details and flight track colours; Figure S2: Aerial view of the three emission source locations (black dots) and flight tracks (coloured dots); Figure S3a: A histogram of methane concentration measurements obtained throughout flight 1 from the sample data; Figure S3b: A histogram of the lowest concentration measurements from flight 1 (red bars) with 15-bin-averages (green line) and log-normal function fit (blue line) plotted against the left hand axis; Table S3: Calculated background concentrations for each flight from the sample data; Figure S4: σ y (x) plotted for various stability classes (red lines) using the equations in Turner (1994); Figure S5: σ z (x) plotted for various stability classes (red lines) using the equations in Turner (1994); Table S4: y c , σ y (x) and σ z (x) for each UAV flight from the sample data and the average parallel distance of the sampling plane from the emission source; Figure S6: An example of a virtually-generated random walk to replicate flight 8 from the sample data, with the duration of the random walk set to the equivalent time period (21 min); Figure S7: F e,v, i values, simulating flight 3 from the sample data, as a function of t (blue crosses), plotted against left-hand axis, with F t (black line) also shown; Table S5: Flux convergence parameters required to approximate the fundamental threshold uncertainty associated with flux measurements as a function of sampling time, which are unique for each flight in the sample data; Table S6: Average virtual test F e values for each flight in the sample data with eight hours of random sampling of a static plume; Table S7: Initial flux values and corresponding fractional lower and upper uncertainty bounds, for each flight from the sample data; Table S8: A comparison between ∆F − me and ∆F − mo for each flight from the sample data, as a fraction of ∆F − ; Table S9: A comparison between ∆F − me,U(z) , ∆F − me,ρ and ∆F − me,ib for each flight from the sample data, as a fraction of ∆F − me . single-storey office facilities. The site has a maximal length of~600 m and a maximal width of~400 m. During each flight, fluxes were released at one of two prescribed emission rates from three locations around the site from a~6.2 m tower, selected for optimal downwind UAV sampling for the prevailing wind direction at the time of each flight experiment (see Table S1 for release details). The MRU was ideal for the acquisition of sample data, as it was distant (>5 km) from any expected significant sources of off-site (and therefore interfering) methane such as large towns, cities or landfill sites. Four 10 Hz sonic 3D anemometers (HS-50 Horizontal-Head Ultrasonic Anemometer, Gill Instruments Limited) were positioned on the site at heights of 2, 10, 25, and 50 m above ground level. The average wind speed and wind direction for the duration of each UAV flight was interpolated (using a shape-preserving piecewise cubic interpolation) as a function of z to derive a continuous wind profiles as a function of height for each flight. Thus, θ took into account wind shear by interpolating wind speed measurements from the four anemometer towers. The variability in wind speed measured by each anemometer was interpolated as a function of altitude to produce a continuous uncertainty profile.

Author
During each controlled release, a modified DJI Spreading Wings S900 hexacopter UAV ( Figure A1) was manually piloted by experienced personnel on the ground and guided for sampling by an assistant, acting as a lookout. The pilot had a qualification in UK airspace legislation and was experienced and competent in flying small UAVs. As the UAV was used for scientific (rather than commercial) purposes, a pilot qualification from the Civil Aviation Authority (CAA) was not required, in accordance with UK airspace legislation. The UAV flew at an average speed of (0.4 ± 0.1) m s −1 . It was directed mid-flight to capture as much of the plume as possible based on mid-flight knowledge of the plume position, using the release of smoke flares. Though the UAV was limited to a physical height of 150 m (i.e., the length of the tubing), the UAV only flew up to (24 ± 3) m to ensure extensive vertical sampling of the emission plume; the maximum height of the plume was determined using a visible plume observed from the release of smoke flares. Prior to each flight, the CRF release point and UAV sampling area were defined in a planning and safety briefing, using forecasts of wind direction to optimise the available operational sampling area around on-site buildings and obstructions (such as the anemometer masts) and to ensure minimum UAV safety distances (50 m) from personnel and property beyond the pilot's control, as defined by CAA regulations [77]. Each flight track was roughly aligned in a vertical sampling plane perpendicular to θ by manual piloting, (90 ± 30) m from the source. However, there was a clear inclination between the sampling plane and the plane perpendicular to θ, which was accounted for in the NGI method. UAV longitude, latitude, and altitude were recorded at a 5 Hz frequency by an on-board Pixhawk computer (independent of the on-board flight controller), connected to a Cirocomm 00BY satellite geolocation antenna, which was secured to the top of the UAV. For a continuous stationary 60 s sampling period, x, y, and x varied with a standard deviation of no more than ±0.15 m. Thus, the precision of the sensor was suitable for our flux quantification purposes. single-storey office facilities. The site has a maximal length of ~600 m and a maximal width of ~400 m. During each flight, fluxes were released at one of two prescribed emission rates from three locations around the site from a ~6.2 m tower, selected for optimal downwind UAV sampling for the prevailing wind direction at the time of each flight experiment (see Table S1 for release details). The MRU was ideal for the acquisition of sample data, as it was distant (>5 km) from any expected significant sources of off-site (and therefore interfering) methane such as large towns, cities or landfill sites. Four 10 Hz sonic 3D anemometers (HS-50 Horizontal-Head Ultrasonic Anemometer, Gill Instruments Limited) were positioned on the site at heights of 2, 10, 25, and 50 m above ground level. The average wind speed and wind direction for the duration of each UAV flight was interpolated (using a shape-preserving piecewise cubic interpolation) as a function of z to derive a continuous wind profiles as a function of height for each flight. Thus, θ took into account wind shear by interpolating wind speed measurements from the four anemometer towers. The variability in wind speed measured by each anemometer was interpolated as a function of altitude to produce a continuous uncertainty profile.
During each controlled release, a modified DJI Spreading Wings S900 hexacopter UAV ( Figure  A1) was manually piloted by experienced personnel on the ground and guided for sampling by an assistant, acting as a lookout. The pilot had a qualification in UK airspace legislation and was experienced and competent in flying small UAVs. As the UAV was used for scientific (rather than commercial) purposes, a pilot qualification from the Civil Aviation Authority (CAA) was not required, in accordance with UK airspace legislation. The UAV flew at an average speed of (0.4 ± 0.1) m s −1 . It was directed mid-flight to capture as much of the plume as possible based on mid-flight knowledge of the plume position, using the release of smoke flares. Though the UAV was limited to a physical height of 150 m (i.e., the length of the tubing), the UAV only flew up to (24 ± 3) m to ensure extensive vertical sampling of the emission plume; the maximum height of the plume was determined using a visible plume observed from the release of smoke flares. Prior to each flight, the CRF release point and UAV sampling area were defined in a planning and safety briefing, using forecasts of wind direction to optimise the available operational sampling area around on-site buildings and obstructions (such as the anemometer masts) and to ensure minimum UAV safety distances (50 m) from personnel and property beyond the pilot's control, as defined by CAA regulations [77]. Each flight track was roughly aligned in a vertical sampling plane perpendicular to θ by manual piloting, (90 ± 30) m from the source. However, there was a clear inclination between the sampling plane and the plane perpendicular to θ, which was accounted for in the NGI method. UAV longitude, latitude, and altitude were recorded at a 5 Hz frequency by an on-board Pixhawk computer (independent of the on-board flight controller), connected to a Cirocomm 00BY satellite geolocation antenna, which was secured to the top of the UAV. For a continuous stationary 60 s sampling period, x, y, and x varied with a standard deviation of no more than ±0.15 m. Thus, the precision of the sensor was suitable for our flux quantification purposes. [CH 4 ] was recorded at~1 Hz by securing a 150 m long, 4.76 mm inner diameter perfluoroalkoxy alkane (PFA) tube to the UAV, through which air was pumped (at a continuous flow rate of 0.5 litres per minute) to a commercially available 15 kg UGGA methane sensor ( Figure A3) which was fixed on the ground, operating at a fixed pressure of 140 Torr. The UGGA uses off-axis integrated cavity output spectroscopy to derive concentration measurements (see Baer et al. [78] and Paul et al. [79] for full spectroscopic details): A 1650 nm tuneable laser is reflected between two mirrors in a high finesse cavity, with the output spectrum integrated as a function of wavelength, to derive concentration using the Beer-Lambert law. Live concentration measurements from the UGGA could be viewed on a computer screen on the ground, where data was saved for subsequent analysis. The UGGA was powered in the field using a 12 V lead-acid battery. Though the ground-based sampling used here allows for high-precision heavier instrumentation (like the UGGA) to be used, it also suffers a number of drawbacks, such as lag time through the tubing (which we account for) and the kinking of the tube, resulting in air failing to reach the sensor on the ground (which we omitted during flux analysis). A small loop at a fixed position on the tubing was fastened to the same point beneath the UAV, prior to flight, to ensure consistency in the length of the tubing between flights. In this work, the tip of the inlet was positioned~20 cm above the plane of the rotor blades to ensure that sampled air was not disturbed by propwash or air flow through the rotors (see Zhou et al. [66] for example). Computational fluid dynamics (CFD) model simulations (not shown) demonstrated that the inlet was positioned sufficiently far enough away from the rotors to exclude this. The CFD simulations showed that there was no significant vertical disturbance to air a sufficient distance above the plane of the propellers. [CH4] was recorded at ~1 Hz by securing a 150 m long, 4.76 mm inner diameter perfluoroalkoxy alkane (PFA) tube to the UAV, through which air was pumped (at a continuous flow rate of 0.5 litres per minute) to a commercially available 15 kg UGGA methane sensor ( Figure A3) which was fixed on the ground, operating at a fixed pressure of 140 Torr. The UGGA uses off-axis integrated cavity output spectroscopy to derive concentration measurements (see Baer et al. [78] and Paul et al. [79] for full spectroscopic details): A 1650 nm tuneable laser is reflected between two mirrors in a high finesse cavity, with the output spectrum integrated as a function of wavelength, to derive concentration using the Beer-Lambert law. Live concentration measurements from the UGGA could be viewed on a computer screen on the ground, where data was saved for subsequent analysis. The UGGA was powered in the field using a 12 V lead-acid battery. Though the ground-based sampling used here allows for high-precision heavier instrumentation (like the UGGA) to be used, it also suffers a number of drawbacks, such as lag time through the tubing (which we account for) and the kinking of the tube, resulting in air failing to reach the sensor on the ground (which we omitted during flux analysis). A small loop at a fixed position on the tubing was fastened to the same point beneath the UAV, prior to flight, to ensure consistency in the length of the tubing between flights. In this work, the tip of the inlet was positioned ~20 cm above the plane of the rotor blades to ensure that sampled air was not disturbed by propwash or air flow through the rotors (see Zhou et al. [66] for example). Computational fluid dynamics (CFD) model simulations (not shown) demonstrated that the inlet was positioned sufficiently far enough away from the rotors to exclude this. The CFD simulations showed that there was no significant vertical disturbance to air a sufficient distance above the plane of the propellers. The accuracy of UGGA concentration measurements was enhanced using a post-flight water vapour correction (see O'Shea et al. [80] for details). Measurements were then further calibrated using a two term linear correction (derived from an average of eight calibration procedures, described in detail by Pitt et al. [81]), such that dry methane concentration measurements were robustly traceable to the World Meteorological Organization greenhouse gas calibration scale [82]. Equation (A1) was used to calculate σi([CH4]) values specific to the UGGA (required to calculate of uncertainty in qme in Equation (3)). The accuracy of UGGA concentration measurements was enhanced using a post-flight water vapour correction (see O'Shea et al. [80] for details). Measurements were then further calibrated using a two term linear correction (derived from an average of eight calibration procedures, described in detail by Pitt et al. [81]), such that dry methane concentration measurements were robustly traceable to the World Meteorological Organization greenhouse gas calibration scale [82]. Equation (A1) was used to calculate σ i ([CH 4 ]) values specific to the UGGA (required to calculate of uncertainty in q me in Equation (3)). (A1) GF was the average of eight previously calculated gain factors, and σ GF was its corresponding standard deviation. The ratio between σ GF and GF was 0.00282. WC([CH 4 ]) is a water correction factor which was directly applied to raw concentration measurements and is a function of water concentration. σ WC was the uncertainty associated with the water correction.
Each longitude, latitude, and altitude measurement from the UAV was interpolated onto the lower frequency UGGA time grid (avoiding extrapolation of [CH 4 ] measurements). UGGA measurements were also corrected to the UAV sampling time by calculating the lag time (τ lag ), in order to geospatially reference UGGA measurements (by subtracting τ lag from the UGGA time-stamp), to account for the time taken for air to travel through the tubing. τ lag was tested prior to and immediately following each flight (by blowing into the inlet and timing the carbon dioxide concentration response). τ lag varied flight-to-flight with a mean (and maximal range from flight-to-flight) of (300 ± 10) s. τ lag was measured before and after each flight and found to be constant. However, it gradually varied throughout the day. This small variation may have been due to changes in ambient pressure, which slightly alter the mass flow rate through the pump in the UGGA. τ lag was assumed to remain constant throughout each flight, as was it controlled by the strength of the pump (operating at a controlled pressure), with ambient pressure changes assumed to be minimal over the short period of each UAV flight.

Appendix B
Geospatially mapped [CH 4 ] me measurements for each UAV flight have been plotted in Figure A3, where z has been plotted on the y-axis, and y (obtained by rotating the distance from the Prime Meridian and the Equator, respectively, by θ) has been plotted on the x-axis. For each of the eight flight tracks shown in Figure A3, the CRF emission point was from a tower~6.2 m above ground level. For each flight bar flight 4, the highest concentration measurements were typically detected by the UAV to the left of the position of the emission source projected by θ onto the sampling plane (set to 0 m in the y direction in Figure A3). The CRF emission flux for flights 1-6 was (3.0 ± 0.2) g s −1 , while the emission flux for flights 7 and 8 was (1.5 ± 0.1) g s −1 . Methane was released at these flux rates throughout the period of sampling.
Atmosphere 2019, 10, x FOR PEER REVIEW 20 of 25 Figure A3. In-situ [CH4]me measurements obtained by the UAV for the sample data, plotted according to geolocation corrected for time lag. Flights are numbered according to assigned flight numbers given in Table S2 in the supplement. A logarithmic color scale has been used. z is plotted on the vertical axis and y is plotted on the horizontal axis.
The highest recorded concentrations during flight 1 were observed between ~−10 m and ~−35 m horizontally, where the UAV clearly intersected the centre of the advecting emission plume. The highest concentration measurements from flight 2 were made nearer to ground level from between ~−25 m and ~+15 m horizontally. These concentration measurements were the highest experienced during all eight flights (~25 ppm). High concentration measurements were also detected during transects ~10 m above ground level, with some to the very right (~+30 m) of the extent of UAV Figure A3. In-situ [CH 4 ] me measurements obtained by the UAV for the sample data, plotted according to geolocation corrected for time lag. Flights are numbered according to assigned flight numbers given in Table S2 in the supplement. A logarithmic color scale has been used. z is plotted on the vertical axis and y is plotted on the horizontal axis.

5.
The penultimate F e value must be less than 0.98 times the maximum constraining F e value (which ensures that F e has stopped increasing).
The first three conditions ensured that τ z has stabilised. The fourth and fifth conditions ensured that the maximum constraining F e value was sufficiently large to avoid it constraining F e . stopped increasing). 2. The final τz value must be less than 0.98 times τz, max for that modelled run (which ensures that τz, max has sufficiently exceeded τz). 3. The penultimate τz value must be less than 0.98 times τz, max for that modelled run (which ensures that τz, max has sufficiently exceeded τz). 4. The final Fe value must be less than 0.98 times the maximum constraining Fe value (which ensures that Fe has stopped increasing). 5. The penultimate Fe value must be less than 0.98 times the maximum constraining Fe value (which ensures that Fe has stopped increasing).
The first three conditions ensured that τz has stabilised. The fourth and fifth conditions ensured that the maximum constraining Fe value was sufficiently large to avoid it constraining Fe. Figure A4. The variation in the modelled output value of Fe (green dots plotted against the left-hand axis) and τz (red dots plotted against the right-hand axis) as a function of τz, max, for flight 7 from the sample data. The maximum value of τz is plotted as a green line against the right-hand axis. Figure A4. The variation in the modelled output value of F e (green dots plotted against the left-hand axis) and τ z (red dots plotted against the right-hand axis) as a function of τ z,max , for flight 7 from the sample data. The maximum value of τ z is plotted as a green line against the right-hand axis.