Currently in the petroleum industry, operators often flare the produced gas instead of commodifying it. The flaring magnitudes are large in some states, which constitute problems with energy waste and CO2 emissions. In North Dakota, operators are required to estimate and report the volume flared. The questions are, how good is the quality of this reporting, and what insights can be drawn from it?

Apart from the company-reported statistics, which are available from the North Dakota Industrial Commission (NDIC), flared volumes can be estimated via satellite remote sensing, serving as an unbiased benchmark. Since interpretation of the Landsat 8 imagery is hindered by artifacts due to glow, the estimated volumes based on the Visible Infrared Imaging Radiometer Suite (VIIRS) are used. Reverse geocoding is performed for comparing and contrasting the NDIC and VIIRS data at different levels, such as county and oilfield.

With all the data gathered and preprocessed, Bayesian learning implemented by Markov chain Monte Carlo methods is performed to address three problems: county level model development, flaring time series analytics, and distribution estimation. First, there is heterogeneity among the different counties, in the associations between the NDIC and VIIRS volumes. In light of such, models are developed for each county by exploiting hierarchical models. Second, the flaring time series, albeit noisy, contains information regarding trends and patterns, which provide some insights into operator approaches. Gaussian processes are found to be effective in many different pattern recognition scenarios. Third, distributional insights are obtained through unsupervised learning. The negative binomial and Gaussian mixture models are found to effectively describe the oilfield flare count and flared volume distributions, respectively. Finally, a nearest-neighbor-based approach for operator level monitoring and analytics is introduced.

Vectors and matrices are in bold type. A subscript asterisk, such as in $y_*$, indicates reference to a test set quantity or a prediction.

*General Nomenclature
*Probability and Statistics
*Probability Distributions
*Gaussian Processes
*Vectors and Matrices
*Standard Functions

[General Nomenclature]Data set: $\mathcal{D}=\{(\mathbf{x}_i, y_i) \mid i=1, \ldots , n\}$$\mathcal{D}$
[General Nomenclature]Prediction for a test input $\mathbf{x}_*$$y_*$
[General Nomenclature]Proportional to; e.g., $p(x \mid y) \propto p(x,y)$ means that $p(x\mid y)$ is
equal to $p(x,y)$ times a factor which is independent of $x$$\propto$
[General Nomenclature]$a$ is defined as $b$$a \coloneqq b$
[General Nomenclature]$b$ is defined as $a$$a \eqqcolon b$
[General Nomenclature]$\ell_2$ norm$\norm{\cdot}$
[General Nomenclature]$n$-dimensional vector space of real numbers$\mathds{R}^{n}$
[General Nomenclature]Natural numbers with zero$\mathds{N}_{0}$
[General Nomenclature]Real numbers$\mathds{R}$
[General Nomenclature]Pi (italic) representing a variable$\pi$
[General Nomenclature]Pi (upright) denoting the transcendental constant ($3.14159\dotsm$)$\uppi$
[General Nomenclature]Universal quantifier: for all $x$$\forall x$

[Probability and Statistics]Conditional probability density function$p(\cdot \mid \cdot)$
[Probability and Statistics]Probability density function$p(\cdot)$
[Probability and Statistics]Probability mass function$P(\cdot)$
[Probability and Statistics]Random variable $X$ is distributed according to $p$$X \sim p$
[Probability and Statistics]Expectation; expectation of $g(x)$ when $x \sim p$$\mathds{E}$ or $\mathds{E}_{p}[g(x)]$
[Probability and Statistics]Variance; variance of $g(x)$ when $x \sim p$$\mathds{V}$ or $\mathds{V}_{p}[g(x)]$

[Probability Distributions]Univariate Gaussian distribution with parameters $\mu, \, \sigma$$\mathcal{N}(\mu, \, \sigma)$
[Probability Distributions]Multivariate Gaussian distribution with parameters $\bm{\mu}, \, \bm{\Sigma}$$\operatorname{MVNormal}(\bm{\mu}, \, \bm{\Sigma})$
[Probability Distributions]Poisson distribution with parameter $\lambda$$\operatorname{Poisson}(\lambda)$
[Probability Distributions]Exponential distribution with parameter $\lambda$$\operatorname{Exponential}(\lambda)$
[Probability Distributions]Binomial distribution with parameters $n, \, p$$\operatorname{Binomial}(n, \, p)$
[Probability Distributions]Gamma distribution with parameters $\alpha, \, \beta$$\operatorname{Gamma}(\alpha, \, \beta)$
[Probability Distributions]Half-Normal distribution with parameter $\sigma$$\operatorname{Half-Normal}(\sigma)$
[Probability Distributions]Half-Cauchy distribution with parameter $\gamma$$\operatorname{Half-Cauchy}(\gamma)$
[Probability Distributions]LKJ distribution with parameter $\eta$$\operatorname{LKJcorr}(\eta)$
[Probability Distributions]Distribution over Cholesky decomposed covariance
matrices with parameters $\eta, \, \bm{\sigma}$$\operatorname{LKJCholeskyCov}(\eta, \, \bm{\sigma})$
[Probability Distributions]Student's $t$-distribution with parameters $\nu, \, \mu, \, \sigma$$\stut (\nu, \, \mu, \, \sigma)$
[Probability Distributions]Dirichlet distribution with parameter $\bm{\alpha}$$\operatorname{Dirichlet}(\bm{\alpha})$
[Probability Distributions]Categorical distribution with parameter $\mathbf{p}$$\operatorname{Categorical}(\mathbf{p})$
[Probability Distributions]Negative binomial distribution with parameters $\mu, \phi$$\operatorname{NegBinomial}(\mu, \phi)$
[Probability Distributions]Continuous uniform distribution with parameters $a, b$$\operatorname{Uniform}(a, b)$

[Gaussian Processes]Gaussian process: $f \sim \mathcal{GP}(m(\mathbf{x}), k(\mathbf{x}, \mathbf{x}'))$, the function $f$ is
distributed as a Gaussian process$\mathcal{GP}$
[Gaussian Processes]Mean function evaluated at $\mathbf{x}$$m(\mathbf{x})$
[Gaussian Processes]Covariance function evaluated at $\mathbf{x}$ and $\mathbf{x}'$$k(\mathbf{x}, \mathbf{x}')$
[Gaussian Processes]Vector of latent function values, $\mathbf{f} = (f(\mathbf{x}_1), \dots, f(\mathbf{x}_n))^\top$$\mathbf{f}$

[Vectors and Matrices]Transpose of matrix $\mathbf{L}$$\mathbf{L}^\top$
[Vectors and Matrices]Vector of parameters$\bm\theta$
[Vectors and Matrices]Cholesky decomposition: $\mathbf{L}$ is a lower triangular matrix
such that $\mathbf{L} \cdot \mathbf{L}^\top=\mathbf{K}$$\operatorname{Cholesky}(\mathbf{K})$
[Vectors and Matrices]Vector of all $0$'s of length $n$$\mathbf{0}_n$
[Vectors and Matrices]Vector of all $1$'s of length $n$$\mathds{1}_n$
[Vectors and Matrices]Identity matrix of size $n \times n$$\mathbf{I}_n$

[Standard Functions]Natural exponential function$\exp(\cdot)$
[Standard Functions]Natural logarithm function$\log(\cdot)$
[Standard Functions]Inverse-logit function$\operatorname{logit}^{-1}(\cdot)$

[Units]thousand cubic feet
[Units]million cubic feet
[Units]billion cubic meter
[Units]barrels of oil equivalent
[Units]barrels per day
[Units]dollars per barrel
[Units]dollars per million Btu
[Units]thousand cubic feet per barrel

North Dakota Industrial CommissionNDIC
Landsat 8L8
Deep learningDL
Visible Infrared Imaging Radiometer SuiteVIIRS
VIIRS NightfireVNF

National Oceanic and Atmospheric AdministrationNOAA
National Aeronautics and Space AdministrationNASA
Markov chain Monte CarloMCMC
Hamiltonian Monte CarloHMC
Independent and identically distributedi.i.d.
Radiant heatRH
Short-wave infraredSWIR

Exempli gratia (Latin: for example)e.g.
Id est (Latin: that is)i.e.
Gaussian processGP
Gaussian mixture modelGMM
Autoregressive integrated moving averageARIMA
Long short-term memoryLSTM
World Geodetic SystemWGS
North American DatumNAD
Maximum likelihood estimationMLE
Maximum a posterioriMAP
Hierarchical Density-Based Spatial Clustering
of Applications with NoiseHDBSCAN
Density-Based Spatial Clustering of Applications
with NoiseDBSCAN
Gas oil ratioGOR
Credible intervalCI
Prediction intervalPI
Right-hand sideRHS
Highest density intervalHDI
Interquartile rangeIQR
Energy Information AdministrationEIA
Zero-inflated PoissonZIP
Zero-inflated negative binomialZINB
Kernel density estimationKDE

In the very first place, I want to express my deepest appreciation to my advisor, Dr. Jennifer L. Miskimins. Dr. Miskimins has been my MS/PhD advisor and mentor since 2011. Since I returned to Mines to start my PhD in 2017, Dr. Miskimins has been providing me with the best guidance, the greatest support, and the most opportunities that I could imagine. During the first semester, I worked as a lab assistant in the High Bay; in a later semester, I worked as a teaching assistant in her well stimulation course; ever since I started to become interested in machine learning, she has provided me with a huge number of opportunities to connect with different groups of people, for brainstorming and pursuing my research interest. To a certain extent, I feel like I finally become a “qualified” FAST student member, thanks to all of these precious experience. What I have achieved, including this dissertation, would have never been possible without the guidance and support from Dr. Miskimins. Her world-class technical expertise, attitudes toward work/life, and art of managing different teams at various levels are what I hope I can learn from in my career and personal life.

I am deeply grateful to my dissertation committee members: Dr. Soutir Bandyopadhyay, Dr. Alfred W. Eustes III, Dr. Yilin Fan, and Prof. Jim Crompton. My competency in my research field, as well as the shape of this dissertation are built with the help of those fruitful discussions and insightful comments from them.

I am indebted to my mentors, colleagues, and friends from the Payne Institute for Public Policy. Especially, I want to thank Dr. Mikhail N. Zhizhin, Dr. Christopher D. Elvidge, and Dr. Morgan D. Bazilian. It is such an eye-opening experience for me to work with these world-class experts in remote sensing and satellite imagery. I would particularly like to thank Dr. Zhizhin for his help, insights, and time.

I am really grateful to Dr. Bandyopadhyay and Dr. Luis Tenorio from the AMS Department, for their fantastic teaching, knowing me personally and motivating me to work hard. Looking back at what I have learned in machine learning which makes this dissertation possible, taking their classes are definitely the most important resources for myself (excuse me for not being a probabilist at this moment). By taking their statistical methods classes, I started to appreciate what really is machine learning, and falling in love with mathematics, more specifically, probability theory and statistical modeling.

The TA experience at Mines makes me a better PhD student. What I have learned, technical or non-technical, made their way into this dissertation. I want to thank Prof. Crompton, Dr. Eustes, Dr. Mark G. Miller, Dr. Linda A. Battalora, and Dr. Miskimins for providing me with those valuable TA opportunities. I am grateful to all of my students for their support and feedback.

I want to thank Dr. Yu-Shu Wu, Dr. Xiaolong Yin, and Dr. Yilin Fan for their care, support, and encouragement throughout my PhD study. I would like to thank Denise Winn-Bower, Rachel McDonald, and Joe Chen for their help.

I really appreciate the feedback from the FAST member companies' representatives. A lot of the discussions and the reflections following those were incorporated into this dissertation. Especially, I want to thank Ty Woodworth for his time and help, in the process of collecting plunger lift data for me. I got very warm welcomes every time I visited their Windsor office in Northern Colorado. Ty kindly introduced me to the team he led, and I got the great opportunities to ask questions and discuss with many field experts in different areas. Those discussions helped me tremendously.

Special thanks go to the open source community. In the process of conducting this research and typesetting this dissertation, I benefited a lot from the ecosystems around Linux/GNU, /, and Python. Especially, I want to thank the people behind PyMC3, a probabilistic programming language that this dissertation is heavily dependent upon.

Last but not least, I would like to thank my family and friends. Thank you to my beloved wife Xiaodan, for all her love, support, and delicious dishes. I also want to thank my parents and parents-in-law for their support, encouragement, and understanding.

I dedicate this work to my mother, Dr. Lingying Ni, and my father, Mr. Honggang Lu.

CHAPTER: INTRODUCTION

Currently in the petroleum industry, for wells which produce both crude oil and natural gas, operators often choose to flare the produced gas instead of commodifying it. The rationales behind such decisions are multifold. Variations in natural gas price can be an important factor, especially when the processing and transportation cost is higher than the value of gas Srivastava2019DistributionalPolicy. The amount of gas being flared each year on a national level is huge, and an increasing trend can be observed for the top flaring countries (<ref>).

[!htbp]top30_flare_nationfigures/NewrankingTop30flaringcountries20142018Top 30 countries ranked by flared gas volume in 2018. United States ranks No. 4 and has a large increase from 2017 to 2018 WorldBank2019GlobalData.

Due to the boom of unconventional resources (e.g., shale gas reservoirs) development in the recent decade, the United States has been among the top flaring countries in terms of total volume flared. This is backed by the data from the U.S. Energy Information Administration (EIA)  showing North Dakota, which is underlain by the Bakken Formation, and Texas, which houses the Permian Basin and the Eagle Ford Shale, are the top two flaring states since 2013. The two states' annual flaring volume time series are shown in <ref>. Some flaring sites can be clearly identified from Google Earth's imagery (<ref>).

[!htbp]nd_vs_txfigures/nd_vs_txThe time series show the trend of gas flaring for the top two states in the United States EIA2019NaturalFlared. Texas regained the lead in 2015.

[!htbp]goog_earth_day_imgfigures/googleEarth_bakken_flare_aThis Google Earth imagery shows gas flaring being conducted on a well location in North Dakota GoogleEarth2019SatelliteDakota.

Natural gas flaring constitutes a problem of energy waste and CO2 emissions. In recent years, various organizations and government agencies have advocated reducing or eliminating routine gas flaring. For example, the  (NDIC) introduced a gas flaring regulatory policy (Order 24665) in 2014, with goals of reducing flaring in different aspects (e.g., volume of gas flared). The World Bank launched the “Zero Routine Flaring by 2030” initiative in 2015. To monitor and benchmark flaring activity's magnitude, a precise and accurate method to obtain quantitative flaring information is desirable. However, in certain situations, this information is only available through self-reporting mechanisms. Inaccuracies might be introduced either intentionally or unintentionally.

Satellite remote sensing is one unbiased approach for solving this problem. It can help detect active flares especially during nighttime and can be used to calibrate the estimation for flared gas volume. For this work, two different types of sensors are considered, including the Landsat 8 (L8)'s Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS), as well as the Visible Infrared Imaging Radiometer Suite (VIIRS) that is on the Suomi National Polar-orbiting Partnership (NPP) and NOAA-20 satellites. In the remainder of this dissertation, they are referred to as L8 and VIIRS, respectively. An example of detecting flaring with VIIRS low light imaging data is shown in <ref>.

[!htbp]eog_viirs_poster_crop_usfigures/eog_viirs_poster_crop_usPart of the original poster EarthObservationGroupatPayneInstitute2019EOGFile which uses one year accumulation of VIIRS low light imaging data to showcase human activities, e.g., gas flaring, fishing, and city lights. As annotated, North Dakota's flaring activities are very visible from space at night.

§.§ Research Goal

This research is undertaken to achieve the following goals:

* Evaluate the methodology for estimating flared gas volume leveraging satellite imagery; and,

* Find insights into operators' gas flaring behavior.

§.§ Dissertation Objectives

To achieve the goals outlined in Section <ref>, more specific objectives are listed below:

* Compare and contrast the flaring data from VIIRS and NDIC.

* Compare the VIIRS flared volumes to the NDIC, using the NDIC as a benchmark.

* Evaluate the effectiveness of using Landsat 8 nighttime images to improve flare detection and volume estimation.

* Determine the detection limits of Landsat 8 and compare it with VIIRS' capabilities.

* Investigate operator approaches for gas flaring.

* Determine the correlation between gas price / oil price / oil production and flared gas volume.

* Evaluate if the North Dakota regulatory policy (Order 24665) achieved its goals.

* Develop a model that can predict flared gas volume at a state level.

* Find any hidden structure/clusters from all the producing entities.

§.§ Outline and Contributions

The main contribution of this dissertation is demonstrating that Bayesian learning implemented by Markov chain Monte Carlo methods is very effective in flaring data analytics. A series of parametric and nonparametric machine learning models are developed for various analytics goals and granularities, providing direct guidance for future modeling endeavors. To demonstrate the effectiveness and robustness, they are all tested with real data. The superiority of this approach is based on the fact that the inference stage is entirely probabilistic, in that the parametric uncertainties arising from probable models as well as the stochastic uncertainties arising from noisy observations are all properly characterized and quantified. It makes the extracted insights robust and interpretable for decision- and policy-making by, for example, a state government.

In Chapter <ref>, a literature review is given for the state of the art in satellite imagery processing, Bayesian inference, Markov chain Monte Carlo methods, and machine learning.

In Chapter <ref>, the data gathering processes are discussed. Results from some exploratory data analysis are presented.

In Chapter <ref>, county level models are built to study the correlations between VIIRS and NDIC, and to explore the heterogeneity among the counties in North Dakota.

In Chapter <ref>, flaring time series analytics is presented for the purposes of revealing trends and patterns at different levels.

In Chapter <ref>, unsupervised learning is applied on flaring data to characterize the latent structures.

In Chapter <ref>, a method of operator level monitoring and analytics is introduced, and some discussions about applying Bayesian learning are given.

In Chapter <ref>, major conclusions drawn are presented. Recommendations based on this work are given. A number of future research areas are outlined.

CHAPTER: LITERATURE REVIEW

In the 1990s, the World Bank started gathering nighttime satellite images, from which big cities and oilfields were both bright and needed to be sorted using extra information. The situation changed in 2012 when infrared data became available from VIIRS Rassenfoss2018EP2018. One of the data products, VIIRS Nightfire (VNF) specializes in natural gas flaring observation and is even able to distinguish between biomass burning and gas flaring Elvidge2017SupportingViirs.

VNF's development was based upon VIIRS imagery. To improve the performance of flare detection and gas volume estimation, other sources of information, such as L8 imagery, can be leveraged. <ref> presents a comparison of L8 and VIIRS spatial and temporal resolutions NASA2019Landsat8, Wikipedia2019VisibleSuite. <ref> illustrates L8's spatial resolution. In addition, L8 collects data in 11 different spectral bands of the electromagnetic spectrum. VIIRS has 22 bands. Both L8 and VIIRS are in near-polar orbits of the earth and can reveal rich features in the landscape. Therefore, L8 should be able to identify smaller gas flares compared to VIIRS' capability, although its longer satellite revisit time poses a challenge to identify less persistent flares. More details on the processing steps of VNF are discussed in Section <ref>, the essence of which will be applied to L8.

Resolutions of Landsat 8 and VIIRS

table-text-alignment = right

2cResolution Type
Spatial []     Temporal []

Landsat 8     15100$\dagger$     16.0$\ddagger$

VIIRS     375750$\dagger$     0.5

$\dagger$ Depends on the band of the electromagnetic spectrum

$\ddagger$ For daytime mode

[!htbp]L8pixelsizesfigures/L8pixelsizesLandsat 8's spatial resolution NASA2020LandsatOverview. Each Landsat pixel (30 by 30 meter area) is roughly the size of a baseball diamond.

Nowadays, one resource which is more than abundant is data. For a certain discipline or research field, new sources of data bring in new dimensions of information, such as satellite images are now playing a role in gas flaring analytics. How to analyze data effectively and intelligently to gain insights is a central problem. In the petroleum engineering domain, for example, data driven approaches have been proposed to analyze stimulation treatments Kazakov2011ApplicationSystems and predict screenouts Yu2020AWells. Machine learning is a powerful tool for this purpose. It is at the core of artificial intelligence and data science, and lies at the intersection of statistics and computer science Jordan2015MachineProspects. Frameworks in computational learning theory, such as the PAC learning proposed by pac, help provide a theoretical backbone for some learning algorithms.

One subset of machine learning, deep learning (DL), had its debut in 2006 when  introduced Deep Belief Networks (DBN), but it did not gain wide acceptance until 2012 when AlexNet showed the breakthrough performance on classification accuracy in the ImageNet competition NIPS2012_4824. AlexNet is a DL-based model (more specifically a convolutional neural network) and achieved an error rate of 15.3, which is more than 10 lower than the runner-up. DL dominated the competition thereafter, and DL-based models finally surpassed human performance on the classification data set in 2015 He_2015_ICCV.

Although neural network-based models have gained much success in recent years, it should be noted that no one type of model can always be the best candidate for all problems. This has been formally shown by Wolpert1996TheAlgorithms, and is usually referred to as the “no free lunch” (NFL) theorem. More recently, 2017arXiv170805070O empirically assessed 13 classification algorithms on 165 different problem sets, and the results aligned with the theorem: even the union of the top five best performing algorithms cannot dominate all of the problem sets.

In the following sections, a detailed review is given for the aspects below, which serve as the foundation and inspiration for this work:

* Satellite image processing

* Bayesian inference

* Markov chain Monte Carlo

* Machine learning

* Analytics toolset

§.§ Satellite Image Processing

Satellite images are utilized to estimate flared gas volume. The fire detection algorithm based on Planck curve fitting and physical laws, known as VIIRS Nightfire (VNF) due to Elvidge2013VIIRSNight, serves as a starting point for analyzing L8 images in this research. The method consists of several major steps:

* Detection of hot pixels

During nighttime, the sensors mainly record instrument noise which approximately follows a Gaussian distribution, except for the few pixels that contain an infrared emitter such as a gas flare. Therefore hot pixels can be identified by setting a cutoff on the tail of the distribution, e.g., those pixels with digital numbers exceeding the mean plus four standard deviations.

* Noise filtering

Hot pixels that are detected in only one spectral band are treated as noise and filtered out.

* Atmospheric correction

Losses in radiance due to scattering and absorption effects can be corrected. MODTRAN ® 5 modtran, parameterized with atmospheric water vapor and temperature profiles, is used to derive the correction coefficients for each spectral band.

* Planck curve fitting

Planck curves are modeled for gas flares, which appear as gray bodies because they are sub-pixel sources. Therefore the output of the fitting is an estimate of the temperature and an emission scaling factor (the emissivity term in the Planck function). The latter is used subsequently to estimate the source area.

* Calculation of source area

The source area $S$ is calculated using

\begin{equation}
S = \varepsilon A\,,
\end{equation}

where $\varepsilon$ is the emission scaling factor and $A$ is the size of the pixel footprint.

* Calculation of radiant heat

The radiant heat is calculated using the Stefan–Boltzmann law:

\begin{equation}
\mathrm{RH} = \sigma T^4 S\,,
\end{equation}

where $\mathrm{RH}$ is the radiant heat in , $\sigma$ is the Stefan–Boltzmann constant, $T$ is the temperature in , and $S$ is the source area in □.

Once $\mathrm{RH}$ is obtained, previous work by Elvidge2015MethodsData developed a calibration for estimating flared gas volume, utilizing nation-level flaring reporting provided by Cedigaz2015AnnualCountry and state-level reporting from Texas and North Dakota. The developed calibration can then be applied to each individual flaring site worldwide for estimation of flared gas volume, etc.

§.§ Bayesian Inference

Bayesian inference leverages conditional probability theory to establish a formal procedure for learning from data Betancourt2018ConditionalEngineers. Bayesian models provide full joint probability distributions $p(\mathcal{D},\bm\theta)$ over observable data $\mathcal{D}$ and unobservable model parameters $\bm\theta$. The essence of Bayesian analysis is to obtain the posterior distribution $p(\bm\theta \mid \mathcal{D})$, which characterizes the conditional probability of parameters $\bm\theta$ given some data $\mathcal{D}$. It can be derived through Bayes' theorem:

\begin{align}
p(\bm\theta \mid \mathcal{D}) &= \dfrac{p(\mathcal{D} \mid \bm\theta) \, p(\bm\theta)}{p(\mathcal{D})} \label{eq:bayes_thm} \\[0.5ex]
&= \dfrac{p(\mathcal{D} \mid \bm\theta) \, p(\bm\theta)}{\int p(\mathcal{D} \mid \bm\theta') \, p(\bm\theta')\diff\bm\theta'} \label{eq:posterior-total} \\[0.5ex]
&\propto p(\mathcal{D} \mid \bm\theta) \, p(\bm\theta)\,,
\end{align}

where $p(\mathcal{D} \mid \bm\theta)$ is the likelihood (also referred to as the observation model) which denotes how likely the data is given a certain set of parameters, and $p(\bm\theta)$ is the prior which models the probability of the parameters before observing any data. The prior encodes domain expertise. Once some observations are given, it is updated into a posterior which quantifies how consistent the model configurations are with both the domain knowledge and the observed data Betancourt2018ConditionalEngineers. After the posterior is obtained, most if not all inferential questions can then be answered with posterior expectation values of certain functions Betancourt2019ProbabilisticInference:

\begin{equation}
\mathds{E}_p [ g(\bm\theta) ]
\int g(\bm\theta) \, p(\bm\theta \mid \mathcal{D}) \diff\bm\theta \,, \label{eq:infer-question-posterior-expect}
\end{equation}

where $g(\bm{\theta})$ is the function encoding some inferential question (e.g., where in the model configuration space the posterior concentrates).

Predictions can be made in the form of a posterior predictive distribution:

\begin{equation}
p(y_* \mid \mathbf{x}_*, \mathcal{D}) = \int p(y_* \mid \bm\theta, \mathbf{x}_*) \, p(\bm\theta \mid \mathcal{D}) \diff\bm\theta \,, \label{eq:posterior-pred}
\end{equation}

where $y_*$ is the predictions based on the training set $\mathcal{D}$ for a test input $\mathbf{x}_*$. Essentially this is integrating the prediction $p(y_* \mid \bm\theta, \mathbf{x}_*)$ over the posterior distribution of parameters gpml. Note that by giving the final results in terms of a probability distribution, richer information and more reliable inferences are accessed compared to merely giving a point estimate through MLE or MAP (as some machine learning models do under the frequentist framework). This is achieved by incorporating into the inference process the uncertainty in the posterior parameter estimate. Other benefits include posterior predictive checks, which are conducted by checking for auto-consistency between generated data ($\mathbf{y}_*$) and observed data ($\mathbf{y}$).

§.§ Markov Chain Monte Carlo

Many of the integration problems central to Bayesian statistics, including those in <Ref>, are analytically intractable. A class of sampling algorithms, known as Markov chain Monte Carlo (MCMC), can be applied to approximate these Andrieu2003AnLearning. Suppose for some function of interest $f(x)$, the objective is to obtain its integral, with respect to a non-standard target distribution $p(x)$ from which samples cannot be drawn directly:

\begin{equation}
I(f) = \int f(x) \, p(x) \diff x \,.
\end{equation}

By constructing Markov chains that have $p(x)$ as the invariant distribution, MCMC samplers, while traversing the sample space $\mathcal{X}$, are able to generate samples $x^{(i)}$ that mimic samples drawn directly from the target distribution $p(x)$. In other words, this mechanism makes it possible to draw a set of samples $\{ x^{(i)} \}^{N}_{i=1}$ from $p(x)$.

Then, by the Monte Carlo principle, the integral $I(f)$ can be approximated with a sum $I_N(f)$:

\begin{equation}
I_N(f) = \dfrac{1}{N} \sum\limits_{i=1}^N f(x^{(i)}) \xrightarrow[N\xrightarrow{}\infty]{\text{a.s.}} I(f) = \int f(x) \, p(x) \diff x \,.
\end{equation}

That is, the estimate $I_N(f)$ is unbiased and by the strong law of large numbers, it will converge almost surely (a.s.) to $I(f)$. That's why MCMC is a powerful tool in Bayesian analysis. In practice, the Metropolis-Hastings (MH) algorithm and Gibbs sampling have been popular MCMC methods Andrieu2003AnLearning, but only when the parameter space is not too high-dimensional mcelreath2020rethink.

Due to limited computing resources, it is impossible to run Markov chains infinitely long. In other words, inference has to be made based on finitely many draws. One approach, which is effectively leveraged in this research, is to run multiple chains in parallel and monitor various statistics for diagnosing non-convergence. Besides the effective sample size per transition of the Markov chain, the Gelman-Rubin statistic Gelman1992InferenceSequences, denoted by $\hat{R}$, is used in this dissertation. The $\hat{R}$ statistic quantifies whether the ensemble of Markov chains initialized from diffuse points in parameter space finally converge to the same equilibrium phase Betancourt2017RobustRStan. When $\hat{R}$ is sufficiently close to $1$ (for example $\hat{R}<1.05$), convergence is declared to be achieved. As an example, <ref> presents how four chains are started in different corners but approach stationarity and convergence after a certain number of iterations.

[!htbp]mcmc_converge_demofigures/mcmc_converge_demoThe evolution of four random walk Metropolis Markov chains carpenter-prob-stats, each started in a different location. The target density is a bivariate normal with unit variance and correlation 0.9. After $M=5000$ iterations, the four chains have mixed well and explored most of the target density.

For many of the problems in practice, including the models in this dissertation, the parameter space is very high-dimensional and involves highly curving regions. The Metropolis-Hastings algorithm and Gibbs sampling are far from efficient in these situations. Hamiltonian Monte Carlo (HMC), originally proposed by Duane1987HybridCarlo, really outshines the other algorithms at this point and is the main sampling strategy adopted in this dissertation. Specifically, No-U-Turn Sampler (NUTS) introduced by NUTS, which is an extension to HMC, is employed for sampling from posterior distributions.

§.§ Machine Learning

Machine learning was defined by mitchellMLbook as computers improving automatically through experience. It can also be viewed as a function estimation problem Vapnik2000TheTheory, or as the process of extracting important patterns and trends from data hastie2009elements.

In terms of tasks, common types of learning consist of supervised, unsupervised, semi-supervised, and reinforcement burkov2019hundred_page_ml_book. Let $\mathbf{x}_i \in \mathcal{X} \subseteq \mathds{R}^d$ represent input, and $y_i \in \mathcal{Y}$ represent target, then the goals of the first two types are:

* Supervised learning aims to use the dataset, consisting of $\mathbf{X}=\{ \mathbf{x}_i \}_{i=1}^{n}$ and $\mathbf{y}=\{ y_i \}_{i=1}^{n}$, to produce a model that is able to predict an output ($y_j$) given some new/unseen input ($\mathbf{x}_j$), i.e., learning the underlying mapping $f\colon \mathcal{X} \to \mathcal{Y}$.

* Unsupervised learning is used to find the hidden patterns in $\mathbf{X}$; in this case there does not exist any labels ($\mathbf{y}$) or predefined targets.

Another variation of learning is online learning, in which case training data is fed to the algorithm continuously or one example at a time aml_book. In other words, streaming data is available that the algorithm has to process on the run. This is different from batch learning, where data is provided beforehand and “frozen” during the learning process. Online learning can be applied to the different tasks as discussed above (supervised and others).

In terms of model characteristics, machine learning models can be categorized into parametric and nonparametric models. Parametric models are characterized by a fixed number of parameters, whereas nonparametric models have an infinite-dimensional parameter space. For example, in the latter case the parameter space can be the set of continuous functions in a regression setting Orbanz2010. In this dissertation, supervised and unsupervised learning are leveraged while exploiting both parametric and nonparametric models.

From Bayesian's perspective, machine learning is essentially computing the posterior deFreitas2013UBCLearning, which is then used for inference and prediction tasks. This is conducted exactly through <Ref>. In practice, machine learning conducted under Bayesian's framework follows a principled workflow (<ref>), which is adapted for the modeling in this dissertation.

[!htbp]bayes_workflowfigures/workflow_flowchartThe flowchart adapted from Betancourt2020TowardsWorkflow shows a principled Bayesian workflow.

§.§ Analytics Toolset

For the past five to ten years, prosperity in contributions and progress in the open source community has been witnessed. Ecosystems around Python, R, and Julia have been prototyped, tested, and deployed in production environments in various industries. Powerful probabilistic programming languages (PPL), for example Stan mc_stan and PyMC3 Salvatier2016ProbabilisticPyMC3, have become the workhorse for Bayesian machine learning.

The majority of this work is implemented in Python. Specifically, Bayesian learning is performed by leveraging PyMC3. Some analytic visualizations are produced employing ArviZ Kumar2019ArviZPython. Geospatial operations are performed with the help of GeoPandas geopandas. Satellite imagery is processed and analyzed in MATLAB, with implementations mainly following Elvidge2013VIIRSNight.

CHAPTER: DATA PREPROCESSING AND EXPLORATORY DATA ANALYSIS

In this chapter, an overview of the flaring data is given. Some other variables which might be correlated with the flaring statistics are also considered. Exploratory data analysis is performed for choosing the subset of the variables as the focus in this dissertation. A state level model is developed in the end which motivates the work in the next two chapters.

§.§ Data Gathering

Four sources of data, L8 satellite images, VIIRS estimated flared volumes, NDIC monthly production reports, and county/oilfield shapefiles for North Dakota were gathered for the analysis used in this research.

§.§.§ Landsat 8 Images

In total, 167 images (since 2013) were downloaded from Google Cloud using the criteria below:

* From five Path/Row's: 126/216, 126/217, 126/218, 127/216, and 127/217.

According to the Worldwide Reference System (WRS), the satellite imagery of any portion of the world can be queried using Path and Row numbers. These five Path/Row's cover the majority of the areas in North Dakota that have production and flaring activities.

* Nighttime images.

Only nocturnal Landsat 8 imagery are used for the purpose of flare detection.

* Cloud cover less than 10.

Images with low cloud cover percentages reveal more clearly land features including gas flares, and thus are ideal for validating the developed methodologies.

* GeoTIFF Data Product.

Both the georeferencing information and the raw images of all the spectral bands are preserved through the GeoTIFF format, which are necessary for the analysis.

§.§.§ VIIRS Estimated Volumes

The VIIRS flare inventory and estimated volume dataset obtained from Mikhail N. Zhizhin (personal communication) are used in this dissertation. This dataset includes monthly flare detection records in North America from March 2012 to December 2018 (both inclusive) with their associated:

* Timestamps giving the specific month

* Latitudes and longitudes in WGS 84 coordinates

* Flared volume estimations in

§.§.§ NDIC Monthly Production Reports

All the monthly production reports from May 2015 to April 2020 (both inclusive) which have flaring information have been downloaded from NDIC. There is one Excel spreadsheet per month; each row corresponds to a well (that was active in that month), and columns are for various types of information, including flared gas volume (estimated and reported by operator), oilfield, oil production, etc. A screenshot of the top ${\sim}50$ rows in one of the spreadsheets is displayed in <ref>.

A screenshot of the top ${\sim}50$ rows in the October 2018 production report. Each row corresponds to a well. There are in total 17135 rows in this spreadsheet, with the first row being the header.

§.§.§ NDIC Shapefiles

The shapefiles for the counties and oilfields in North Dakota are downloaded from the NDIC GIS Map Server. All the polygons are described in NAD 27 coordinates. The shapefiles are for reverse geocoding the satellite detection locations to readable addresses, specifically which county and oilfield is a flare located in.

§.§ Satellite Image Processing

As discussed in Section <ref>, all the available L8 images have been downloaded. They are processed in batch, following the workflow as outlined in Section <ref>. To compare and contrast with VIIRS' performance, specifically the nighttime combustion source detection limits, all the flares detected from all of the L8 images are gathered and used to generate the source area versus temperature scattergram shown in <ref>.

[VIIRS performance Elvidge2019ExtendingData]

[L8 performance; figure provided by Mikhail N. Zhizhin (personal communication)]

The nighttime combustion source detection limits of VIIRS (top) and L8 (bottom). For natural gas flaring whose temperature is generally greater than 1500, L8 detected flares show source areas (around d-2□) orders of magnitude less than that of VIIRS (around d0□).

Although it is expected that L8 would pick up smaller flares than VIIRS (which is capable of detecting flares around the size of a whole cooktop area), the majority of the detections as indicated on the scattergram are too small for natural gas flaring. To verify if some hot pixels are clustered together and actually representing a single flare or flaring site, HDBSCAN hdbscan with an implementation due to McInnes2017 is executed on every L8 detection map to find out if large blobs of hot pixels are present. HDBSCAN is a density-based clustering algorithm which keeps all the advantages of the original DBSCAN dbscan, for example the capacity of finding clusters of arbitrary shapes. It also outperforms DBSCAN by being able to build clusters of varying density burkov2019hundred_page_ml_book. Further, to get the most accurate results in this case, haversine metric is chosen to handle the great-circle distances between the hot pixels; leaf clustering is used instead of the default Excess of Mass method to produce more fine grained clusters. The clustering results are illustrated in <ref>.
[!htbp]L8_cluster_count_pltfigures/L8_cluster_count_pltA count plot showing the distribution of cluster sizes: clearly there are a certain number of large clusters (as shown by the tail to the right). For example, there exists 2 clusters each of which contains 120 hot pixels and there is one cluster with 84 hot pixels.

To verify whether these clusters are really single flares or they are actually a large number of neighboring wells (in which case each hot pixel still represents an individual flare), they are tracked down by looking further into each detection map (KMZ file). It is found that some large blobs of hot pixels are clustered and indeed represent single (huge) flares. One of the examples is shown in <ref>. This poses a challenge to situations where an accurate estimate of the flare count is needed.

[Band 6 (SWIR)]

[KMZ view]

A large flare consisting of many hot pixels (detections), which is found by running the nightfire algorithm on L8 images. Both the Band 6 (grayscale image) and the KMZ view are shown and provided by Christopher D. Elvidge (personal communication).

The reason for this processing artifact is that, for large flares, there is glow surrounding the flare that was treated as many individual combustion sources. There are potential approaches to mitigate this to make the interpretation and estimation out of L8 more accurate. In this work, the flares detected from VIIRS and the gas volumes estimated out of those are the focus for analytics.

§.§ Reverse Geocoding

By reverse geocoding, the county information of every VIIRS flare that is in North Dakota can be retrieved. For most of the flares, the oilfield information is also retrievable. Thereafter, the flaring statistics from VIIRS and NDIC can be compared and contrasted at different levels, for a certain point or period of time.

Shapefiles as discussed in Section <ref> are used. With the help of GeoPandas, the procedures for extracting counties and oilfields are the same:

* Read the VIIRS records into a geospatial data object, with their original coordinates in WGS 84.

* Read the shapefile into a geospatial data object, with its original coordinates in NAD 27.

* Transform all the geometries in the shapefile to WGS 84 coordinates.

* Perform a spatial join of the two data objects to get the county or oilfield information for each flare, if a specific county/oilfield's polygon and the flare intersect, i.e., having any boundary or interior point in common.

§.§ Correlational Analysis

To study the correlations between oil/gas prices, flaring statistics, and production performance, various time series are extracted for May 2015 to December 2018 (both inclusive). The below list describes all the variables used with their associated labels:
[ ]NDIC flaring well count

VIIRS flared vol monthly flared gas volume from VIIRS

NDIC flared vol monthly flared gas volume from NDIC

WTI oil price WTI crude oil price given by EIA2020PetroleumPrices

Henry Hub gas price Henry Hub natural gas price given by EIA2020NaturalPrices

NDIC oil prod monthly oil production from NDIC

NDIC gas prod monthly gas production from NDIC

VIIRS flare count monthly flare detections count from VIIRS

NDIC flaring well count monthly wells count which conduct flaring from NDIC

NDIC GOR ratio of the NDIC gas production to the NDIC oil production

First, the monthly observations are extracted from each time series, and Spearman's $\rho$ is employed to measure the statistical dependence between the variables. Spearman's $\rho$ is a rank correlation, which quantifies the correlation between the rankings of two variables. Compared to Pearson's $r$, it assesses monotonic relationships which can be nonlinear and is more robust to outliers, therefore is used in this section. The pairwise correlations between the variables are presented in <ref>. Since a correlation matrix is always symmetric with unit diagonals, only the lower triangular part without the diagonal is plotted to minimize the information redundancy.
[!htbp]oil_pr_corrfigures/oil_pr_corrA heat map showing the pairwise Spearman correlations between the original time series' monthly observations. The values are annotated in each cell, the corresponding variables of which can be obtained by reading off the tick labels from the vertical and horizontal axes.

It can be observed that most pairs show positive correlations. Financial factors (i.e., the oil and gas prices) are not among any of the highly correlated pairs (e.g., above 0.80). Nevertheless, it is indicated that the NDIC and VIIRS reportings have a positive correlation, and oil production is positively correlated with flared gas volume.

In this analysis, due to the nature of the procedure (i.e., extract the monthly data and then measure the rank correlations), all the information on the time scale is neglected. To explore the correlations in the context of time series, the first differences (i.e., lag-$1$ differences) are taken for each variable

\begin{equation}
y'_t = y_t - y_{t-1},
\end{equation}

and then pairwise Spearman's $\rho$ is evaluated and visualized in <ref>. In this case, there aren't many pairs of variables which are highly correlated, except the oil and gas production are shown to be monotonically related on the lag-$1$ differences, which is unsurprising. In the remainder of this dissertation, the focus is put on flaring and production related statistics instead of the financial factors.
[!htbp]oil_pr_corr_lag1figures/oil_pr_corr_lag1A heat map showing the pairwise Spearman correlations between the time series after applying the first differences. The values are annotated in each cell, the corresponding variables of which can be obtained by reading off the tick labels from the vertical and horizontal axes.

§.§ State Level Flaring Model

In this section, a regression model is built for the purpose of investigating the statistical relationships between the NDIC and VIIRS reportings. Data from both sources are visualized in <ref>, which demonstrate a positive correlation.
[!htbp]eda_state_data_visfigures/eda_state_data_visVisualizations of both the NDIC and VIIRS reportings. Left figure shows the time series. Right figure presents the scatterplot using the data points of each month.

Assuming a Gaussian observation model for the NDIC reporting with the location parameter encoding VIIRS' information, the model is specified through eda_state_begineda_state_end:

\begin{align}
\alpha &\sim \operatorname{Half-Normal}(0.2) \label{eda_state_begin} \\[0.5ex]
\beta &\sim \operatorname{Gamma}(2, 2) \\[0.5ex]
\sigma &\sim \operatorname{Half-Cauchy}(0.1) \\[0.5ex]
\mu_i &= \alpha + \beta \times \mathrm{VIIRS}_i \label{eq:state_lin_mod}\\[0.5ex]
\mathrm{NDIC}_i &\sim \mathcal{N}(\mu_i, \sigma) \label{eda_state_end}
\end{align}

where $\alpha$ is the intercept and $\beta$ is the slope, both of which are constrained to be non-negative based on the nature of flaring volume; $\sigma$ is the standard deviation in the Gaussian likelihood function, which has to be non-negative as well; $\mu_i$ is the expected NDIC reporting of month $i$, while $\mathrm{NDIC}_i$ and $\mathrm{VIIRS}_i$ are the observed data (i.e., reported volumes) from NDIC and VIIRS in month $i$, respectively. The notation used in defining this model communicates the data generating process unambiguously and is adopted throughout this dissertation. Priors and hyperpriors are on the top while the observation model is at the bottom. The prior distributions for this model and all the others in this dissertation are chosen following the principles below:

* Prefer weakly informative priors, i.e., choose the priors based on the domain expertise at hand before observing any data. They should be strong enough to reflect the domain expertise and be weak enough to “let the data speak”, i.e., let the likelihood dominate when there is a decent amount of data. For example, a prior of a gamma distribution with mean $\mathds{E}\beta=2/2=1$ is placed on $\beta$, reflecting the assumption that the satellite interpretation workflow gives the same flared volume as the NDIC reporting, before one observes any data.

* Prefer priors with soft constraints as opposed to hard constraints, i.e., follow Cromwell's rule. For example, $\alpha$, $\beta$ and $\sigma$ all have prior distributions with support on $\mathds{R}_{>0}$ or $\mathds{R}_{\geq 0}$. Counterexamples include using a triangular distribution or a continuous uniform distribution as the prior for such quantities, for which the author does not recommend.

* Prefer maximum entropy distributions, i.e., make the most conservative assumptions based on all the information at hand (obeying all the known constraints). For example, the Gaussian and the binomial distributions are maximum entropy distributions and used in this dissertation, the fact of which can be formally shown leveraging the definition of Kullback–Leibler (KL) divergence.

Once the priors and likelihood are established, four Markov chains of Hamiltonian Monte Carlo are run in parallel to sample from the posterior. The parameter estimates are reported in <ref>, and the posterior distributions and trace plots are presented in <ref>. The four chains are plotted separately with different colors. The $x$-axis of the trace plot shows the number of iterations. This layout is used consistently for the remainder of this dissertation.

Parameter Estimates of State Level Flaring Model

Parameter     Variable     Point Estimate     90 Credible Interval
$\alpha$     Intercept     0.061     (0.044, 0.079)

$\beta$     Slope     0.535     (0.482, 0.590)

$\sigma$     Reporting variability     0.030     (0.024, 0.035)

[!htbp]eda_state_tracepltfigures/eda_state_tracepltPosterior distributions (left column) and trace plots (right column) for the state level flaring model. Well mixing and convergence of the Markov chains have been achieved as shown by the trace plots.

Utilizing the model and the trace, posterior predictive samples are generated to construct the intervals (<ref>). Point estimates and point predictions are easy to obtain for a certain machine learning model, however it is the properly constructed intervals that will provide insights into the uncertainty for decision making. The author would like to emphasize the importance of quantifying uncertainties when using machine learning, no matter for inference, prediction, or building intermediate models for integration into physics-based models. This is unfortunately neglected or ignored in some of the applications/publications in the petroleum engineering domain. The importance of properly quantifying the uncertainties will also be stressed in the following chapters.

[!htbp]eda_state_ppcfigures/eda_state_ppcIntervals are constructed using posterior predictive samples. In both figures, the line shows the “best” fit using point estimates (posterior means) of $\alpha$ and $\beta$. Shaded area in the left figure presents the 90 credible interval (CI) of the regression mean. Shaded area in the right figure demonstrates the 90 prediction interval for the future NDIC reporting, for which most of the existing observations fall within.

Whenever only one model specification is needed for making point predictions, it can be recovered by the parameter estimates from <ref>:

\begin{equation}\label{eq:state_point_est_mod}
\mathrm{NDIC}_i = 0.061 + 0.535 \times \mathrm{VIIRS}_i \, ,
\end{equation}

where $\mathrm{NDIC}_i$ and $\mathrm{VIIRS}_i$ are flared volumes in bcm of month $i$. The model also provides clear interpretations for the NDIC reporting regression mean, on the whole state level:

* The intercept indicates on average there is 90 probability that 0.040.08 reported volume per month will not be captured by the current VIIRS processing workflow. The posterior mean is 0.061 ($\approx \SI{2150}{\mmcf}$).

* The slope indicates on average when satellite estimated volume increases by one unit, under 90 probability the NDIC reporting will increase by 0.48 unit to 0.59 unit. The posterior mean is 0.535 unit.

This model, while serving as a decent calibration and estimation tool for NDIC reporting on the state level, makes the assumption that the heterogeneity within the state (e.g., among different counties) is negligible and all the monthly observations are conditionally independent and identically distributed (i.i.d.). For the scenarios in which these assumptions do not hold, other types of models can be built and are discussed in Chapter <ref> and Chapter <ref>.

CHAPTER: COUNTY LEVEL FLARING MODEL

“Multilevel regression deserves to be the default form of regression.”— mcelreath2016rethink

§.§ Learning the Heterogeneity

In this chapter, the author explores the heterogeneity in correlations between the state-reported and satellite-detected flaring statistics, among different counties in North Dakota. The motivations are threefold:

* Provide more granular insights than merely investigating the whole state's flaring statistics.

* Compare and contrast different counties' reporting consistencies with the baseline (i.e., the satellite detections).

* Develop a dedicated model for each county for calibration and prediction purposes.

§.§ Hierarchical Model

A common problem in learning from data is modeling individuals or units of a population. For example, building models for different counties in a state, or for different well pads in an oilfield. Usually from domain expertise, it is expected that the units would demonstrate some differences, however they do not necessarily represent completely independent data generating processes. In other words, the units are different in some ways, while being similar in others. Unfortunately, the following two common modeling approaches are extreme and not ideal:

* Complete pooling

* This ignores heterogeneity and assumes that the observations from all the units are generated/described by the exact same process. One set of parameters is learned for the whole population. In this situation, the variance might be smaller, however the bias could be huge.

* No pooling

* This lets each unit learn its own set of parameters from its own data. The assumption is that the information from each unit tells one nothing about any other unit. In this situation, the bias might be smaller, however the variance could be huge.

In practice, neither of these approaches will be able to generalize well for insight extraction or prediction tasks, due to the total generalization error being large. In fact, these two extremes can be compromised by explicitly modeling the entire population of units. That is, in order to investigate the correlations among the individual units, an explicit model is introduced for the population. In the learning phase, the individual posteriors are used to fit some population distribution, while the information of the population is then fed back to the individuals. What happens in this case is that the individuals with diffuse likelihood functions (e.g. with less data) are dragged more towards the population distribution, whereas the individuals which are well informed by their data will have their posteriors mostly unchanged. In this process, dynamic regularization is achieved, i.e., the total generalization error is much smaller by partially pooling the data and balancing between the bias and variance.

In the context of county level model development, the question is now how might one model the population. To motivate the choice of a particular class of models, some characteristics of the counties have to be examined. In this work, the counties are considered to be exchangeable, i.e., the joint probability $p(\bm{\theta}_1, \bm{\theta}_2, \dots, \bm{\theta}_n)$ is invariant to permutation of the indices, where $\bm{\theta}_i$, $i=1,2,\dots,n$ is the parameters for the $i$-th county. That is, for any permutation $\pi$,

\begin{equation}
p(\bm{\theta}_1, \bm{\theta}_2, \dots, \bm{\theta}_n) = p(\bm{\theta}_{\pi_1}, \bm{\theta}_{\pi_2}, \dots, \bm{\theta}_{\pi_n}) \, .
\end{equation}

Furthermore, the list of counties can grow, i.e., although one might only look at a few counties at this point, in the future new counties in terms of flaring activities might be considered.

If a population being modeled is exchangeable, and the population can grow arbitrarily large, de Finetti's theorem shows that the only distribution that respects exchangeability is a hierarchical distribution:

\begin{equation}
p(\bm{\theta}_1, \bm{\theta}_2, \dots, \bm{\theta}_n) = \int \Biggl[ \prod_{i=1}^{n} p(\bm{\theta}_i \mid \phi) \Biggr] p(\phi) \diff\phi \, ,
\end{equation}

where $\phi$ is a population parameter (which can be generalized to multiple population parameters) and $p(\phi)$ is a population prior. It asserts an important fact that if exchangeable data is used for analytics, there must exist a population model Jordan2010UCBTheorem. This provides guidance for the development of the county level flaring models in this chapter.

Equivalently, the individual and population parameters can be fitted jointly, achieving a dynamic pooling of the data:

\begin{equation}\label{eq:de_finetti_joint_mod}
p(\bm{\theta}_1, \bm{\theta}_2, \dots, \bm{\theta}_n, \phi) = \Biggl[ \prod_{i=1}^{n} p(\bm{\theta}_i \mid \phi) \Biggr] p(\phi) \, ,
\end{equation}

in which process not only the $\bm{\theta}$'s but also $\phi$ are learned. After adding the observations component ($\mathcal{D}=\{(\mathbf{x}_j, y_j) \mid j=1, \ldots , m\}$) to it, the joint model becomes:

\begin{equation}\label{eq:hier_gen_form}
p\!\left( \left\{ y_j, \mathbf{x}_j, \bm{\theta}_{\mathsmaller{\textsc{county}} [j]}, \psi_j \right\}_{j=1}^{m}, \phi \right) = \Biggl[ \prod_{j=1}^{m} p(y_j \mid \mathbf{x}_j, \bm{\theta}_{\mathsmaller{\textsc{county}} [j]}, \psi_j) \, p(\bm{\theta}_{\mathsmaller{\textsc{county}} [j]} \mid \phi) \Biggr] p(\phi) \, ,
\end{equation}

where $\bm{\theta}_{\mathsmaller{\textsc{county}} [j]}$ stands for the parameters for the $j$-th observation based on its county assignment, and $\psi$ are some other parameters in the likelihood function that are not necessarily distributed according to a population model. <Ref> characterizes a hierarchical model that fits nicely into the Bayesian framework and is exploited for building the models in this chapter.

As a fundamental approach to model heterogeneity, hierarchical models have been depended upon routinely in various fields including ecological science bolker2008ecological, political science gelman2006Hierarchical, and biological science mcelreath2016rethink. The author believes that they should be widely accepted and utilized in the petroleum engineering domain as well, where the dataset is usually presented in hierarchies. For example, the shale gas wells in a given basin were completed by different oilfield service companies. The information can then be pooled among the service companies. A further discussion is given in Section <ref>. One caveat, though, is that de Finetti's theorem is based on the assumption that the population (of units) is exchangeable and can grow arbitrarily large. Just like every other assumption in machine learning, it should not be taken for granted and does not always hold. In the context of county level flaring model development, one might argue that there are currently 53 counties in North Dakota and there might not be many new counties (as administrative divisions) in any finite amount of time. In that regard, the author agrees with the claim of Box2009StatisticalAdjustment that, since assumptions “are never exactly true”, what shall be sought is the useful models as opposed to the correct ones. That is the goal for applying the hierarchical models in this chapter.

It is worth noting that the terminologies are not consistent when referring to these types of models: some argue that hierarchical model and multilevel model are different names for the same modeling technique bolker2008ecological,mcelreath2016rethink, while others tried to differentiate them Carpenter2019AllModeling. In this dissertation, the model assumptions are communicated via the mathematical structures instead of the terminologies, by writing out the full model definitions whenever possible.

§.§ Data Description

After performing the reverse geocoding as outlined in Section <ref>, there are twelve counties found to have reported flaring activities from both VIIRS and NDIC. For each county's historical data from May 2015 to December 2018 (both inclusive), only the months that have reported volumes from both sources are extracted. A scatterplot for each of the 12 counties is presented in <ref>, where the county abbreviations follow the convention from the NDIC monthly production reports. <ref> lists the full county names associated with each abbreviation.
[!htbp]county_scatterfigures/county_scatterScatterplots of NDIC and VIIRS reportings for different counties. Both the $x$- and $y$-axis are shared among all the subplots. The $x$-axis is the monthly VIIRS reporting of the flared volume in , and the $y$-axis is for the NDIC reporting in the same unit.

North Dakota County Abbreviations

Abbreviation     County
MCK              McKenzie County

DUN              Dunn County

WIL              Williams County

MTL              Mountrail County

BOW              Bowman County

DIV              Divide County

BRK              Burke County

MCL              McLean County

BIL              Billings County

STK              Stark County

SLP              Slope County

GV               Golden Valley County

It can be seen that the flaring magnitudes in terms of the flared volumes are quite diverse for the different counties. To better visualize all of them, a zoomed-in view for each county is shown in <ref>. It becomes clear that most of the counties except SLP and GV have more than ${\sim}12$ data points; however, only the four counties in the top row (i.e., MCK, DUN, WIL and MTL) have the largest amount of data and indicate stronger positive correlations between VIIRS and NDIC.
[!htbp]county_scatter_ALLzoomedfigures/county_scatter_ALLzoomedScatterplots of NDIC and VIIRS reportings for different counties, without sharing neither $x$- nor $y$-axis for all the subplots. Within each subplot, equal scaling and limits are set for $x$- and $y$-axis. The axes' meanings are the same as in <ref>.

For the purpose of building county level models and investigating the heterogeneity among the counties, the no pooling option discussed in the previous section will fail. Especially with counties SLP (which has 3 observations) and GV (which has 2 observations), if a linear model such as <Ref> is fitted, the learned slope parameters $\beta_{\mathsmaller{\textsc{county}}}$ will have point estimates $\hat{\beta}_{\textsc{slp}} \approx 0$ and $\hat{\beta}_{\textsc{gv}} \gg 0$ with their associated samples. The interpretation of the slope parameter (which was discussed right after <Ref>) implies that such inferences are never possible. Some other counties, even with more data points (e.g., MCL), suffer from the noise levels in their observations. Using their own dataset will frustrate accurate inferences. Therefore, in order to build models robustly at a county level, the hierarchical model discussed in the previous section is exploited.

§.§ Model Specification

Motivated by the discussions in Section <ref>, partial pooling is performed by explicitly modeling the entire population of counties. In this way, the counties such as MCL can leverage the information from other counties to learn their own parameters. Counties with “strong data” (i.e., very informative data which makes the likelihood dominate the structure of the posterior), such as those in the top row of <ref>, indicate a positive correlation between VIIRS and NDIC. Therefore, a similar strategy as in Model <ref> is adopted for the counties, i.e., one set of slope and intercept is learned for each county.

Since the slope and intercept are very interpretable, the meanings of which were discussed right after <Ref>, partial pooling is also enabled across parameter types (i.e., intercepts and slopes). In other words, knowing how much flared volume is missed from VIIRS (i.e., the information carried by the intercept) might improve learning how VIIRS and NDIC will covary (i.e., the information carried by the slope). Specifically, a population model with a multivariate normal density is used for the different counties' parameters.

The hierarchical model is specified through county_cp_begincounty_cp_end:

\begin{align}
\mu_\alpha &\sim \operatorname{Half-Normal}(0.1) \label{county_cp_begin} \\[0.5ex]
\mu_\beta &\sim \operatorname{Gamma}(2, 2) \\[0.5ex]
\sigma_\alpha &\sim \operatorname{Half-Normal}(0.1) \\[0.5ex]
\sigma_\beta &\sim \operatorname{Half-Normal}(0.1) \\[0.5ex]
\sigma &\sim \operatorname{Half-Normal}(0.05) \\[0.5ex]
\mathbf{R} &\sim \operatorname{LKJcorr}(2) \label{exp:cp_core_beg} \\[0.5ex]
\bm{\Sigma} &= \begin{pmatrix}
\sigma_\alpha & 0 \\
0 & \sigma_\beta
\end{pmatrix} \cdot \mathbf{R} \cdot \begin{pmatrix}
\sigma_\alpha & 0 \\
0 & \sigma_\beta
\end{pmatrix} \\[0.5ex]
\begin{bmatrix} \alpha_{\mathsmaller{\textsc{county}}} \\ \beta_{\mathsmaller{\textsc{county}}} \end{bmatrix} &\sim \operatorname{MVNormal} \!\left( \begin{bmatrix} \mu_\alpha \\ \mu_\beta \end{bmatrix}, \bm{\Sigma} \right) \label{exp:mvn_cp} \\[0.5ex]
\mu_j &= \alpha_{\mathsmaller{\textsc{county}} [j]} + \beta_{\mathsmaller{\textsc{county}} [j]} \times \mathrm{VIIRS}_j \\[0.5ex]
\mathrm{NDIC}_j &\sim \mathcal{N}(\mu_j, \sigma) \label{county_cp_end}
\end{align}

μ_α is the average intercept for all the counties;

μ_β is the average slope for all the counties;

σ_α is the standard deviation among different counties' intercepts;

σ_β is the standard deviation among different counties' slopes;

σ is the the standard deviation in NDIC reporting within the counties;

𝐑 is the correlation matrix distributed according to an LKJ distribution. It is 2-by-2 in size and encodes the correlation between the intercepts and slopes;

Σ is the covariance matrix for the population model, which is constructed by multiplying the correlation matrix from both sides by a diagonal matrix of standard deviations;

α_county and β_county are the intercept and slope for each county, whose prior distributions are defined by a two-dimensional Gaussian population model;

countyj  (in the subscript) denotes the county index, i.e., $\footnotesize{\textsc{county}}\lbrack j \rbrack \in \{k\in \mathds{N}_{0} \mid k\leq 11\}$, such that $\alpha_{\mathsmaller{\textsc{county}}[j]}\text{ and }\beta_{\mathsmaller{\textsc{county}}[j]}$ are the intercept and slope for the $j$-th observation based on its county assignment;

VIIRS_j is the VIIRS reported volume of the $j$-th observation;

μ_j denotes the underlying flared volume of the $j$-th observation;

NDIC_j is the NDIC reported volume of the $j$-th observation.

The LKJ distribution due to maxnames99Lewandowski2009GeneratingMethod is a distribution over positive-definite symmetric matrices with unit diagonals, i.e., correlation matrices. In the model specification above, it directly influences the prior for the covariance matrix. Before it was introduced and when HMC was not widely applicable, the usual choices for modeling covariance matrices were Wishart or inverse-Wishart distributions, due to their nice conjugacy properties. However, LKJ is better suited for modern Bayesian computational settings Betancourt2015CommentArguments,lambert2018bayes and therefore employed in this work.

LKJ has a single parameter $\eta$, which can be interpreted as the shape parameter of a symmetric beta distribution bda3. As $\eta$ gets larger, the prior is more skeptical of large correlations in the matrix, i.e., providing regularizing effects. The probability density of LKJ with a few $\eta$ values are displayed in <ref>. In this work, $\operatorname{LKJcorr}(\eta=2)$ is chosen to define a weakly informative and regularizing prior.
[!htbp]lkj_densityfigures/lkj_density0.425$\operatorname{LKJcorr}(\eta=\mathtt{eta})$ probability density. As $\eta$ increases, larger correlations become less plausible.

Model <ref>, while being expressive in the data generating process, is a centered parameterization of the hierarchical structure Papaspiliopoulos2007AModels. In this parameterization, the hierarchical parameters (such as $\beta_{\mathsmaller{\textsc{county}}}$) and the lower-level parameters in the prior (e.g., $\mu_{\beta}$ and $\sigma_{\beta}$) are tightly coupled, and they are highly correlated in the posterior. Since this model involves complex geometries and interactions in the posterior, HMC is leveraged for sampling. When there is not a lot of data (which is the case for the current NDIC and VIIRS reportings), this parameterization leads to very inefficient sampling and non-convergences stan_man. The noncentered parameterization is preferable in these cases and therefore employed for building the county level models.

§.§ Model Reparameterization

Reparameterization of hierarchical models can be applied to any distribution in the location-scale family, for which the normal distribution is a good candidate. In the case of reparameterizing a multivariate normal prior, suppose the prior for $\bm{\theta}$ is a multivariate normal with mean vector $\bm{\mu}$ and covariance matrix $\bm{\Sigma}$ (such as <ref>), then a noncentered parameterization is given by:

\begin{align}
\widetilde{\bm{\theta}} &\sim \operatorname{MVNormal}(\mathbf{0}_n, \, \mathbf{I}_n) \\[0.5ex]
\bm{\varphi} &= \bm{\mu} + \mathbf{L} \cdot \widetilde{\bm{\theta}}
\end{align}

where $\widetilde{\bm{\theta}}$ has the same dimensions as $\bm{\theta}$ and all of its elements i.i.d. according to $\mathcal{N}(0,1)$, $\mathbf{L}$ satisfies $\mathbf{L} \cdot \mathbf{L}^\top=\bm{\Sigma}$, and $\bm{\varphi}$ recovers the exact same prior distribution for $\bm{\theta}$. This reparameterization leads to more efficient sampling by reducing the dependence between $\bm{\mu}$, $\mathbf{L}$, and $\widetilde{\bm{\theta}}$. One choice for $\mathbf{L}$ is the Cholesky factor of $\bm{\Sigma}$, which provides implementation convenience for the multivariate normal cases stan_man and is adopted in this work.

The noncentered county level model is specified through county_ncp_begincounty_ncp_end, with the reparameterized part (corresponding to Model <ref>) highlighted in blue:

\begin{align}
\mu_\alpha &\sim \operatorname{Half-Normal}(0.1) \label{county_ncp_begin} \\[0.5ex]
\mu_\beta &\sim \operatorname{Gamma}(2, 2) \\[0.5ex]
\sigma_\alpha &\sim \operatorname{Half-Normal}(0.1) \\[0.5ex]
\sigma_\beta &\sim \operatorname{Half-Normal}(0.1) \\[0.5ex]
\sigma &\sim \operatorname{Half-Normal}(0.05) \\[0.5ex]
\color{Blue}\mathbf{L} &\color{Blue}\sim \operatorname{LKJCholeskyCov}\!\left(\eta=2, \begin{bmatrix} \sigma_\alpha & \sigma_\beta \end{bmatrix}^\intercal \right) \label{exp:ncp_core_beg} \\[0.5ex]
\color{Blue}\begin{bmatrix} z_\alpha \\ z_\beta \end{bmatrix} &\color{Blue}\sim \operatorname{MVNormal} \!\left( \begin{bmatrix} 0 \\ 0 \end{bmatrix}, \begin{pmatrix}
1 & 0 \\
0 & 1
\end{pmatrix} \right) \\[0.5ex]
\color{Blue}\begin{bmatrix} \alpha_{\mathsmaller{\textsc{county}}} \\ \beta_{\mathsmaller{\textsc{county}}} \end{bmatrix} &\color{Blue}= \begin{bmatrix} \mu_\alpha \\ \mu_\beta \end{bmatrix} + \mathbf{L} \cdot \begin{bmatrix} z_\alpha \\ z_\beta \end{bmatrix} \label{exp:ncp_core_end} \\[0.5ex]
\mu_j &= \alpha_{\mathsmaller{\textsc{county}} [j]} + \beta_{\mathsmaller{\textsc{county}} [j]} \times \mathrm{VIIRS}_j \\[0.5ex]
\mathrm{NDIC}_j &\sim \mathcal{N}(\mu_j, \sigma) \label{county_ncp_end}
\end{align}

𝐋 is the Cholesky factor of the covariance matrix which has LKJ distributed correlations;

z_α and z_β are the standardized intercept and slope for each county.

The rest of the symbols have the same meaning as in Model <ref>. The noncentered model imposes the exact same probabilistic structure as in Model <ref>, and is implemented for making inference on each county's parameters.

§.§ Model Fitting

Four chains are sampled from the posterior distributions. The posterior distributions and trace plots for the slopes and intercepts are presented in <ref> and <ref>, respectively. Well mixing and convergence have been achieved as shown by the trace plots.

Posterior distributions and trace plots of the slopes for each county.

Posterior distributions and trace plots of the intercepts for each county.

To better compare and contrast the different counties' parameters, the forest plots of 90 highest density intervals (HDI) for the slopes and intercepts are given in <ref> and <ref>, respectively. In both figures, counties are ordered by the VIIRS reported volumes, and those with the least amount of estimated volumes (such as SLP and GV) are at the bottom. The thin lines present the 90 HDI's and the thicker line segments stand for the interquartile ranges (IQR). The points represent the posterior means.
[!htbp]county_slope_forestpltfigures/county_slope_forestplt0.5A forest plot showing the uncertainties around each county's slope estimate. The counties at the bottom have insufficient or noisy datasets, therefore their estimates are largely pulled towards the partially-pooled mean.
[!htbp]county_intercept_forestpltfigures/county_intercept_forestplt0.5A forest plot showing the uncertainties around each county's intercept estimate. The dotted line labels the zero intercept, for which some counties' estimates are not significantly different from.

In the case of the slopes (<ref>), it can be seen the top four counties are quite diverse. MTL has the largest point estimate in the entire population ($\hat{\beta}_{\textsc{mtl}} > 0.6$) while DUN has the smallest one ($\hat{\beta}_{\textsc{dun}} < 0.5$). Furthermore, the HDI's for DUN and MTL rarely overlap, indicating that it is almost certain that MTL has a larger slope than DUN. The counties with fewer observations (remaining eight counties) have greater uncertainties in their parameter estimates, while all of their point estimates are pulled towards the partially-pooled mean which is between 0.5 and 0.6. When there is not enough data for some counties, the hierarchical model strives to reinforce information sharing among different counties, thus providing more sensible results and also quantifying the uncertainties in such processes. From domain expertise, these results make more physical sense than the no-pooling estimates discussed in Section <ref> (i.e., $\hat{\beta}_{\textsc{slp}} \approx 0$ and $\hat{\beta}_{\textsc{gv}} \gg 0$).

In the case of the intercepts (<ref>), there is also heterogeneity among the counties. In particular, by plotting a dotted line labeling the zero intercept, some counties are found to likely have zero intercept (e.g., zero is covered by the IQR or HDI) while others have intercepts that are significantly different from zero. It might not be surprising to get close-to-zero intercepts and greater uncertainties for those counties with less data (such as SLP and GV), however it is interesting to obtain the HDI for MTL that covers zero. Recall that the intercept parameter can be interpreted as the NDIC reported volume which is not captured by VIIRS. This finding for MTL, along with the fact that MTL has the largest slope point estimate (where a larger slope denotes closer proximity to the satellite estimation), convinces the author that MTL used to have persistent and stronger gas flares. They kept VIIRS from missing the flaring events in general, and lead to the reported volumes from NDIC and VIIRS being closer to each other. On the contrary, DUN's smaller slope and larger intercept characterize its flares as sporadic and weaker. One thing worth mentioning is that, with the current interpretation of the intercept, it does not make much physical sense to have negative intercepts. Although every county has positive point estimates for their intercepts, some counties' HDI's show coverage over the negative values. This is a limitation of choosing a 2D Gaussian population model for the intercepts and slopes. Since the 2D Gaussian is supported on $\mathds{R}^2$, in the context of some counties having “weak data”, negative values make an appearance in their HDI's.

The discussions above naturally lead to the question of whether the slopes and intercepts are correlated. It turns out that, by partially pooling the different types of parameters, a probable negative correlation between the slopes and intercepts is revealed (<ref>). The correlation is learned from the heterogeneity in flare characteristics among the counties:

* Persistent flares yield smaller intercepts and larger slopes.

* Sporadic flares yield larger intercepts and smaller slopes.

In other words, intercepts and slopes covary in the entire population of counties. By pooling information across parameter types, what the model learns in the intercept can improve learning about slopes, and vice versa. With this “experience” or “knowledge”, the hierarchical model will be able to quickly update its expectation for any new counties' parameters even with just a few observations in the beginning. It should be noted that there is also some probability mass for the positive correlation values, i.e., the negative correlation is not very strong. This could be due to that some counties do not have a lot of data at this time. The posterior will be updated as more data is brought in.
[!htbp]rho_prior_vs_postfigures/rho_prior_vs_post0.5Correlation between the intercepts and slopes. Blue: Posterior distribution of the correlation, the mode of which is below zero. Dashed: Prior distribution, the $\operatorname{LKJcorr}(2)$ density.

Finally, the parameter estimates are reported in <ref>, from which the parametric model for each county can be recovered, and then deployed in calibration and prediction usage scenarios.

Parameter Estimates of County Level Flaring Model

Parameter                      Variable                        County     Point Estimate     90  CI
12*$\alpha_{\mathsmaller{\textsc{county}}}$     12*Intercept     MCK        0.019                (0.015, 0.023)

DUN        0.008                (0.004, 0.013)

WIL        0.010                (0.007, 0.013)

MTL        0.002                (-0.001, 0.006)

BOW        0.015                (0.013, 0.017)

DIV        0.003                (0.001, 0.005)

BRK        0.003                (0.001, 0.004)

MCL        0.000                (-0.001, 0.002)

BIL        0.001                (-0.001, 0.003)

STK        0.001                (-0.003, 0.004)

SLP        0.001                (-0.005, 0.007)

GV         0.002                (-0.005, 0.009)
black!3012*$\beta_{\mathsmaller{\textsc{county}}}$      12*Slope         MCK        0.519                (0.493, 0.542)

DUN        0.464                (0.385, 0.547)

WIL        0.549                (0.495, 0.605)

MTL        0.623                (0.553, 0.693)

BOW        0.516                (0.370, 0.677)

DIV        0.554                (0.395, 0.719)

BRK        0.556                (0.389, 0.715)

MCL        0.563                (0.391, 0.730)

BIL        0.560                (0.395, 0.727)

STK        0.562                (0.393, 0.729)

SLP        0.561                (0.406, 0.752)

GV         0.560                (0.398, 0.731)

§.§ Model Extensibility

Looking back at the hierarchical model and the reparameterization strategy from the previous sections, there are four potential deployment scenarios that are worth discussing. They demonstrate the extensibility and flexibility of the chosen approach in the context of flaring data analytics:

* New counties are present in terms of the reported flaring statistics from both VIIRS and NDIC.

At this time, there are 12 counties that have reported flaring statistics from both VIIRS and NDIC. If flaring data becomes available for some other counties in the future, the hierarchical model allows the population to be immediately expanded to accommodate the new counties. This can be seen from the conditional structure in <Ref>: by taking a model for $n+1$ counties

\begin{equation}
p(\bm{\theta}_1, \bm{\theta}_2, \dots, \bm{\theta}_n, \bm{\theta}_{n+1},\phi) = \Biggl[ \prod_{i=1}^{n+1} p(\bm{\theta}_i \mid \phi) \Biggr] p(\phi) \, ,
\end{equation}

then pulling out the term for the $(n+1)$-th county from the right-hand side (RHS)

\begin{equation}
p(\bm{\theta}_1, \bm{\theta}_2, \dots, \bm{\theta}_n, \bm{\theta}_{n+1},\phi) = p(\bm{\theta}_{n+1} \mid \phi) \Biggl[ \prod_{i=1}^{n} p(\bm{\theta}_i \mid \phi) \Biggr] p(\phi) \, ,
\end{equation}

it can be recognized that the remaining part on the RHS is the hierarchical model for $n$ counties

\begin{equation}
p(\bm{\theta}_1, \bm{\theta}_2, \dots, \bm{\theta}_n, \bm{\theta}_{n+1},\phi) = p(\bm{\theta}_{n+1} \mid \phi) \, p(\bm{\theta}_1, \bm{\theta}_2, \dots, \bm{\theta}_n, \phi) \, .
\end{equation}

This indicates the newly introduced counties will only depend on the population parameters $\phi$, i.e., how the new counties interact with the existing ones (from the initial dataset) is not explicitly specified but being mediated through $\phi$. This mechanism allows the population (of counties) to be expanded arbitrarily. In practice, without any modification, Model <ref> can be re-fitted with the new dataset as a whole.

* More data are available for those counties which used to have very few observations.

In the event of more data becoming available for those counties with wide HDI's such as SLP and GV, the posteriors will be updated according to that information. Their HDI's would become narrower and narrower as more and more data are available, and since the hierarchical model pools information among the counties, these counties will contribute to updating the population model's and other counties' parameters. Similar to Item 1 above, Model <ref> does not need modifying and can be re-fitted with the new data.

* Sample sizes among counties become more unbalanced.

In general, when there is a lot of data for each county, the centered parameterization (Model <ref>) is more efficient. When the sample size is not large, which is the case for the current VIIRS and NDIC reportings, the noncentered parameterization (Model <ref>) is better. However, the parameterization for hierarchical models is not a monolithic tactic. If the reported flaring data becomes very unbalanced across counties, e.g., some counties have a huge amount of data whereas others have very little data, then each county can be parameterized differently. More specifically,

* For the counties that have strong data such that their likelihood functions dominate, centered parameterization can be applied through exp:cp_core_begexp:mvn_cp.

* For the counties that have weak data such that their prior models dominate, noncentered parameterization can be applied through exp:ncp_core_begexp:ncp_core_end.

All in all, this is still one hierarchical model which defines the exact same probabilistic structure as Model <ref> or Model <ref>, but avoids inefficiencies and non-convergences in the sampling from posteriors.

* Oilfield level heterogeneity needs to be examined.

Under the assumptions that the oilfields in North Dakota are exchangeable and the population of oilfields (which conduct flaring) can grow, the hierarchical model developed in this chapter can be directly applied to investigate the heterogeneity in different oilfields' parameters. Following the reverse geocoding as discussed in Section <ref>, there are 258 oilfields that have both NDIC and VIIRS reportings for the same study period as in this chapter. Some oilfields have very few observations and can benefit from the hierarchical model through pooling information among the entire population of oilfields. Furthermore, due to the number of oilfields being relatively large, the population model could be learned with more ease (because more information is available for the population). In the case of the county level model developed in this chapter, since there are only 12 individuals (counties) in the population, some uncertainties about the population are inevitably present and reflected through the posteriors.

The models developed in this chapter, while capturing the heterogeneity among the different counties in North Dakota, rely on the assumption that all the monthly observations within a certain county are conditionally i.i.d. For situations where the temporal structure has to be taken into consideration, other types of models can be built and are discussed in the next chapter.

CHAPTER: FLARING TIME SERIES ANALYTICS

“Were neural networks over-hyped, or have we underestimated the power of smoothing methods?

I think both these propositions are true.”— mackay_IT

§.§ Learning the Flaring Pattern and Behavior

In this chapter, the author develops a generic framework for revealing flaring patterns and behaviors. The main challenges are fourfold:

* Observed data are noisy.

* Companies estimate the flaring volumes and conduct self-reporting. Satellites could miss some events. However, having knowledge about the underlying process is vital in lots of situations including when the state and local governments need to make key decisions based on the data. In the meantime, understanding the underlying process helps with anomaly detection by differentiating between true anomalies in reporting and ordinary noise or stochasticity.

* A probabilistic approach is desirable to be adopted.

* A set of most probable functions (characterizing the underlying process) are preferable over one single best fit function.

* The observations of a certain entity are time series.

* The temporal structure is intrinsic to the dataset and thus must be harnessed.

* The framework should be generic enough for automated insights extraction.

* There are more than 200 operators and 500 oilfields operating in North Dakota. Choosing a specific parametric form of model (e.g., ARIMA or LSTM) for each entity and then fitting the model to the data is not only time consuming, but also prevents easy integration into automation pipelines (for extracting insights for example).

It is striking that the elegant properties of Gaussian process make it a natural choice to tackle all of these challenges and is therefore employed in this chapter.

§.§ Gaussian Process

A Gaussian process (GP) can be viewed as a distribution over infinite-dimensional Hilbert space of functions. It is formally defined as “a collection of random variables, any finite number of which have a joint Gaussian distribution” gpml. Gaussian processes are extremely powerful nonparametric learning techniques, which provide a composite of flexibility and interpretability. They are well suited to problems which necessitate principled handling of uncertainty and interpretation, in the presence of noisy and dynamic datasets. Such scenarios include smoothing Deisenroth2012RobustProcesses and time series modeling Roberts2013GaussianModellingb. They are also well established in different fields under various names, for example kriging in geostatistics and Kalman filters both correspond to Gaussian processes mackay1998introductionGP.

In this work, the motivation is to develop a generic framework for recognizing the underlying unknown processes $f(\mathbf{x})$ which reflect flaring strategies and behaviors. Thus inference is conducted directly in the function space employing GP as a prior. A Gaussian process is completely specified by its mean function $m(\mathbf{x})$ and covariance function $k(\mathbf{x}, \mathbf{x}')$ Bandyopadhyay2018CSMStatistics, which are defined as:

\begin{align}
m(\mathbf{x}) &= \mathds{E} [f(\mathbf{x})] \,, \\[0.5ex]
k(\mathbf{x}, \mathbf{x}') &= \mathds{E} [ (f(\mathbf{x}) - m(\mathbf{x})) (f(\mathbf{x}') - m(\mathbf{x}'))] \,,
\end{align}

and the function distributed as a Gaussian process is denoted by

\begin{equation}
f(\mathbf{x}) \sim \mathcal{GP}\!\left(m(\mathbf{x}), \, k(\mathbf{x}, \mathbf{x}') \right).
\end{equation}

§.§.§ Mean Function

In this work, the mean functions are always chosen to be zero, since there is no prior knowledge on the mean of the latent processes. In the meantime, for GPs with a zero mean function, the mean of the posterior process is not confined to be zero gpml. All the latent functions modeled with a GP prior in this dissertation follow

\begin{equation}
f(\mathbf{x}) \sim \mathcal{GP}\!\left(0, \, k(\mathbf{x}, \mathbf{x}') \right),
\end{equation}

where $k$ is some covariance function.

§.§.§ Covariance Function

Covariance function, also known as kernel, is the crucial ingredient in a GP, as it encodes one's assumptions about how the function should behave by defining similarity. The fundamental assumption is that data points with inputs $\mathbf{x}$ which are close would have similar target values $y$. This assumption is usually very reasonable in areas including time series modeling, and it is theoretically backed by Tobler's first law of geography. The covariance functions used in this dissertation include:

* The Matérn class of covariance functions, which is given by:

\begin{equation}
k_{\nu }(r)={\frac {2^{1-\nu }}{\Gamma (\nu )}}{\Bigg (}{\sqrt {2\nu }}{\frac {r}{\ell }}{\Bigg )}^{\!\nu} K_{\nu }{\Bigg (}{\sqrt {2\nu }}{\frac {r}{\ell }}{\Bigg )} \, ,
\end{equation}

where $\Gamma(\cdot)$ is the gamma function, $K_{\nu }$ is a modified Bessel function of the second kind of order $\nu$, $r=\norm{\mathbf{x}-\mathbf{x}'}$, and $\ell$ is the lengthscale controlling the smoothness from one perspective: large $\ell$ characterizes functions which change slowly and can be reliably extrapolated further away.

The Matérn covariance functions can be written as a product of an exponential and a polynomial of order $p$, when $\nu$ is half-integer: $\nu=p+1/2, \, p \in \mathds{N}_{0}$. The hyperparameter $\nu$ controls the smoothness from another perspective: when $\nu=1/2$, the Matérn kernel becomes the exponential kernel (continuous but not differentiable); as $\nu \to \infty$, it becomes the exponentiated quadratic kernel (infinitely differentiable). gpml argued that the most interesting cases for machine learning would be $\nu=3/2$ and $\nu=5/2$.

For gas flaring time series, as operators might change flaring strategy at any given time due to policy changes, gas processing facility deployment, gas price fluctuation, etc., the latent process might not be as smooth as infinitely differentiable. Instead the Matérn kernel is harnessed which is capable of inducing non-smooth function realizations to handle those discontinuities. Specifically the Matérn kernel with $\nu=5/2$ is chosen for this dissertation with the input space $\mathcal{X} \subseteq \mathds{R}^1$:

\begin{equation}
k_{\textsc{mat\'ern52}}(x, x'; \ell) \coloneqq \left(1 + \frac{\sqrt{5(x - x')^2}}{\ell} +
\frac{5(x-x')^2}{3\ell^2}\right)
\exp\!\left[ - \frac{\sqrt{5(x - x')^2}}{\ell} \right],
\end{equation}

where $x$ vary over the time domain.

* The standard periodic kernel due to mackay1998introductionGP:

\begin{equation}
k_{\textsc{periodic}}(x, x'; T, \ell) \coloneqq \exp\!\left( -\frac{\sin^2(\uppi |x-x'| \frac{1}{T})}{2\ell^2} \right),
\end{equation}

where $T$ denotes the period. This kernel is used for modeling seasonal behaviors.

* The white noise kernel, which is given by:

\begin{equation}
k_{\mathsmaller{\textsc{WhiteNoise}}}(x, x'; \delta) \coloneqq \delta^2 \mathbf{I}_n,
\end{equation}

where $\delta^2$ is the variance of the noise. In this dissertation, the usage of the white noise kernel is for stabilizing the computation of the covariance matrix. Adding a small value of diagonal shift will try to guarantee the resulting covariance matrix is always positive semi-definite.

A nice property is that the sum and product of the established kernels are still valid kernels. This fact is also exploited in the model building process in this work.

§.§.§ Inference and Model Reparameterization

In practice, one always works with a dataset of finite size. In such situations, a multivariate normal prior distribution is placed on the vector of function values $\mathbf{f}$,

\begin{equation}
\mathbf{f} \sim \operatorname{MVNormal}(\mathbf{m}_{x},\, \mathbf{K}_{xx}) \, , \label{eq:gp_latent_cp}
\end{equation}

where the vector $\mathbf{m}_x$ and the matrix $\mathbf{K}_{xx}$ are the mean function and covariance function evaluated over the inputs $\mathbf{x}$.

A key question which has significant impact on the inference is how to learn the hyperparameters from data. A natural (and popular) approach is to conduct maximum likelihood estimation, i.e., generating point estimates leveraging the data. However, as Betancourt2017RobustStan showed with experiment results, both regularized and unregularized maximum marginal likelihood have limited performance in terms of fitting robustly and recovering the true data generating process. Technically, given a particular kernel with particular hyperparameters, a GP does not support an entire Hilbert space but only a slice through that space; changing the hyperparameters by an infinitesimal amount yields a different slice which has no overlap with the original one. Therefore in this dissertation, a full Bayesian approach is taken for the GP inference, i.e., the entire Hilbert space of functions is considered by taking into account all of the possible hyperparameters for a specific kernel.

For the class of problems which have Gaussian observation models, GP has nice closed-form posterior results. However, for the situations which do not have Gaussian observation models, for examples the ones in this dissertation which employ $\stut$ or Poisson likelihood, there does not exist analytical solutions. HMC as discussed in Section <ref> is used to sample from the posteriors.

Specifically, the noncentered parameterization of the latent multivariate Gaussian is exploited. The reparameterized model is

\begin{align}
\widetilde{\mathbf{f}} &\sim \operatorname{MVNormal}(\mathbf{0}_n, \, \mathbf{I}_n) \label{eq:gp_repar_begin} \\[0.5ex]
\mathbf{L} &= \operatorname{Cholesky}(\mathbf{K}_{xx}) \\[0.5ex]
\mathbf{f} &= \mathbf{m}_{x} + \mathbf{L} \cdot \widetilde{\mathbf{f}} \label{eq:gp_repar_end}
\end{align}

which defines the same distribution as <ref> but induces a nicer posterior geometry for HMC to explore and sample from Betancourt2017RobustStan.

Once the learning on hyperparameters is done, posterior predictive distribution of the latent function values which are not part of the original dataset is obtained by

\begin{equation}\label{eq:gp-cond}
\mathbf{f}_{*} \mid \mathbf{f} \sim \operatorname{MVNormal} \!\left(
\mathbf{m}_{*} + \mathbf{K}_{x*}^{\top} \mathbf{K}_{xx}^{-1} (\mathbf{f}-\mathbf{m}_{x}), \,
\mathbf{K}_{**} - \mathbf{K}_{x*}^{\top} \mathbf{K}_{xx}^{-1}\mathbf{K}_{x*} \right),
\end{equation}

where $\mathbf{m}_{*}$ is the mean function evaluated at the new inputs, $\mathbf{K}_{**}$ is the covariance between the new inputs, and $\mathbf{K}_{x*}$ is the covariance between the original inputs and the new inputs.

§.§ Suite of Models for Pattern Recognition

This section presents models built from various angles, with the goal of providing a coherent framework for learning the flaring pattern and behavior in a principled manner. Each model is tested on real flaring data from North Dakota. Whenever more granular analytics capabilities are demonstrated through investigations at oilfield level or operator level, the data from a major producing field, the Blue Buttes Oilfield Alexeyev2017WellFormation, and one operator, denoted by `Operator A' are used.

§.§.§ Modeling Proportion of Gas Flared

The proportion of gas production that is flared is an indicator of flaring intensity and energy efficiency. It is interesting to investigate whether the proportion has changed over a period of time for certain operators and oilfields. The model is specified through gas_p_begingas_p_end:

\begin{align}
\ell &\sim \operatorname{Gamma}(2, 1) \label{gas_p_begin} \\[0.5ex]
\eta &\sim \operatorname{Half-Cauchy}(5) \\[0.5ex]
\nu &\sim \operatorname{Gamma}(2, 0.1) \\[0.5ex]
\hat{\sigma}^2 &\sim \operatorname{Half-Cauchy}(5) \\[0.5ex]
k &= \eta^2 \times k_{\textsc{mat\'ern52}}(x, x'; \ell) \\[0.5ex]
f &\sim \mathcal{GP}(0, k) \\[0.5ex]
\pi_i &= \operatorname{logit}^{-1}\!\left(f(x_i) \right) \\[0.5ex]
\mu_i &= \pi_i \times G_i \\[0.5ex]
F_i &\sim \stut(\nu, \mu_i, 1/\hat{\sigma}^2) \label{gas_p_end}
\end{align}

ℓ is the lengthscale for the Matérn kernel;

η is the marginal deviation parameter controlling how strongly the latent functions vary in the output space;

ν is the degrees of freedom for the $\stut$ likelihood;

σ̂^2 controls the inverse scaling parameter of the $\stut$ likelihood (analogous to the precision of a Gaussian distribution);

k is the covariance function for the GP;

f denotes the latent process, which is distributed according to the GP;

π_i is the underlying flaring gas proportion of month $i$. Since proportion is bounded between 0 and 1, the inverse-logit function is applied to the latent process;

G_i is the total gas production of month $i$;

μ_i denotes the underlying flared volume of month $i$;

F_i is the reported flared volume, which is modeled using a $\stut$ observation model.

The reasoning behind choosing a $\stut$ observation model is to make the model specification be able to generalize to as many entities as possible and be robust to (potentially many) outliers and noisy data points. This is due to the fact that at this time, operators have to estimate the flared volume by their own procedures and conduct reporting, in which case inaccuracies are introduced unintentionally or intentionally. The heavier tail of Student's $t$-distribution is a natural decision in modeling to deal with those phenomena. This line of thought, i.e., design models that are generic and robust, is indeed reflected in choosing the half-Cauchy priors (which are heavy-tailed and very weakly informative) and GP as a nonparametric regression technique.

To demonstrate this model's capability on real data, both the Blue Buttes Oilfield and Operator A are tested. The production and flared volumes coming from NDIC are used. For the oilfield, the posterior distributions and trace plots of the hyperparameters are presented in <ref>. The posterior predictive samples for the underlying process of gas flaring proportions ($\pi_i$) are demonstrated in <ref>, which depict the trend very clearly. The colored bands have the below coverage for the posterior samples:

* The darkest colored band (in the center at a certain $x$ location) represents the 49 percentile to 51 percentile;

* The lightest colored band (characterized by the widest interval at a certain $x$ location) represents the 1 percentile to 99 percentile.

Additionally, 30 random samples are drawn from the GP posterior and plotted on the same figure, showing as thin lines. The latent functions do not go through all the observed data points, in which case the model would have been overfitted; instead they present the possible functions which are most compatible with the data as well as the assumptions inherent in the model. On one hand, the insights are already obtained, i.e., the underlying process is inferred. On the other hand, this serves as an anomaly detection tool. For example, the state government might be interested to look into that observed data in the second half of 2019 which deviated quite a lot from the “true” process, e.g. to audit the reporting for that month or to investigate what had happened that led to a sudden huge drop in flaring in just one month.
[!htbp]field_gas_proportion_tracefigures/field_gas_proportion_tracePosterior distributions and trace plots for the Blue Buttes Oilfield gas flaring proportion model. Well mixing and convergence have been achieved.
[!htbp]field_gas_proportion_ppcfigures/field_gas_proportion_ppcPosterior predictive samples showing the gas flaring proportion variations at the Blue Buttes Oilfield. Blue points are the observed data while red lines present the posterior samples.

With the exact same model specification, the model is also run with the operator's data. The posterior distributions and trace plots of the hyperparameters are presented in <ref>. The posterior predictive samples for the underlying process of gas flaring proportions ($\pi_i$) are demonstrated in <ref>. It can be seen this operator's flaring proportion time series is more jagged than the Blue Buttes Oilfield (which is operated by more than five companies). A operator can change flaring strategies more swiftly which can be captured as well. Nevertheless the long-term trend is also available. Comparing <ref> and <ref>, it can be seen the posterior distributions are very different. However the priors for them were specified in the exact same way. This showcases the power of Bayesian approach. Taking $\ell$ as an example, a $\operatorname{Gamma}(2,1)$ prior is placed on it. However, after conditioning on the data, the operator model reports smaller lengthscale values on average (indicating jagged processes), whereas the oilfield model reports larger lengthscale values (suggesting smoother processes).
[!htbp]operator_gas_proportion_tracefigures/operator_gas_proportion_tracePosterior distributions and trace plots for the Operator A gas flaring proportion model. Well mixing and convergence have been achieved. Notice the differences between these inference results and those in <ref>, both of which are based on exactly the same priors and likelihood, demonstrating the model specification's wide applicability.
[!htbp]operator_gas_proportion_ppcfigures/operator_gas_proportion_ppcPosterior predictive samples showing the gas flaring proportion variations of Operator A. Blue points are the observed data while red lines present the posterior samples.

Order 24665, which is established by the , defines the gas capture percentage $p_{\textsc{cap}}$ as

\begin{equation}
p_{\textsc{cap}} = \frac{G_{\textsc{sold}} + G_{\textsc{used}} + G_{\textsc{proc}}}{G_{\textsc{prod}}} \, ,
\end{equation}

G_sold is the monthly gas sold;

G_used is the monthly gas used on lease;

G_proc is the monthly gas processed;

G_prod is the monthly gas produced.

Since North Dakota bans the venting of natural gas U.S.DepartmentofEnergy2019NorthSheet, it is obvious the model developed in this section provides a powerful tool for NDIC to evaluate compliance with the gas capture goals: at a given month $i$, $p_{\textsc{cap}}=1-\pi_i$. Furthermore, when looking at the model specification, there is nothing special that encodes the data sources and location information. A user of this model is free to use satellite estimation as the observed data or apply it to the Permian Basin, and conduct inference on the flaring proportion. This is a benefit from using nonparametric and interpretable models as opposed to black box models (such as the neural networks, in which case the learned weights and bias inside the network provide little or no domain insights). The author hopes this section provides a comprehensive view in terms of how and why to use GP, with real data. Models built and presented in later sections follow a similar flow.

§.§.§ Modeling Proportion of Wells Flaring

The proportion of wells that conduct flaring in a month can reflect a company's flaring strategy and is an indicator of flaring magnitude. It is interesting to investigate how this indicator varies for a certain entity in a certain time period. The model is specified through well_p_beginwell_p_end:

\begin{align}
\ell &\sim \operatorname{Gamma}(2, 1) \label{well_p_begin} \\[0.5ex]
\eta &\sim \operatorname{Half-Cauchy}(5) \\[0.5ex]
k &= \eta^2 \times k_{\textsc{mat\'ern52}}(x, x'; \ell) \\[0.5ex]
f &\sim \mathcal{GP}(0, k) \\[0.5ex]
p_i &= \operatorname{logit}^{-1}\!\left(f(x_i) \right) \\[0.5ex]
W_i &\sim \operatorname{Binomial}(N_i, p_i) \label{well_p_end}
\end{align}

where $p_i$ is the unobserved “true” proportion of wells that conduct flaring in month $i$, $N_i$ is the total number of active wells in month $i$, and $W_i$ is the observed (i.e., estimated and reported by company) number of wells that conduct flaring in month $i$. The rest of the symbols have the same meaning as in Model <ref>.

To demonstrate this model's capability on actual data, both the Blue Buttes Oilfield and Operator A are tested. For the oilfield, the posterior distributions and trace plots of the hyperparameters are presented in <ref>. The posterior predictive samples for the underlying process of well flaring proportion ($p_i$) are demonstrated in <ref>. The visualization strategy (different colors represent different percentiles, etc.) is the same as in Section <ref>.
[!htbp]field_well_proportion_tracefigures/field_well_proportion_tracePosterior distributions and trace plots for the Blue Buttes Oilfield well flaring proportion model. Well mixing and convergence have been achieved.
[!htbp]field_well_proportion_ppcfigures/field_well_proportion_ppcPosterior predictive samples showing the well flaring proportion variations at the Blue Buttes Oilfield. Blue points are the observed data while red lines present the posterior samples.

With the exact same model specification, this model is also tested with the operator's data. The posterior distributions and trace plots of the hyperparameters are presented in <ref>. The posterior predictive samples for the underlying process of well flaring proportion ($p_i$) are demonstrated in <ref>. Comparing the two sets of figures from the oilfield and the operator, it can be seen:

* With the same prior placed on the lengthscale $\ell$, the oilfield model learns from the data and gives a posterior mode around 5.5, whereas the operator model gives a posterior mode around 10.0. This is also reflected in the posterior samples time series plot: the oilfield experienced some well flaring proportion changes in relative shorter time periods, whereas the operator underwent changes on a longer time span.

* The oilfield's posterior samples time series show narrower percentile bands while the operator's show wider percentile bands. This is due to the fact that the operator chosen here had smaller number of wells than the oilfield. Since the binomial observation model is used for each month's flaring well count, this naturally represents and quantifies the uncertainties (i.e., binary data contains less information especially when the sample size is small), as well as aligns with the expectation that when there is more data, there should be less uncertainties; when there is less data, there should be more uncertainties.

[!htbp]operator_well_proportion_tracefigures/operator_well_proportion_tracePosterior distributions and trace plots for the Operator A well flaring proportion model. Well mixing and convergence have been achieved. Notice the differences between these inference results and those in <ref>, both of which are based on exactly the same priors and likelihood, demonstrating the model specification's wide applicability.
[!htbp]operator_well_proportion_ppcfigures/operator_well_proportion_ppcPosterior predictive samples showing the well flaring proportion variations of Operator A. Blue points are the observed data while red lines present the posterior samples.

This really showcases how and why to encode domain expertise in flaring data analytics while exploiting machine learning models, which is also the reason to choose the Bayesian approach. One could fit a black box model either with target values $W_i \in \mathds{R}$, or without any probabilistic view (e.g., to optimize for the best deterministic function mapping in the hypothesis space). But either of those would be fundamentally flawed. Domain expertise indicates the well count has to be a non-positive integer, i.e., $W_i \in \mathds{N}_{0}$. Furthermore, neither the NDIC reporting nor the satellite estimation is ever produced in a noise-free environment, and therefore probabilistic modeling is a must. Compared to frequentist machine learning, Bayesian learning is entirely probabilistic and gives one the capability and freedom to encode his/her domain expertise.

§.§.§ Modeling Flare Detection Count

Satellite detected flare count provides an unbiased indicator of flaring intensity. How this indicator varies in a certain time period for a certain entity is valuable information to obtain. The model is specified through cox_begincox_end. Essentially the latent process is modeled as a Gaussian Cox process cox_process, where the Poisson process has varying intensity across time domain and a GP prior is placed on this intensity.

\begin{align}
\ell &\sim \operatorname{Gamma}(2, 1) \label{cox_begin} \\[0.5ex]
\eta &\sim \operatorname{Half-Cauchy}(5) \\[0.5ex]
k &= \eta^2 \times k_{\textsc{mat\'ern52}}(x, x'; \ell) \\[0.5ex]
f &\sim \mathcal{GP}(0, k) \\[0.5ex]
\lambda_i &= \exp\!\left(f(x_i) \right) \\[0.5ex]
C_i &\sim \operatorname{Poisson}(\lambda_i) \label{cox_end}
\end{align}

where $\lambda_i$ is the unobserved flaring intensity (“true” count) in month $i$ and $C_i$ is the reported VIIRS detection count in month $i$. Since $\lambda_i$ is bounded to be positive, the natural exponential function is applied to the latent process. The rest of the symbols have the same meaning as in Model <ref>.

For the task of flaring pattern recognition, the author believes this approach (leveraging a Gaussian Cox process) is a nicer surrogate than a popular change point model presented in bmh,Salvatier2016ProbabilisticPyMC3,stan_man, which is specified by:

\begin{align}
e   &\sim  \operatorname{Exponential}(r_e) \\[0.5ex]
l   &\sim  \operatorname{Exponential}(r_l) \\[0.5ex]
s   &\sim  \operatorname{Uniform}(1, T) \\[0.5ex]
C_i &\sim  \operatorname{Poisson}(i < s \; ? \; e \; : \; l)
\end{align}

where $e$ and $l$ are the early and late rates respectively, $r_e$ and $r_l$ controls the priors for the early and late rates, $s$ is the change point, $T$ is the total time period, and the rate in the Poisson likelihood is decided through a ternary conditional operator (?:). The reason is that, although this model could be generalized to more than one change point, its usage is restricted by the assumption that any period between two adjacent change points has a constant rate. This limitation becomes obvious when analyzing the actual flaring data in the discussions below, and is a major disadvantage of the change point model.

The Gaussian Cox process model is tested with the Blue Buttes Oilfield's data. Since only VIIRS data is used, the whole time series is analyzed beginning in 2012. The posterior distributions and trace plots of the hyperparameters are presented in <ref>. The posterior predictive samples for the underlying process of flare count ($C_i$) are demonstrated in <ref>. The visualization strategy (different colors represent different percentiles, etc.) is the same as in Section <ref>. From the time series plot, it can be seen the observations from 2014 to 2017 can possibly be described by a change point model (with late 2015 being a potential change point), but the steady growth before and after that time span will frustrate accurate inference with such a model.
[!htbp]field_coxProcess_tracefigures/field_coxProcess_tracePosterior distributions and trace plots for the Blue Buttes Oilfield flare count model. Well mixing and convergence have been achieved.
[!htbp]field_coxProcess_ppcfigures/field_coxProcess_ppcPosterior predictive samples showing the flare count variations at the Blue Buttes Oilfield. Blue points are the observed data while red lines present the posterior samples.

This model's inference results serve as a type of confirmation, if not evidence, in terms of whether or not an entity achieves the goal/target in reducing the number of wells flaring, when the detection count is used as a surrogate for the number of wells flaring. In practice, reducing the number of wells flaring is exactly the second goal of the regulatory policy introduced by the  in 2014. If the state government is interested in this order's effectiveness from a macroscopic standpoint, the model can also be used to conduct inferences with the state level data. In this case, the posterior distributions and trace plots of the hyperparameters are presented in <ref>. The posterior predictive samples for the underlying process of flare count ($C_i$) are demonstrated in <ref>.
[!htbp]stateAllYr_coxProcess_tracefigures/stateAllYr_coxProcess_tracePosterior distributions and trace plots for the North Dakota flare count model. Well mixing and convergence have been achieved. Notice the differences between these inference results and those in <ref>, both of which are based on exactly the same priors and likelihood, demonstrating the model specification's wide applicability.
[!htbp]stateAllYr_coxProcess_ppcfigures/stateAllYr_coxProcess_ppcPosterior predictive samples showing the flare count variations in North Dakota. Blue points are the observed data while red lines present the posterior samples.

The percentile bands in this case are quite narrow, which indicate greater confidence in the inferences about the data generating process given the model assumptions. By not (over)fitting to each and every observation, interesting patterns are discovered, for example in every year there is one and only one peak that happened around June. It is worth pointing out that there is no model that can tell the modeler if his/her assumptions are good, only domain expertise might. This model employing a Poisson observation model could be considered “rigid” due to the fact that a Poisson likelihood has only one parameter $\lambda$ (to control both the mean and variance) and, furthermore, when $\lambda$ is large as in this scenario, a Poisson distribution is well approximated by a normal distribution. Whenever the state government believes that overdispersion might exist, other observation models such as the negative binomial distribution could be considered. In such cases, only <ref> needs to be changed to the negative binomial likelihood, with a prior added for the overdispersion parameter. The specific parameterization is given by <Ref> in Section <ref>. This really showcases both the flexibility and interpretability of taking a Bayesian approach for high-stakes decision making areas including flaring data analytics.

§.§.§ Modeling Proportion of Oil Flared

As crude oil (as opposed to natural gas) is the main commodity at this time, the amount of gas in a barrel of oil equivalent () that is flared provides an indicator of production efficiency due to flaring. In this work, the normalized quantity, proportion of oil production being flared, is used such that the model specification is generic for large and small entities. The model is specified through oil_p_beginoil_p_end:

\begin{align}
\ell &\sim \operatorname{Gamma}(2, 1) \label{oil_p_begin} \\[0.5ex]
\eta &\sim \operatorname{Half-Cauchy}(5) \\[0.5ex]
\nu &\sim \operatorname{Gamma}(2, 0.1) \\[0.5ex]
\hat{\sigma}^2 &\sim \operatorname{Half-Cauchy}(5) \\[0.5ex]
k &= \eta^2 \times k_{\textsc{mat\'ern52}}(x, x'; \ell) \\[0.5ex]
f &\sim \mathcal{GP}(0, k) \\[0.5ex]
\pi_i &= \operatorname{logit}^{-1}\!\left(f(x_i) \right) \\[0.5ex]
\mu_i &= \pi_i \times O_i \\[0.5ex]
c &= \frac{\SI{6}{\mcf}}{\SI{1}{\boe}} \\[0.5ex]
F_i/c \eqqcolon E_i &\sim \stut(\nu, \mu_i, 1/\hat{\sigma}^2) \label{oil_p_end}
\end{align}

π_i is the underlying flaring  proportion of month $i$;

O_i is the total oil production of month $i$;

μ_i denotes the “true” flared  of month $i$;

c denotes the conversion factor that 6 equals 1, given by the usgs_boe;

E_i is the reported flared , which is modeled using a $\stut$ observation model.

The rest of the symbols have the same meaning as in Model <ref>. To test this model's performance on real data, both the Blue Buttes Oilfield and Operator A are used. For the oilfield, the posterior distributions and trace plots of the hyperparameters are presented in <ref>. The posterior predictive samples for the underlying process of  flaring proportion ($\pi_i$) are demonstrated in <ref>. The visualization strategy (different colors represent different percentiles, etc.) is the same as in Section <ref>.
[!htbp]field_BOE_proportion_tracefigures/field_BOE_proportion_tracePosterior distributions and trace plots for the Blue Buttes Oilfield  flaring proportion model. Well mixing and convergence have been achieved.
[!htbp]field_BOE_proportion_ppcfigures/field_BOE_proportion_ppcPosterior predictive samples showing the  flaring proportion variations at the Blue Buttes Oilfield. Blue points are the observed data while red lines present the posterior samples.

With the exact same model specification, this model is also tested with the operator's data. The posterior distributions and trace plots of the hyperparameters are presented in <ref>. The posterior predictive samples for the underlying process of  flaring proportion ($\pi_i$) are demonstrated in <ref>.
[!htbp]operator_BOE_proportion_tracefigures/operator_BOE_proportion_tracePosterior distributions and trace plots of the  flaring proportion model for Operator A. Well mixing and convergence have been achieved. Notice the differences between these inference results and those in <ref>, both of which are based on exactly the same priors and likelihood, demonstrating the model specification's wide applicability.
[!htbp]operator_BOE_proportion_ppcfigures/operator_BOE_proportion_ppcPosterior predictive samples showing the  flaring proportion variations of Operator A. Blue points are the observed data while red lines present the posterior samples.

Comparing the two sets of figures from the oilfield and the operator, it can be observed:

* With the same prior placed on the lengthscale $\ell$, which has a mean of 2 (months), both models have updated the posterior to move away from this mean, reflecting a long range variation. The oilfield has a posterior mode about 1 year while the operator has a mode around 15 months. The operator has much larger reporting variability, shown by the parameter $\hat{\sigma}^2$.

* With a $\stut$ likelihood, both models demonstrate robustness to outliers and overfitting. This can be seen from the oilfield's late 2019 observations and the operator's early 2016 observations. For the posterior function samples, shown as the thin lines, some of them are indeed pulled towards those “outliers”. However, the percentile plots (shown as the colored bands) are not impacted and those really can be interpreted as the trend which is most compatible with the data and the assumptions. This built-in Occam's razor of the Bayesian approach when choosing appropriate priors is very impressive. In many of the frequentist machine learning methods, if the regularization strategy is not implemented well especially when the sample size is not huge enough for the asymptotic properties to kick in, outliers become “influential observations” that will have a huge undesirable effect on the inference results.

§.§.§ Modeling Scale Factor between VIIRS and NDIC

Both NDIC and VIIRS reporting give (estimated) flared gas volume. The scale factor between the two sources provides insights into whether NDIC reporting is consistent:

* for different entities (e.g., among a group of operators), and

* for one entity when looking at a certain time period.

This is based on the fact that the satellite detection processing algorithm is unbiased and consistent. Item 2 is particularly interesting in terms of time series analytics. The model is specified through sf_beginsf_end:

\begin{align}
\ell_{\textsc{mat}} &\sim \operatorname{Gamma}(8, 2) \label{sf_begin} \\[0.5ex]
\eta_{\textsc{mat}} &\sim \operatorname{Half-Cauchy}(5) \\[0.5ex]
T &\sim \mathcal{N}(12, 1) \\[0.5ex]
\ell_{\textsc{per}} &\sim \operatorname{Gamma}(4, 3) \\[0.5ex]
\eta_{\textsc{per}} &\sim \operatorname{Half-Cauchy}(5) \\[0.5ex]
\nu &\sim \operatorname{Gamma}(2, 0.1) \\[0.5ex]
\hat{\sigma}^2 &\sim \operatorname{Half-Cauchy}(5) \\[0.5ex]
k_{\textsc{mat}} &= \eta_{\textsc{mat}}^2 \times k_{\textsc{mat\'ern52}}(x, x'; \ell_{\textsc{mat}}) \\[0.5ex]
k_{\textsc{per}} &= \eta_{\textsc{per}}^2 \times k_{\textsc{periodic}}(x, x'; T, \ell_{\textsc{per}}) \\[0.5ex]
k_{\textsc{wn}} &= k_{\mathsmaller{\textsc{WhiteNoise}}}(x, x'; \delta=\num[output-exponent-marker=\ensuremath{\mathrm{e}}]{1e-6}) \\[0.5ex]
f &\sim \mathcal{GP}(0, k_{\textsc{mat}} + k_{\textsc{per}} + k_{\textsc{wn}}) \\[0.5ex]
\beta_i &= \exp\!\left(f(x_i) \right) \\[0.5ex]
\mu_i &= \beta_i \times \mathrm{VIIRS}_i \\[0.5ex]
\mathrm{NDIC}_i &\sim \stut(\nu, \mu_i, 1/\hat{\sigma}^2) \label{sf_end}
\end{align}

ℓ_mat is the lengthscale for the Matérn kernel;

η_mat is the marginal deviation for the Matérn kernel;

T is the period for the periodic kernel;

ℓ_per is the lengthscale for the periodic kernel;

η_per is the marginal deviation for the periodic kernel;

k_mat is the Matérn kernel (component);

k_per is the periodic kernel (component);

k_wn is the white noise kernel (component);

f denotes the latent process, which is distributed according to a GP whose covariance function is the sum of 3 kernels;

β_i is the underlying scale factor between VIIRS and NDIC of month $i$. Since this scale factor is bounded to be positive, the natural exponential function is applied to the latent process;

VIIRS_i is the VIIRS reported volume of month $i$;

μ_i denotes the underlying flared volume of month $i$;

NDIC_i is the NDIC reported volume of month $i$, which is modeled using a $\stut$ observation model.

The rest of the symbols have the same meaning as in Model <ref>. The reason for adding a periodic kernel is to investigate if there are any seasonal patterns. Maintaining a proper Bayesian workflow lets the data speak for itself, i.e., whether there exists seasonal behaviors or not, as shown by the two case studies in this section.

The model is first fitted with the state level data to investigate the macroscopic reporting consistency. The posterior distributions and trace plots of the hyperparameters are presented in <ref>. The posterior predictive samples for the underlying process of the scale factor variations ($\beta_i$) are demonstrated in <ref>. The visualization strategy (different colors represent different percentiles, etc.) is the same as in Section <ref>. From the posterior time series plot, it can be seen in general the volumes from NDIC reporting is smaller than that of VIIRS reporting, except for the times when the total flaring magnitude was small (indicated by the smaller points). More importantly, within each and every year from 2015 to 2018, there is a decreasing trend in the values of the scale factor ($\beta_i$) around midyear. Each year's latent process from Q2 to Q3 can be viewed as a “seesaw”, with July being the middle pivot point and the months after July always going down. Note that within each year, the NDIC reporting of flared volumes might increase steadily or a lot (which was actually happening from the time series plot in <ref>), however this scale factor declining trends indicate the satellites observed much greater flaring activities than what was reported by the companies! This finding suggests that the NDIC reporting is very likely not consistent throughout the year, and the state government should be concerned that some companies might underreport their flared volumes especially in the second half of the year.
[!htbp]nd_SF_tracefigures/nd_SF_tracePosterior distributions and trace plots for the North Dakota VIIRS-NDIC scale factor model. Well mixing and convergence have been achieved.
[!htbp]nd_SF_ppcfigures/nd_SF_ppcPosterior predictive samples showing the scale factor variations of North Dakota. Blue points are the observed data while red lines present the posterior samples. Larger points indicate greater flaring magnitude as observed from VIIRS.

A interesting question arises: is this seasonal behavior universal across all the entities? The answer is unfortunately no, which indicates some operators likely reported their flared volume in an inconsistent manner throughout the entire year. In fact, if the Blue Buttes Oilfield data is used to fit the model, rather consistent behavior is observed. In this case, the posterior distributions and trace plots of the hyperparameters are presented in <ref>. The posterior predictive samples for the underlying process of the scale factor variations ($\beta_i$) are demonstrated in <ref>. With the exact same model specification incorporating the periodic kernel, no apparent seasonal behaviors are discerned by the inference process. There are much uncertainties around the time of early 2016, where the point sizes indicate the overall flaring magnitudes were small as observed from VIIRS, and the NDIC reported volumes were actually larger than that of VIIRS. This could be due to the truncation effects instead of the reporting inconsistencies, i.e., when the flares are sporadic and weaker, they are not easily captured by the satellites, resulting in a truncated sample for the VIIRS processing workflow. By applying this model and workflow to the other major producing fields, it will likely pick up the ones who have the “seesaw” behaviors in their reporting.
[!htbp]oilfield_SF_tracefigures/oilfield_SF_tracePosterior distributions and trace plots for the Blue Buttes Oilfield VIIRS-NDIC scale factor model. Well mixing and convergence have been achieved. Notice the differences between these inference results and those in <ref>, both of which are based on exactly the same priors and likelihood, demonstrating the model specification's wide applicability.
[!htbp]oilfield_SF_ppcfigures/oilfield_SF_ppcPosterior predictive samples showing the scale factor variations in the Blue Buttes Oilfield. Blue points are the observed data while red lines present the posterior samples. Larger points indicate greater flaring magnitude as observed from VIIRS.

§.§.§ Predicting NDIC Flared Volume

GP is not only fully capable of making predictions once the model hyperparameters are learned, but it can provide rigorously constructed intervals quantifying uncertainties as well through <ref>, for which many of the frequentist machine learning methods fail to do. The author chooses to present one particular prediction case study, that is to predict NDIC reported volume based on the projected scale factor between VIIRS and NDIC. This will be a particular interesting deployment scenario once fast satellite detection/estimation is available, which takes less time than waiting on company reports followed by compiling everything into an analytics-ready format.