Methane Emissions: Remote Mapping and Source Quantification Using an Open‐Path Laser Dispersion Spectrometer

Reducing man‐made greenhouse gas emissions depends on the effective detection and location of sources. We present a new method that remotely detects, locates, and quantifies gas emission rates by sequentially steering an optical beam between multiple retro‐reflectors. The novel open‐path laser gas sensor uses Laser Dispersion Spectroscopy (LDS), with seven beams up to 98 meters long deployed across open, flat terrain. LDS offers high precision (10–20 ppb), high dynamic range and linearity, enhanced immunity to atmospheric perturbations, with fast response to probe an area in 3 s. Simultaneous wind and concentration data were collected for four calibrated methane gas release schemes with emission rates of ~1.3 kg/hr. The resulting data were processed using a Bayesian, Markov chain Monte‐Carlo inverse solver to locate the sources and quantify their mass emission rates and uncertainty bounds. All the sources were located to within a few meters and mass emission rates established within the associated confidence bounds.


Introduction
Effective action to reduce climate forcing caused by anthropogenic emissions of greenhouse gases (Balcombe et al., 2018) is impeded by the technical difficulty of remotely detecting, locating, and quantifying individual source emission rates. To be broadly applicable such monitoring would need to be automated. Numerous advanced technologies have been developed and deployed to detect leaks; however, without localization and emission rate quantification, remediation efforts cannot be focused on the relatively few substantial emitters that are found to be responsible for most of the gas emitted (Brandt et al., 2014). Among those studies that have systematically evaluated source location and quantification performance (Feitz et al., 2018) (Ravikumar et al., 2019), reliable emission rate quantification remains a major challenge, even for manned operations. Until source emission rates are readily measurable, it will remain impossible to demonstrate effective progress towards reducing greenhouse gas emissions. The incentives proposed to drive emissions reductions, taxes, fines, and emissions trading schemes, are all dependent on the availability of measured and auditable emissions reporting at local and national scales.
We present results from a field-trial in which a series of calibrated gas releases are used to evaluate a new optical remote sensing approach to quantifying the mass emission rates of sources distributed across a large area. Two novel aspects underpin the approach: (1) a prototype multibeam open-path optical gas sensor based on laser dispersion spectroscopy (LDS), which measures high precision path-averaged concentrations (PAC) along multiple separate paths within the area of interest. (2) From these measurements, combined ©2020. The Authors. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. with wind velocity data, median estimates and confidence bounds for the location(s) and mass emission rates of the source(s) were obtained by using a Markov chain Monte-Carlo (MCMC) algorithm to sample the Bayesian posterior.

10.1029/2019GL086725
The experimental evaluation is based on 17 gas releases from diffuse area ground level sources distributed across a flat open area and organized in four distinct release schemes. Experimental data comprise gas PAC along each of seven optical beams. A full spatial measurement cycle lasts~3 s. The wind velocity vector is measured simultaneously and used to derive turbulence intensities in the horizontal and vertical planes that drive a simple atmospheric eddy gas dispersion model. We present the inversion method and results from four different multirelease test configurations of single and double sources, situated inside and outside the beam array. The results include confidence bounds for the inferred source locations and median mass emission rates.

Inferring Source Locations and Their Mass Emission Rates
While simple gas sensors estimate concentration at a single point location, open-path sensors provide a measurement averaged along the optical beam and report this as a PAC, typically quoted in parts per million (ppm). PAC is a more stable (i.e., low temporal variance) measurement than point data and better characterizes how much gas is present within an area of interest (Alden et al., 2019). The open-path beam samples concentration more extensively and is more likely to encounter a plume of gas dispersing from an unknown source location. By sequentially steering the optical beam along multiple different paths, spatially extensive, high-precision PAC data are collected. When combined with appropriately sampled wind data and a suitable gas dispersion model, the locations and mass emission rates of the sources responsible for the measured data can be inferred. Figure 1 illustrates the key elements of the experimental method.
Using Bayesian analysis we can solve the inverse problem of where sources are and what the mass emission rates need to be to match the measured PAC and wind data. The inference does not seek the best fit possible but fits that are consistent within the limited precision of the data. The area of interest is 120 × 120 m 2 ; it is divided into a mesh of 22 × 22 cells with a 2 × 2 m 2 area source of unknown-but positive-emission rate at the center of each cell. A Gaussian plume atmospheric eddy dispersion model (Hirst et al., 2013) serves as the forward model in which the horizontal and vertical opening angles are driven by the respective turbulence intensities (Draxler, 1976). The dispersion model provides predictions of the ensemble averaged PAC values that would result for each of the mass emission rate source patterns proposed; it does this for the meteorological conditions prevailing for each individual measurement. The low concentration gas (so no buoyancy effects) is passively transported by advection and diluted by turbulence in the atmosphere, which is a stochastic process. Effective gas dispersion prediction is fundamentally dependent on representative measurements of wind direction, wind speed, and turbulence intensities and the use of an appropriate gas dispersion model. To be representative of the actual steady-state plume, such a simple model requires the measured data to be temporally averaged over a time-scale representative of the gas transit time from source to detection.
Knowing the locations of the instrument and the retroreflectors, each beam's PAC c i (t,s j ) at a given time can be calculated from the ensemble average concentration map provided by the gas dispersion model, for proposed source strength s j . The model is simple and computationally light but assumes the following: (1) The terrain is flat and open; (2) data are time averaged to reflect the effect of atmospheric turbulence; (3) the source locations are fixed and emission rates stable and positive; (4) the gas is conserved during transport from source to measurement location; and (5) the height of the Atmospheric Boundary Layer (ABL) is much greater than the dispersion distances considered.
In the absence of any methane sources discernible by the gas sensor, the measurements of concentration within the area of interest would comprise a slowly varying well-mixed background concentration b(t); this is assumed to be spatially uniform at the scale of these tests. In addition, to represent measurement biases and noise, each PAC has been assigned a constant offset d i and measurement noise ε i (t). As a result, the PAC y i measured at beam i is written as equation 1.
The model can be rewritten using matrix notation and algebra as described in the supporting information.
Solving the inverse problem consists of finding combinations of source locations and emission rates s j that adequately explain the observed temporal evolution of PAC data given the wind data and all associated experimental uncertainties.
In Bayesian inversion, candidate solutions are obtained through exploration of the posterior p(s, b, d, λ|y) given by equation 2 and further described in the supporting information.
Evaluating the posterior using direct numerical integration would be too computationally onerous (Diaconis & Ylvisaker, 1979). Instead, the posterior distribution is sampled using a MCMC scheme (Robert, 2015). The conditional probability distributions of each of the variables are sampled in turn, given fixed values of all the others. The source posterior distribution is sampled using a Metropolis Hastings type MCMC scheme, with gradient-based proposals (Girolami & Calderhead, 2011), while the other parameter distributions are characterized using a Gibbs sampler (Geman & Geman, 1984).

Laser Dispersion Open-Path Spectroscopy
LDS principles and demonstration have been described in the literature (Wysocki & Weidmann, 2010), as well as a multibeam open-path deployment (Plant et al., 2015). The open-path LDS sensor deployed for this experiment has been discussed in detail previously (Daghestani et al., 2014), with the exception that here the operating wavelength was 3.33 μm rather than 7.8 μm. To that end, an intraband cascade laser (ICL) instead of a quantum cascade laser was used to trade off optical power for lower power consumption.
To summarize, unlike optical absorption spectrometers, the LDS instrument additionally measures the molecular optical dispersion occurring in the vicinity of a molecular resonance. At the macroscopic level, optical dispersion can be described as a slight change of refractive index occurring in the gas medium owing to the absorption processes. LDS measures a differential change in refractive index based on phase changes Figure 1. Illustration of the method for the case of source location 2. As the wind changes direction, speed, and turbulence intensity, the plume from a persistent gas source changes orientation and shape and intersects subsets of the multiple measurement beams in ways that allows the source location and median mass emission rate to be inferred.
in the propagating electromagnetic fields from which the molecular number density can be deduced. Common phase variations induced in the propagating medium (such as atmospheric scintillation for instance) are suppressed. Owing to the linearity of the optical dispersion, improved selectivity and large measurement dynamic ranges are obtained, which is well suited to measuring a wide range of potential emission rates and ranges to source in a complex mixture. With optimized settings, the measured signal has no offset, which facilitates optimization of the analogue to digital converter dynamic range and ensures baseline stability.
The LDS sensor is deployed in a 7-beam open-path configuration as illustrated in Figure 1. The ICL operates in the ν 3 ro-vibrational band (stretching mode) of methane, targeting an intense transition (a doublet at 2,998.9940 and 2,999.0604 cm −1 with intensities of 4.57×10 −20 and 3.05×10 −20 cm −1 /(mole/cm −2 ), respectively). These transitions are more than 100 times stronger than those from the near infrared harmonic bands (1.6 microns). From the ICL, the beam is collimated and directed to an acousto-optical modulator in order to produce two fully coherent laser beams separated by a known frequency difference. These beams are then superimposed and directed to distant retro-reflectors. To provide rapid, multiple, open-path measurements in the horizontal plane, the transmitted/received beams are steered by a one-dimensional galvo mirror, which in this implementation is limited to a ±20°angular amplitude. All seven retro-reflectors are interrogated in~3 s, and the spectral scan acquisition for each open-path lasts 200 ms for a 30 spectra average.
The LDS instrument delivers molecular absorption and dispersion spectra from which the PAC needs to be deduced. To that end, a full physics model describing the instrument outputs and associated noise is used to calculate synthetic spectra and fit them to the experimental ones to derive concentration and associated errors (Daghestani et al., 2014). The single point precision of the PAC measurement is between 10 and 20 ppb yielding a normalized precision of~1 ppm.m/√Hz. The sensor is particularly suited to long-path length measurements (up to~1 km) of high precision PAC.

Controlled Methane Release Tests
The controlled release tests were designed to evaluate the combined performance of the multiple open-path LDS and the inverse solver algorithm for locating and quantifying the source emission rates. The tests were arranged to match the assumptions underpinning the method as well as possible: flat terrain without any major flow perturbations, nevertheless representative of monitoring large, flat, open areas for signs of gas emissions from the ground's surface. The tests were done at the Chilbolton Observatory, Hampshire, UK; the program lasted for 5 days. Instrument and retro-reflectors were installed as illustrated in Figure 1. The longest beam path was 98 m, the shortest 34 m.
The gas sources, laid flat on the ground, comprised 2 × 2 m aluminum frames with a plastic cover sheet perforated with a 1 cm spaced matrix of small holes. This ensures the released gas is distributed evenly and enters the atmosphere with minimal momentum: as would a low flux rate of gas seeping from the ground. Hoses of equal length and flow resistance carried the gas from the cylinders to the release frame(s). The mass emission rates were obtained by logging the gas cylinder weight during the releases (see supporting information).
A 3-D ultrasonic anemometer (circle in Figure 1) was positioned at the beam height (1.6 m) and aligned to true north to a~1°precision. It measured the windspeeds, directions, and turbulence intensities in the horizontal and vertical planes at 20 Hz. It was located close to the center of the beam array, to give a good representation of the average velocity and turbulence that drives the gas dispersion within the test area. All equipment locations were measured to~2 cm precision using a surveying GPS.
We collected data for four different source configurations (squares in Figure 1): a single source within the beam array, a single source at two separate locations outside the beam array, and two simultaneous sources outside the beam array. Given time-and gas-constraints, it was impractical to release gas continuously for one case and wait for the wind to provide a sufficient variety of wind directions and speeds. Instead, we intermittently collected data for each of the four different configurations as the wind conditions evolved during the week. It is crucial to the analysis method that we collect data from different wind directions, including if/ when the wind's direction would result in the plumes missing all the beams. Slow variabilities in background atmospheric concentration are handled in the data analysis code by ensuring its "flexibility" is linked to the elapsed time between successive data points.
An individual test commenced by measuring atmospheric background concentration for about 10 minutes before starting the release. We then allowed 2-3 minutes for the plume to establish "steady state" across the area of interest. The gas release was then sustained for up to an hour. After stopping the release, we continued collecting data as the residual plume dispersed and the PAC returned to background levels. The total release measurement duration for the week was 13h25min.

Temporal Concentration Series and Data Analysis
The experimental data comprises temporal series of methane PAC, measured over the seven beam paths, as well as associated 3-s averaged wind vectors. Figure 2 shows the data set for release 3. The high temporal resolution data resolves the stochastic fluctuations of CH 4 PAC along each of the beams.
The data analysis starts by concatenating the experimental data sets for each of the four source configurations. Wind and PAC data were post-processed to obtain 1-minute averages; this is representative of the air transit time across the test area and is used for averaging of the Gaussian dispersion plume model. Those data with average windspeeds of <2 m/s were excluded because the gas dispersion model becomes less reliable at lower windspeeds.
Each cell from the spatial gridding of the test area is considered to have an unknown strength diffuse gas source of 2 m width at its center; its emission rate is permitted to be positive, up to a defined maximum value, which in this analysis was set to 3 kg/hr. The area of interest and the mesh size were chosen to ensure all sources could be captured (given transport conditions) while providing a reasonable trade-off between MCMC sampling time and mesh resolution. The maximum source emission rate is derived from a conservative estimate of expected rates; inappropriate choices become evident in the resulting solutions and associated coverage plots.
To infer a realistic solution, the temporal evolution of the atmospheric background concentration needs to be well represented-so the contributions due to the sources are precisely identified. This local background methane concentration is driven by changes in the ABL during the day. Shallow and poorly mixed in the morning, solar insolation drives mixing and the growth of the ABL, reducing the local CH 4 background concentrations until the early afternoon minimum is reached. In addition to this slow diurnal variation of a few ppb/hr, regional transport can also introduce well mixed higher methane concentration air from very remote sources. A detectable, localized, and persistent source will produce high frequency (spatial and temporal) methane concentration signals in contrast to the slowly varying background component. The rate of variation of the background concentration is represented by the background covariance matrix and the associated flexibility parameter η, which is critical in conditioning the inversion. Again, inspection of graphical output from preliminary runs can be used to ensure the background flexibility is appropriately set.
The source sparsity setting γ is the second input parameter critical to the inversion, and it interacts with the background flexibility η. If the background concentration is too "flexible" source-related PAC will wrongly be attributed to the background. Conversely, too "stiff" a background flexibility will produce concentration artifacts that need to be explained by weak (and spurious) sources. Setting these two parameters appropriately is nontrivial. The method was (1) set a high source sparsity parameter and find the background flexibility parameter that provides a convincingly smooth background concentration behavior, as evaluated by early MCMC runs; (2) gradually relax the source sparsity until a stable source pattern appears and persists, providing a stable burn-in result; (3) last, the background flexibility is reviewed. In the case presented here, the atmospheric background was very stable, and the main adjustments needed were on the sparsity parameter to determine the minimum total emissions required to explain the data.
In each run of the MCMC algorithm, a large number (~100,000) of "burn-in" iterations are run to allow the chain to converge to stability. Convergence is assessed through inspection of "trace plots" of parameter values versus iteration number: Any trends in parameter values indicate non-convergence, suggesting the need for more burn-in iterations or re-evaluation of run settings. The acceptance rate of the source emission rate parameter proposals is also used to assess performance of the algorithm: Too low (<25%) an acceptance rate may indicate too high a proposal step size (giving poor proposals), while the opposite case (acceptance rate > 90%) may indicate too low a step size, which gives poor exploration of the solution space.

Mapping and Quantification Results
The results for one of the four release cases are shown in Figure 3, where the actual source was located within the beam array. Figure 3a shows the mapping and quantification results. The location was obtained with high confidence, as estimated from Figure 3b that shows the associated interquartile range (IQR), which reflects how tightly grouped the post burn-in estimates of source emission rate were for each source location. Unlike other inference methods, MCMC provides a group of solutions from which uncertainty bounds are derived and are expressed here as IQR. The cell coinciding with the actual location (pink "×") was found to emit 0.85 kg/hr (with the lowest IQR of 0.13 kg/hr) whereas the satellite sources were 0.1 kg/hr (IQR 0.22 kg/hr) and 0.45 kg/hr (IQR 0.32 kg/hr). An artifactual "spread" of the source is observed, while the total emission rate inferred is close to the actual value (1.40 kg/hr inferred for an actual of 1.38 ± 0.16 kg/hr). Figure 3c shows the source detectability map derived from the coupling parameter A ij (t). It provides an indicator of source detectability as a function of location, given the wind data collected and detector sensitivity. The white area-also present in Figure 3a-reveals that no <3 kg/hr sources would be detectable there due to the absence of suitable winds to bring dispersing plumes from those locations across any of the beams.
The data used for this result came from a cumulative release duration of 6h03min. The corresponding wind variability was 3.12 m/s for velocity and 161 degrees for direction (variability given as the difference between the 10th and the 90th percentile). The quality of the solution was also assessed by inspecting the concentration residuals of the 1-minute averaged PAC measured for each of the beams. Except for one beam, the statistical distributions of the residuals are very close to a normal distribution suggesting a fully evolved solution with limited data conflicts.
The MCMC analyses were carried out similarly for the three other release configurations. The solutions are shown in Figure 4 and summarized in Table S3 in the supporting information. For the case in Figure 4a, the solution source is located in exactly the correct cell. The inferred median emission rate is 0.72 kg/hr with an IQR of 0.50 kg/hr, while the averaged experimental release rate was 1.36 kg/hr. The cumulative duration of the release was 2h19min, made up from three individual tests; the horizontal wind direction and velocity variability were 100 degrees and 1.92 m/s, respectively. Figure 4b corresponds to another case of a single source outside the beam array, again located in the correct cell. The inferred median emission rate is 1.45 kg/hr with an IQR of 0.40 kg/hr. The experimental emission rate was 1.47 kg/hr. The cumulative release duration was 57 minutes with a wind direction and speed variability of 58 degrees and 2.14 m/s. Lastly,   Figure 4c shows the case for two sources emitting simultaneously, both located outside the beam array. The inferred locations are consistently offset by 6 and 12 m from the true ones, which suggests a systematic bias.
The inferred median emission rates were 1.59 and 1.53 kg/hr with IQRs of 0.20 and 0.50 kg/hr, respectively. The actual emission rates were 1.23 kg/hr per source over a series of four release tests. The cumulative release duration was 3h03min, and the wind speed and direction spread were 2.86 m/s and 99 degrees.

Conclusions
Controlled methane release sources of~1.3 kg/hr within a wide (120 × 120 m 2 ) open flat area have been successfully located remotely and their mass emission rates quantified. Single and double sources have been successfully detected and located, including sources outside of the beam array. The mapping and quantification method demonstrates the benefits of combining (1) a multibeam open-path laser sensor using middle infrared laser dispersion spectroscopy together with (2) a Markov chain Monte Carlo Bayesian inference approach to solve the associated inverse problem. For these tests the simple Gaussian plume gas dispersion model proved adequate.