# Simultaneous quantification and changepoint detection of point source gas
emissions using recursive Bayesian inference

Amir Montazeri
Sibley School of Mechanical and Aerospace Engineering
Cornell University
Ithaca, NY, USA
<EMAIL_ADDRESS>

& Xiaochi Zhou
School of Civil and Environmental Engineering
Cornell University
Ithaca, NY, USA
<EMAIL_ADDRESS>
&John D. Albertson
School of Civil and Environmental Engineering
Cornell University
Ithaca, NY, USA
<EMAIL_ADDRESS>

Corresponding author.

###### Abstract

Recent findings suggest that abnormal operating conditions of equipment in the
oil and gas supply chain represent a large fraction of anthropogenic methane
emissions. Thus, effective mitigation of emissions necessitates rapid
identification and repair of sources caused by faulty equipment. In addition
to advances in sensing technology that allow for more frequent surveillance,
prompt and cost-effective identification of sources requires computational
frameworks that provide automatic fault detection. Here, we present a
changepoint detection algorithm based on a recursive Bayesian scheme that
allows for simultaneous emission rate estimation and fault detection. The
proposed algorithm is tested on a series of near-field controlled release
mobile experiments, with promising results demonstrating successful detection
(>90% success rate) of changes in the leak rate when the emission rate is
tripled after an abrupt change. Moreover, we show that the statistics of the
measurements, such as the coefficient of variation and range are good
predictors of the performance of the algorithm. Finally, we describe how this
methodology can be easily adapted to suit time-averaged concentration data
measured by stationary sensors, thus showcasing its flexibility.

_K_ eywords changepoint detection $\cdot$ Bayesian inference $\cdot$ Methane
$\cdot$ environmental sensing

##  1 Introduction

Methane is a potent greenhouse gas (GHG) with a global warming potential (GWP)
that is approximately 84 and 28 times greater than carbon dioxide (CO2) on 20
and 100 year time scales, respectively [1]. Methane emissions from the oil and
gas industry are among the largest anthropogenic sources of methane in the
United States, accounting for approximately 30% of total emissions in 2019
[2]. Recent studies have found that emissions from almost all subsectors of
the oil and gas supply chain demonstrate a “fat-tail” distribution, such that
a relatively small number of large emitters are responsible a large fraction
of the emissions [3, 4, 5, 6, 7, 8, 9]. While the presence of these large
emitters is concerning, it offers the potential of expedient reduction in GHG
emissions and costs if the largest emitters are rapidly identified and
repaired.

Measurements of methane emissions are often classified as either top-down or
bottom-up [10]. Top-down studies rely on ambient methane measurements using
aircraft, satellites, or tower networks to estimate aggregate emissions from
all contributing sources across large geographies [10]. On the other hand,
bottom-up methods aggregate and extrapolate emissions from individual pieces
of equipment, operations, or facilities, using measurements made directly at
the emission point or, in the case of facilities, directly downwind [10, 11,
12, 13, 14]. Recent integrated research efforts have found that while emission
estimates from facility-based bottom-up approaches and top-down approaches are
in agreement, these estimates are significantly higher than component-based
estimates (i.e., when emissions are extrapolated from individual pieces of
equipment) [3-4]. Detailed investigation of the discrepancy between component-
based aggregates and estimates from other approaches (i.e., facility-based
bottom-up and top-down methods) suggests that component-based methods miss
high emissions caused by abnormal operating conditions (e.g., malfunctions).
Such abnormal conditions are the defining attribute of large emitters that
contribute the majority of emissions in the oil and gas sector. Therefore,
prompt identification of abnormal process conditions (i.e., fault detection)
can lead to substantial reductions in emission and costs for operators.

The abnormal operating conditions observed in the largest emitters are
spatially and temporally variable [15, 16, 17, 4]. For example, emissions can
significantly increase for a production site due to malfunctions at a certain
point in time. Hence, emission reduction requires monitoring approaches that
enable efficient and timely responses to the appearance of abnormal process
conditions. Continuous monitoring offers the capability to rapidly detect
faulty behavior that is necessary for reducing emissions from the largest
emitters. Moreover, recent efforts to develop innovative technologies and
algorithms to mitigate methane emissions have also highlighted the advantages
of continuous monitoring over “snapshot-in-time” approaches in rapid
identification of large emission sources [18, 19].

While continuous monitoring of oil and gas facilities are not commonplace yet,
changes in the near future are likely due to the following: 1) Advancements in
sensor technology and wireless communications allow for continuous
measurements to be made and stored in the cloud [18], and 2) monitoring
mandates at the state and federal level, such as the U.S. Environmental
Protection Agency (USEPA) fence-line monitoring program for early detection of
benzene emissions [20] and continuous monitoring of air quality during pre-
production and early-production of drilling operations producing gas and
liquid hydrocarbons required by the Colorado Department of Public Health and
environment [21]. Meanwhile, innovations in atmospheric inversion modeling are
required to enable automatic emission estimation and fault detection using
continuous measurements. It is worth noting that efficient and rapid fault
detection can serve as an incentive for operators to employ continuous
surveillance systems, as it can lead to significant reductions in unwanted
emissions and costs.

Here, we propose a recursive Bayesian inference model that utilizes
measurements from continuous surveillance systems (e.g., network of fixed
sensors, and mobile sensors on unmanned vehicles) for simultaneous estimation
and changepoint (fault) detection in point-source emission rates. Note that we
use changepoint or fault to refer to a sudden increase in point-source
emission rate as expected under abnormal operating conditions described
earlier. The Bayesian inference model is useful in this context for multiple
reasons: 1) It is equipped to deal with noisy data which in this case are
caused by the stochastic nature of turbulence that drives the emitted gas and
other measurement uncertainties, 2) it permits the determination of the
uncertainty in estimated emission rates [22] and 3) it allows for “online”
changepoint detection and emission estimates, i.e., the emission rate
estimates and probability of detecting changes in emission rates are updated
with every new measurement that arrives incrementally [23].

In this chapter, we first introduce an instantaneous plume dispersion
formulation that is used to guide our Bayesian analysis. Then, the
mathematical framework of the Bayesian inference and its application in point-
source estimation and changepoint detection is described. Next, we apply the
proposed Bayesian framework to near-field (with source-to-sensor distances
$\leq$ 30m) mobile measurements of a controlled point-source emission and
evaluate its performance in changepoint detection. Although the Bayesian
framework in this chapter is tailored towards mobile measurements, it can
easily be adapted for continuous measurements made by networks of fixed
sensors.

##  2 Theory

Point-source characterization approaches often rely on time-averaged
measurements. However, our analysis of the mobile sensor data requires a
formulation for the instantaneous plume that is adapted from [24] and
described in section 2.1. The presented formulation is applicable to passive
scalars, which are diffusive contaminants in low concentrations such that they
have no dynamical effect on the motion of the surrounding flowing fluid [25].
In subsequent sections, the following assumptions are made: 1) The emission
rate from the point source is constant until an abrupt change causes the leak
rate to increase to a new constant value. 2) The emission rates before and
after the changepoint, and therefore mass concentrations before and after the
changepoint are independent.

###  2.1 Instantaneous view of plume transport

Consider a steady-state point source located at the origin $O$ of a local
coordinate system (Figure 1). The wind velocity components in $x,y,$ and $z$
directions are defined as $u,v,$ and $w$, respectively. We define a control
volume starting at the origin $O$ to a downwind vertical plane located at
$x_{m}$, extending from $y_{min}$ to $y_{max}$ laterally, and from $z_{min}$
to $z_{max}$ vertically, such that the control volume encompasses the entire
plume upwind of the mobile sensor. For this control volume, the application of
conservation of mass yields the following expression for the emission rate
(mass per time), $Q_{0}$:

$Q_{0}=F(x_{m},t)+\frac{dS(t)}{dt},$ (1)

where $S(t)$ is total mass of the emitted gas in the control volume, $t$ is
time, and $F(x_{m},t)$ is the mass flow rate exiting the control volume
through the vertical plane at $x_{m}$. The control volume is defined such that
no mass exits anywhere other than the vertical plane at $x_{m}$, therefore the
mass flow rate can be expressed as

$F(x_{m},t)=\int_{z_{min}}^{z_{max}}\int_{y_{min}}^{y_{max}}c(x_{m},y,z,t)u(x_{m},y,z,t)dydz,$
(2)

where $c$ is the mass concentration of the passive scalar. It is worth noting
that in the atmospheric boundary layer, large Reynolds numbers are typically
observed and the flow is highly turbulent, therefore, molecular diffusion is
ignored [26].

Figure 1: A control volume containing a point emission source (located at the
origin, $O$) with a mass flow rate of $Q_{0}$, and a cross-plane view of the
mass flow rate at downwind distance, $x_{m}$.

It is useful to define a normalized distribution of the mass concentration,
labeled $D$, and a plume-weighted advection velocity, labeled $u_{e}$, at the
exit plane of the control volume as

$D(x_{m},y,z,t)=\frac{c(x_{m},y,z,t)}{\int_{z_{min}}^{z_{max}}\int_{y_{min}}^{y_{max}}c(x_{m},y,z,t)dydz}$
(3)
$u_{e}(x_{m},t)=\int_{z_{min}}^{z_{max}}\int_{y_{min}}^{y_{max}}D(x_{m},y,z,t)u(x_{m},y,z,t)dydz,$
(4)

We can rewrite $u_{e}(x_{m},t)$ by substituting equation (3) into (4) and
applying equation (2) as:

$\displaystyle u_{e}(x_{m},t)$
$\displaystyle=\frac{\int_{z_{min}}^{z_{max}}\int_{y_{min}}^{y_{max}}c(x_{m},y,z,t)u(x_{m},y,z,t)dydz}{\int_{z_{min}}^{z_{max}}\int_{y_{min}}^{y_{max}}c(x_{m},y,z,t)dydz}$
$\displaystyle=\frac{F(x_{m},t)}{\int_{z_{min}}^{z_{max}}\int_{y_{min}}^{y_{max}}c(x_{m},y,z,t)dydz}.$
(5)

Equation (2.1) can then be used in conjunction with equation (1) to relate the
mass concentration trajectory when traversing the plume, $c(x_{m},y,z,t)$, to
the other relevant variables as

$c(x_{m},y,z,t)=\frac{Q_{0}-dS(t)/dt}{u_{e}(x_{m},t)}D(x_{m},y,z,t).$ (6)

In practice, $u_{e}(x_{m},t)$ can be approximated using nearby meteorological
measurements. The vertical scaling of the wind profile based on the Monin-
Obukhov similarity theory (MOST) can then be applied to adjust these
meteorological measurements by height difference as needed [27], following the
methodology previously described in [24]. Accordingly, $u_{e}(x_{m},t)$ is
replaced by $u_{e}^{M}(x_{m},t)\delta_{u}(x_{m},t)$, where
$\delta_{u}(x_{m},t)$ accounts for the ratio between the actual and
approximated plume-weighted advection velocities, $u_{e}(x_{m},t)$ and
$u_{e}^{M}(x_{m},t)$, respectively. The superscript $M$ is used to highlight
model estimated quantities. Similarly, we introduce
$\delta_{S}(x_{m},t)=1-\frac{dS(t)/dt}{Q_{0}}$ to represent the non-steadiness
in the total mass stored in the control volume, normalized by $Q_{0}$.
Therefore, equation 6 can be expressed as follows

$c(x_{m},y,z,t)=\frac{Q_{0}}{u_{e}^{M}(x_{m},t)}\left(\frac{\delta_{S}(x_{m},t)}{\delta_{u}(x_{m},t)}\right)D(x_{m},y,z,t),$
(7)

which describes an instantaneous view of plume transport while having the same
underlying form as commonly used models based on an ensemble-averaged view
[28]. The key differences between this instantaneous view and common ensemble-
averaged models are the time dependence of the distribution $D$ that
represents the stochastic nature of the turbulent plume, and the presence of
$\delta_{u}(x_{m},t)$ and $\delta_{S}(x_{m},t)$ that accounts for non-
stationarity in wind speed and mass storage in the control volume.

$D(x_{m},y,z,t)$ is a random variable that captures the plume movement in time
and the lateral and vertical directions as it responds to the instantaneous
turbulent velocity components in these directions. It is therefore expected
that $D(x_{m},y,z,t)$ scales with the standard deviations of the velocity
components in the $y$ and $z$ directions ($\sigma_{v}$ and $\sigma_{w}$,
respectively). It is well-understood that local scaling approaches based on
MOST often provide an acceptable description of $\sigma_{w}$ [29]. However,
$\sigma_{v}$ is affected by random large scale motions in the atmosphere that
cannot be described accurately by local scaling laws [30]. Therefore, a
greater degree of randomness is expected in $y$ than in $z$ directions for
$D(x_{m},y,z,t)$.

This observation highlights the benefits of integrating both sides of equation
(7) with respect to $y$ (i.e. across the plume), since the uncertainty
associated with the lateral plume dispersion can be effectively removed [31]:

$c^{y}(x_{m},z,t)=\frac{Q_{0}}{u_{e}^{M}(x_{m},t)}\left(\frac{\delta_{S}(x_{m},t)}{\delta_{u}(x_{m},t)}\right)D_{z}(x_{m},z,t)$
(8)

where $c^{y}(x_{m},z,t)=\int_{y_{min}}^{y_{max}}c(x_{m},y,z,t)dy$ is the
cross-plume integrated concentration, and the reflection of the vertical
profile of mass conservation at $x_{m}$ is described by
$D_{z}=\int_{y_{min}}^{y_{max}}D(x_{m},y,z,t)dy$. Note that $D_{z}(x_{m},z,t)$
is a random variable that is mainly driven by the stochastic nature of the
vertical transport dynamics in the turbulent flow. In field applications, the
sensor path is typically constrained by adjacent roadways, which are not
always perpendicular to the wind direction. When the road segments are at a
significant angle to the wind direction, $c^{y}(x_{m},z,t)$ can be estimated
by numerical integration of the mass concentration along the path as [31]:

$c^{y}(x_{m},z,t)=\sum_{y_{min}}^{y_{max}}c(x_{m},y,z,t)\Delta
tV\sin(\theta_{r})$ (9)

where $\Delta t$ is the sensor acquisition time step, $V$ is the vehicle
velocity and $theta_{r}$ is the acute angle between the road segment and the
wind direction.

We can account for all the stochasticity in $c^{y}(x_{m},z,t)$ by introducing
a fluctuating variable. To this end, we define
$D_{z,e}=\frac{\delta_{S}(x_{m},t)}{\delta_{u}(x_{m},t)}D_{z}(x_{m},z,t)$ to
include the stochastic nature of $\delta_{S}(x_{m},t),\delta_{u}(x_{m},t)$ and
$D_{z}(x_{m},z,t)$. This new fluctuating variable is helpful in empirical
analysis of the cross-plume integrated mass concentration and can be used to
rewrite equation (8) as follows:

$c^{y}(x_{m},z,t)=\frac{Q_{0}}{u_{e}^{M}(x_{m},t)}D_{z,e}(x_{m},z,t).$ (10)

Equation (10) can be used in a forward manner to estimate the downwind cross-
plume integrated mass concentration $c^{y}(x_{m},z,t)$ for a given emission
rate $Q_{0}$, as well as in the inverse problem of inferring $Q_{0}$ given
downwind measurements of $c^{y}(x_{m},z,t)$. Most dispersion models only offer
an approximation of the ensemble-averaged $D_{z,e}(x_{m},z,t)$, therefore, we
apply Bayesian inference to account for the fluctuation of the instantaneous
$D_{z,e}(x_{m},z,t)$ from its ensemble averaged as detailed in section 2.2.
For simplicity of notation the independent variables, $x_{m},z,$ and $t$ will
be dropped hereafter.

###  2.2 Bayesian inference for source estimation

Following Bayes’ rule, and the notation introduced by Arumpalam et al. for
recursive Bayesian inference [32], the posterior probability distribution of
the emission rate $Q$ based on the measurements of $c^{y}$ at time step $k$
(or after the $k$’th sensor pass) is [31], [22]

$p\left(Q_{k}\big{|}c^{y}_{1:k}\right)=\frac{p\left(Q_{k}\big{|}c^{y}_{1:k-1}\right)p\left(c^{y}_{k}\big{|}Q_{k}\right)}{p\left(c^{y}_{1:k}\right)},$
(11)

where $p\left(Q_{k}\big{|}c^{y}_{1:k}\right)$,
$p\left(Q_{k}\big{|}c^{y}_{1:k-1}\right)$, and $p\left(c^{y}_{1:k}\right)$ are
probability density functions (PDFs).
$p\left(Q_{k}\big{|}c^{y}_{1:k-1}\right)$ is the prior that is being updated
through the recursion, $p\left(c^{y}_{k}\big{|}Q_{k}\right)$ is the likelihood
function, and $p\left(c^{y}_{1:k}\right)$ is the evidence term that ensures
$p\left(Q_{k}\big{|}c^{y}_{1:k}\right)$ integrates to unity. Note that the
notation $c^{y}_{a:b}$ refers to the contiguous set of measurements between
time (sensor pass) $a$ and $b$ inclusive.

In practical applications, past measurements of similar facilities [33] may be
used to formulate the prior probability distribution at the first time step,
i.e. before any sampling activities. Under the assumption that the only prior
knowledge of $Q$ is its lower and upper bounds, a uniform prior distribution
can be adopted [22, 34]. This uniform distribution is often considered as
sufficiently uninformative based on the principle of maximum entropy [35] and
can be expressed as follows

$p(Q_{1})=\frac{1}{Q_{max}-Q_{min}},$ (12)

where $Q_{max}$ and $Q_{min}$ are the prescribed upper and lower bounds of the
emission rate, respectively.

The likelihood function, $p\left(c^{y}_{k}\big{|}Q_{k}\right)$, describes the
probability of observing $c^{y}$ given $Q$ at the $k$’th sensor pass and
encodes all the information provided by the mass concentration measurements
about the unknown emission rate [22]. Since the underlying distribution of the
concentration measurements are unknown, the principle of maximum entropy
supports the application of a Gaussian distribution with a prescribed error
scale [36]. This choice for the likelihood function has proven useful in
previous studies [37, 34, 38, 31], and is therefore adopted here. Furthermore,
the dataset investigated in this study has been previously tested for leak
estimation using the Gaussian likelihood function with satisfactory results
[24] (more details on the dataset are provided in section 3) . Consequently,
The Gaussian likelihood function in this study is expressed as

$p\left(c^{y}_{k}\big{|}Q_{k}\right)=\frac{1}{\sigma_{e}\sqrt{2\pi}}\exp\left[-\frac{1}{2}\left(\frac{c^{y}_{k}-c^{y,M}_{k}(Q_{k})}{\sigma_{e}}\right)^{2}\right],$
(13)

where $c^{y,M}(Q)=\frac{Q}{u_{e}^{M}}D_{z}^{M}$ is the cross-plume integral of
a modeled concentration for a given candidate value of $Q$. $D_{z}^{M}$ is the
estimated value of $D_{z,e}$ based on a Lagrangian Stochastic Model (LSM). The
LSM is used to describe plume dispersion in a turbulent flow by modeling paths
of fluid particles, which are driven by the random velocity field modeled by
the generalized Langevin equation [39]. In this study, we impose the so-called
well-mixed conditions and adopt Thomson’s simplest solution for statistically
stationary and horizontally homogeneous turbulence [40]. The LSM takes
meteorological measurements (friction velocity, surface roughness, standard
deviation of $u$ and $w$ and Obukhov Length) and the estimated distance
between the emission source and the sensor as input parameters. The LSM is
previously described in further detail in [24].

In equation (13), $\sigma_{e}$ is the uncertainty scale parameter, which can
be estimated for observed data as

$\sigma_{e}=\sqrt{\frac{1}{N-1}\sum_{k=1}^{N}\left(c^{y}_{k}-c^{y,M}(Q_{0})\right)^{2}},$
(14)

where $N$ is the number of passes per experiment. In this study, $\sigma_{e}$
is estimated from the controlled release experiments and is known prior to the
application of the Bayesian inference approach. In cases where $\sigma_{e}$
cannot be estimated from prior measurements, it can be estimated using the
error propagation method [41, 42]. Briefly, $\sigma_{e}$ is due to the
following error scale parameters: 1) error due to the stochastic nature of
atmospheric plume dispersion, 2) error due to the plume dispersion model and
3) measurements error including errors from the model input data. The
parameterization of each of these error parameters is dependent upon the local
meteorological conditions, the dispersion model used, and the quality of
measurements.

The recursive Bayesian formulation of equation (11) used in conjunction with
the uniform prior of equation does not have an analytical solution and should
be solved numerically. For the numerical solution, $Q$ is discretized from
$Q_{min}$ to $Q_{max}$ with a uniform step size of $\Delta Q$ to form a vector
of candidate values for $Q$ to be considered. For each measurement of $c^{y}$,
the likelihood function is evaluated at all candidate $Q$ values using
equation (13) and multiplied by the prior probability distribution.
Subsequently, the evidence term after the $k$’th sensor pass can be calculated
through numerical integration as follows:

$p\left(c^{y}_{1:k}\right)=\sum_{Q_{min}}^{Q_{max}}p\left(Q_{k}\big{|}c^{y}_{1:k-1}\right)p\left(c^{y}_{k}\big{|}Q_{k}\right)\Delta
Q.$ (15)

After calculating the evidence term, equation (11) can be applied to evaluate
the posterior distribution of the emission rate at each candidate $Q$ value.
The same procedure is repeated after each mobile sensor pass and the posterior
distribution is updated using a new prior (posterior at previous sensor pass)
and a newly calculated likelihood function with the most recent $c^{y}$
measurement.

After each sensor pass, the posterior PDF can be used to estimate the emission
rate and the associated uncertainty. For instance, the emission rate after
sensor pass $k$, can be calculated as the mean, median or mode of of the
posterior PDF $p\left(Q_{k}\big{|}c^{y}_{1:k}\right)$. Given that an
uninformative, uniform prior distribution was adopted, we expect that the
median and the mean to be heavily affected by the prior in the early stages of
analysis. Therefore, the mode of the posterior PDF
$p\left(Q_{k}\big{|}c^{y}_{1:k}\right)$ is used as the estimated emission
rate, i.e.:

$E_{k}^{Q}=\operatorname*{arg\,max}_{Q_{k}}\left[p\left(Q_{k}\big{|}c^{y}_{1:k}\right)\right].$
(16)

As number of mobile passes are increased, the effects of the prior
distribution are reduced and the mode, median and mean of the posterior PDF
grow closer in value. Furthermore, in cases where an informative prior can be
derived prior from past experiments, the mean or median of the posterior PDF
may be better candidates for the emission rate, as they better incorporate the
prior information than the mode.

Finally, the associated uncertainty of the emission rate estimation using
Bayesian inference is often calculated as the standard deviation,
$\sigma_{k}^{Q}$, of the posterior PDF:

$\left(\sigma_{k}^{Q}\right)^{2}=\int_{Q_{min}}^{Q_{max}}\left(Q-\bar{Q_{k}}\right)^{2}\times
p\left(Q_{k}\big{|}c^{y}_{1:k}\right)dQ,$ (17)

where $Q_{k}$ is the expectation of the posterior PDF evaluated as follows:

$\bar{Q_{k}}=\int_{Q_{min}}^{Q_{max}}Q\times
p\left(Q_{k}\big{|}c^{y}_{1:k}\right)dQ.$ (18)

###  2.3 Bayesian inference for changepoint detection

In order to detect a change in the source emission rate, we apply the Bayesian
Online Changepoint Detection (BOCD) methodology [23]. In this approach,
changepoints are found by first estimating the posterior distribution over the
run length, i.e. the time (or in this case, number of sensor passes) since the
last changepoint, given the data observed so far. Denoting the length of the
current run after sensor pass $k$ with $r_{k}$, and applying the definition of
conditional probability the run length posterior distribution can be expressed
as

$p\left(r_{k}\big{|}c^{y}_{1:k}\right)=\frac{p\left(r_{k},c^{y}_{1:k}\right)}{p(c^{y}_{1:k})},$
(19)

where the run length evidence term is calculated using
$p(c^{y}_{1:k})=\sum_{r_{k}}p\left(r_{k},c^{y}_{1:k}\right)$.

After every sensor pass there are two possibilities regarding the changepoint:
1) No changes occur after the sensor pass and therefore the run length is
increase by 1 or 2) change occurs and run length is reset to 0. The
probability that no changepoint occurs is referred to as the "growth
probability" as it indicates that the run length is growing by 1 compared to
the previous sensor pass. Similarly, we refer to the probability that a
changepoint occurs as the "changepoint probability". We denote the growth
probability such that the run length reaches $i\in\\{0,1,2,...,k-1\\}$ after
$k$ sensor passes by $\alpha_{k}(i)$ and derive a recursive estimate as
follows

$\displaystyle\alpha_{k}(i)$ $\displaystyle=p\left(r_{k}=i,c^{y}_{1:k}\right)$
$\displaystyle=p\left(r_{k}=i,r_{k-1}=i-1,c^{y}_{k},c^{y}_{k-1}\right)$
$\displaystyle=p\left(r_{k}=i,c^{y}_{k}\big{|}r_{k-1}=i-1,c^{y}_{1:k-1}\right)p\left(r_{k-1}=i-1,c^{y}_{1:k-1}\right)$
$\displaystyle=p\left(r_{k}=i,c^{y}_{k}\big{|}r_{k-1}=i-1,c^{y}_{1:k-1}\right)\alpha_{k-1}(i-1).$
(20)

The first step in equation (2.3) is due to the fact that the run length can
only increase by 1 after each sensor pass and the second step follows from the
chain rule in probability theory. Furthermore, the recursive estimate in
equation (2.3) requires the evaluation of
$p\left(r_{k}=i,c^{y}_{k}\big{|}r_{k-1}=i-1,c^{y}_{1:k-1}\right)$ which can be
achieved by employing the chain rule again:

$\displaystyle p\left(r_{k},c^{y}_{k}\big{|}r_{k-1},c^{y}_{1:k-1}\right)$
$\displaystyle=p\left(r_{k}\big{|}r_{k-1},c^{y}_{1:k-1},c^{y}_{k}\right)p\left(c^{y}_{k}\big{|}r_{k-1},c^{y}_{1:k-1}\right)$
$\displaystyle=p\left(r_{k}|r_{k-1}\right)p\left(c^{y}_{k}\big{|}c^{y}_{k-i:k-1}\right),$
(21)

where $r_{k}=i$ and $r_{k-1}=i-1$ are implied and not explicitly written for
simplicity of notation. The final step of equation (2.3) follows from the
independence of the measured mass concentrations before and after a
changepoint, which is true based on our assumption of independence of leak
rates before and after a changepoint. Further, the condition on $r_{k}$ in the
second term on right hand side of the equation is absorbed by only limiting
the conditional probability on measurements since the last changepoint leading
to the subscript $k-i:k-1$. Equation (2.3) suggests that the growth
probability can be computed based on two calculations: 1) The prior over
$r_{k}$ given $r_{k-1}$ (also referred to as the changepoint prior) and 2) The
predictive distribution over the new measurement after sensor pass $k$, given
the data since the last changepoint.

The changepoint prior in equation (2.3) has nonzero mass at only two outcomes,
because after each sensor pass the run length either continues to grow such
that $r_{k}=r_{k-1}+1$ or a changepoint occurs and $r_{k}=0$. Furthermore, the
hazard function can be used to quantify each of these outcomes, since by
definition the hazard function quantifies the probability that a changepoint
occurs at a given time step conditioned that no changepoint has occurred prior
to that time step [23, 43]. Therefore, the changepoint prior can be expressed
as

$p\left(r_{k}|r_{k-1}\right)=\begin{cases}H\left(r_{k-1}+1\right)&\text{if}\;r_{k}=0\\\
1-H\left(r_{k-1}+1\right)&\text{if}\;r_{k}=r_{k-1}+1\\\
0&\text{otherwise},\end{cases}$ (22)

where $H\left(r_{k-1}\right)$ is the hazard function. The hazard function
depends on the discrete a priori probability distribution over the interval
between changepoints. However, in this study, we consider the special case
where the a priori probability distribution is a discrete geometric
distribution with timescale $\lambda$, i.e. the run length distribution is due
to a memoryless process and the hazard function is constant at $1/\lambda$.
The timescale $\lambda$ can be set through prior knowledge, for instance, the
average number of passes completed before a change in the emission rate
occurs. In this study, $\lambda$ is set to 15 based on the conducted
experiments described in section 3.1.

The predictive distribution $p\left(c^{y}_{k}\big{|}c^{y}_{k-i:k-1}\right)$
can be described by an equivalent distribution by utilizing the plume
transport model of equation (10) leading a one-to-one correspondence between
the predictive distribution $p\left(c^{y}_{k}\big{|}c^{y}_{k-i:k-1}\right)$
and $p\left(Q_{k}\big{|}c^{y}_{k-i:k-1}\right)$. In this case, the predictive
distribution $p\left(c^{y}_{k}\big{|}c^{y}_{k-i:k-1}\right)$ can be found
through scaling of $p\left(Q_{k}\big{|}c^{y}_{k-i:k-1}\right)$. We note that
$p\left(Q_{k}\big{|}c^{y}_{k-i:k-1}\right)$ is the prior distribution in
equation (11), given the data since the last changepoint. Therefore,
$p\left(Q_{k}\big{|}c^{y}_{k-i:k-1}\right)$ and consequently
$p\left(c^{y}_{k}\big{|}c^{y}_{k-i:k-1}\right)$ can be estimated using the
recursive Bayesian approach detailed in section 2.2, leading to simultaneous
estimation of leak rate and changepoint detection. The computational details
of calculating the growth probability term are provided in section 2.4.

The changepoint probability is evaluated in a similar manner to the growth
probability. By noting that with the occurrence of a changepoint, the run
length $r_{k}$ drops to 0, we can derive a recursive estimate for the
changepoint probability as

$\displaystyle\alpha_{k}(0)$ $\displaystyle=p\left(r_{k}=0,c^{y}_{1:k}\right)$
$\displaystyle=\sum_{j=1}^{k-2}p\left(r_{k}=0,r_{k-1}=j,c^{y}_{k},c^{y}_{k-1}\right)$
$\displaystyle=\sum_{j=1}^{k-2}p\left(r_{k}=0|r_{k-1}=j\right)p\left(c^{y}_{k}\big{|}c^{y}_{k-j-1:k-1}\right)\alpha_{k-1}(j)$
(23)

where $\alpha_{k}(0)$ is the changepoint probability. The summation after the
first step of equation (2.3) appears due to marginalization over the run
length at sensor pass $k-1$. The intermediate steps in the derivation are
omitted as they are identical to the derivation of the growth probability as
outlined in equations (2.3) and (2.3). The estimation of each term on the
right hand side of equation (2.3) follows the same procedure as equation
(2.3).

Figure 2 illustrates the algorithm used to estimate the posterior distribution
over the run length as described by equations (2.3) and (2.3). In this
diagram, the solid blue lines correspond to growth probability calculations,
while the red dashed lines are associated with changepoint probability
evaluations. We note that multiple dashed lines arriving at a node in Figure
2, correspond to the marginalization over the run length at the previous
sensor pass.

Figure 2: Visual description of the message passing algorithm used to estimate
the run length distribution over the observed data. The circles represent run
length hypotheses and the lines between the circles show recursive transfer of
mass between sensor pass (or time step). Solid lines indicate that probability
mass is being passed upwards, causing the run length to grow at after the next
sensor pass and dashed lines indicate that the current run is truncated, and
the run length drops to zero.

In practice, to automatically detect a changepoint and alert the system that a
change in leak rate has occurred a prescribed condition on the changepoint
probability should be put into place. The appropriate detection condition can
be chosen based on the application and the available ancillary information
regarding the emission conditions. In this study, we use a changepoint
probability threshold such that when the changepoint probability is above this
threshold the system automatically registers a change in leak rate, the
posterior probability over the emission rate at the previous sensor pass is
retained and the prior distribution in equation (11) is reset to the uniform
prior of equation (12).

###  2.4 Computational details of growth probability estimation

The calculation of the growth probability distribution after sensor pass $k$
relies on knowledge of the probability distribution of the emission rate
conditioned on the observed data from the last changepoint until the previous
sensor pass at $k-1$, i.e., $p\left(Q_{k}\big{|}c^{y}_{k-i:k-1}\right)$ as
shown by the derivation in equation (2.3). Furthermore, the growth probability
is calculated for every value of run length $i\in\\{0,1,2,...,k-1\\}$, which
requires the recursive calculation of the posterior probability of emission
rate for every $i$ through equation (11). Therefore, the growth probability
calculation can be computationally expensive for large values of $k$. To
overcome this challenge, we note that the likelihood function in equation (11)
is independent of $i$, hence, the calculation of the posterior probability of
emission rate for every $i$ should completed in a vectorized or Single
Instruction, Multiple Data (SIMD) manner that allows for parallel processing
of data and significantly improves the run time of the probability estimation
[44].

##  3 Materials and Methods

In this section we first describe the conducted field experiments that yielded
the data used to examine the Bayesian framework. Throughout these experiments,
the leak rate was kept constant, therefore, a data synthesis procedure is
implemented to simulate a step change in the leak rate before the application
of the changepoint detection algorithm. A series of performance measures are
then defined to quantify the performance of the changepoint detection
algorithm.

###  3.1 Field experiments

Experiments of controlled releases of methane were conducted at the McGovern
soccer training field of Cornell University (Game Farm Rd, Ithaca, NY. USA) in
early August 2016. During the experiments, the site was covered with short
grass ($\sim$5cm) and located in a relatively open field with a 500m distance
from a residential area in the West, 150m distance from small forest in the
North and approximately 400m (500m) distance from roads on the east (south)
side. Point-source emission of methane (99.9% pure gas) was controlled by a
mass flow controller (SmarTrak 100 from Sierra Instruments Inc., Monterey, CA,
USA), with a mass flow accuracy of $\pm 1\%$. The height of the methane
release was similar to the height of the grass at 5cm. On the west side of the
field, a 3D sonic anemometer (CSAT-3, Campbell Scientific Inc., Logan, UT,
USA) was installed on a small tower to measure local meteorological
conditions. The height of the tower was 2.31m with the sonic anemometer
measuring the three components of wind velocity and air temperature at a
frequency of 10 Hz. To mimic emissions from a source surrounded by other low-
level structures, a 1.4m barrier (windbreak) was established in a circle
around the emission source. This setup can for instance, approximate a well
head located in densely organized well pad or a pipeline within a natural gas
metering station.

A mobile measurement platform (MMP) was configured with a precise GPS unit
(Trimble Geo 7X handheld from Trimble Inc., Sunnyvale, CA, USA) to track its
position at a sampling frequency of 1 Hz. The accuracy of the GPS unit was
approximately 5-15 cm for $>$97% of the measured data points after post
processing. The MMP was equipped with a LI-COR LI-7700 open-path methane
analyzer (LI-COR Biosciences, Lincoln, NE, USA), that outputs methane mixing
ratios in the unit of parts per million (ppm). The operating frequency of the
analyzer was set to 10 Hz, and it was positioned at a height of 1.3m.
Furthermore, the analyzer was calibrated by the manufacturer less than a month
before the experiment, and is designed with an open-path configuration for
long term monitoring without regular re-calibration.

Regarding the measurements, a conversion factor is applied to translate above-
ambient mixing ratios ($c_{a}$, in ppm) into mass concentrations ($c$, in
g/m3) in equation (2). This conversion factor is dependent on the molecular
weight of the released gas (16.04 g/mol for methane) and the ambient
temperature which affects the molar of the gas. The above-ambient mixing
ratios are found by subtracting the ambient methane mixing ratios from the raw
methane mixing ratios measured by the open-path analyzer, with the ambient
mixing ratio calculated as the 5th percentile of the ranked time series of raw
mixing ratio measurements [33, 45, 24, 31]. The estimated ambient mixing ratio
was compared to methane mixing ratios measured prior to the experiments with
minimal differences found ($<$2%), suggesting that the ambient mixing ratio
was determined robustly.

To ensure perpendicular sensor passes with respect to the wind direction,
stake flags were placed in three circles centered at the emission source with
radii of 10, 20 and 30m, and repeated passes were made along each of the
circles. The average sensor speed was very low during the experiments
(approximately 2 m/s) to better capture the plume structure. In addition, the
start of a pass took place approximately one minute after the end of the
previous pass, warranting the independence of the measurements of each pass
from those of previous passes. Data was aggregated within 30 minute intervals,
during which the measurements from the meteorological tower was used to
estimate meteorological parameters such as the Obukhov length ($L$) and the
friction velocity ($u_{*}$). Each set of passes completed in a 30-minute
period is considered and analyzed as a single experiment.

###  3.2 Data synthesis

A total of 18 experiments were conducted with the 1.4m barrier present in the
field, with six experiments at each source-to-sensor distance (i.e. $x_{m}$ of
10, 20 and 30m). Two selection requirements were established to filter out
experiments performed under unacceptable conditions. First, experiments under
stable atmospheric conditions were excluded, and only experiments conducted
under neutral or unstable conditions were retained. Second, experiments
conducted under low wind ($\bar{u}<$ 1.0 m/s) and high turbulent intensity
($I_{u}>$0.5) conditions were discarded. As a result of this selection
criteria, 14 experiments were retained for further analysis, the details of
which are summarized in Table 1 alongside the meteorological conditions
reported by the meteorological tower.

Table 1: Summary of experimental conditions, including experiment identification number (ID), approximate source-to-sensor distance ($x_{m}$), number of sensor passes for the experiment ($N$), and sampling day of year (DOY), and meteorological conditions as measured by the nearby meteorological tower, including mean streamwise velocity ($\bar{u}$), standard deviation of streamwise velocity ($\sigma_{u}$), turbulent intensity ($I_{u}$), friction velocity ($u_{*}$), mean wind direction ($\theta_{m}$) clockwise from north, sensible heat flux ($H$), and atmospheric stability ($z/L$). The meteorological variables are derived from data collected during the corresponding experiments ( 30 min). ID |  $x_{m}$ (m) | $N$ | DOY |  $\bar{u}$ (m/s) |  $\sigma_{u}$ (m/s) |  $I_{u}$ (-) |  $u_{*}$ (m/s) |  $\theta_{m}$ (deg) |  $H$ (W/m2) |  $z/L$ (-)
---|---|---|---|---|---|---|---|---|---|---
1 | 30 | 14 | 217 | 2.94 | 0.98 | 0.28 | 0.18 | 147 | 77.67 | -0.32
2 | 20 | 16 | 217 | 2.72 | 1.06 | 0.31 | 0.23 | 149 | 148.97 | -0.31
3 | 10 | 15 | 217 | 2.49 | 0.98 | 0.33 | 0.29 | 172 | 138.13 | -0.14
4 | 30 | 12 | 217 | 2.72 | 1.15 | 0.31 | 0.24 | 152 | 161.21 | -0.30
5 | 20 | 16 | 217 | 2.95 | 1.08 | 0.29 | 0.22 | 146 | 171.01 | -0.41
6 | 10 | 16 | 217 | 2.48 | 0.97 | 0.33 | 0.21 | 159 | 159.46 | -0.46
7 | 30 | 14 | 218 | 3.41 | 1.34 | 0.28 | 0.37 | 204 | 219.77 | -0.11
8 | 20 | 14 | 218 | 3.70 | 1.27 | 0.26 | 0.36 | 211 | 223.92 | -0.12
9 | 10 | 13 | 218 | 3.82 | 1.33 | 0.26 | 0.38 | 212 | 207.86 | -0.10
10 | 30 | 16 | 218 | 4.18 | 1.31 | 0.24 | 0.36 | 194 | 171.36 | -0.09
11 | 20 | 13 | 218 | 4.31 | 1.31 | 0.24 | 0.37 | 184 | 173.08 | -0.08
12 | 10 | 13 | 218 | 3.99 | 1.23 | 0.25 | 0.40 | 179 | 129.53 | -0.05
13 | 20 | 12 | 219 | 2.75 | 1.04 | 0.30 | 0.16 | 318 | 125.83 | -0.77
14 | 10 | 13 | 219 | 2.29 | 1.04 | 0.35 | 0.20 | 315 | 140.29 | -0.47

In all the conducted experiments, the emission rate was kept constant at a
rate of $Q_{0}$ = 0.083 g/s. Therefore, a data synthesis procedure was
established to artificially simulate a change in the emission rate. The
proposed data synthesis approach relies on the observation that the order of
the measurements of the cross-plume integrated mass concentrations through
each pass of the MMP can affect the performance of the changepoint detection
algorithm. For instance, consider two permutations of the measurements of the
same experiment (experiment ID 14 in Table 1) as shown in Figure 3. In the
first permutation, series of measurements leading to multiple low $c^{y}$
values are followed by a comparatively high measurement (at sensor pass 9),
therefore the algorithm is likely to detect a changepoint while in reality the
emission rate has been constant. However, in the second permutation, the
difference in consecutive $c^{y}$ measurements are smaller than the first
permutation and the probability of detecting a changepoint is lowered.

Figure 3: Two instances (a) and (b) of the same experiment (experiment ID 4)
created through random shuffling of the measurements. The measurements are
cross-plume integrated above-ambient mass concentrations of methane calculated
over sensor passes.

This observation motivates the use of random shuffling in synthesizing
experiments consisting of a change in the leak rate after a few passes of the
MMP around the point source. For each experiment, the data synthesis steps are
as follows (Experiment ID 4 in Table 1 used in Figure 4):

  1. 1.

For each pass, $c^{y}$ is calculated to create a time series for the
experiment, which will be referred to as the "original signal" as depicted in
Figure 4a.

  2. 2.

The original signal is duplicated and scaled by a given constant that is
chosen based on the ratio of the leak rate before and after a simulated
changepoint resulting in a "scaled signal" as presented in Figure 4b.

  3. 3.

Both original and scaled signals are shuffled to create a random permutation
of each respective signal with the results shown in Figure 4c and 4d

  4. 4.

The two shuffled signals are concatenated such that the first measurement of
the shuffled scaled signal follows the last measurement of the shuffled
original signal, creating a signal that consists of a changepoint as
illustrated in Figure 4e.

Figure 4: Summary of the data synthesis procedure for one instance of an
experiment (ID 4), where the (a) original measurements from the field
experiment are (b) scaled by multiplying all measurements by a prescribed
constant. Two sets of new measurements are created by random shuffling leading
to a (c) shuffling of the original measurements and a (d) shuffling of the
scaled measurements. The two shuffled sets are then (e) concatenated to create
one instance of synthesized measurements consisting of a step change.

Steps 2-4 of the above procedure are repeated 1000 times, to create a total of
1000 signals with changepoints for each experiment. The changepoint detection
algorithm is then applied to these signals and the performance of the
algorithm is evaluated according to the performance measures described in
section 3.3.

###  3.3 Performance measures

The performance of changepoint detection methods are often evaluated through a
series of commonly used measures. The importance of each measure is dependent
upon the application of the changepoint detection system. Here, we introduce
four different performance measures, that are slightly altered with respect to
common definitions to better suit the context of changepoint detection in
emission rates. In the following description of the performance measures, each
unique signal consisting of a changepoint is referred to as an "instance" of
an experiment.

#### Recall

This refers to the portion of the changepoints that are detected after exactly
one sensor pass following the change in emission rate. We label these
successful changepoint detections as "True Positive" instances denoted by TP.
Instances where the changepoints are detected with a delay, i.e. where
changepoints are detected after at least two sensor passes after the change in
emission rate are labeled "Delayed True Positive" (DTP), and instances where
changepoints were not detected are referred to as "False Negative" (FN)
instances. Therefore, recall which is a measure of how effective the
changepoint algorithm is in detecting changepoints as soon as they occur is
expressed as

$Recall=\frac{TP}{TP+DTP+FN}.$ (24)

#### Detection Recall

This refers to the portion of changepoints that are detected any time after
the change in emission rate has occurred. Employing the labels introduced
earlier, detection recall as a measure of how effective the changepoint
algorithm is in detecting the changepoints is expressed as

$Detection\;Recall=\frac{TP+DTP}{TP+DTP+FN}.$ (25)

#### Detection Delay

This measures the average number of passes that it takes to detect the
changepoint after the emission rate has changed. This measure is evaluated
only for experiments where the changepoints were detected for all instances
(with or without delay) and is evaluated as

$Detection\;Delay=\frac{\sum_{i=1}^{TP+DTP}Predicted\;CP-Actual\;CP}{TP+DTP},$
(26)

where $Predicted\;CP$ refers to the sensor pass after which the changepoint is
detected and $Actual\;CP$ refers to the sensor pass after which the change in
leak rate has occurred.

#### False Positive Rate

This refers to the ratio of number of instances where changepoints are
detected prior to the change in emission rate to to total number of instances.
Here "False Positive" (FP) refers to instances where data points that are not
changepoints are recognized as changepoints. The False Positive Rate is a
measure that reflects how many false alarms would be generated by the
changepoint detection algorithm and is expressed as follows

$False\;Positive\;Rate=\frac{FP}{All\;instances}=\frac{FP}{1000},$ (27)

noting that $All\;instances$ refers to the 1000 signals with changepoints
created for each experiment.

To account for the variability introduced through random shuffling during the
data synthesis stage, for each experiment the data synthesis procedure is
repeated 100 times and therefore 100 different estimates of each performance
measure are computed. These 100 values of the performance measures are then
collected and used in a bootstrapping significance test analysis to establish
95% confidence intervals for the computed performance measures [46].

##  4 Results and Discussion

Before presenting the results related to the performance of the changepoint
detection algorithms across all experiments, we explore one instance of an
experiment. For this instance, it is shown how the changepoint detection
algorithm is coupled with the emission estimation procedure to approximate the
leak rates before and after a change in the emission rate.

###  4.1 Leak estimation and changepoint detection

We first show the changepoint detection procedure for one instance of an
experiment (Experiment ID 4 in Table 1). In this instance, 12 sensor passes
are made before the leak rate is significantly increased from $Q_{1}$ = 0.083
g/s to $Q_{2}$ = 0.332 g/s which is four times as large as $Q_{1}$. The
$c^{y}$ measurements including the changepoint are shown in Figure 5a. After
each sensor pass, the changepoint probability is evaluated through the
procedure described in section 2.3 with the values presented in Figure 5b. A
changepoint is detected when the changepoint probability surpasses a
prescribed probability threshold of 0.8, after which the prior to the
recursive Bayesian inference of equation (11) is reset to the uniform prior
(equation(12)), so that the new emission rate can be approximated.
Furthermore, the figure shows an increase in the changepoint probability after
the 19th sensor pass which can be attributed to the high value of the $c^{y}$
measurement (in comparison to previous measurements) corresponding to this
sensor pass. In this study, the changepoint probability threshold is chosen
through trial and error to lower the false alarm rate of the detection
algorithm, and its effect on false positive rate is investigated in section
4.2.

Figure 5: Application of the changepoint algorithm to an instance of an
experiment (ID 4). (a) Synthesized instance with a step change in leak rate
from 0.083 g/s to 0.332 g/s, where the vertical dashed line indicates the
first sensor pass after the change. (b) The changepoint probability plotted
after every sensor pass, where the horizontal dashed line represents the
changepoint threshold probability, above which the algorithm registers a
changepoint and resets the recursive Bayesian inference for leak estimation.

For the same experiment instance as above, Figure 6 illustrates the evolution
of the posterior PDF of the emission rate after each sensor pass before and
after the changepoint. In this case, the lower and upper bounds of the
emission rate, denoted by $Q_{min}$ and $Q_{max}$ are specified as 0 and 5.0
g/s. The choice for $Q_{min}$ is trivial as the emission rate can only take
positive values. $Q_{max}$ is determined through trial and error such that the
tail of the derived posterior PDF of the emission rate is close to zero. Using
a larger $Q_{max}$ does not affect the accuracy of the Bayesian inference
procedure, however it is deemed unnecessary as it increases the computational
cost of the recursive Bayesian inference scheme. Figure 6b shows that
posterior PDF is fairly small at $Q$ = 2.0 g/s, suggesting that the choice of
$Q_{max}$ = 5.0 g/s is effective. Starting from a relatively broad posterior
PDF, suggesting a large uncertainty in the emission rate, the posterior PDF
tends to approach a more narrow shape with additional sensor passes. It is
worth noting that after the change in emission rate, the variation in the
$c^{y}$ measurements are much larger compared to measurements at the original
emission rate. Therefore, it is necessary to use a new estimate for
$\sigma_{e}$ in equation (13) for approximating the emission rate after the
changepoint. In practice, the emission rate after the change is not known,
hence, a larger and more conservative choice for $\sigma_{e}$ can be used to
accommodate this lack of information. In the example shown in Figure 6, we
employ $\sigma_{e,2}=10\times\sigma_{e,1}$ where $\sigma_{e,1}$ and
$\sigma_{e,2}$ refer to the error scale parameters before and after the
changepoint, respectively, which is a conservative choice given that the error
scale after the change is four times the error scale prior to the change. This
conservative choice leads to a higher projected uncertainty when it comes to
estimating the emission rate after the changepoint.

Figure 6: The evolution of the posterior probability $p\left(Q|c^{y}\right)$,
of the emission rate $Q$ after each sensor pass (a) before and (b) after the
change in leak rate as detected by the changepoint detection algorithm for one
instance of an experiment (ID 4). The posterior probability obtained after the
final sensor pass before the changepoint in (a) and after the overall final
pass are presented with a solid red line. The vertical dashed lines indicate
the actual emission rate of (a) 0.083 g/s and (b) 0.332 g/s.

###  4.2 Changepoint detection performance

We investigate the performance of the changepoint detection method using the
measures introduced in section 3.3 by systematically varying the magnitude of
the change in leak rate when synthesizing the data. To this end, Figure 7
shows recall for varying values of "jump-to-noise ratio" (JNR), where JNR is
the ratio of the absolute difference in the average $c^{y}$ before and after
the change (i.e., the jump) to the standard deviation of $c^{y}$ before the
change in leak rate (i.e., the noise). As expected, when the change in leak
rate is of the order of the noise in the measurements, or in other words JNR
is of the order of 1, changepoints are difficult to detect and therefore
recall is low for all experiments. Further, as JNR is increased a monotonic
rise in performance is observed across all experiments, with similar recall
values observed in almost all cases. This similarity of recall values across
experiments for each JNR motivates the idea of grouping all experiments based
on the source-to-sensor distance. Figure 8 presents the recall averaged across
all experiments within each group as a function of JNR, where the vertical
bars indicate the 95% confidence intervals. The trends observed in Figure 7
and 8 suggest that JNR can solely predict the recall for the changepoint
algorithm irrespective of the source-to-sensor distance and the measurement
noise.

Figure 7: Evolution of recall for a series of jump to noise ratios varying
between 1.5 and 15.5 for source-to-sensor distances, $x_{m}$ of (a) 10m, (b)
20m and (c) 30m. ID refers to the experiment ID as seen in Table 1. Vertical
bars represent the 95% confidence intervals. Figure 8: The evolution of recall
for a series of jump to noise ratios varying between 1.5 and 15.5 after
grouping experiments based on source-to-sensor distance, $x_{m}$. For each
$x_{m}$, the set of jump to noise ratios are identical, however they are
plotted in an offset to improve visibility. Vertical bars represent the 95%
confidence intervals.

In practice, it is more constructive to predict the performance of the
changepoint algorithm based on the ratio of the emission rate before and after
the changepoint. To this end, Figure 9 illustrates recall as a function of
increasing leak rate ratio (LRR) for all experiments, where leak rate ratio is
the ratio of the leak rate after the change to the leak rate before the
change. In this figure, for each source-to-sensor distance, the experiments
are sorted based on coefficient of variation (CV) of $c^{y}$. For each
experiment, CV is calculated as the ratio of the standard deviation of $c^{y}$
measurements to the average $c^{y}$ measurements in the original signal (e.g.,
Figure 4a). It can be seen that a higher CV is a predictor for lower recall as
an indicator for the performance of the changepoint algorithm. This
relationship between recall and CV can be explained through a comparison
between LRR and JNR.

Figure 9: Evolution of recall for a series of leak rate ratios varying between
1.5 and 7.5 for source-to-sensor distances, $x_{m}$ of (a) 10m, (b) 20m and
(c) 30m. ID refers to the experiment ID as seen in Table 1, and CV refers to
coefficient of variation of $c^{y}$ measurements calculated for each
experiment. Vertical bars represent the 95% confidence intervals.

According to equation (10), $c^{y}$ is directly proportional to the leak rate
$Q$, therefore in our synthesized data, the leak rate ratio is the same as the
ratio of the mean $c^{y}$ after and before the change. Therefore we can write

$LRR=\frac{Q_{2}}{Q_{1}}=\frac{\mu_{c_{2}}}{\mu_{c_{1}}},$ (28)

where $Q_{2}$ and $Q_{1}$ are the emission rate after and before the change,
and $\mu_{c_{2}}$ and $\mu_{c_{1}}$ are the average $c^{y}$ measurements after
and before the change, respectively. With this definition, we can relate JNR
and LRR as follows

$JNR=\frac{\mu_{c_{2}}-\mu_{c_{1}}}{\sigma_{c}}=\frac{LRR-1}{CV},$ (29)

where $\sigma_{c}$ is the standard deviation of $c^{y}$ measurements in the
original signal of an experiment. Based on equation (29), for a constant LRR,
a higher value of CV corresponds to a smaller JNR, which according to Figure 8
points to a lower recall.

In most practical applications, delayed detection of the changepoint is
acceptable. Therefore, Figure 10 depicts the detection recall as a function of
increasing LRR. In this figure, for each source-to-sensor distance, the
experiments are sorted based on $\sigma_{c}$ as a measure of noise in the
measurements. It can be seen that experiments with higher values of
$\sigma_{c}$ correspond to higher detection recalls. This behaviour is
expected due to our data synthesis procedure, where high values of
$\sigma_{c}$ lead to significantly large $c^{y}$ measurements after the change
in emission rate which are easily detected by the changepoint detection
algorithm. Moreover, above a leak rate ratio of 3, the changepoints are rarely
missed if we account for delayed detection, therefore showing the
effectiveness of the algorithm in raising the alarm when a substantial change
in the emission rate occurs. The significance of this result can be
highlighted by noting that in a recent study of natural gas well pads in
California, it was shown that well pads for which facility-based emission
estimates were at least 3 times the component-based emission estimations were
responsible for 80% the total measured emissions [9].

Given that in some applications the change in emission rate can be
intermittent, it is also important to quantify the delay in changepoint
detection. Therefore, Figure 11 illustrates the detection delay for
experiments where the changepoints are successfully detected across all
instances (i.e., Detection Recall = 1) against increasing leak rate ratio. In
this case, there is no clear trend between the noise in the measurements and
the detection delay. However, the delay in changepoint detection monotonically
decreases with increasing LRR as expected. It is worth noting that even at the
lowest LRR where all changepoints are detected, the detection delay is less
than one sensor pass, showcasing the speed of the changepoint detection
algorithm.

Figure 10: Evolution of detection recall for a series of leak rate ratios
varying between 1.5 and 7.5 for source-to-sensor distances, $x_{m}$ of (a)
10m, (b) 20m and (c) 30m. ID refers to the experiment ID as seen in Table 1,
and $\sigma_{c}$ refers to the standard deviation of $c^{y}$ measurements
before the changepoint that is calculated for each experiment. Vertical bars
represent the 95% confidence intervals. Figure 11: Evolution of detection
delay (with units of number of passes) for a series of leak rate ratios
varying between 4.5 and 7.5 for source-to-sensor distances, $x_{m}$ of (a)
10m, (b) 20m and (c) 30m. ID refers to the experiment ID as seen in Table 1,
and $\sigma_{c}$ refers to the standard deviation of $c^{y}$ measurements
before the changepoint that is calculated for each experiment. Vertical bars
represent the 95% confidence intervals.

Next, we investigate the sensitivity of the changepoint detection algorithm to
the changepoint probability threshold, which in earlier results was set to a
value of 0.8. Figure 12 presents the false positive rate when varying the
changepoint probability threshold from 0.5 to 0.95. It can be seen that even
at the lowest chosen threshold the false positive rate is less than 12%,
highlighting the robustness of the changepoint detection algorithm.
Furthermore, while the probability of false alarms is generally higher for
experiments with larger noise, noise is not the sole predictor of the false
positive rate in experiments. Range of the measurement distribution in each
experiment, i.e., the difference between the maximum and minimum $c^{y}$
measurements in each experiment, seems to be a better predictor than standard
deviation of measurements for the false positive rate. Consequently, the
changepoint detection algorithm can be adversely affected by the presence of
outliers in the data. There are multiple possible solutions for alleviating
the sensitivity of the changepoint detection algorithm to outliers. One
possible solution is to modify the changepoint detection condition first
introduced in section 2.3. For example, the condition can be adapted such that
a changepoint is retained only if the changepoint probability is above a
threshold for multiple measurements over the next few sensor passes. This
requires the algorithm to delay resetting the Bayesian inference of equation
(11) until the detection condition is satisfied. The downside of this solution
is potential poor changepoint detection when the change in emission rate is
intermittent and temporary.

Figure 12: Evolution of false positive rate for changepoint probability
thresholds varying between 0.5 and 0.95 for source-to-sensor distances,
$x_{m}$ of (a) 10m, (b) 20m and (c) 30m. ID refers to the experiment ID as
seen in Table 1, $\sigma_{c}$ refers to the standard deviation of $c^{y}$
measurements before the changepoint, and $R$ refers to range of $c^{y}$
measurements before the changepoint for each experiment. Vertical bars
represent the 95% confidence intervals.

##  5 Conclusions

In this study, we addressed the problem of detecting changes in the emission
rate of a point-source by developing a recursive Bayesian scheme. This
methodology directly builds on a recursive Bayesian framework that was
previously used to estimate the emission rate from point sources. As a result,
the introduced recursive Bayesian methodology has the ability to
simultaneously detect changepoints in the emission rate and estimate the
emission rate before and after changepoints. In addition, we applied our
changepoint detection algorithm to a series of controlled release experiments,
where a mobile sensor traversed cross-sections of the plume emitting from a
point-source at different downwind distances, in the presence of an obstacle
close to the source. Several measures were used to evaluate the performance of
the changepoint detection methodology noting that the importance of each
performance measure depends on the practical application at hand. We found
that the changepoint algorithm is extremely effective ($>$90% success rate) in
identifying changes when the emission rate is tripled. This level of success
is significant given recent findings suggesting that majority of emissions
from the oil and gas sector can be caused by abnormal operations that
drastically increase the emission rate [9, 16]. Further, the results showed
that the statistics of the cross-plume mass concentration measurements such as
the mean, standard deviation and the range can be used as predictors of the
performance of the changepoint detection algorithm. Particularly, it was found
that at a given leak rate ratio, lower values of coefficient of variation
correspond to higher recall values which translates to higher effectiveness of
the algorithm in detecting changes immediately after they occur. Moreover, it
was shown that the false positive rate of the changepoint detection algorithm
was less than 2% when using a prescribed changepoint probability threshold of
0.8 for all the controlled release experiments.

Although the changepoint detection algorithm was applied to mobile sensor
measurements in the near field, the methodology can be easily adapted for
fenceline monitoring applications using networks of fixed sensors or far-field
measurements using a single stationary sensor. In these examples, mass
concentrations and meteorological conditions are often averaged over 30-minute
periods. By treating each 30-minute interval similar to a single sensor pass
in the experiments described in the current study, changes in the emission
rate can be found using the detection algorithm.

With the changepoint detection methodology presented here applied to
synthesized data from a single emission source, future work will be focused on
evaluating the performance of the algorithm under more real-world scenarios
such as intermittent faulty operation, and multiple emission rates caused by
various operating conditions. Moreover, for practical settings, it is
necessary to investigate the training time required to learn all the baseline
parameters for the Bayesian inference scheme, most importantly the range of
values used in the prior and the uncertainty term ($\sigma_{e}$) in the
likelihood function of equation (13) before a change occurs. More studies on
fault detection using advance sensing and measurement technologies will be
beneficial in effective and rapid identification of large emitters which can
lead to significant reductions in methane emissions from the oil and gas
industry.

## Declaration of competing interest

The authors declare they have no actual or potential competing financial
interests.

## Acknowledgments

This project was supported by David R. Atkinson Center for a Sustainable
Future (ACSF) at Cornell University, and DOE ARPA-E’s Methane Observation
Networks with Innovative Technology to Obtain Reductions (MONITOR) program
under grant DE-AR0000749.

## References

  * [1] T. F. Stocker, Climate Change 2013: The Physical Science Basis: Working Group I Contribution to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change.  Cambridge University Press, 2014.  Google-Books-ID: o4gaBQAAQBAJ.
  * [2] “U.S. Environmental Protection Agency. Inventory of U.S. Greenhouse Gas Emissions and Sinks: 1990-2019.” Available online: https://www.epa.gov/ghgemissions/inventory-us-greenhouse-gas-emissions-and-sinks-1990-2019, 2021\.  (Accessed on 30 June 2021).
  * [3] G. Pétron, A. Karion, C. Sweeney, B. R. Miller, S. A. Montzka, G. J. Frost, M. Trainer, P. Tans, A. Andrews, J. Kofler, D. Helmig, D. Guenther, E. Dlugokencky, P. Lang, T. Newberger, S. Wolter, B. Hall, P. Novelli, A. Brewer, S. Conley, M. Hardesty, R. Banta, A. White, D. Noone, D. Wolfe, and R. Schnell, “A new look at methane and nonmethane hydrocarbon emissions from oil and natural gas operations in the Colorado Denver-Julesburg Basin,” Journal of Geophysical Research: Atmospheres, vol. 119, no. 11, pp. 6836–6852, 2014.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/2013JD021272.
  * [4] A. R. Brandt, G. A. Heath, E. A. Kort, F. O’Sullivan, G. Pétron, S. M. Jordaan, P. Tans, J. Wilcox, A. M. Gopstein, D. Arent, S. Wofsy, N. J. Brown, R. Bradley, G. D. Stucky, D. Eardley, and R. Harriss, “Methane Leaks from North American Natural Gas Systems,” Science, vol. 343, pp. 733–735, Feb. 2014.  Publisher: American Association for the Advancement of Science Section: Policy Forum.
  * [5] C. Frankenberg, A. K. Thorpe, D. R. Thompson, G. Hulley, E. A. Kort, N. Vance, J. Borchardt, T. Krings, K. Gerilowski, C. Sweeney, S. Conley, B. D. Bue, A. D. Aubrey, S. Hook, and R. O. Green, “Airborne methane remote measurements reveal heavy-tail flux distribution in Four Corners region,” Proceedings of the National Academy of Sciences, vol. 113, pp. 9734–9739, Aug. 2016.
  * [6] A. P. Ravikumar, J. Wang, M. McGuire, C. S. Bell, D. Zimmerle, and A. R. Brandt, ““Good versus Good Enough?” Empirical Tests of Methane Leak Detection Sensitivity of a Commercial Infrared Camera,” Environmental Science & Technology, vol. 52, pp. 2368–2374, Feb. 2018.  Publisher: American Chemical Society.
  * [7] D. Zavala-Araiza, D. R. Lyon, R. A. Alvarez, K. J. Davis, R. Harriss, S. C. Herndon, A. Karion, E. A. Kort, B. K. Lamb, X. Lan, A. J. Marchese, S. W. Pacala, A. L. Robinson, P. B. Shepson, C. Sweeney, R. Talbot, A. Townsend-Small, T. I. Yacovitch, D. J. Zimmerle, and S. P. Hamburg, “Reconciling divergent estimates of oil and gas methane emissions,” Proceedings of the National Academy of Sciences, vol. 112, pp. 15597–15602, Dec. 2015.
  * [8] D. Zavala-Araiza, D. Lyon, R. A. Alvarez, V. Palacios, R. Harriss, X. Lan, R. Talbot, and S. P. Hamburg, “Toward a Functional Definition of Methane Super-Emitters: Application to Natural Gas Production Sites,” Environmental Science & Technology, vol. 49, pp. 8167–8174, July 2015.  Publisher: American Chemical Society.
  * [9] X. Zhou, S. Yoon, S. Mara, M. Falk, T. Kuwayama, T. Tran, L. Cheadle, J. Nyarady, B. Croes, E. Scheehle, J. D. Herner, and A. Vijayan, “Mobile sampling of methane emissions from natural gas well pads in California,” Atmospheric Environment, vol. 244, p. 117930, Jan. 2021.
  * [10] R. A. Alvarez, D. Zavala-Araiza, D. R. Lyon, D. T. Allen, Z. R. Barkley, A. R. Brandt, K. J. Davis, S. C. Herndon, D. J. Jacob, A. Karion, E. A. Kort, B. K. Lamb, T. Lauvaux, J. D. Maasakkers, A. J. Marchese, M. Omara, S. W. Pacala, J. Peischl, A. L. Robinson, P. B. Shepson, C. Sweeney, A. Townsend-Small, S. C. Wofsy, and S. P. Hamburg, “Assessment of methane emissions from the U.S. oil and gas supply chain,” Science, vol. 361, pp. 186–188, July 2018.  Publisher: American Association for the Advancement of Science Section: Report.
  * [11] C. W. Rella, T. R. Tsai, C. G. Botkin, E. R. Crosson, and D. Steele, “Measuring Emissions from Oil and Natural Gas Well Pads Using the Mobile Flux Plane Technique,” Environmental Science & Technology, vol. 49, pp. 4742–4748, Apr. 2015.  Publisher: American Chemical Society.
  * [12] M. Omara, M. R. Sullivan, X. Li, R. Subramanian, A. L. Robinson, and A. A. Presto, “Methane Emissions from Conventional and Unconventional Natural Gas Production Sites in the Marcellus Shale Basin,” Environmental Science & Technology, vol. 50, pp. 2099–2107, Feb. 2016\.  Publisher: American Chemical Society.
  * [13] A. M. Robertson, R. Edie, D. Snare, J. Soltis, R. A. Field, M. D. Burkhart, C. S. Bell, D. Zimmerle, and S. M. Murphy, “Variation in Methane Emission Rates from Well Pads in Four Oil and Gas Basins with Contrasting Production Volumes and Compositions,” Environmental Science & Technology, vol. 51, pp. 8832–8840, Aug. 2017.  Publisher: American Chemical Society.
  * [14] A. R. Brandt, G. A. Heath, and D. Cooley, “Methane Leaks from Natural Gas Systems Follow Extreme Distributions,” Environmental Science & Technology, vol. 50, pp. 12512–12520, Nov. 2016.  Publisher: American Chemical Society.
  * [15] T. L. Vaughn, C. S. Bell, C. K. Pickering, S. Schwietzke, G. A. Heath, G. Pétron, D. J. Zimmerle, R. C. Schnell, and D. Nummedal, “Temporal variability largely explains top-down/bottom-up difference in methane emission estimates from a natural gas production region,” Proceedings of the National Academy of Sciences, vol. 115, pp. 11712–11717, Nov. 2018.  ISBN: 9781805687115 Publisher: National Academy of Sciences Section: Physical Sciences.
  * [16] D. Zavala-Araiza, R. A. Alvarez, D. R. Lyon, D. T. Allen, A. J. Marchese, D. J. Zimmerle, and S. P. Hamburg, “Super-emitters in natural gas infrastructure are caused by abnormal process conditions,” Nature Communications, vol. 8, p. 14012, Jan. 2017.  Number: 1 Publisher: Nature Publishing Group.
  * [17] R. M. Duren, A. K. Thorpe, K. T. Foster, T. Rafiq, F. M. Hopkins, V. Yadav, B. D. Bue, D. R. Thompson, S. Conley, N. K. Colombi, C. Frankenberg, I. B. McCubbin, M. L. Eastwood, M. Falk, J. D. Herner, B. E. Croes, R. O. Green, and C. E. Miller, “California’s methane super-emitters,” Nature, vol. 575, pp. 180–184, Nov. 2019.  Bandiera_abtest: a Cg_type: Nature Research Journals Number: 7781 Primary_atype: Research Publisher: Nature Publishing Group Subject_term: Atmospheric science;Carbon cycle Subject_term_id: atmospheric-science;carbon-cycle.
  * [18] S. Coburn, C. B. Alden, R. Wright, K. Cossel, E. Baumann, G. Truong, F. Giorgetta, C. Sweeney, N. R. Newbury, K. Prasad, I. Coddington, and G. B. Rieker, “Regional trace-gas source attribution using a field-deployed dual frequency comb spectrometer,” Optica, vol. 5, pp. 320–327, Apr. 2018.  Publisher: Optical Society of America.
  * [19] ARPA-E, “Methane Observation Networks with Innovative Technology to Obtain Reductions.” Available online: https://arpa-e.energy.gov/technologies/programs/monitor, 2015.  (Accessed on 30 June 2021).
  * [20] “U.S. Environmental Protection Agency. Petroleum Refinery Sector Rule (Risk and Technology Review and New Source Performance Standards).” Available online: https://www.epa.gov/stationary-sources-air-pollution/petroleum-refinery-sector-rule-risk-and-technology-review-and-new, 2020\.  (Accessed on 30 June 2021).
  * [21] Colorado Department of Public Health and Environment, “Regulation 7: Control of Ozone via Ozone Precursors and Control of Hydrocarbons via Oil and Gas Emissions (Emissions of Volatile Organic Compounds and Nitrogen Oxides).” Available online: https://cdphe.colorado.gov/aqcc-regulations, 2021.  (Accessed on 14 July 2021).
  * [22] E. Yee, “Bayesian probabilistic approach for inverse source determination from limited and noisy chemical or biological sensor concentration measurements,” in Chemical and Biological Sensing VIII, vol. 6554, p. 65540W, International Society for Optics and Photonics, Apr. 2007.
  * [23] R. P. Adams and D. J. C. MacKay, “Bayesian Online Changepoint Detection,” arXiv:0710.3742 [stat], Oct. 2007.  arXiv: 0710.3742.
  * [24] X. Zhou, A. Montazeri, and J. D. Albertson, “Mobile sensing of point-source gas emissions using Bayesian inference: An empirical examination of the likelihood function,” Atmospheric Environment, vol. 218, p. 116981, Dec. 2019.
  * [25] Z. Warhaft, “Passive Scalars in Turbulent Flows,” Annual Review of Fluid Mechanics, vol. 32, pp. 203–240, Jan. 2000.  Publisher: Annual Reviews.
  * [26] R. B. Stull, An Introduction to Boundary Layer Meteorology.  Dordrecht, Netherlands: Kluwer Academic Publishers, first ed., 1988.
  * [27] A. M. Obukhov, “Turbulence in an atmosphere with a non-uniform temperature,” Boundary-Layer Meteorology, vol. 2, pp. 7–29, Mar. 1971.
  * [28] T. W. Horst and J. C. Weil, “Footprint estimation for scalar flux measurements in the atmospheric surface layer,” Boundary-Layer Meteorology, vol. 59, pp. 279–296, May 1992.
  * [29] J. C. Kaimal, Atmospheric Boundary Layer Flows: Their Structure and Measurement.  Oxford, UK: Oxford Universtiy Press, 1994.
  * [30] J. L. Lumley and H. A. Panofsky, The structure of atmospheric turbulence (Interscience monographs and texts in physics and astronomy).  New York, NY, USA: Wiley, 1964.
  * [31] J. D. Albertson, T. Harvey, G. Foderaro, P. Zhu, X. Zhou, S. Ferrari, M. S. Amin, M. Modrak, H. Brantley, and E. D. Thoma, “A Mobile Sensing Approach for Regional Surveillance of Fugitive Methane Emissions in Oil and Gas Production,” Environmental Science & Technology, vol. 50, pp. 2487–2497, Mar. 2016.  Publisher: American Chemical Society.
  * [32] M. Arulampalam, S. Maskell, N. Gordon, and T. Clapp, “A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking,” IEEE Transactions on Signal Processing, vol. 50, pp. 174–188, Feb. 2002.  Conference Name: IEEE Transactions on Signal Processing.
  * [33] H. L. Brantley, E. D. Thoma, W. C. Squier, B. B. Guven, and D. Lyon, “Assessment of Methane Emissions from Oil and Gas Production Pads using Mobile Measurements,” Environmental Science & Technology, vol. 48, pp. 14508–14515, Dec. 2014.  Publisher: American Chemical Society.
  * [34] E. Yee, “Theory for Reconstruction of an Unknown Number of Contaminant Sources using Probabilistic Inference,” Boundary-Layer Meteorology, vol. 127, pp. 359–394, June 2008.
  * [35] E. T. Jaynes, “Prior Probabilities,” IEEE Transactions on Systems Science and Cybernetics, vol. 4, pp. 227–241, Sept. 1968.  Conference Name: IEEE Transactions on Systems Science and Cybernetics.
  * [36] E. T. Jaynes, Probability Theory: the Logic of Science.  Cambridge university press, 2003.
  * [37] A. Keats, E. Yee, and F. Lien, “Bayesian inference for source determination with applications to a complex urban environment,” Atmospheric Environment, vol. 41, pp. 465–479, Jan. 2007.
  * [38] E. Yee and T. K. Flesch, “Inference of emission rates from multiple sources using Bayesian probability theory,” Journal of Environmental Monitoring, vol. 12, pp. 622–634, Mar. 2010.  Publisher: The Royal Society of Chemistry.
  * [39] J. D. Wilson and B. L. Sawford, “Review of Lagrangian stochastic models for trajectories in the turbulent atmosphere,” Boundary-Layer Meteorology, vol. 78, pp. 191–210, Feb. 1996.
  * [40] D. J. Thomson, “Criteria for the selection of stochastic models of particle trajectories in turbulent flows,” Journal of Fluid Mechanics, vol. 180, pp. 529–556, July 1987.  Publisher: Cambridge University Press.
  * [41] K. S. Rao, “Uncertainty Analysis in Atmospheric Dispersion Modeling,” pure and applied geophysics, vol. 162, pp. 1893–1917, Oct. 2005.
  * [42] X. Zhou, F. H. Passow, J. Rudek, J. C. von Fisher, S. P. Hamburg, and J. D. Albertson, “Estimation of methane emissions from the U.S. ammonia fertilizer industry using a mobile sensing approach,” Elementa: Science of the Anthropocene, vol. 7, May 2019.
  * [43] C. Forbes, M. Evans, N. Hastings, and B. Peacock, Statistical Distributions.  Wiley, 4th ed., 2010.
  * [44] M. J. Flynn, “Some Computer Organizations and Their Effectiveness,” IEEE Transactions on Computers, vol. C-21, pp. 948–960, Sept. 1972.  Conference Name: IEEE Transactions on Computers.
  * [45] T. A. Foster-Wittig, E. D. Thoma, and J. D. Albertson, “Estimation of point source fugitive emission rates from a single sensor time series: A conditionally-sampled Gaussian plume reconstruction,” Atmospheric Environment, vol. 115, pp. 101–109, Aug. 2015.
  * [46] B. Efron and R. Tibshirani, “Bootstrap Methods for Standard Errors, Confidence Intervals, and Other Measures of Statistical Accuracy,” Statistical Science, vol. 1, pp. 54–75, Feb. 1986.  Publisher: Institute of Mathematical Statistics.