Quantifying methane emissions from coal mining ventilation shafts using an unmanned aerial vehicle (UAV)-based active AirCore system

CH 4 mole fractions downwind of a single coal mining ventilation shaft. Besides CH 4 , we also measured CO 2 , CO, atmospheric temperature, pressure, and relative humidity. Wind-speed and wind-direction measurements were made using a lightweight balloon-borne radiosonde. Fifteen UAV flights were performed flying perpendicular to the wind direction at several altitude levels, to effectively build a ‘curtain ’ of CH 4 mole fractions in a two-dimensional plane at a distance between 150 and 350 m downwind of a single ventilation shaft. Furthermore, we have developed an inverse Gaussian approach for quantifying CH 4 emissions from a point source using the UAV-based observations, and have applied it as well as the mass balance approach to both simulated data and actual flight data to quantify CH 4 emissions. The simulated data experiments revealed the importance of having multiple transects at different altitudes, appropriate vertical spacing between the in- dividual transects, and proper distance between the center height of the plume and the center flight transect. They also showed that the inverse Gaussian approach performed better than the mass balance approach. Our estimate of the CH 4 emission rates from the sampled shaft ranges from 0.5 to 14.5 kt/year using a mass balance approach, and between 1.1 and 9.0 kt/year using an inverse Gaussian method. The average difference between the mass balance and the inverse Gaussian approach was 2.3 kt/year. Based on the observed correlation between CO 2 and CH 4 (R-squared > 0.69), the CO 2 emissions from the shaft were estimated to be between 0.3 and 9.8 kt/ year. This study demonstrates that the UAV-based active AirCore system provides an effective way of quantifying coal mining shaft emissions of CH 4 and CO 2.


Introduction
Atmospheric methane (CH 4 ) is the second most abundant anthropogenic greenhouse gas (GHG), and its global warming potential is 28-32 times greater than that of carbon dioxide (CO 2 ) for a 100-year horizon (Etminan et al., 2016). The mole fraction of CH 4 has been increasing over the last 150 years, with more than a doubling of its concentration since the 1850s (IPCC, 2014). CH 4 is produced naturally when organic matter decomposes in an environment containing low oxygen. This process is known as natural CH 4 production and comes from sources like wetlands, fresh waters, and oceans. Currently, most CH 4 is emitted through human-made (anthropogenic) operations, where sources include agriculture and landfills, coal mining, oil and gas industry and biomass burning (Kirschke et al., 2013;Saunois et al., 2020). A significant part of the anthropogenic CH 4 emissions come from fossil exploitation (Kirschke et al., 2013). However, the fossil fuel related emissions of CH 4 are often not well characterized, carrying significant uncertainties ranging from 20 to 30% for both bottom-up and top-down estimates, which is especially true for the substantial contribution from coal mining where top-down estimates are lacking (Kirschke et al., 2013;Saunois et al., 2020). Both the production (extracting from the coalbed) and processing of the coal (crushing and pulverization) releases large amounts of CH 4 , where underground mines are responsible for the majority of coal-related emissions. The CH 4 from underground mines are diluted in air and emitted to the atmosphere via ventilation shafts located at the surface (Zazzeri et al., 2016). The CH 4 embedded in the coalbeds comes from the carbonification of vegetable substances (Swolkień, 2020). The amount of CH 4 emitted to the atmosphere is poorly determined as it is removed through the ventilation shaft due to non-precise measurements of the flowrate and CH 4 mole fraction, leaving large uncertainties in the total CH 4 emission. The Silesia area in Poland houses a large number of coal mining facilities. It was, according to the European Pollutant Release and Transfer Register (E-PRTR), responsible for emitting 447 kiloton CH 4 in 2017 (E-RPTR, 2017). The E-PRTR inventories are based on reported CH 4 activity determined by occasionally taking airbag samples and analyzing them using a flame ionization detector in a laboratory. However, verification of emission estimates from individual shafts is missing.
Unmanned aerial vehicles (UAVs) have recently become an emerging platform for atmospheric GHG measurements. They are versatile in their use and comparatively cheap to purchase and maintain (Villa et al., 2016). There are currently three principal approaches used to obtain GHG measurements using UAVs: continuous measurements using sensors aboard (Berman et al., 2012;Khan et al., 2012;Nathan et al., 2015;Golston et al., 2017;Kunz et al., 2018;Allen et al., 2019;Tuzson et al., 2020), continuous measurements with an analyzer on the ground analyzing air sampled through a tethered tube attached to a flying UAV (Turnbull et al., 2014;Brosy et al., 2017;Wolf et al., 2017;Shah et al., 2019Shah et al., , 2020, and finally, discrete or semi-continuous measurements by on-board sampling of air for later analysis (Chang et al., 2016;Greatwood et al., 2017;Andersen et al., 2018).
There have been several studies using UAV measurements to estimate CH 4 emissions. Nathan et al. (2015) flew a fixed-wing model aircraft with an onboard laser-based open-path CH 4 sensor, where they quantified CH 4 emitted from a compressor station using a mass balance approach and an inverse Gaussian approach. They obtained fluxes with uncertainties of about 55%. In Shah et al. (2020), a UAV tethered to a Los Gatos analyzer sampled air downwind of a controlled release of CH 4 . They quantified the CH 4 flux with uncertainties ranging from 71 to 147% using a near-field Gaussian method. Allen et al. (2019) inferred CH 4 fluxes from a landfill using proxy measurements of CO 2 concentration measured via a UAV, in combination with a mass balance approach. They obtained fluxes with uncertainties ranging between 54 and 61% at 1σ.
With this study, we aim to quantify the emitted CH 4 from a single coal mine ventilation shaft located in Silesia, Poland, using an active AirCore system. According to Polish mining regulation, a CH 4 drainage system must be placed if the anticipated total methane-bearing capacity is above 10 m 3 /min (~3.8 kton CH 4 /yr). The captured CH 4 from the drainage system is transported via dedicated pipelines to a drainage station, either combusted for the purpose of heating/cooling onsite or is sold to nearby power plants. However, not all captured CH 4 is combusted or sold due to capacity issues, and the fraction of CH 4 utilized by coal mines varies seasonally, and there may be occasions when unused CH 4 is released from the drainage station to the atmosphere. It is important to notice that the occasional release of CH 4 takes place at the drainage station, and is not vented through ventilation shafts. For the particular coal mine we studied, there are three ventilation shafts, and the CH 4 drainage station is in the mine facility and is ~3.5 km from the studied ventilation shaft. To this end, CH 4 emissions from the drainage station should be regarded as a separate and occasional point source, but has limited influence on our sampling due to the long distance. On the other hand, CH 4 emissions from the ventilation shaft are continuous, and are one of three dominant sources from the coal mine. According to the data from State Mining Authority (WUG), the average percentage of CH 4 emitted to the atmosphere via ventilation shafts vs. captured but vented to the atmosphere from drainage stations in 2020 is 82% vs. 18%, and this percentage may differ for different coal mines. Also, reported CH 4 emissions from coal mines are registered as a sum of CH 4 emitted from the ventilation shafts and released from the drainage system. Therefore, quantification of the ventilation shaft CH 4 emissions allows us to obtain an independent estimate of a major part of CH 4 emissions from the coal mine and to verify the reported CH 4 emissions. The active AirCore system is an air sampling system consisting of a long coiled stainless-steel tubing, a small pinhole orifice, and a small pump that drags air through the tube (Andersen et al., 2018). The high mobility and versatility of the UAV, together with the active AirCore's ability to capture and preserve an atmospheric air sample along the trajectory of the UAV, creates an opportunity to make both horizontally and vertically semi-continuous atmospheric measurements of CH 4 , CO 2 , and CO downwind of a point source. This downwind trajectory creates a two-dimensional map of accurate CH 4 mole-fractions of the point source, which is subsequently used to quantify the CH 4 emission rate. The main objective of the work is to develop a methodology for sampling and quantifying CH 4 point source emissions using the active AirCore system.
The field work was performed in the summer of 2017, as part of a joint pre-campaign of CoMet (Carbon dioxide and Methane mission) in the Silesia coal mining region. The overall aim of CoMet was to improve the methodology for measuring the efficiency of gas streams (CO 2 and CH 4 ) at a local and regional scale (Fix et al., 2018;Swolkień, 2020). A variety of different approaches have been implemented during the CoMet campaign: quantification of CH 4 emissions from single shafts and a subregion of the Upper Silesia Coal Basin (USCB) using a ground-based, portable, sun-viewing Fourier transform spectrometer (Luther et al., 2019); estimation of the emissions of CH 4 , CO 2 , and CO of the entire USCB using a mass balance approach and dispersion modeling based on in situ aircraft measurements, respectively (Fiehn et al., 2020;Kostinek et al., 2021); determination of CH 4 emissions from small clusters of ventilation shafts with airborne passive remote sensing measurements of CH 4 column anomalies (Krautwurst et al., 2020). The advantage of using a UAV-based active AirCore is that individual plumes from ventilation shafts can be accurately mapped with in situ concentration measurements of CH 4 , CO 2 , and CO on a low-cost and flexible platform, which can be used to quantify ventilation shaft emissions and validate model dispersion modeling.
The instrumentation, flight planning, simulated data experiment, and quantification method are presented in the methodology (section 2). The results (section 3) show CH 4 emission quantifications from the simulated data experiments and flights made during the campaign, using the two methodologies: correlation data between CO 2 and CH 4 ; CO 2 emission rate quantification; and an estimate of the ventilation shaft's hourly CH 4 emission data. Section 3 also discusses the results, the effect of atmospheric variability, and the uncertainty of the methods. Finally, section 4 presents the conclusions.

Active AirCore
The active AirCore system is capable of continuously collecting air samples to obtain atmospheric observations of CO 2 , CH 4 , and CO mole fractions along the flight track of a UAV flight. The details of the system were described in Andersen et al. (2018), and we have made further improvements to the system since then. Previously, the on-board relative humidity sensor, and all three of the temperature sensors, were located inside the carbon fiber flight-box along with the active AirCore. Prior to the CoMet 0.5 campaign, the relative humidity sensor, and one of the three temperature sensors, were moved outside the flight box, and fastened underneath. This added ambient atmospheric measurements of relative humidity and temperature during the flights. In addition, the dryer tube was extended from 7.5 cm to 17.5 cm. This reduced the mole fraction of H 2 O in the active AirCore samples from ~0.30% to below 0.06% when the mole fraction of ambient H 2 O is ~1%. The inlet, previously located at the bottom of the flight-box, was repositioned to the side. This makes it easier to attach the dryer to the active AirCore, and to remove the cap of the dryer prior to the flights, and put on the cap of the dryer after the flights had been completed. Three active AirCores were used for this study. Table 1 below shows the characteristics of these active AirCores.

Analysis system
The collected air samples were analyzed by a cavity ring-down spectrometer (CRDS, Picarro Inc. CA, USA, model G2401-m) through an inlet system that is similar to the one described in Andersen et al. (2018), with some minor changes to the analysis system, that among others make it easier to adjust the excess flow rate of the fill and reference gas independently. The analysis system is now also capable of additional sub-sampling of AirCore samples. Additionally, the active AirCore is connected to the analysis box using 1/8 in. ultra-Torr fittings (SS-2-UT-6) instead of the previously used 1/8 in. Swagelok fittings with stainless-steel nuts and Teflon ferrules. These additions have two benefits: one being increased ease-of-use to connect the active AirCore to the analysis box (no tools required to make a leak tight connection); the other being less metal-to-metal friction, which reduces the CH 4 generated during the connection of the active AirCore to the analysis box (Higaki et al., 2006).

2.1.2.1.
The sub-sampling of AirCore samples. A section has been added downstream of the pump used for the CRDS. This part allows for the filling of collected atmospheric air samples into Supel multi-layer foil gas sampling bags (Supelco, product no. 30279-U). These bags have a volume of 0.6 L. This added element provides an option of collecting an integrated atmospheric air sample, of roughly the volume of the active AirCore, for laboratory analyses of isotopic compositions of 13 C and 2 H in CH 4 . The isotopic composition measurements will not be shown in this paper, but in a separate one. A three-way valve is located downstream of the exhaust of the vacuum pump. One side leads to a multilayer foil gas sampling bag, while the other side leads to an Alicat MB-100SCCM-D/5M flowmeter. A needle valve is located between the CRDS and the vacuum pump to have an adjustable analysis flow rate. During the sample analysis, the analysis flow rate was kept constant at ~22-23 [ccm], and was not affected by the sub-sampling downstream of the pump.

2.1.2.2.
The cavity ring-down spectrometer. The CRDS was configured to a high-range CH 4 mode during the sample analysis to be able to measure largely enhanced CH 4 mole fractions. In the normal mode, the CH 4 spectral absorption lines were saturated around 20 ppm while the sampled mole fraction of CH 4 in this study was up to 125 ppm during some flights.
A two-point calibration was conducted using a known WMO-scale gas mixture around ambient CH 4 concentration (WMO X2007, X2004A, and X2014A scales for CO 2 , CH 4 , and CO, respectively), and a certified mole-fraction gas mixture from the Dutch National Metrology Institute (VSL) containing a high mole-fraction CH 4 (see Table 2).

Meteorological measurements
Radiosonde launches were performed a few minutes before and during the UAV active AirCore flights. The radiosonde (Sparv Embedded AB, Sweden, model S1H2-R) measures ambient temperature, pressure, relative humidity, and most importantly, wind speed, and wind direction.
The temperature sensor has an operating range of − 40 ~ +80 • C, with an accuracy of 0.3 • C. The ambient pressure sensor has a range of 300-1100 mbar, with an accuracy of 1 mbar. The relative humidity sensor has a range between 0 and 100%, and measures with an accuracy of approximately 2%. The wind speed and direction measurements are calculated based on the temporal change of the position of the radiosonde during a launch. The radiosonde weighs ~12 g, and was launched with a 9 g latex balloon filled with ~30 L helium gas. The balloon was tethered to a thin fishing line so that both the balloon and the radiosonde could be easily recycled. During the days of operation, the communication between the radiosonde and its positional satellites was good. It is thus assumed that the uncertainty in the wind direction is low. The wind speed can be estimated within a range of 0-150 m/s, with an uncertainty of approximately 5% or 0.2 m/s at a wind speed of less than 4 m/s.

The spatial resolution of AirCore measurements
The spatial resolution of the active AirCore flights were discussed in Andersen et al. (2018), and is comprised of three contributing parts; uncertainty in the GPS positioning, uncertainty due to diffusion and Taylor dispersion, and uncertainty due to mixing of air in the cavity of the analyzer. The dominant effect is the mixing of air in the analyzer cavity, and causes the analyzed sample to be smoothed, removing larger peaks and troughs in the signal. To understand the effect of the cavity mixing on the active AirCore profiles, we attempt to de-tangle the analyzed signal with the actual signal, using a moving average with a specific averaging kernel as an approximate solution.
The inverse Gaussian approach uses an averaging kernel which specifies the number of seconds used during a moving average. In the determination of the averaging kernel for the moving average, a deconvolution of the measured CRDS signal with the response function of the analyzer was done to obtain the original in-situ signal. Using a simulated flight, we could compare the deconvoluted signal with that of a moving average to find which moving average kernel best coincide with the retrieved signal. The Gaussian function given in eq. (3) was used to generate a virtual flight profile with elevated mole fractions of CH 4 at three altitude levels.
The dominating factor in the determination of the spatial resolution of the active AirCore flights comes from mixing of air in the cavity of the analyzer. To further study this effect, the response function of the analyzer is required. Using the Picarro, two cylinders were used to  switch from a low mole fraction of CH 4 to a high fraction of CH 4 , forming a plug flow. The flow rate was held at a constant 21-22 sccm, which is the same flowrate that was used to measure the active AirCore data during the study. This allowed us to fit an exponential function to the plug flow and obtain the response function of the analyzer. As explained in Winderlich et al. (2010), the measured concentration c(t) through a well-mixed volume can be calculated by convolving the in-situ concentration s(t) with the response function g(t) of the volume, such that: In this case, c(t) is the integrated concentration measured by the CRDS analyzer. Taking the derivative of equation (1) yields a solution as an approximation for the highly variable in-situ data: We found the best match between the convoluted signal and the moving average using an averaging kernel of 33-34 s. This was thus used for all flights during the processing of the data.

UAV flights
During this study, 15 active AirCore flights were performed by a DJI Inspire 1 UAV, using three different active AirCore systems. The drone is a rotary UAV which weighs 2.9 kg and has a maximum flight time of approximately 15 min. With the active AirCore system the weight was ~4.0 kg and was able to make ~12 min flights.

Flight track
All the flights were performed downwind of a single coal mine ventilation shaft, located in the USCB in Poland (Coordinates: 49.975459 • N, 18.735420 • E). To ensure that each flight was executed downwind of the shaft, the wind direction was monitored prior to the launch of each flight by launching a radiosonde. The flight patterns were then planned accordingly. We aimed to obtain flight tracks that were perpendicular to the wind direction, starting at background (i.e., not plume affected) mole fractions of CH 4 . The drone would take off and ascend to a predetermined altitude. It would then fly in a horizontal line, perpendicular to the wind direction while maintaining the same altitude, attempting to transect the plume. Once the drone had passed the assumed plume location, it would ascend to a higher altitude to perform another horizontal transect. This pattern continued until the drone battery dropped to 20% of its capacity, which prompted the drone to return to its home position. The flight time ranged from 9 to 11 min. The horizontal distance from the coal mining ventilation shafts and the transect plane of the drone flights ranged from 150 m-350 m. Fig. 1 illustrates these flights.

Flight selection criteria
In this study, we have selected the flights that meet the following two requirements for further analysis: 1. The wind speed is larger than 2 m/s. At lower wind speeds, the wind direction is ill-defined, meaning the plume is less likely to follow a stable wind direction. Besides this, the relative uncertainty of the wind speed, which is proportional to the uncertainty of the emission rate estimates, is larger when the wind speed is small. 2. The flight trajectory is perpendicular to the wind direction (within 15 • ). When the flight plane deviates significantly from the direction that is perpendicular to the wind direction, the chance of capturing the plume decreases, causing a sampling bias. In addition, the uncertainty of the quantified emission rate increases with the deviation of the angle between the wind direction and the flight plane from 90 • .
Out of the 15 flights, 8 fulfilled these requirements. The main cause for flights to have failed the requirements was not being perpendicular to the wind direction, with 6 out of 15 flights failing due to this. Three out of 15 failed the requirement of the wind speed being too low (see Table A.2). The successful flights are presented in Table 3.  Table 3 The wind speed and its uncertainty, as well the angle on the wind direction for each successful flight and date at which each flight was performed.

Emission rate quantification approaches
In this study, two methodologies were used to quantify the emission rates from the coal mining ventilation shafts, namely an inverse Gaussian approach and a mass balance approach. Commonly for both approaches, a linear fit is applied to the positional coordinates of the flight to ensure that the flight is not deviating strongly from a line (this is only done for flights where the drone had been flying seemingly along a line). This is done to simplify the estimate of the mean deviation of the flight plane from the direction that is perpendicular to the wind, and the calculation of the distance between flight points along the flight track. It also allows us to see whether the flight is deviating much from a straight line, which would complicate the assumption that the flight track is perpendicular to the wind direction (see Fig. 2). For example, the angle between the wind direction and the perpendicular line on the linear fit was 9 • , 8 • , and 9 • for flights #6, #12, and #18. The positional coordinates, namely the longitude and latitude, were then projected onto this fitted line. The average distances that the positional coordinates are moved when projected onto the line are 1.9 m, 2.0 m, and 7.8 m for flights #6, #12, and #18, respectively, with the standard deviations (1 σ) of the moved distances of 1.6 m, 1.6 m, and 8.1 m for the same flights.
Once the points are all projected onto the line, the distance between the take-off position and each individual point is calculated. The mean distance between the points are 2.7 m, 2.7 m, and 2.4 m, with standard deviations of 2.3 m, 1.4 m, and 1.4 m for flights #6, #12, and #18, respectively.

Inverse Gaussian approach
The inverse Gaussian approach uses the sampled three-dimensional mole fraction data from the active AirCore together with a Gaussian dispersion model for point sources from Lagzi et al. (2013). This model is given by equation (3): Where is the dry mole fraction at a given position x, y, and z, which are the projected positional coordinates (see previous section) downwind of the plume, across the plume horizontally, and across the plume vertically. The units of C ′ (x, y, z) in mol mol − 1 , and the units of x, y, and z are given in m. The emission rate Q is given in kg s − 1 , the wind speed u⋅ Cos(θ) in m s − 1 (with θ being the angle between the wind direction and the perpendicular line of the flight trajectory), and the stack height h is given in m. The parameters σ y and σ z describe the mixing of the pollutants in the horizontal-and vertical direction, respectively, and have units of m. V is the dry molar volume in m 3 mol − 1 . The molar volume is expressed by equation (4): where R is the universal gas constant of 8.314 J mol − 1 K − 1 , T is the ambient temperature given in K, P is the ambient pressure in Pa, and P H2O is the water vapor pressure in Pa. Lastly, M CH4 is the molar mass of CH 4 , given in. kg mol − 1 . The AirCore flight data (Y) is paired with the Gaussian dispersion model, where the emission rate (Q), the mixing of pollutants in the horizontal-and vertical direction (σ y and σ z ), and the center height of the plume (h) are parameters to be optimized. A best fit for eq. (3) to the data can be found for these four parameters by minimizing the cost function J(x) using a standard square error (SSE) approach (Rawlings et al., 1998), such that where x is a vector of the 4 unknowns (Q, σ y , σ z , and H), C(x) is the modelled mole fractions, and Y is the measured mole fraction. The minimization of J(x) is done using the nonlinear solver fmincon optimizer from MATLAB (The MathWorks, 2019). The solver puts constraints on the unknowns given in x, which are pre-selected prior to the optimization.
To ensure that not only a local minimum is found, the optimization is run several times while the starting points for the optimization parameters are randomized between the upper and lower boundary conditions of its respective parameter. These boundary conditions are determined by looking at the sampled active AirCore data. The three transects of the flight cut through the plume at different heights, typically with the middle transect having a larger CH 4 mole-fraction than the lower and upper transect. This allows us to constrain the center height of the plume h. As if drawing lines through a circle at different places, the different transects will have varying mole-fractions of CH 4 . This will allow us to not only constrain σ z , in the vertical direction, but also σ y in the horizontal direction by obtaining the best fit between the inverse Gaussian model and the different transects. Constraints can then be put on both σ y and σ z to ensure that the dispersion in both y-and z-direction is accounted for. This provides us with the upper bounds of σ y , σ z , and h, which will change for each flight. The lower bound is always set to 1 m. That leaves Q, which we set to a lower bound of 0. The upper bound is chosen as a large emitter in comparison to other known sources of this type, in our case coal mine ventilation shafts, and has been set to 31 [kt/ year]. The simulated data was used to determine how many times (N) the optimization was run.
Running the optimization 1000 (N) times, while randomly drawing parameters between the parameter boundaries, ensured that the parameters consistently optimized to within 5% of the known simulated data parameter values. All flights were optimized with N ≥ 1000 for consistency.
2.3.1.1. Uncertainty. The uncertainty of the obtained inverse Gaussian emission rate comes from drawing random starting points for the optimizer within the constrained bounds. The inverse Gaussian is run N times, each time drawing a new random starting point for the optimizer. A normal distribution is fitted to the output of the optimized parameters, and the uncertainty is given as the standard deviation of the normal distribution.

Mass balance approach
The second emission rate quantification method uses a mass balance approach. The first step with this approach is to take the sampled 'curtain' of mole fractions obtained from the active AirCore, and inter-and extrapolate the data into a full two-dimensional grid. This was done using the statistical kriging method. This method estimates, in our case, CH 4 mole fractions at locations where no measurements were taken; mainly in the altitude gaps between the flight transects. Based on this mole fraction curtain, an emission rate corresponding to each flight was estimated by applying a mass balance approach. This is a method that was similarly used by Karion et al. (2013); Krautwurst et al. (2017); Krings et al. (2011Krings et al. ( , 2013Krings et al. ( , 2018 where the output of the emission rate Q is in kg s − 1 , u⋅ Cos(θ) is the wind speed in m s − 1 (with θ being the angle between the wind direction and the perpendicular line of the flight trajectory) and assumed to be constant throughout the flights, k i is the number of horizontal grid boxes in the kriged plane, k j is the number of vertical grid boxes in the kriged plane, M CH4 is the molecular mass of CH 4 in kg mol − 1 , C i,j is the mole fraction of grid box i, j in mol mol − 1 , ΔX is the area of each grid box in m 2 , R is the universal gas constant in J K − 1 mol − 1 , T is the temperature in K, and P i,j is the pressure at each grid box in Pa. The term H 2 O is the mole fraction of water vapor in the air sample in %. As mentioned in section 2.1.3., the wind speed was measured from radiosondes launched during the flights and assumed to be constant during the time duration of the flights (8-12 min).

Uncertainty mass balance.
The uncertainty of the mass balance approach stems mainly from one factor; the wind. The uncertainty of the wind speed is a product of the uncertainty of the wind speed measurements, the atmospheric variability of the wind speed, as well as the wind direction. The wind speed measurement has an uncertainty of approximately 5%, and the standard deviations of the wind speed ranged from 24 to 59% (see Table 3). As the wind speed is a product of the wind speed and wind direction, these are added and summed linearly. The relative uncertainty of wind speed equates to 24-60%. The relative uncertainty of the CH 4 mole fraction measurements is small due to the large signal-to-noise ratio, and is negligible compared to the relative uncertainty of the wind speed. The relative uncertainty that stems from the kriging interpolation is also a small term (Vinkovich et al., in prep). The mass balance uncertainty is therefore only a product of the wind speed, and yields a relative uncertainty for the mass balance approach ranging from 24 to 60%.

Simulated data experiment
The objective of the simulated data experiment is to assess the performance of the two emission rate quantification approaches while knowing the real emission rate of the simulated data and the mean wind field, and to explore and determine the optimum sampling parameters for different plume scenarios. The simulated data is generated based on flight tracks and concentration mappings similar to those of the sampled active AirCore flights. The assumption is that the coal mining ventilation shafts emit as a point source, with a plume following a Gaussian dispersion model as described by eq. (3). The emissions can be regarded as a point source because the distance of the observations away from the shafts are 150-300 m, which is significantly larger than the diameters of the shafts of ~ 5m.

Generating flight tracks
The flight track is generated from user-specified parameters, starting from the point source. Wind directions between 0 and 360 • and the distance from the source in meters are used to project the plume center downwind of the source. A line is then drawn perpendicular to the wind direction to simulate the flight track. The length of the transects and the center of the plume are then specified to fully generate the flight tracks. The pass at each altitude level, along with the incremental increase in altitude for each transect, is then given to derive a full flight profile.

Generating the concentration mapping
Having generated the flight track, CH 4 mole fractions along the flight track are added following the Gaussian distribution model for point sources given by eq. (3). Necessary parameters are provided to generate the concentration mapping. These include the wind speed in m/s, the background dry mole fraction in ppb, the emission rate Q in kg/s, and the center of the plume in m. If the track length is 300 m, the plume center is specified as a number between 0 and 300. CH 4 mole fractions can be generated using the Pasquill stability classes from A -F (Pasquill, 1961). The variability in y-and z-direction, σ y and σ z , along with plume rise h, are calculated from the specified stability class, ambient temperature in K, source temperature in K, and ambient pressure in Pa. CH 4 mole fractions can also be generated by bypassing the Pasquill stability classes. This is done by directly specifying σ y , σ z , and the center height of the plume, H.
Random noise is added to the mole fraction data (up to 100 ppb) and positional coordinates; latitude, longitude, and altitude (up to 1 m). Finally, a moving average is applied to the mole fraction mapping to smooth the data, with a window length equivalent to 34 s, equivalent to a drone speed of 2.13 m/s. The smoothing is applied to the data to simulate the diffusion processes that occur while the active AirCore is sampling and being analyzed. The determination of the averaging kernel has been described in a previous section.
Three simulated data sets were produced to mimic similar characteristics to that of the sampled active AirCore flights. Similar magnitudes of wind speed, down-wind distance, background mole fractions of CH 4 , plume width, plume height, and plume center height to that of the sampled flights were selected. A table showing the parametric values for the generation of each data set (from 1 to 3) can be found in the supplement information, and shown in table A1. The data is shown in Fig. 3 a, d, and g.

Emission rates derived from simulated data
The simulated-data experiment shows that under various realistic sampling conditions and observation-based estimates of the errors, our quantification approaches are capable of retrieving the correct combination of mean emission estimates and dispersion parameters. Table 4 provides the estimated source strength Q from the two quantification approaches described in section 2.3.1 and 2.3.2 and the estimated other parameters of the inverse Gaussian approach. The estimated source strength using the inverse Gaussian approach has a difference of 2.7%, 3.2%, and 3.6% relative to the true value for simulated flights #1, #2, and #3, respectively, with uncertainties ranging between 8.2 and 16.0% of the estimated source strength. For the mass balance approach, this relative differences for the same flights are 9.1%, 1.6%, and 10.0%, respectively, with uncertainties ranging between 25.9 and 50.8%. Fig. 3 shows the generated simulated data sets, the kriging extrapolations, and the calculated Gaussian plume for all simulated data sets. The parameters used to create each of the simulated data flights can be found in appendix A.1. Both the shape and the magnitude of the retrieved CH 4 plumes are similar for the two methodologies. The results suggest that both approaches can accurately retrieve the true source strength, and the inverse Gaussian approach is advantageous over the mass balance approach, especially concerning the significantly reduced uncertainties.
A possible reason for the good agreement between the inverse Gaussian approach and the simulated data is the fact that they both follow the same model, and the wind speed used for the emission rate estimation is the same for both approaches, and set to be the same as for the generated data sets. The generation of the simulated data follows equation (3), although noise has been added, the signal has been smoothed, and the spatial peak position has been shifted. However, given the short averaging time, the Gaussian distribution assumption may not be accurate, and that differences between the assumed Gaussian shape and observed average turbulent plume would negatively affect this result (further discussed in section 3.2.1). Therefore, for real sampled flight data, the estimated emission rate will likely have a larger uncertainty.

Emission rate estimates using individual transects
Using the simulated data, we investigated the benefit of using all transects of a flight simultaneously to estimate one source strength (Q), compared to using one transect at a time to estimate multiple source strengths. The inverse gaussian approach was applied to one of the three transects that have transected the plume for a single flight at a time, and finally compared with the results obtained with the methodology applied to the full flight profile. For this test, simulated data set #1 was used. As can be seen in Table 5, the estimated parameters using all transects compare well with the parameters that were used to create the simulated data, with relative differences between 0.3 and 5.3%.
The likely reason why the methodology using the full flight profile (all transects) estimates parameters closer, and more reliably, to the true Fig. 3. The simulated data (a, d, and g), kriged extrapolation mapping (b, e, and h), and inverse Gaussian curtain (c, f, and i) for the simulated flights 1 (a-c), 2 (d-f), and 3 (g-i) that resemble three selected actual flights.

Table 4
The estimated source strength Q (kt/yr) from both the mass balance approach and the inverse Gaussian approach, and the retrieved model parameters from the inverse Gaussian approach in a simulated data experiment with realistic conditions and uncertainties of the parameters. value than when only using one transect at a time, is due to the data constraining certain parameters in the model. If the top transect only hits a portion of the plume, but the model treats this transect as the center of the plume, it will estimate a smaller horizontal dispersion (σ y ) than what was truly set. This also applies to the center height H, where having several transects will constrain the center to a certain height. This does not happen when only one transect is used, forcing the center height of the plume to placed somewhere along the height of the transect, as seen in the estimated plume center heights in Table 5. The vertical dispersion (σ z ) will be constrained within several transects, but will run free when only one transect is used. Essentially, when using only one transect at a time, σ y , σ z , and H will not be properly constrained, causing unreliable estimations of the emission source strength Q. The use of multiple vertical transects to constrain plume parameters, such as height and plume width, has previously been done for manned aircraft with a primary focus on the use of the mass balance approach Karion et al., 2015;Pitt et al., 2019;Fiehn et al., 2020;France et al., 2020), but has also for UAV-based emission rate quantification been done by Allen et al. (2019) and Shah et al. (2020).

Emission rate estimates using different transect spacing
Using the simulated data, we studied how different flight patterns would affect the emission rate estimates and found that while planning potential flights, the spacing between the transects of the flights, and the altitude at which the flight is planned, should aim to stay within 2.5 times the vertical distribution parameter (σ z ) of the plume. Fig. 4a shows the relative difference between optimized variables (Q, σ y, σ z , and H) and their true value for variable plume center heights. The ratio indicated in Fig. 4a is the summed up relative difference between these four variables as such: Where X i is one of the previously mentioned inverse Gaussian parameters. This was done for 5 different transect spacings, in multiples of σ z .
What can be seen is that the spacing between the individual transects of more than 2.5 times the plume's vertical dispersion will greatly increase the uncertainty of the estimated parameters. Similarly, ensuring that the center transect of the flight is no more than 2.5 times the plumes vertical dispersion away from the center height of the plume will also greatly decrease the uncertainty of the estimated parameters (see Fig. 4b). These results provide a useful guidance in the fight pattern planning to quantify point source emissions. In practice, the magnitude of σ z can be estimated according to Thoma and Squier (2014) prior to a flight. The atmospheric stability indicator (ASI) is an indicator of the atmospheric conditions and the stability class, and can be obtained by having a 3D-sonic anemometer to measure the turbulence intensity (TI) and standard deviation of the wind direction (σθ). From Thoma and Squier (2014), the ASI ranges from 1 (TI > 0.205, σθ > 27.5 • ) to 7 (TI < 0.08, σθ < 7.5 • ). This corresponds roughly to A -G of the Pasquill stability classes, from which σ y and σ z can be estimated.

Emission rates derived from flight data
The quantified CH 4 emission rates from the eight flights fulfilling the criteria posed in section 2.2.2 are shown in Fig. 5a. The flight profiles, as well as the kriged extrapolation plane for the mass balance approach and the inverse Gaussian created concentration mapping are shown in Figs. 6 and 7. Most of the flights suggest a wide horizontal, but narrow vertical spread plume. The estimated CH 4 emission rate on day 1 is 4.9 ± 1.4 [kt/ year] for the inverse Gaussian approach and 3.7 ± 2.3 [kt/year] for the mass balance approach. For day 2, the estimated CH 4 emission rate is 4.3 ± 3.1 [kt/year] for the inverse Gaussian approach, and 6.9 ± 5.0 [kt/year] for the mass balance approach. Comparing the quantified CH 4 emission rates between the inverse Gaussian approach and the mass balance approach seen in Fig. 5a, we get a mean difference of 2.3 [kt/ year] between the two methodologies over the two days. The difference between the quantification methods is 1.2 [kt/year] during the first day, and 2.6 [kt/year] during the second day. There is a significant variation in the estimated emission rates, indicating that the shaft emissions are highly variable, or that the estimates from the quantification methods are not precise for individual flights if there is no or little variation in the actual emissions. Hourly CH 4 emission data from the ventilation shaft demonstrates highly variable emissions. The provided CH 4 emission data comes from hourly shaft data, and is calculated from raw CH 4 concentration measurements and air flux measurements from the shafts. The mean and standard deviation of the CH 4 emission data estimates using hourly data for day 1 and day 2 are shown in Fig. 5a, and are 14.4 ± 4.9 [kt/year] and 8.2 ± 2.9 [kt/year], respectively. A discrepancy is seen when comparing daily averaged quantified emission rates to daily inventories, where the average quantified emission rates on day 1 are lower than the facility emission data, while more in agreement with the facility emission data on day 2. A possible explanation to this could be due to the snapshot nature of the quantified flights. The large variability of the shaft emissions could mean that we sampled during a time of low emission while the facility emission data averages out over the whole day to a higher number. As mentioned in the introduction, the drainage station is located ~3.5 km away from the ventilation shaft, and the occasionally vented CH 4 to the atmosphere from the drainage station has limited influence on our observations. Although the criteria posed in section 2.2.2 have been fulfilled for all these flights, the quantification methodologies do not always work well. The uncertainty of flights #4, #8, #14, #15, and #17 for the inverse Gaussian is of the same order as the quantified emission rates, which indicates that the inverse Gaussian methodology struggles to obtain a good estimate. A possible cause for this might be that there are enhancements of CH 4 in several areas of the 2D concentration mapping, which causes the optimization of the Gaussian parameters to not be able to find a minimum. The mass balance approach can be applied successfully to more of the flights than the inverse Gaussian, however, may face issues with overestimation of the emission rate if the spacing between the transects is too large, as in flight #12. Underestimation of the emission rate can also happen if the enhanced CH 4 is on the outer edge of the concentration mapping plane. This leads to underestimation of the plume, since part of the plume may be effectively missing from the kriged concentration mapping plane which the mass balance approach calculates the emission rate from. The determined emission rate is most likely dominated by the shaft emissions. The closest two ventilation shafts are located 3.66 km and 4.3 km west of the ventilation shaft. To the north the closest ventilation shaft was 8.6 km away. To the east the closest shaft is 20 km away, and no shafts within 20 km are located to the south. As the wind was primarily coming from the south, south-west on day 1 and from the east, northeast on day 2, we do not anticipate any significant influence on the observed CH 4 plumes by other nearby shafts. However, we were unable to rule out potential influence from other parts of the coal mine infrastructure.
As mentioned, the emissions from the ventilation shafts seem to vary throughout the day, reflected in both hourly facility emission data emissions and in the quantified emission rates. During the first day, there is a large difference between the facility emission data and the quantified emission rates. This could be because the flights are 10-min snapshots of the emission rates, whereas the facility emission data is an hourly averaged data that has been averaged into a mean emission for a particular day. There is better agreement on the second day, especially with the two flights where the inverse Gaussian approach performs the best, #12 and #18.
Quantified CO 2 emission rates are shown in Fig. 5b, and show similar variabilities as seen in the CH 4 emission estimates. The CO 2 emission rates are further discussed in section 3.2.2.

Atmospheric turbulence
For the inverse Gaussian approach, we follow the assumption that  T. Andersen et al. the plumes observed at a certain distance to the coal mining ventilation shafts are Gaussian shaped. The Gaussian plume distribution is a result of sufficient averaging, and the required averaging timescales vary depending on the distance to the source (Fritz et al., 2005). The required averaging timescales can be compensated by multiple transects in mobile van measurements (Caulton et al., 2018). Since we use instantaneously observed plumes at different altitudes by a UAV, sufficient averaging timescale at a certain altitude is lacking, which adds significant uncertainty to the estimation. However, the observations from multiple transects at different heights provide additional constraints on the location and the shape of the plume, especially concerning the location of the center of the plume, as is demonstrated in the simulated data experiments in section 3.1.1. As shown in Dinger et al. (2018) during their observation of turbulent dispersion of SO 2 puffs, under cloudy conditions, the direction and the height of the plume can change significantly within a short time span. Although our UAV flights were made under clear sky conditions and under relatively steady wind conditions, we have observed significant variabilities in the wind speed and the wind direction. An example of how much the plume can spread is seen in flight #15, where elevated mole fractions of CH 4 are all over the concentration mapping plane. Profiles such as this make it difficult to use inverse Gaussian approach. Flights #4, #8, #14, #15, and #17 end up with a large uncertainty (close to 100%) using the Inverse Gaussian approach. This is an indication that these profiles are not well optimized for the use of the inverse Gaussian approach, and that the inverse Gaussian approach plume requires a welldefined plume, a steady wind speed, and wind direction. Table 6 shows the estimated dispersion coefficients in y-and z-direction for all flights. In all flights, the horizontal dispersion coefficients (σ y ) are significantly larger than the vertical dispersion coefficients (σ z ), as is also seen in Figs. 6 and 7 a, d, g, and j, the captured plumes from the flights often have wide horizontal distribution and narrow vertical distribution. These stretched plumes are seen in the measurements possibly due to the variability in the wind directions at different altitudes along the plume trajectory, but also in part due to the non-instantaneous sampling of the plume, i.e. sequential transects at different altitudes. It took approximately one and a half minutes for the UAV to pass through a plume twice at adjacent transects. In this time, the turbulent atmosphere could cause the plume to meander in both horizontal and vertical directions. This may lead to the observed plume to look narrow and elongated. From the analytical point of view, from the moment at which

Table 6
The estimated dispersion in y-and z-direction for the flights. sampling with the AirCore starts, to the moment the sampled air has all been analyzed, 20-30 min have passed. Throughout this time, molecular diffusion is present within the AirCore. This leads to smoothed sampled mole-fractions, as described in Andersen et al. (2018). This smoothing leads to the measured plumes looking wider, however, the smoothing effect has been taken into account in the inverse Gaussian approach, leading to unbiased plumes in Figs. 6 and 7 c, f, i, and l.

Correlation between CO 2 and CH 4
The measured CO 2 and CH 4 for all flights are correlated, having an Rsquared value of above 0.69 with N > 124 (see Fig. 8). The extraction of coal deposits is accompanied by CO 2 as well as CH 4 (Swolkien, 2020), which could explain the large correlation. The calculated slope ranges from 3.54 ± 0.25 [ppm CH4 /ppm CO2 ] to 5.27 ± 0.07 [ppm CH4 /ppm CO2 ]. In addition, a flask sample was taken directly from the coal mine ventilation shaft during the span of this campaign. This sample was too high in CH 4 concentration to measure directly on a CRDS analyzer, and was diluted by a factor of 500 to measure the CH 4 concentration. The concentration of CH 4 was 0.15% ± 0.05%. The ratio between CH 4 and CO 2 of this sample was determined to be 4.1 ± 0.9 [ppm CH 4 /ppm CO 2 ], which falls well in the slope range calculated from the individual flights.
Using this range for the slope, and the CH 4 emission rate for each flight, we can calculate the CO 2 emission rate using equation (8) below: Where Q CH4 is the CH 4 emission rate, M CO2 is the molar mass of CO 2 , and M CH4 is the molar mass of CH 4 . Using the CH 4 emission rates found using both the inverse Gaussian approach and the mass balance approach, the calculated CO 2 emission rates for all flights are shown in Fig. 5

Conclusions
A large quantity of CH 4 is emitted to the atmosphere via ventilation shafts of coal mines. To understand and help reduce these emissions, it is important to independently estimate the magnitude of the emissions and verify the reported facility emission data. Using the UAV-based active AirCore system, we have made atmospheric measurements of CH 4 and CO 2 mole fractions downwind of ventilation shafts of a coal mine in the USCB, Poland. Furthermore, we have developed an inverse Gaussian approach for quantifying CH 4 emissions from a point source using the UAV-based observations, and have applied the inverse Gaussian approach as well as the mass balance approach to both simulated data and actual flight data to quantify CH 4 emissions.
The simulated data experiment indicated that both the inverse Gaussian approach and the mass balance approach were capable of estimating CH 4 point source emissions by flying downwind of the source while transecting the plume at multiple altitude levels perpendicular to the wind direction. Simulated data tests showed that the inverse Gaussian approach performed better than the mass balance approach. The retrieved source strengths were within 3.6% and 10% of the true source strength, with the uncertainty ranges of 8.2-16% and 25.9-50.8% for the inverse Gaussian approach and the mass balance approach, respectively. Furthermore, the simulated data experiment showed the importance of multiple transects at different altitudes, which provides sufficient constraints on the parameters of the inverse Gaussian function. Moreover, we found that vertical spacing between the individual transects and the distance between the center height of the plume and the center flight transect were important, and that the inverse Gaussian approach performed the best when the vertical spacing and the distance were smaller than 2.5 times the vertical distribution (σ z ) of the plume.
For the actual flight data, the estimated CH 4 emission rate on day 1 is 4.9 ± 1.4 [kt/year] for the inverse Gaussian approach and 3.7 ± 2.3 [kt/ year] for the mass balance approach. For day 2, the estimated CH 4 emission rate is 4.3 ± 3.1 [kt/year] for the inverse Gaussian approach, and 6.9 ± 5.0 [kt/year] for the mass balance approach. These results show a large variability in the emission estimates, both day-to-day and flight-to-flight, which is likely caused by actual emission variations as is reflected in the hourly facility emission data. All eight flights had a large correlation (R-squared > 0.69) between measured CH 4 and CO 2 , with the slope of the linear fit ranging from 3. Having unstable wind is a limitation for the use of these quantification strategies. The UAV active AirCore is not an instantaneous sampling method, which means that unstable winds can cause the plume to meander and be potentially picked up more than one time in one flight, or perhaps not at all. The inverse Gaussian approach assumes a perfect Gaussian shape of the plume, which is the case with sufficient averaging time. Since we use single-pass transects at different altitudes, sufficient averaging time at a certain altitude is lacking, which adds significant uncertainty to the estimation. It would be of great advantage if more than one UAV could be organized to fly at different altitudes downwind of a point sources, to sample the plume simultaneously and to possibly reach sufficient averaging time. We emphasize that structured flight planning is important, ideally planning the flights with spacings between transects below 2.5 times the vertical distribution of the plume that can be estimated prior to the flight as described by Thoma and Squier (2014). Furthermore, when such measurements are applied to a number of ventilation shafts in a region, e.g. the USCB, the regional budget of CH 4 emissions of the entire region could be estimated.
The UAV-based active AirCore measurement in combination with an inverse Gaussian approach and a mass balance approach has been demonstrated as an effective way of quantifying the source strength of a point source. Due to the flexibility and the mobility of the UAV system, it is capable of quantifying the emissions of point sources that could otherwise be inaccessible. This can be a valuable tool to obtain a better understanding of the emissions from poorly verified ventilation shafts, and to verify and reduce uncertainties of the estimate of greenhouse gas emissions from point sources.

CRediT authorship contribution statement
Truls Andersen: Performed the field work, Analyzed the data and Wrote the manuscript. Katarina Vinkovic: Analyzed the data and Wrote the manuscript. Marcel de Vries: Performed the field work, Analyzed the data and Wrote the manuscript. Bert Kers: Performed the field work, Analyzed the data and Wrote the manuscript. Jaroslaw Necki: Analyzed the data and Wrote the manuscript. Justyna Swolkien: Analyzed the data and Wrote the manuscript. Anke Roiger: Analyzed the data and Wrote the manuscript. Wouter Peters: Designed the research, Analyzed the data and Wrote the manuscript. Huilin Chen: Designed the research, Performed the field work, Analyzed the data and Wrote the manuscript.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.