Journal cover Journal topic
Atmospheric Measurement Techniques An interactive open-access journal of the European Geosciences Union
Journal topic
Atmos. Meas. Tech., 12, 4659–4676, 2019
https://doi.org/10.5194/amt-12-4659-2019
Atmos. Meas. Tech., 12, 4659–4676, 2019
https://doi.org/10.5194/amt-12-4659-2019

Research article 02 Sep 2019

Research article | 02 Sep 2019

# Bayesian atmospheric tomography for detection and quantification of methane emissions: application to data from the 2015 Ginninderra release experiment

Bayesian atmospheric tomography for detection and quantification of methane emissions: application to data from the 2015 Ginninderra release experiment
Laura Cartwright1, Andrew Zammit-Mangion1, Sangeeta Bhatia2,a, Ivan Schroder3, Frances Phillips4, Trevor Coates5, Karita Negandhi6,b, Travis Naylor4, Martin Kennedy6, Steve Zegelin7, Nick Wokker8, Nicholas M. Deutscher4, and Andrew Feitz3 Laura Cartwright et al.
• 1School of Mathematics and Applied Statistics, University of Wollongong, Wollongong, Australia
• 2Centre for Research in Mathematics, Western Sydney University, Parramatta, Australia
• 3Geoscience Australia, Canberra, Australia
• 4Centre for Atmospheric Chemistry, School of Earth, Atmospheric and Life Sciences, University of Wollongong, Wollongong, Australia
• 5School of Agriculture and Food, University of Melbourne, Melbourne, Australia
• 6Department of Earth and Planetary Sciences, Macquarie University, Sydney, Australia
• 7CSIRO Oceans and Atmosphere, Canberra, Australia
• 8Department of Industry, Innovation and Science, Canberra, Australia
• acurrently at: School of Public Health, Imperial College London, London, UK
• bcurrently at: Office of Environment and Heritage, Parramatta, Australia

Correspondence: Laura Cartwright (lcartwri@uow.edu.au)

Abstract

Detection and quantification of greenhouse-gas emissions is important for both compliance and environment conservation. However, despite several decades of active research, it remains predominantly an open problem, largely due to model errors and assumptions that appear at each stage of the inversion processing chain. In 2015, a controlled-release experiment headed by Geoscience Australia was carried out at the Ginninderra Controlled Release Facility, and a variety of instruments and methods were employed for quantifying the release rates of methane and carbon dioxide from a point source. This paper proposes a fully Bayesian approach to atmospheric tomography for inferring the methane emission rate of this point source using data collected during the experiment from both point- and path-sampling instruments. The Bayesian framework is designed to account for uncertainty in the parameterisations of measurements, the meteorological data, and the atmospheric model itself when performing inversion using Markov chain Monte Carlo (MCMC). We apply our framework to all instrument groups using measurements from two release-rate periods. We show that the inversion framework is robust to instrument type and meteorological conditions. From all the inversions we conducted across the different instrument groups and release-rate periods, our worst-case median emission rate estimate was within 36 % of the true emission rate. Further, in the worst case, the closest limit of the 95 % credible interval to the true emission rate was within 11 % of this true value.

1 Introduction

Methane (CH4) is an important transition fuel for decarbonisation of the global energy system . As countries increase the renewable energy mix into their existing electricity networks, CH4 can firm up network stability and supply . Utilisation of biogas or natural gas with carbon capture and storage offers a lower cost pathway to achieve deep decarbonisation targets . One of the disadvantages of CH4, however, is that its global warming potential is much greater than that of carbon dioxide (CO2), so that only a few percent of losses of CH4 into the atmosphere can negate any climate-change mitigation advantages from reducing conventional coal-fired power production . For this reason, it is critical that losses of CH4 along the supply chain are accurately accounted for to ensure public confidence in climate-change mitigation benefits of switching to natural gas. Unfortunately, while several types of instrumentation are available to aid the detection and estimation of fugitive emissions, harnessing acquired data for reliable emission detection and quantification remains a notoriously difficult problem.

Several controlled-release experiments of CH4 and CO2 have been conducted in order to improve techniques for estimating greenhouse-gas emissions . Building on this body of work, in 2015 a CH4 and CO2 controlled-release experiment was held at the Ginninderra Controlled Release Facility in Canberra, Australia . This large multidisciplinary, multi-institutional blind-release trial (i.e. the participants did not know the true release rate) simultaneously assessed eight different CH4 emission-rate estimation techniques, using data from both mobile and stationary instrumentation. These eight techniques included tracer ratio techniques, backwards Lagrangian stochastic modelling, forward Lagrangian stochastic modelling, Lagrangian stochastic footprint modelling, and atmospheric tomography techniques. A full description of the methods and results is given in .

Every group involved in the analysis presented in used a unique combination of instrumentation and estimation technique when carrying out the analysis, making it hard to establish the respective merits (or otherwise) of the employed techniques from the inversion results. Nonetheless, an interesting observation from the study is that none of the eight techniques deployed during the blind-release trial had a leakage uncertainty range (95 % interval) that included the true emission rate, while some estimates (including one obtained using atmospheric tomography) were factors of 2 or more off from the true value. Given that atmospheric methane concentration and meteorological instrument measurement uncertainty is generally low for each of the different approaches, it suggests that the techniques that were used did not adequately account for the variability of atmospheric measurements or the uncertainty introduced through parameterisation of atmospheric mixing conditions (e.g. Monin–Obukhov lengths and/or Pasquill stability classes; see Sect. 3.1) and atmospheric dispersion/transport model uncertainty.

A number of studies have highlighted the importance of atmospheric-model error in estimating emission rates or fluxes (e.g. Chevallier et al.2010; Basu et al.2018). For example, showed that flux estimates are sensitive to the chosen spatio-temporal resolution of the fluxes and the chosen transport model. Uncertainty in the meteorological fields driving the transport model is also known to play a big role (e.g. Miller et al.2015). While ensemble inversions are frequently used to highlight the sensitivity of the results to atmospheric models and meteorological fields, learning unknown parameters associated with transport concurrently with the emission rate is not often done. This is largely due to the computational implications of such an approach. Key here are the use of surrogate models (or emulators) to obtain simplified transport representations. For example, use decision/regression trees as a surrogate for FLEXPART-WRF (a sophisticated atmospheric transport model integrating weather research and forecasting into a Lagrangian particle dispersion model), which allows for quick simulation at various parameter settings that can in turn be used to make inference. For the Ginninderra data we employ the more traditional Gaussian plume model, which can be seen as a surrogate for a full-blown transport model. While known to work well in the small domain (an area of approximately 100 m×100 m) setting we consider (e.g. Riddick et al.2017), importantly this plume model is quick to simulate from, giving us the opportunity to calibrate it while estimating the emission release rate (e.g. Borysiewicz et al.2012). As we see in our sensitivity analysis of our results in Sect. 6, online plume-model calibration is crucial for obtaining accurate emission-rate estimates with our data.

The transport model plays an important role in inverse modelling. Calibration of the transport model from observations can be done within the classic inverse theory framework of . This framework is in turn seated within a Bayesian paradigm, which underpins several of the inversion systems in place today (see, for example, ). Inference in such cases is often done using sampling techniques such as Markov chain Monte Carlo (MCMC) or importance sampling (e.g. Rajaona et al.2015). Quick evaluation of the transport/dispersion model (or surrogate) is crucial when repeatedly evaluating it within an MCMC framework; the Gaussian plume model is hence a popular choice in these frameworks (e.g. Jones et al.2016; Wang et al.2017). MCMC is also our method of choice for Bayesian atmospheric tomography, because it allows relatively easy computation of posterior distributions of parameters that are deeply nested within a hierarchical model. It is also ideally suited for the case of point-source emissions, where the dimensionality of the latent space is low (unlike, for example, when performing regional emission quantification).

Atmospheric tomography, a term inspired from medical imaging, combines data from a collection of measurement sites with Bayesian inversion to detect and quantify emissions. The primary contribution of this article is an extension of the atmospheric tomography technique described in Sect. 2.4.2 of . In , atmospheric tomography was only used on one type of instrument and did not account for uncertainty in the transport model. The technique we propose accounts for uncertainty in our data, in our process models, and in our parameters; is applicable to both point- and path-sampling instruments; and takes into account instrument-specific bias. Inference is made on all unknown parameters using MCMC, and uncertainty in the transport-model parameters are propagated to our posterior inferences on the release rate. We demonstrate the efficacy and utility of the unifying Bayesian framework on data from point- and path-sampling instruments used in the Ginninderra experiment. A secondary contribution is the curated provision of a data set containing a large portion of the Ginninderra data at a 5 min resolution, which we hope will serve as a resource for other researchers to validate their own emission-rate estimation techniques on. The data and scripts required to reproduce the results in this article are available from https://github.com/Lcartwright94/BayesianAT (last access: 1 August 2019).

The remainder of the article is organised as follows. Section 2 gives an overview of the experimental setup and the data collected during the 2015 Ginninderra experiment. Section 3 describes the atmospheric transport model used, while Sect. 4 details the hierarchical model we employ and the Bayesian methodology we develop for emission-rate estimation. Section 5 gives the results from application of our Bayesian atmospheric tomography technique on the Ginninderra data. Section 6 examines how our results would change if certain components in our model (e.g. relating to the plume model) are (erroneously) assumed fixed and known. Section 7 concludes.

2 The 2015 Ginninderra release experiment

A full description of the experimental setup, measurement techniques, and quantification methods used in the 2015 Ginninderra release experiment are given in . Briefly, CH4 (together with CO2 and nitrous oxide) was released from a small chamber located in a fallow agricultural field from 23 April to 12 June 2015 and from 23 to 24 June 2015. A variety of CH4 sensors were placed around the release chamber. The measurement data considered in this study were obtained from two Picarro G2201-i analysers (positioned in the predominant upwind (NW) and downwind (SE) location of the release chamber, labelled Picarro.West and Picarro.East, respectively), four eddy covariance (EC) towers equipped with Li-COR 7700 open-path CH4 sensors (labelled EC.A, EC.C, EC.D, and EC.E, respectively), two scanning Fourier-transform infrared (FTIR) spectrometers with four retro-reflectors terminating six measurement paths (labelled P1 to P6, respectively), and a scanning GasFinder 2 Boreal laser with seven reflectors forming seven measurement paths (labelled R1 to R7, respectively); see the left panel of Fig. 1. Meteorological data were collected from EC.A equipped with a Vaisala HMP50 relative humidity and temperature sensor, a CSI EC150 CO2H2O sensor, a Li-COR 7700 CH4 sensor, a Kipp and Zonen CNR4 radiometer, a CSI CSAT3 sonic anemometer, and a Gill WindSonic anemometer. Wind speed and wind direction were measured by the CSAT3 sonic anemometer and the Gill WindSonic anemometer. As part of data quality control, horizontal wind speed and wind direction data from the two instruments were compared, with no arising issues. Both sonic anemometers were using factory calibration. Wind directions were determined by manually aligning the sonic anemometers so that the reference direction was true north. Data from CSAT3 sonic anemometer were logged at 10 Hz and data from the Gill WindSonic anemometer at 1 Hz.

The gases were released at a height of 0.3 m, and the standard CH4 release rate was 5.8 g min−1, limited mostly to daylight hours. On brief occasions, the CH4 release rate was varied between 2.9 and 20 g min−1 to enable testing of mobile CH4 sensor platforms. Towards the end of the experiment (8–12 June), the CH4 release rate was decreased from 5.8 to 5.0 g min−1, and the setup for the Boreal laser measurements was modified with the number of retro-reflectors and paths reduced to six (labelled R8 to R13, respectively; see the right panel of Fig. 1). The location of all other CH4 sensors did not change over the duration of the experiment. The Picarro analysers were not deployed until 21 May, and the CH4 release rate on 23 and 24 June was constantly varied. Hence, in this article we only consider data between 21 May and 12 June, excluding 26 and 27 May where the release rate was also constantly varied.

Figure 1(a) Layout of instruments in the 2015 Ginninderra release experiment between 21 May and 7 June 2015. (b) Layout of instruments between 8 and 12 June 2015. R1 to R13 are the paths formed between the Boreal laser and reflectors; P1 to P6 are the paths formed between the FTIR spectrometers and retro-reflectors; EC.A to EC.E are the EC towers; and Picarro.East and Picarro.West are the Picarro analysers. All coordinates are relative to EC.A, which is situated at the origin.

The data set used to obtain the results presented in Sect. 5 was compiled by pooling together the separate meteorological and concentration data sets used in the Ginninderra experiment. A common resolution of 5 min was chosen; that is, all measurements of concentration and meteorological variables were averaged over a regular set of 5 min intervals. Measured CH4 concentrations were then matched with corresponding meteorological measurements by time and placed into long-table format, with each row corresponding to a unique data point. For path measurements, two extra columns were used to denote the end-point coordinates of the paths.

Initial preprocessing was carried out to provide a complete data set without outliers. First, data containing missing values considered critical for emission-rate estimation (in particular, air temperature, air pressure, wind speed, and wind direction) were removed from the data set. Second, data points corresponding to upwind measurements that were more than three median absolute deviations away from the instrument's median upwind measured concentration were determined to be outliers and hence removed. A point measurement was classified as upwind if the angle subtended from the source by a line joining the instrument location to the plume centreline was more than 45. A path measurement was classified as upwind if the angles subtended at every point along the path were more than 45.

3 Transport modelling

In this section we detail the plume model employed and how it is used to supply model-predicted concentrations for the path measurements.

## 3.1 Gaussian plume dispersion modelling

As outlined in Sect. 1, we use a transport model that is simply parameterised, and easy to evaluate, so that it can be calibrated online. One of the simplest models that works well on the short distances we consider is the Gaussian plume dispersion model (e.g. Wark et al.1998, chap. 4). Here the true emission rate is denoted by Q in grams per second (g s−1), the height of the CH4 point source by H in metres (m), and the total number of observations by N. The classic Gaussian plume model is given by

$\begin{array}{}\text{(1)}& \begin{array}{rl}& C\left({x}_{i},{y}_{i},{z}_{i},Q,{U}_{i},H,{\mathbit{\theta }}_{{k}_{i}}\right)\\ & =\frac{Q}{\mathrm{2}\mathit{\pi }{U}_{i}{\mathit{\sigma }}_{{y}_{i},{k}_{i}}{\mathit{\sigma }}_{{z}_{i},{k}_{i}}}\mathrm{exp}\left(-\frac{{y}_{i}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{{y}_{i},{k}_{i}}^{\mathrm{2}}}\right)\\ & \left[\mathrm{exp}\left(-\frac{\left({z}_{i}-H{\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{{z}_{i},{k}_{i}}^{\mathrm{2}}}\right)+\mathrm{exp}\left(-\frac{\left({z}_{i}+H{\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{{z}_{i},{k}_{i}}^{\mathrm{2}}}\right)\right],\end{array}\end{array}$

where C is the model-predicted concentration in grams per cubic metre (g m−3) of CH4 at a single spatial point $\left({x}_{i},{y}_{i},{z}_{i}\right)$ in metres (m) along the direction of the plume corresponding to the ith measurement, Ui is the wind speed associated with the ith measurement in metres per second (m s−1), ${k}_{i}\in \mathit{\left\{}A,B,C,D,E,F\mathit{\right\}}$ represents the Pasquill stability class (a categorisation reflective of the expected level of horizontal and/or vertical spread of the atmospheric particles after emission; see Pasquill1961) associated with the ith measurement, and ${\mathbit{\theta }}_{{k}_{i}}$ represents plume-specific parameters used to construct the standard deviations ${\mathit{\sigma }}_{{z}_{i},{k}_{i}}$ and ${\mathit{\sigma }}_{{y}_{i},{k}_{i}}$. These standard deviations of the plume in the vertical and horizontal directions are given by

$\begin{array}{rl}{\mathit{\sigma }}_{{z}_{i},{k}_{i}}& ={a}_{{k}_{i}}{x}_{i}^{{b}_{{k}_{i}}},\\ {\mathit{\sigma }}_{{y}_{i},{k}_{i}}& =\mathrm{0.4651}{x}_{i}\mathrm{tan}\left({\mathit{\nu }}_{i}\right),\end{array}$

respectively, where ${\mathit{\nu }}_{i}=\mathrm{0.01745}\left({c}_{{k}_{i}}-{d}_{{k}_{i}}\mathrm{ln}\left({x}_{i}/\mathrm{1000}\right)\right)$. Note that the coefficients ${a}_{{k}_{i}},{b}_{{k}_{i}},{c}_{{k}_{i}},{d}_{{k}_{i}}$ correspond to the ith measurement and depend on the stability class associated with that measurement, ${k}_{i}\in \mathit{\left\{}A,B,C,D,E,F\mathit{\right\}}$. Values for these coefficients by stability class are given in Wark et al. (1998, chap. 4) and shown here in Table 1 for completeness. We collect the plume-specific parameters in ${\mathbit{\theta }}_{{k}_{i}}\equiv \left({a}_{{k}_{i}},{b}_{{k}_{i}},{c}_{{k}_{i}},{d}_{{k}_{i}}{\right)}^{\prime }$, where denotes the transpose operator.

Table 1Stability classes to which observations within the Ginninderra experiment are allocated, and the corresponding values of ${a}_{{k}_{i}},{b}_{{k}_{i}},{c}_{{k}_{i}},$ and ${d}_{{k}_{i}}$ used to construct the horizontal $\left({\mathit{\sigma }}_{{y}_{i},{k}_{i}}\right)$ and vertical $\left({\mathit{\sigma }}_{{z}_{i},{k}_{i}}\right)$ standard deviations of the plume when xi is in metres (m).

The stability class to which an observation is allocated is classically based on (i) the Monin–Obukhov length (the theoretical height at which turbulence is produced by buoyancy and mechanical forces in equal amounts; see , chap. 16) and (ii) an effective roughness length. The Monin–Obukhov length (L value) is given by $L=-{u}_{*}^{\mathrm{3}}{\stackrel{\mathrm{‾}}{\mathit{\xi }}}_{v}/\left(qg\left(\stackrel{\mathrm{‾}}{{w}^{*}{\mathit{\xi }}_{v}^{*}}{\right)}_{s}\right)$ (Jacobson2005, chap. 8), where u* is the frictional velocity, ${\stackrel{\mathrm{‾}}{\mathit{\xi }}}_{v}$ is the mean virtual potential temperature, $\left(\stackrel{\mathrm{‾}}{{w}^{*}{\mathit{\xi }}_{v}^{*}}{\right)}_{s}$ is the surface virtual potential temperature flux, q is the von Kármán constant, and g is the acceleration due to gravity. In our case we used WindTrax (http://www.thunderbeachscientific.com/windtrax.html, last access: 27 March 2019) to determine the L value for each observation; we provide the L values with the compiled data. We set the effective roughness length z0=0.01 m, corresponding to a relatively flat area, with short or no grass, and minimal buildings/trees/other obstacles; see Sienfeld and Pandis (2006, chap. 16) and World Meteorological Organisation (2008, chap. 5). This is a suitable choice for the Ginninderra site. We used the results of to allocate a stability class to each observation based on the L values provided by WindTrax and z0=0.01 m.

The coefficients typically used for each stability class could be off by a factor of 2 or more (Wark et al.1998, chap. 4). To show that this is also the case with our categorisation scheme, in Fig. 2 we show the Gaussian-plume-model-predicted outputs together with the observed data enhancements (see Sect. 4.1) at one of our measurement locations (namely, EC.A) between 21 May and 7 June, when scaling ${\mathit{\sigma }}_{{y}_{i},{k}_{i}}$ by 1, 2.5, and 4, respectively. Clearly, with no scaling the predicted plume is too narrow, while with a scaling of 4 it is too broad. A scaling of 2.5 gives good agreement. Importantly, since in Eq. (1) Q only serves to scale the predicted concentrations (i.e. make them larger or smaller by a constant factor), it is apparent that this plume-scaling factor is identifiable, in the sense that we can learn it from the data while estimating the emission rate (provided the source is active). Online plume-model calibration fits naturally within the MCMC framework discussed in Sect. 4.

Figure 2Predicted (blue) and observed (red) enhancements in parts per million (ppm) at EC.A between 21 May and 7 June 2015 when scaling ${\mathit{\sigma }}_{{y}_{i},{k}_{i}}$ by 1, 2.5, and 4, respectively. The mean-squared errors (MSE) between the observed and the predicted enhancements are also shown. Of the three, the best agreement between predicted and observed values occurs when ${\mathit{\sigma }}_{{y}_{i},{k}_{i}}$ is scaled by 2.5.

## 3.2 Low wind speeds

It is well known that the Gaussian plume model is less accurate for low wind speeds (e.g. Turner1994, chap. 2). One reason for this is that the wind speed Ui is in the denominator of the scaling coefficient of Eq. (1); hence, the plume model prediction becomes very sensitive to Ui as it tends towards zero. This is problematic as Ui, although often assumed known, is an average calculated from noisy measurements taken over some time span (in our case 5 min) and is thus itself noisy. From Eq. (1) we see that, when conditioned on all other parameters, the variance of Ci is proportional to the variance of the inverse of Ui, which can be very large for small Ui. Instead of removing data at low wind speeds as is often done (e.g. Feitz et al.2018), we analyse the theoretical relationship between the variance of the inverse wind speed and Ci. We then use this relationship to discount low-wind-speed model predictions in the Bayesian framework in a principled manner. While the analyst still needs to choose a cutoff below which to model this relationship, in separate studies we found that our inferences are not particularly sensitive to the chosen cutoff. Moreover, we found that downweighting instead of excluding was necessary for making inference when not many observations associated with high wind speeds were available.

Each wind speed Ui is an average of a number of wind speeds (say ${n}_{{U}_{i}}$) recorded over 5 min. Therefore Ui is a sample mean and thus an unbiased estimator of the true (population) mean wind speed, say ${\mathit{\mu }}_{{U}_{i}}$, over this time interval. By the central limit theorem, $\sqrt{{n}_{{U}_{i}}}\left({U}_{i}-{\mathit{\mu }}_{{U}_{i}}\right)\stackrel{\mathcal{D}}{⟶}\text{Gau}\left(\mathrm{0},{\mathit{\sigma }}_{{U}_{i}}^{\mathrm{2}}\right),$ where 𝒟 implies convergence in distribution, Gau(μ,σ2) denotes the Gaussian distribution with mean μ and variance σ2, and ${\mathit{\sigma }}_{{U}_{i}}^{\mathrm{2}}$ is the variance of the wind speeds over the ith time interval, which was derived from the raw (disaggregated) data. We can then use the delta method (e.g. Casella and Berger2002, chap. 5) to deduce that

$\sqrt{{n}_{{U}_{i}}}\left(\frac{\mathrm{1}}{{U}_{i}}-\frac{\mathrm{1}}{{\mathit{\mu }}_{{U}_{i}}}\right)\stackrel{\mathcal{D}}{⟶}\text{Gau}\left(\mathrm{0},{\left(\frac{d}{\text{d}{\mathit{\mu }}_{{U}_{i}}}\left(\frac{\mathrm{1}}{{\mathit{\mu }}_{{U}_{i}}}\right)\right)}^{\mathrm{2}}{\mathit{\sigma }}_{{U}_{i}}^{\mathrm{2}}\right).$

Hence, the variance of 1∕Ui is approximately

${\left(\frac{d}{\text{d}{\mathit{\mu }}_{{U}_{i}}}\left(\frac{\mathrm{1}}{{\mathit{\mu }}_{{U}_{i}}}\right)\right)}^{\mathrm{2}}{\mathit{\sigma }}_{{U}_{i}}^{\mathrm{2}}=\frac{\mathrm{1}}{{\mathit{\mu }}_{{U}_{i}}^{\mathrm{4}}}{\mathit{\sigma }}_{{U}_{i}}^{\mathrm{2}}\propto \frac{\mathrm{1}}{{\mathit{\mu }}_{{U}_{i}}^{\mathrm{4}}}.$

Therefore, conditional on all other terms in Eq. (1), the variance of the model-predicted concentrations increases as a quartic of the true inverse wind speed. This is important, as it means that model predictions at low wind speeds, say less than 1 m s−1, could be highly uncertain; we show a way of handling this uncertainty when we detail the Bayesian inversion model in Sect. 4.

## 3.3 Predicted concentrations for point and path measurements

The plume model given by Eq. (1) sets the x axis as its centreline and the CH4 source at the origin. The predicted plume-model concentration at a physical location $\left({\stackrel{\mathrm{̃}}{x}}_{i},{\stackrel{\mathrm{̃}}{y}}_{i},{\stackrel{\mathrm{̃}}{z}}_{i}\right)$ is thus found by first applying a spatial shift and time-dependent rotation (by wind direction) to $\left({\stackrel{\mathrm{̃}}{x}}_{i},{\stackrel{\mathrm{̃}}{y}}_{i},{\stackrel{\mathrm{̃}}{z}}_{i}\right)$ in order to obtain $\left({x}_{i},{y}_{i},{z}_{i}\right)$, which is then used to compute a model-predicted concentration (conditional on Q, Ui, H, and ${\mathbit{\theta }}_{{k}_{i}}$). Conversion to parts per million (ppm) is done via the ideal gas law.

Let Ci be a model-predicted concentration ($i=\mathrm{1},\mathrm{2},\mathrm{\dots },N$). If Ci corresponds to a point measurement, then one needs only to evaluate Eq. (1) at the transformed point-measurement location to obtain a predicted concentration. If Ci corresponds to a path measurement, however, it represents an average of concentrations along the path. Denote the transformed end points of the straight-line path in the horizontal plane as $\left({x}_{i,\mathrm{1}},{y}_{i,\mathrm{1}}\right)$ and $\left({x}_{i,\mathrm{2}},{y}_{i,\mathrm{2}}\right)$, respectively. The line between the given points in the horizontal plane can be parameterised by ${\mathbit{\rho }}_{i}\left(s\right)=\left({\mathit{\rho }}_{x,i}\left(s\right),{\mathit{\rho }}_{y,i}\left(s\right){\right)}^{\prime },$ where

$\begin{array}{r}{\mathit{\rho }}_{x,i}\left(s\right)=s{x}_{i,\mathrm{2}}+\left(\mathrm{1}-s\right){x}_{i,\mathrm{1}},\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}s\in \left[\mathrm{0},\mathrm{1}\right],\\ {\mathit{\rho }}_{y,i}\left(s\right)=s{y}_{i,\mathrm{2}}+\left(\mathrm{1}-s\right){y}_{i,\mathrm{1}},\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}s\in \left[\mathrm{0},\mathrm{1}\right],\end{array}$

so that

$\begin{array}{}\text{(2)}& \begin{array}{rl}{C}_{i}=& \frac{\mathrm{1}}{{T}_{i}}\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}C\left({\mathit{\rho }}_{x,i}\left(s\right),{\mathit{\rho }}_{y,i}\left(s\right),{z}_{i},Q,{U}_{i},H,{\mathbit{\theta }}_{{k}_{i}}\right)\\ & ∥\frac{\text{d}{\mathbit{\rho }}_{i}\left(s\right)}{\text{d}s}∥\phantom{\rule{0.125em}{0ex}}\text{d}s,\end{array}\end{array}$

where ${T}_{i}\in {\mathbb{R}}^{+}$ is the path length and $\parallel \cdot \parallel$ is the standard Euclidean norm. In our case, ${T}_{i}=∥\frac{\text{d}{\mathbit{\rho }}_{i}\left(s\right)}{\text{d}s}∥$ is not a function of s, and so Eq. (2) simplifies to

${C}_{i}=\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}C\left({\mathit{\rho }}_{x,i}\left(s\right),{\mathit{\rho }}_{y,i}\left(s\right),{z}_{i},Q,{U}_{i},H,{\mathbit{\theta }}_{{k}_{i}}\right)\phantom{\rule{0.125em}{0ex}}\text{d}s.$

This integral can be approximated numerically over a fine partitioning of J segments $P=\mathit{\left\{}\left[{s}_{\mathrm{0}},{s}_{\mathrm{1}}\right],\left[{s}_{\mathrm{1}},{s}_{\mathrm{2}}\right],\mathrm{\dots },\left[{s}_{J-\mathrm{1}},{s}_{J}\right]\mathit{\right\}}$, where $\mathrm{0}={s}_{\mathrm{0}}<{s}_{\mathrm{1}}<\mathrm{\dots }<{s}_{J-\mathrm{1}}<{s}_{J}=\mathrm{1}$. Then

${C}_{i}\approx \sum _{j=\mathrm{1}}^{J}C\left({\mathit{\rho }}_{x,i}\left({s}_{j}^{*}\right),{\mathit{\rho }}_{y,i}\left({s}_{j}^{*}\right),{z}_{i},Q,{U}_{i},H,{\mathbit{\theta }}_{{k}_{i}}\right){\mathrm{\Delta }}_{s},$

where ${\mathrm{\Delta }}_{s}\equiv {s}_{j}-{s}_{j-\mathrm{1}}=\mathrm{1}/J$ and ${s}_{j}^{*}=\frac{{s}_{j}+{s}_{j-\mathrm{1}}}{\mathrm{2}}$. In our experiments we set J=100.

4 Bayesian atmospheric tomography

We are ultimately interested in obtaining a range of plausible values for the emission rate, Q, a posteriori, (i.e. after we have observed some data). In this section we present a hierarchical statistical model that relates Q to the observed concentrations via the Gaussian plume model. Although Q itself is univariate, the model contains several other unknown parameters that capture our uncertainty about the physical and the measurement processes; inferences on these parameters and Q are made simultaneously. For ease of exposition we adopt the terminology of to describe the model, which we also summarise graphically in Fig. 3. The top layer in the hierarchy is the data model (the model for the observations, Y, Sect. 4.1), the middle layer is the process model (the model for Q, Sect. 4.2), and the bottom layer is the parameter model (the unknown parameters not of direct interest, τ, ωy,ωz, Sect. 4.3). In Sect. 4.4 we outline the MCMC strategy we use to make inference with the model.

Figure 3Directed acyclic graph showing the conditional dependence relationships between the data (enhancements) Y and the error components ε (Sect. 4.1); the emission rate Q (Sect. 4.2); and the unknown parameters $\mathbit{\tau },{\mathit{\omega }}_{y},$ and ωz (Sect. 4.3).

## 4.1 The data model

Let $\stackrel{\mathrm{̃}}{\mathbit{Y}}\equiv \left({\stackrel{\mathrm{̃}}{Y}}_{\mathrm{1}},{\stackrel{\mathrm{̃}}{Y}}_{\mathrm{2}},\mathrm{\dots },{\stackrel{\mathrm{̃}}{Y}}_{N}{\right)}^{\prime }$ denote the measured concentrations averaged over 5 min intervals. We model each of these averaged measurements as ${\stackrel{\mathrm{̃}}{Y}}_{i}={C}_{i}+{X}_{i}+{\mathit{\epsilon }}_{i},$ where Ci is the ith Gaussian-plume-predicted concentration, Xi is the sum of the ith CH4 background concentration and instrument-specific bias, and εi denotes the random error associated with the ith observed CH4 concentration. The background concentration and bias can be explicitly modelled and predicted . Here, as in , we estimate Xi as the 5th percentile of all the measurements from the instrument associated with the ith measurement. Figure 4 compares the raw averaged concentrations to those corrected for background and instrument-specific bias, which we term enhancements, when plotted against wind direction (in degrees east of north).

Figure 4(a) Raw averaged concentrations, plotted by instrument and against wind direction. (b) Enhancements obtained by subtracting off the background and instrument-specific bias.

Now, let $\mathbit{Y}\equiv \left({Y}_{\mathrm{1}},{Y}_{\mathrm{2}},\mathrm{\dots },{Y}_{N}{\right)}^{\prime }$ denote the enhancements, and $\mathbit{\epsilon }\equiv \left({\mathit{\epsilon }}_{\mathrm{1}},{\mathit{\epsilon }}_{\mathrm{2}},\mathrm{\dots },{\mathit{\epsilon }}_{N}\right)$. It is straightforward to verify that

$\begin{array}{ll}{Y}_{i}& ={\stackrel{\mathrm{̃}}{Y}}_{i}-{X}_{i}\\ & ={C}_{i}+{\mathit{\epsilon }}_{i},\phantom{\rule{1em}{0ex}}i=\mathrm{1},\mathrm{\dots },N.\end{array}$

Therefore, Yi is made up of two main components of variability: the Gaussian-plume-predicted concentration and a random error term. We assume that the εi terms are Gaussian and independent but that they are not identically distributed. Specifically, εi contains two components of variation, one pertaining to the error characteristics of the instrument and one to the stability class with which we have categorised the measurement. Recall also from Sect. 3.2 that we model the variance of the predicted concentrations to be proportional to a quartic of the true mean inverse wind speed for Ui<1 m s−1.

First, we capture instrument-specific measurement error characteristics and stability-condition-specific variation by introducing an auxiliary variable mi (${m}_{i}=\mathrm{1},\mathrm{2},\mathrm{\dots },M$), where M is the total number of unique combinations of stability class and instrument type, and consider M different precision (i.e. inverse variance) parameters $\mathit{\left\{}{\mathit{\tau }}_{{m}_{i}}\mathit{\right\}}$ that need to be estimated, one for each combination. Second, we take the influence of low wind speeds into account by assuming that the precision of εi is ${\mathit{\tau }}_{{m}_{i}}$ multiplied by ${\stackrel{\mathrm{^}}{U}}_{i}$, where, for Ui>0,

$\begin{array}{}\text{(3)}& {\stackrel{\mathrm{^}}{U}}_{i}=\left\{\begin{array}{ll}{U}_{i}^{\mathrm{4}}& \mathrm{0}<{U}_{i}<\mathrm{1},\\ \mathrm{1}& {U}_{i}\ge \mathrm{1},\end{array}\right\\end{array}$

which encapsulates our prior belief that observed model–measurement mismatch variability at low wind speeds (in this case under 1 m s−1) is dominated by the low wind speed.

Putting these two components together, we have that, conditional on the instrument type and stability class encoded in mi,

${\mathit{\epsilon }}_{i}\mid {m}_{i}\sim \text{Gau}\left(\mathrm{0},\mathrm{1}/\left({\stackrel{\mathrm{^}}{U}}_{i}{\mathit{\tau }}_{{m}_{i}}\right)\right),\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}i=\mathrm{1},\mathrm{\dots },N.$

We detail the prior distribution for ${\mathit{\tau }}_{{m}_{i}}$ in Sect. 4.3.1.

## 4.2 The process model

The process of interest in this application is the emission rate, Q, which we assume is constant. Since in this application Q≥0, we model it using a half-normal prior distribution (a Gaussian distribution with mean zero truncated from below at zero),

$\begin{array}{}\text{(4)}& p\left(Q\right)=\left\{\begin{array}{ll}\frac{\sqrt{\mathrm{2}}}{{\mathit{\sigma }}_{Q}\sqrt{\mathit{\pi }}}\mathrm{exp}\left(-\frac{{Q}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{Q}^{\mathrm{2}}}\right),& Q\phantom{\rule{0.33em}{0ex}}\in \left[\mathrm{0},\mathrm{\infty }\right)\\ \mathrm{0}& \text{otherwise},\end{array}\right\\end{array}$

with a standard deviation parameter, σQ, which is known and fixed. In our case we fixed σQ to 1.5 g s−1 (90 g min−1), which results in a relatively uninformative prior distribution.

While addressing nonnegativity, half-normal priors do not contain a point mass at zero and thus do not encode a prior belief that there is a possibility of having exactly a zero emission rate. As a consequence, a posterior estimate or even a credible interval that includes zero is not possible. A spike-and-slab distribution consisting of a diffuse uniform distribution with a point mass at zero could be alternatively used at the cost of a slightly more complex model.

## 4.3 The parameter model

Our parameter model is divided into two parts: one pertaining to the precision parameters $\mathit{\left\{}{\mathit{\tau }}_{{m}_{i}}\mathit{\right\}}$ in the random-error component in the data model; and the other to the standard deviations in the Gaussian-plume dispersion models which, as shown in Sect. 3.1, are also uncertain.

### 4.3.1 The precision parameters

For conjugacy with the Gaussian likelihood, we model each ${\mathit{\tau }}_{{m}_{i}}$ using a Gamma prior distribution, with shape parameter α and rate parameter β:

$p\left({\mathit{\tau }}_{{m}_{i}}\right)=\frac{{\mathit{\beta }}^{\mathit{\alpha }}}{\mathrm{\Gamma }\left(\mathit{\alpha }\right)}{\mathit{\tau }}_{{m}_{i}}^{\mathit{\alpha }-\mathrm{1}}{e}^{-\mathit{\beta }{\mathit{\tau }}_{{m}_{i}}},\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}i=\mathrm{1},\mathrm{\dots },N.$

In our application we set α=1.058 and β=0.621. These values were chosen through quantile matching, such that the 1st and 99th percentiles of the distribution of $\mathrm{1}/\sqrt{{\mathit{\tau }}_{{m}_{i}}}$ are approximately 0.35 and 6.5 ppm, respectively (giving a mode close to 0.7 ppm). Values for these percentiles were selected based on prior exploratory data analysis of the measurements that were taken upwind of the source.

### 4.3.2 The Gaussian plume model parameters

From separate studies into the reliability of the model values for ${\mathit{\sigma }}_{{y}_{i},{k}_{i}}$ and ${\mathit{\sigma }}_{{z}_{i},{k}_{i}}$, briefly discussed in Sect. 3.1, we concluded that these parameters could indeed be off by factors of 2 or more and that, if they are off, they are so by similar amounts for each stability class. These factors correspond to vertical shifts of the Pasquill stability curves when plotted on a log–log scale (e.g. Wark et al.1998, chap. 4). We thus replaced ${\mathit{\sigma }}_{{y}_{i},{k}_{i}}$ and ${\mathit{\sigma }}_{{z}_{i},{k}_{i}}$ in Eq. (1) with ${\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{{y}_{i},{k}_{i}}$ and ${\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{{z}_{i},{k}_{i}}$, respectively, where

$\begin{array}{rl}& {\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{{y}_{i},{k}_{i}}\equiv {\mathit{\omega }}_{y}{\mathit{\sigma }}_{{y}_{i},{k}_{i}}\\ & \text{and}\\ & {\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{{z}_{i},{k}_{i}}\equiv {\mathit{\omega }}_{z}{\mathit{\sigma }}_{{z}_{i},{k}_{i}},\end{array}$

and ${\mathit{\omega }}_{y},{\mathit{\omega }}_{z}\in {\mathbb{R}}^{+}$ are scaling parameters for ${\mathit{\sigma }}_{{y}_{i},{k}_{i}}$ and ${\mathit{\sigma }}_{{z}_{i},{k}_{i}},$ respectively .

We use Gamma prior distributions for ωy and ωz. In our application we set the shape parameters equal to 1.6084 and the rate parameters equal to 0.7361. These parameters give approximate 1st and 99th percentiles of 0.1 and 8, respectively, and a mode close to 1 (representative of no scalar influence on ${\mathit{\sigma }}_{{y}_{i},{k}_{i}}$ or ${\mathit{\sigma }}_{{z}_{i},{k}_{i}}$). This reflects our prior belief that the standard deviations could be up to an order of magnitude off from those derived using classical Pasquill stability-class theory.

## 4.4 Bayesian inference

Recall $\mathbit{Y}\equiv \left({Y}_{\mathrm{1}},{Y}_{\mathrm{2}},\mathrm{\dots },{Y}_{N}{\right)}^{\prime }$ are the N observed enhancements, and let $\mathbit{U}\equiv \left({U}_{\mathrm{1}},{U}_{\mathrm{2}},\mathrm{\dots },{U}_{N}{\right)}^{\prime }$ and $\mathbf{\Theta }\equiv \left({\mathbit{\theta }}_{{k}_{\mathrm{1}}},{\mathbit{\theta }}_{{k}_{\mathrm{2}}},\mathrm{\dots },{\mathbit{\theta }}_{{k}_{N}}{\right)}^{\prime }$. Further, let $\mathbit{\tau }\equiv \left({\mathit{\tau }}_{\mathrm{1}},{\mathit{\tau }}_{\mathrm{2}},\mathrm{\dots },{\mathit{\tau }}_{M}{\right)}^{\prime }$ be the M parameters associated with each combination of instrument type and stability class. The posterior distribution of the emission rate Q is then given by

$\begin{array}{rl}& p\left(Q\mid \mathbit{Y},\mathbit{U},H,\mathbf{\Theta }\right)\\ \propto & \underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{{\mathbb{R}}^{M+}}{\int }p\left(\mathbit{Y},\mathbit{\tau },{\mathit{\omega }}_{y},{\mathit{\omega }}_{z}\mid Q,\mathbit{U},H,\mathbf{\Theta }\right)p\left(Q\right)\phantom{\rule{0.33em}{0ex}}\text{d}\mathbit{\tau }\phantom{\rule{0.125em}{0ex}}\text{d}{\mathit{\omega }}_{y}\phantom{\rule{0.125em}{0ex}}\text{d}{\mathit{\omega }}_{z}\\ =& p\left(Q\right)\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{{\mathbb{R}}^{M+}}{\int }p\left(\mathbit{Y}\mid Q,\mathbit{\tau },{\mathit{\omega }}_{y},{\mathit{\omega }}_{z},\mathbit{U},H,\mathbf{\Theta }\right)\\ & p\left(\mathbit{\tau }\right)p\left({\mathit{\omega }}_{y}\right)p\left({\mathit{\omega }}_{z}\right)\phantom{\rule{0.125em}{0ex}}\text{d}\mathbit{\tau }\phantom{\rule{0.125em}{0ex}}\text{d}{\mathit{\omega }}_{y}\phantom{\rule{0.125em}{0ex}}\text{d}{\mathit{\omega }}_{z},\end{array}$

where p(Q) is given by Eq. (4) and $p\left(\mathbit{Y}\mid Q,\mathbit{\tau },{\mathit{\omega }}_{y},{\mathit{\omega }}_{z},\mathbit{U},H,\mathbf{\Theta }\right)$ is the likelihood, which is Gaussian.

Computation of the posterior distribution $p\left(Q\mid \mathbit{Y},\mathbit{U},H,\mathbf{\Theta }\right)$ involves a high-dimensional integral that is analytically intractable. We therefore use MCMC, specifically a Gibbs sampler, to obtain samples from the posterior distributions of Q, τ, ωy, and ωz (see Gelman et al.2013, for a comprehensive introduction to MCMC). The Gibbs sampler samples each parameter one at a time from their respective full conditional distributions, where conditioning is done using the most recent samples of all other parameters.

In the case of τ, use of Gamma prior distributions leads to full conditional distributions that are also Gamma. Hence, sampling τ is straightforward. However, the prior distributions on the other parameters are not conjugate priors, and hence the full conditional distributions for each of these are not available in closed form. We therefore use standard Metropolis-within-Gibbs to sample from these conditional distributions, with Gaussian proposals and adaptive scaling during the early stages of the MCMC algorithm. Specifically, for each parameter, the standard deviation of the proposal was increased or decreased as appropriate whenever the acceptance rate fell below 10 % or exceeded 80 %.

5 Results and discussion

## 5.1 Observing system simulation experiment

In this section we discuss results from applying our model to simulated data in an observing system simulation experiment (OSSE). To mimic the conditions in the real experiment, we simulated enhancements using the actual Boreal and EC instrument locations, meteorological observations from the Ginninderra data, and realistic variances for the random-error components. We considered the two release-rate periods separately, using a 6 g min−1 emission rate in the first and a 12 g min−1 emission rate in the second. As in the real experiment, the first Boreal laser/reflector setup (seven paths) was used in the first release-rate period, while the second setup (six paths) was used in the second release-rate period; the EC tower locations were kept constant for both periods. We set the precisions ${\mathit{\tau }}_{{m}_{i}}=\mathrm{4}$ for ${m}_{i}=\mathrm{1},\mathrm{\dots },M$ and the scaling factors ${\mathit{\omega }}_{y}={\mathit{\omega }}_{z}=\mathrm{2}$ to assess the algorithm's ability to calibrate the plume online. Following data simulation, we used MCMC to generate 60 000 samples, left out 20 000 of these as burn-in, and used a thinning factor of 10. Adaptation of the Metropolis samplers was only done during burn-in. Convergence was assessed through visual inspection of the MCMC trace plots.

We made inference on Q, as well as all other parameters in the model, for the Boreal- and EC-simulated data and the two emission rate settings. Table 2 shows the posterior median emission rates, the 95 % posterior credible intervals for the emission rate, and the intervals for the plume standard deviation scaling parameters ωy and ωz. In all cases, we see that the true (simulated) emission rate is captured within our posterior credible intervals and that the median estimates are very close to the true values. Interestingly, we see that while the plume-scaling coefficients have been accurately recovered in most cases, the posterior uncertainty over ωy for the Boreal lasers is very wide. This suggests that ωy might not be identifiable for path measurements, possibly because the averaging effect of the line integral renders the measured concentration insensitive to a specific plume width in the horizontal direction.

Table 2Posterior median emission rates in grams per minute (g min−1), and the posterior 95 % credible intervals of the emission rate in grams per minute, ωy, and ωz from the OSSE. Results shown are from simulated data corresponding to the Boreal lasers (B) and EC towers (E) when the emission rate is 6 g min−1 (E1 and B1) and 12 g min−1 (E2 and B2).

Figure 5(a) Posterior empirical distributions of the emission rate Q in grams per minute (g min−1), for the Boreal lasers (B), FTIR spectrometers (F), EC towers (E), Picarro analysers (P), and the ensemble of all instruments (BFEP), for each release-rate period (1 and 2) during the Ginninderra experiment. The 5.8 g min−1 release-rate period is shown in red (B1, F1, E1, P1, and BFEP1), while the 5.0 g min−1 release-rate period is shown in blue (B2, F2, E2, P2, and BFEP2). The vertical dashed lines denote the respective true emission rates, the black dots denote the median estimates, and the black vertical bars denote the upper and lower limits of the 95 % posterior credible intervals. (b) Same as (a) but showing results obtained using measurements taken when the methane point source was inactive. In both cases, we can recover a reasonable range of estimates for the emission rate, with no 95 % posterior credible interval being far from the true emission rate. Further, we see that the posterior emission rate credible intervals move towards zero when the source is inactive, as desired.

## 5.2 Application to the Ginninderra data set

In this section we discuss results from applying our model to enhancements from the compiled Ginninderra data. We considered several settings. In the first setting, we estimated the emission rate separately for each of the four instrument types and for each release-rate period (5.8 and 5.0 g min−1) when the source was active. In addition, for each release-rate period we estimated the emission rate for all the instruments combined, yielding a total of 10 inversion results. In the second setting we estimated the emission rate for the same 10 cases but for periods when the source was switched off. In the third setting we again considered the same 10 cases but using only measurements that were taken when upwind of the source. These three settings serve to demonstrate how our inferences adapt to the various settings one might encounter in the field. In particular, online plume calibration is almost impossible in the latter two settings, and we expect this to result in large posterior uncertainties on the scaling coefficients, and also the emission rate in the third setting. In the second setting downwind measurements are present. Therefore, while online plume calibration is again almost impossible since there is no active source, the absence of a source (Q=0 g min−1) should be reflected in our posterior inferences (recall, however, that use of a half-normal prior distribution precludes the possibility of a zero emission rate being estimated; see Sect. 4.2).

As in the OSSE, we generated 60 000 MCMC samples, left out 20 000 of these as burn-in, and used a thinning factor of 10. In line with what we observed in the OSSE, our initial results showed that, more often than not, ωy is not identifiable (leading to wide posterior distributions and poor MCMC mixing) when attempting to estimate the emission rate with the source switched on with path measurements. We therefore chose to fix ωy=1 (but not ωz) for path measurements, and this choice is reflected in all the results discussed below.

The left panel of Fig. 5 summarises our results for Q in the first setting (both upwind and downwind measurements with the source switched on); full results are given in the first 10 rows of Table A1. While our posterior inferences are reflective of the true underlying emission rate, unlike in the OSSE we see that with the real data the true values were not always captured within our 95 % posterior credible intervals. This suggests that there are other important factors at play (e.g. with the meteorological data such as ambient temperature or wind direction, which we assume are fixed and known) that are not (or not fully) accounted for in our model. A close inspection of the residuals at EC.A revealed mild deviations from our Gaussianity assumption, while posterior predictive distributions on left-out EC tower data in a reanalysis revealed coverage probabilities (specifically, empirical probabilities computed from the quantity of validation data falling into the 68 % and 95 % prediction intervals, respectively) that are slightly too large. Nevertheless, our worst-case scenario, obtained with the combination of all instruments in the 5.0 g min−1 release-rate period, had an interval limit which was only 0.55 g min−1 (approximately 11 %) off from the true value, while all posterior medians were within 36 % of the true value (within 22 % if one ignored results from the Picarro analysers during the 5.0 g min−1 release-rate period). This is encouraging because a single, common inference method was used to obtain the inferences from data at a common temporal resolution – no manual instrument-specific tuning was carried out. The approach thus seems relatively robust to instrument type; in Sect. 6 we show this is no longer the case once certain components in our model are assumed fixed and known.

The first 10 rows of Table A1 also show the 95 % posterior credible intervals for ωy and ωz. None of the obtained credible intervals for ωy contain 1, and the results corroborate the conclusion from our exploratory data analysis in Sect. 3.1 that a plausible value for ωy is about 2 or 3. This result lends credence to our ability to calibrate the Pasquill stability-class curves corresponding to ${\mathit{\sigma }}_{{y}_{i},{k}_{i}}$ while estimating the emission rate with point measurements. There was less agreement on ωz in the inversions, suggesting that something more complex than a simple scaling is required (or that the model used for ${\mathit{\sigma }}_{{z}_{i},{k}_{i}}$ is, in this case, inappropriate) for calibrating the Pasquill stability-class curves corresponding to ${\mathit{\sigma }}_{{z}_{i},{k}_{i}}$. Nonetheless, in Sect. 6 we show that our emission-rate estimates from point measurements were relatively less sensitive to the assumption ωz=1 than to the assumption ωy=1.

The right panel in Fig. 5 summarises our results for Q in the second setting (both upwind and downwind measurements with the source switched off), while full results are given in the second set of 10 rows in Table A1. Recall from Sect. 4.2 that, due to the choice of prior over Q (a half-normal distribution), it is not possible for the 95 % credible interval to include zero. Clearly, however, the intervals for Q are close to zero and are suggestive of a small emission rate. As expected, the plume standard deviation scaling parameters are not well-constrained in this setting when the source is off: narrow credible intervals on the emission rate here are only possible when the measurement is largely insensitive to the plume shape. This is indeed the case for the Boreal paths, some of which pass very close to the source. With other instrument configurations, uncertainty in the plume scalings dominates. In some cases (FTIR spectrometers and Picarro analysers in the 5.0 g min−1 release-rate period) our MCMC algorithm did not converge after the 60 000 samples; these results are thus omitted from Fig. 5 and Table A1.

The bottom 10 rows in Table A1 give full results in the third setting (upwind measurements only with the source switched on). In this setting the 95 % posterior credible intervals produced for the emission rates are very wide (most with a range of over 100 g min−1), as are those produced for ωy and ωz: our posterior distributions are largely uninformative. This was expected since upwind measurements contain no information on both the emission rate and the plume model parameters. These results from upwind measurements serve as verification and confirm that we are indeed relying on useful information from downwind measurements when making inference on the emission rate and other parameters that appear within our model.

6 Sensitivity of results to model components

As detailed throughout Sect. 4, the Bayesian model we employ contains many parameters that are updated using MCMC. A natural question to ask is whether all these parameters do need to be updated and what the effects on the emission rate inferences are when instead some of these are assumed fixed and known. Specifically, we are interested in seeing what happens when (i) considering only one single precision parameter τ for all of the data regardless of stability class and/or instrument group, (ii) considering one ${\mathit{\tau }}_{{m}_{i}}$ per instrument group only, (iii) not accounting for plume-model variability in low wind speeds (i.e. setting $\stackrel{\mathrm{^}}{\mathbit{U}}=\mathbf{1}\right)$, (iv) not updating ωy when using point measurements, (v) not updating ωz, and (vi) not updating both ωy and ωz when using point measurements. The 95 % credible intervals for Q in grams per minute for all these settings and for each of the 10 groupings considered in Sect. 5 are given in Table A2.

Grouping the precision parameters $\mathit{\left\{}{\mathit{\tau }}_{{m}_{i}}\mathit{\right\}}$ by instrument only (instead of by instrument and stability class) had a slightly negative impact on the emission-rate estimates obtained during the second release-rate period but less so during the first release-rate period. Assuming (and fixing) ωz=1 for both the point and path measurements also did not have a serious impact on the emission-rate estimates. Note that this does not mean that these components are not relevant in the general model – for example, from our estimates of ωz in Table A1 we see ωz=1 would be a plausible choice for this experiment if one opted to fix ωz (while ωy=1 would not be).

On the other hand several components in our model appear to be crucial to obtaining reasonable emission-rate estimates. Using a single precision parameter to capture all observed variability due to measurement error and the stability-class categorisation clearly had a negative impact on our emission-rate estimates. Similarly, assuming the variability of the measurements is independent of wind speed when performing inversion resulted in 95 % posterior credible intervals on the emission rate that are considerably shifted in the negative direction. A similar observation was made by Feitz et al. (2018, p. 207) when analysing data from the Boreal lasers. There, observations with wind speeds below 1.5 m s−1 were removed to mitigate this effect.

The scaling factor ωy is clearly also crucial for obtaining emission-rate estimates of practical significance for point measurements, with the ensuing emission-rate estimates often being off by nearly a factor of 2 when ωy=1 is assumed. As expected, the width of the credible intervals on the emission rate decreased substantially when ${\mathit{\omega }}_{y}={\mathit{\omega }}_{z}=\mathrm{1}$ was assumed, indicating that ωy and ωz play a big role in quantifying uncertainty on the emission rate. Therefore, as noted in other studies discussed in Sect. 1, incorporating uncertainty in the transport model by treating parameters within the model itself as uncertain (note that this is different from adding another component of variability in the data model, as is often done) is likely to have a positive impact on emission-rate estimates and uncertainty quantification.

7 Conclusions

In this article we have proposed a fully Bayesian model for atmospheric tomography that takes into account uncertainty in the data measurement process, the physical processes, and parameters appearing in the transport model, when estimating the emission rate. We see that the model is robust to different instrument types and configurations, and it provides useful inferences on the emission rate and the plume dispersion model used. When applied to the Ginninderra data using a variety of instruments in different release-rate periods, we obtain 95 % posterior credible intervals on the emission rate that either encapsulate the true emission rate or have a limit which is no more than 11 % from the true value.

The methods developed in this study are ideal for quantifying local-scale leaks from industrial facilities or from the subsurface (e.g. well heads, buried pipelines, or gas leakage up geological fractures and faults) where a surface leak has been detected but needs to be quantified. It can be used where physical access to the source location is limited, e.g. gas bubbling from a creek or where measurement is hazardous. Depending on the circumstance, detection of leakage can take many different forms, from visible bubble detection, optical gas imaging, handheld sniffers, noise detection, helicopters equipped with lasers, drones equipped with gas sensors, to monitoring die-off in vegetation using remote sensing techniques. Surface leakage typically expresses as small, concentrated hotspots if sourced from the subsurface , for which the quantification approach outlined in this article is ideally suited. Equipment placement can be optimised around the leakage site (i.e. prevailing upwind/downwind) for optimal quantification.

In most applications neither the number of sources nor the source location is known. As such, the framework we construct should be seen as a foundational building block that needs to be extended appropriately for each specific application. For example, if the source location is not known, then source localisation can be incorporated into the Bayesian framework as discussed by . If there are multiple possible sites, and these locations are not known, then the framework needs to be further extended to incorporate multiple Gaussian plume models (one for each site), and joint localisation–inversion will be required. While these extensions are straightforward both mathematically and computationally, in practice they are unlikely to be effective for detection of leakage over large spatial scales. Gas fields or geological storage sites can cover areas of tens to hundreds of square kilometres. Unless there is a high density of sensors (≈100 m scale, van Leeuwen et al.2013; Jenkins et al.2016), the sensitivity of detection will be poor . It is however relatively straightforward to effectively extend the methodology to when the emission is from an area rather than a point source.

Our work is closely connected to other atmospheric tomography techniques but with some small, significant, differences. used a backward Lagrangian particle model to simulate the trajectories of methane and carbon dioxide backwards in time to localise the source and estimate the emission rates. Their approach yielded good quality estimates for the methane emission rates but highly uncertain estimates for the carbon dioxide emission rates and source location parameters. Twenty-three runs of the Lagrangian model required approximately 1 h of computing time, and therefore their framework becomes problematic with thousands of observations as we have in our study. More pertinently, online calibration of the atmospheric transport model would be virtually impossible without the construction and use of a surrogate model or emulator (e.g. Harvey et al.2018). In the study of , carbon dioxide and nitrous oxide emission rates and source locations were estimated relatively well. We do not consider the localisation problem but otherwise extend their method to handle various instrument types and a number of extra levels of uncertainty. The case in our sensitivity analysis in which we fix ${\mathit{\omega }}_{y}={\mathit{\omega }}_{z}=\mathrm{1}$ yields a model that is structurally very similar to that of ; we see from our results that having this hard constraint is not a tenable assumption in practice. Our work also has close connections with that of where the Pasquill stability class for an observation is chosen from a subset of appropriate stability classes, based on the best fit of model-predicted values to observed values. While this may help fit the Gaussian plume dispersion model to the data, it does not take into account the uncertainty arising from stability-class choice. Further, if all plume model standard deviations are off by a factor of 2 or more, there is a distinct possibility that no stability class yields a good fit. Online calibration of these standard deviations is needed to account for lack-of-fit arising from the inherently simple Gaussian plume model.

Our results provide interesting insights into the design and monitoring of sensor networks for detecting and quantifying methane emissions. For example, our sensitivity analysis in Sect. 6 showed that estimates using the two Picarro analysers were particularly sensitive to assumptions made on the model plume parameters. Moreover, when uncertainty on these parameters was considered, the release-rate estimates from these instruments tended to be uncertain. This is despite the Picarro analysers being among the more accurate and expensive instruments used in the study. Uncertainty in our experiment is, as is often the case, dominated by that in the transport model. Hence, the number of instruments used, the proximity of the instruments to the source, and their configuration around the source appear to be more important design criteria than instrument accuracy when the inferential target is emission-rate quantification of a point source. In particular, having more (less expensive) instruments set up to cover many more possible wind directions is better than having only one or two more expensive instruments with which to monitor emissions. If one is limited to using a small number of instruments, then those giving path measurements are preferable to those giving point measurements, as the former will be able to capture a larger range of wind directions. Our results also provide insight on the transport model used. For example, close inspection of our posterior inferences for τ indicated that, across all instrument groups and for both release-rate periods, the model–data mismatch was much lower for the more neutral stability classes C and D than for the more stable/unstable classes A and F.

The fully Bayesian framework we adopt is adaptable to various scenarios. We envision, for example, that source localisation (e.g. Humphries et al.2012; Hirst et al.2013) could be done in tandem with plume-model calibration within an inversion framework, provided several instruments in suitable configurations (as in the Ginninderra experiment) are available. Future work will also investigate how uncertainty in other meteorological variables such as wind-direction, as well as the stability-class categorisation adopted (possibly via z0), could be incorporated within the model.

Code and data availability
Code and data availability.

Software code and data are available at https://github.com/Lcartwright94/BayesianAT (last access: 1 August 2019 Cartwright2019).

Appendix A: Full results

Table A1Posterior median emission rate in grams per minute (g min−1), and the posterior 95 % credible intervals for the emission rate in grams per minute, ωy, and ωz, for the Boreal lasers (B), FTIR spectrometers (F), EC towers (E), Picarro analysers (P), and an ensemble of all instruments (BFEP), for each release-rate period (5.8 g min−1 (1), and 5.0 g min−1 (2)) under various settings. Dashes correspond to parameters that were not updated via MCMC. Results for which MCMC did not converge are marked as n/a.

Table A2Posterior 95 % credible intervals for the emission rates in grams per minute (g min−1) for the Boreal lasers (B), FTIR spectrometers (F), EC towers (E), Picarro analysers (P), and an ensemble of all instruments (BFEP), for each release-rate period (5.8 g min−1 (1), and 5.0 g min−1 (2)) and for various alterations to the model as detailed in Sect. 6. Dashes correspond to the redundant case (e.g. ωy=1 was assumed for all path measurements in the full model).

Author contributions
Author contributions.

LC compiled the data with the help of all authors and ran all the analyses. LC and AZM conducted the research. AF conceptualised and supervised the study. LC, AZM, and AF wrote the manuscript. SB conducted an initial investigation using a simplified version of the proposed model. IS, FP, TC, KN, TN, MK, SZ, NW, and NMD acquired the field data for the study. All authors discussed the results and commented on the manuscript.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Special issue statement
Special issue statement.

This article is part of the special issue “The 10th International Carbon Dioxide Conference (ICDC10) and the 19th WMO/IAEA Meeting on Carbon Dioxide, other Greenhouse Gases and Related Measurement Techniques (GGMT-2017) (AMT/ACP/BG/CP/ESD inter-journal SI)”. It is a result of the 10th International Carbon Dioxide Conference, Interlaken, Switzerland, 21–25 August 2017.

Acknowledgements
Acknowledgements.

Laura Cartwright acknowledges the support of the Australian Government Research Training Program Scholarship. Laura Cartwright, Andrew Zammit-Mangion, and Andrew Feit would like to acknowledge APR.Intern for facilitating the first 5 months of this modelling study. All authors thank Gareth Davies for reviewing an earlier version of the manuscript. The Ginninderra field site was supported by the Australian Government through the Carbon Capture and Storage – Implementation budget measure. The authors also acknowledge funding for the research provided by the Australian Government through the CRC programme and support from the CO2CRC. The National Geosequestration Laboratory is thanked for making the two Picarro instruments available for the study. We would like to thank Phil Dunbar and his staff (CSIRO Plant Industry) for maintaining the site and Dale Hughes (CSIRO) for his assistance with maintenance of the CSIRO EC tower. The authors also wish to acknowledge the assistance of Field Engineering Services at Geoscience Australia. Geoscience Australia and the Western Sydney University team would like to acknowledge Charles Jenkins (CSIRO) for early discussions about the atmospheric tomography line technique and the Australian Mathematical Sciences Institute. The University of Wollongong wishes to acknowledge Joel Wilson, Maximilien Desservettaz, and Ruhi Humphries for their assistance in the site operations. Andrew Feitz and Ivan Schroder publish with the permission of the CEO of Geoscience Australia.

Financial support
Financial support.

This research has been supported by the Australian Research Council (grant nos. DE180100203 and FT180100327).

Review statement
Review statement.

This paper was edited by Hubertus Fischer and reviewed by two anonymous referees.

References

Ars, S., Broquet, G., Yver Kwok, C., Roustan, Y., Wu, L., Arzoumanian, E., and Bousquet, P.: Statistical atmospheric inversion of local gas emissions by coupling the tracer release technique and local-scale transport modelling: a test case with controlled methane emissions, Atmos. Meas. Tech., 10, 5017–5037, https://doi.org/10.5194/amt-10-5017-2017, 2017. a, b

Basu, S., Baker, D. F., Chevallier, F., Patra, P. K., Liu, J., and Miller, J. B.: The impact of transport model differences on CO2 surface flux estimates from OCO-2 retrievals of column average CO2, Atmos. Chem. Phys., 18, 7189–7215, https://doi.org/10.5194/acp-18-7189-2018, 2018. a

Berliner, L. M.: Hierarchical Bayesian time series models, in: Maximum Entropy and Bayesian Methods, edited by: Hanson, K. M. and Silver, R. N., Springer, New York, NY, 15–22, 1996. a

Borysiewicz, M., Wawrzynczak, A., and Kopka, P.: Stochastic algorithm for estimation of the model's unknown parameters via Bayesian inference, Proceedings of the Federated Conference on Computer Science and Information Systems, Wroclaw, Poland, 501–508, 2012. a, b

Casella, G. and Berger, R. L.: Statistical Inference, 2nd edn., Duxbury Press, Pacific Grove, CA, 2002. a

Cartwright, L.: Bayesian atmospheric tomography with application to data from the 2015 Ginninderra release experiment, available at: https://github.com/Lcartwright94/BayesianAT, last access: 1 August 2019. a

Chevallier, F., Feng, L., Bösch, H., I. Palmer, P., and Rayner, P.: On the impact of transport model errors for the estimation of CO2 surface fluxes from GOSAT observations, Geophys. Res. Lett., 37, L21803, https://doi.org/10.1029/2010GL044652, 2010. a

Etheridge, D., Luhar, A., Loh, Z., Leunning, R., Spencer, D. Steele, P., Zegelin, S., Allison, C., Krummel, P., Leist, M., and van der Schoot, M.: Atmospheric monitoring of the CO2CRC Otway Project and lessons for large scale CO2 storage projects, Energy Proceedia, 4, 3666–3675, 2011. a

Feitz, A., Leamon, G., Jenkins, C., Jones, D. G., Moreira, A., Bressan, L., Melo, C., Dobeck, L. M., Repasky, K., and Spangler, L. H.: Looking for leakage or monitoring for public assurance?, Energy Proceedia, 63, 3881–3890, 2014. a

Feitz, A., Schroder, I., Phillips, F., Coates, T., Negandhi, K., Day, S., Luhar, A., Bhatia, S., Edwards, G., Hrabar, S., Hernandez, E., Wood, B., Naylor, T., Kennedy, M., Hamilton, M., Hatch, M., Malos, J., Kochanek, M., Reid, P., Wilson, J., Deutscher, N., Zegelin, S., Vincent, R., White, S., Ong, C., George, S., Maas, P., Towner, S., and Griffith, D.: The Ginninderra CH4 and CO2 release experiment: an evaluation of gas detection and quantification techniques, Int. J. Greenh. Gas Con., 70, 202–224, 2018. a, b, c, d, e, f, g, h

Flesch, T. K., Wilson, J. D., Harper, L. A., Crenna, B. P., and Sharpe, R. R.: Deducing ground-to-air emissions from observed trace gas concentrations: A field trial, J. Appl. Meteorol., 43, 487–502, 2004. a, b

Forde, O. N., Mayer, K. U., and Hunkeler, D.: Identification, spatial extent and distribution of fugitive gas migration on the well pad scale, Sci. Total Environ., 652, 356–366, 2019. a

Ganesan, A. L., Rigby, M., Zammit-Mangion, A., Manning, A. J., Prinn, R. G., Fraser, P. J., Harth, C. M., Kim, K.-R., Krummel, P. B., Li, S., Mühle, J., O'Doherty, S. J., Park, S., Salameh, P. K., Steele, L. P., and Weiss, R. F.: Characterization of uncertainties in atmospheric trace gas inversions using hierarchical Bayesian methods, Atmos. Chem. Phys., 14, 3855–3864, https://doi.org/10.5194/acp-14-3855-2014, 2014. a

Ganesan, A. L., Manning, A. J., Grant, A., Young, D., Oram, D. E., Sturges, W. T., Moncrieff, J. B., and O'Doherty, S.: Quantifying methane and nitrous oxide emissions from the UK and Ireland using a national-scale monitoring network, Atmos. Chem. Phys., 15, 6393–6406, https://doi.org/10.5194/acp-15-6393-2015, 2015. a

Gelman, A., Stern, H. S., Carlin, J. B., Dunson, D. B., Vehtari, A., and Rubin, D. B.: Bayesian Data Analysis, 3rd edn., Chapman & Hall/CRC Press, Boca Raton, FL, 2013. a

Golder, D.: Relations among stability parameters in the surface layer, Bound.-Lay. Meteorol., 3, 47–58, 1972. a

Harvey, N. J., Huntley, N., Dacre, H. F., Goldstein, M., Thomson, D., and Webster, H.: Multi-level emulation of a volcanic ash transport and dispersion model to quantify sensitivity to uncertain parameters, Nat. Hazards Earth Syst. Sci., 18, 41–63, https://doi.org/10.5194/nhess-18-41-2018, 2018. a

Hirst, B., Jonathan, P., del Cueto, F. G., Randell, D., and Kosut, O.: Locating and quantifying gas emission sources using remotely obtained concentration data, Atmos. Environ., 74, 141–158, 2013. a, b

Houweling, S., Bergamaschi, P., Chevallier, F., Heimann, M., Kaminski, T., Krol, M., Michalak, A. M., and Patra, P.: Global inverse modeling of CH4 sources and sinks: an overview of methods, Atmos. Chem. Phys., 17, 235–256, https://doi.org/10.5194/acp-17-235-2017, 2017. a

Humphries, R., Jenkins, C., Leuning, R., Zegelin, S., Griffith, D., Caldow, C., Berko, H., and Feitz, A.: Atmospheric tomography: a Bayesian inversion technique for determining the rate and location of fugitive emissions, Environ. Sci. Technol., 46, 1739–1746, 2012. a, b, c, d, e, f

International Energy Agency: Energy Technology Perspectives 2017, OECD/IEA, Paris, France, 2017. a, b

Jacobson, M. Z.: Fundamentals of Atmospheric Modeling, 2nd edn., Cambridge University Press, New York, NY, 2005. a

Jenkins, C., Kuske, T., and Zegelin, S.: Simple and effective atmospheric monitoring for CO2 leakage, Int. J. Greenh. Gas Con., 46, 158–174, 2016. a, b

Jenkins, J. D., Luke, M., and Thernstrom, S.: Getting to zero carbon emissions in the electric power sector, Joule, 2, 2487–2510, 2018. a

Jones, M., Goldstein, M., Jonathan, P., and Randell, D.: Bayes linear analysis for Bayesian optimal experimental design, J. Stat. Plan. Infer., 171, 115–129, 2016. a

Kinnon, M. A. M., Brouwer, J., and Samuelsen, S.: The role of natural gas and its infrastructure in mitigating greenhouse gas emissions, improving regional air quality, and renewable resource integration, Prog. Energ. Combust., 64, 62–92, 2018. a

Lewicki, J. L. and Hilley, G. E.: Eddy covariance mapping and quantification of surface CO2 leakage fluxes, Geophys. Res. Lett., 36, L21802, https://doi.org/10.1029/2009GL040775, 2009. a

Loh, Z. M., Leuning, R., Zegelin, S. J., Etheridge, D. M., Bai, M., Naylor, T., and Griffith, D.: Testing Lagrangian atmospheric dispersion modelling to monitor CO2 and CH4 leakage from geosequestration, Atmos. Environ., 43, 2602–2611, 2009. a

Lucas, D. D., Simpson, M., Cameron-Smith, P., and Baskett, R. L.: Bayesian inverse modeling of the atmospheric transport and emissions of a controlled tracer release from a nuclear power plant, Atmos. Chem. Phys., 17, 13521–13543, https://doi.org/10.5194/acp-17-13521-2017, 2017. a

Luhar, A. K., Etheridge, D. M., Leuning, R., Loh, Z. M., Jenkins, C. R., and Yee, E.: Locating and quantifying greenhouse gas emissions at a geological CO2 storage site using atmospheric modeling and measurements, J. Geophys. Res.-Atmos., 119, 10959–10979, 2014. a, b, c, d

Miller, S. M., Hayek, M. N., Andrews, A. E., Fung, I., and Liu, J.: Biases in atmospheric CO2 estimates from correlated meteorology modeling errors, Atmos. Chem. Phys., 15, 2903–2914, https://doi.org/10.5194/acp-15-2903-2015, 2015. a

Mitchell, T. J. and Beauchamp, J. J.: Bayesian variable selection in linear regression, J. Am. Stat. Assoc., 83, 1023–1032, 1998. a

Pasquill, F.: The estimation of the dispersion of wind-borne material, Meteorol. Mag., 90, 33–49, 1961. a

Peylin, P., Baker, D., Sarmiento, J., Ciais, P., and Bousquet, P.: Influence of transport uncertainty on annual mean and seasonal inversions of atmospheric CO2 data, J. Geophys. Res.-Atmos., 107, 4385, https://doi.org/10.1029/2001JD000857, 2002. a

Rajaona, H., Septier, F., Armand, P., Delignon, Y., Olry, C., Albergel, A., and Moussafir, J.: An adaptive Bayesian inference algorithm to estimate the parameters of a hazardous atmospheric release, Atmos. Environ., 122, 748–762, 2015. a

Riddick, S. N., Connors, S., Robinson, A. D., Manning, A. J., Jones, P. S. D., Lowry, D., Nisbet, E., Skelton, R. L., Allen, G., Pitt, J., and Harris, N. R. P.: Estimating the size of a methane emission point source at different scales: from local to landscape, Atmos. Chem. Phys., 17, 7839–7851, https://doi.org/10.5194/acp-17-7839-2017, 2017. a

Sepulveda, N. A., Jenkins, J. D., de Sisternes, F. J., and Lester, R. K.: The role of firm low-carbon electricity resources in deep decarbonisation of power generation, Joule, 2, 2403–2420, 2018. a

Sienfeld, J. H. and Pandis, S. N.: Atmospheric Chemistry and Physics: From Air Pollution to Climate Change, 2nd edn., John Wiley & Sons, Hoboken, NJ, 2006. a, b

Tarantola, A.: Inverse Problem Theory and Methods for Model Parameter Estimation, SIAM, Philadelphia, PA, https://doi.org/10.1137/1.9780898717921, 2005. a

Turner, B.: Workbook of Atmospheric Dispersion Estimates, 2nd edn., Lewis Publishers, Boca Raton, FL, 1994. a

van Leeuwen, C., Hensen, A., and Meijer, H. A. J.: Leak detection of CO2 pipelines with simple atmospheric CO2 sensors for carbon capture and storage, Int. J. Greenh. Gas Con., 19, 420–431, 2013. a, b

Wang, Y., Huang, H., Huang, L., and Ristic, B.: Evaluation of Bayesian source estimation methods: A comparison of likelihood functions and distance measures, Atmos. Environ., 152, 519–530, 2017. a

Wark, K., Warner, C. F., and Davis, W. T.: Air Pollution: Its Origin and Control, Addison Wesley Longman, Menlo Park, CA, 1998. a, b, c, d

White, E. D., Rigby, M., Lunt, M. F., Smallman, T. L., Comyn-Platt, E., Manning, A. J., Ganesan, A. L., O'Doherty, S., Stavert, A. R., Stanley, K., Williams, M., Levy, P., Ramonet, M., Forster, G. L., Manning, A. C., and Palmer, P. I.: Quantifying the UK's carbon dioxide flux: an atmospheric inverse modelling approach using a regional measurement network, Atmos. Chem. Phys., 19, 4345–4365, https://doi.org/10.5194/acp-19-4345-2019, 2019. a

Wilson, P., Feitz, A., Jenkins, C., Berko, H., Loh, Z., Luhar, A., Hibberd, M., Spencer, D., and Etheridge, D.: Sensitivity of CO2 leak detection using a single atmospheric station, Energy Proceedia, 63, 3907–3914, 2014. a

World Meteorological Organisation: Guide to Meteorological Instruments and Methods of Observation, available at: https://library.wmo.int/pmb_ged/wmo_8_en-2012.pdf (last access: 27 March 2019), 2008.  a

Zammit-Mangion, A., Cressie, N., Ganesan, A. L., O'Doherty, S., and Manning, A. J.: Spatio-temporal bivariate statistical models for atmospheric trace-gas inversion, Chemometr. Intell. Lab,. 15, 227–241, 2015. a