Enhanced monitoring of atmospheric methane from space with hierarchical
Bayesian inference

ID Clayton Roberts1∗,  ID Oliver Shorttle1,2,  ID Kaisey Mandel1,3,4, Matthew
Jones5,

Rutger Ijzermans6,  ID Bill Hirst7 and  ID Philip Jonathan5,8

∗ Corresponding author
1 Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge
CB3 0HA, UK
2 Department of Earth Sciences, University of Cambridge, Downing Street,
Cambridge CB2 3EQ, UK
3 Statistical Laboratory, Department of Pure Mathematics and Mathematical
Statistics, University of Cambridge, Wilberforce Road, Cambridge CB3 0WB, UK
4 The Alan Turing Institute, Euston Road, London NW1 2DB, UK
5 Shell Global Solutions (U.K.), London SE1 7NA, UK
6 Shell Global Solutions International B.V., Grasweg 31, 1031 HW Amsterdam,
the Netherlands
7 Atmospheric Monitoring Sciences, Amsterdam, the Netherlands
8 Department of Mathematics and Statistics, Lancaster University, LA1 4YW, UK

## Abstract

Methane is a strong greenhouse gas, with a higher radiative forcing per unit
mass and shorter atmospheric lifetime than carbon dioxide. The remote sensing
of methane in regions of industrial activity is a key step toward the accurate
monitoring of emissions that drive climate change. Whilst the TROPOspheric
Monitoring Instrument (TROPOMI) on board the Sentinal-5P satellite is capable
of providing daily global measurement of methane columns, data are often
compromised by cloud cover. Here, we develop a statistical model which uses
nitrogen dioxide concentration data from TROPOMI to accurately predict values
of methane columns, expanding the average daily spatial coverage of
observations of the Permian Basin from 16% to 88% in the year 2019. The
addition of predicted methane abundances at locations where direct
observations are not available will support inversion methods for estimating
methane emission rates at shorter timescales than is currently possible.

## Introduction

Methane and carbon dioxide are the two dominant anthropogenic greenhouse
gasses responsible for warming Earth’s atmosphere above its pre-industrial era
temperature [1]. The current average atmospheric concentrations of these
gasses are $\sim$1,880 parts per billion by volume (ppbv) for methane [2] (an
increase of over 1,000 ppbv in the past 250 years as a result of human
activity [3]) and $\sim$413 parts per million by volume (ppmv) for carbon
dioxide [4]; the mass of carbon dioxide in the atmosphere now is roughly 600
times the mass of methane. However, methane is a dramatically stronger
absorber of thermal radiation than carbon dioxide, and despite its much lower
concentration is responsible for trapping about 58% of the heat that carbon
dioxide does when the effect of methane’s inevitable decay products are also
included [5].

Satellite-borne remote sensing and monitoring of greenhouse gas emissions is
playing an increasingly important role in assessing mankind’s impact on the
climate [6, 7], as direct measurements can displace or complement inventory
estimates [8]. Measurements of methane column concentrations from space began
in 2003 with the SCanning Imaging Absorption SpectroMeter for Atmospheric
CHartographY (SCIAMACHY) on board ENVISAT, an ESA mission which terminated in
2012 [9]. SCIAMACHY was succeeded by the Greenhouse Gases Observing Satellite
(GOSAT, 2009-present) and GOSAT2 (2018-present) [10, 11], operated by JAXA.
Both GOSAT and GOSAT2 have improved capabilities with pixel resolutions of
10x10 km2 and global coverage in 3 days when compared to SCIAMACHY, which
required 6 days for global coverage with a ground pixel resolution of 30x60
km2. In 2016 the private enterprise GHGSat-D instrument was launched, with a
greatly improved pixel resolution over GOSAT of 50x50 m2 for targeted viewing
of selected methane sources [6]. Next generation instruments such as the
Advanced Hyperspectral Imager (AHSI) launched on board China’s GaoFen5
satellite in 2018 provide methane retrievals with a pixel resolution down to
30 metres [12], but such missions are sporadically operated for targeted areas
and do not provide the same level of spatial coverage as GOSAT [13]. Leading
the field in regional observation is the TROPOspheric Monitoring Instrument
(TROPOMI) on board the ESA’s Sentinal-5 Precursor satellite. Launched in 2017,
TROPOMI delivers daily global coverage, observing full methane columns with a
ground pixel resolution of 7x7 km2 [14, 15]. With early data successfully
compared with GOSAT [16], TROPOMI is the cornerstone of the ESA’s commitment
to monitoring national pledges towards emission reductions under the Paris
Climate Accord [17].

Remote sensing of methane emissions using satellites is attractive because it
can be frequent, recurring and non-intrusive to operations on the ground.
Satellites therefore have the potential to detect intermittent methane sources
or emissions released during abnormal operating conditions [18]. Observations
from SCIAMACHY [19] and GOSAT [20] have successfully provided top-down
estimates of regional methane emission rates. More recently, observations from
TROPOMI have provided top-down estimates of methane emission rates from the
Permian Basin (shown in Fig. 1, extending from western Texas to southeastern
New Mexico) of 2.7 teragrams per year over a 12-month period from March 2018
to March 2019 [21].

Although the TROPOMI instrument covers the entire surface of the Earth at
least once a day, the actual amount of methane concentration data from TROPOMI
is limited by a variety of factors [14]. Cloud cover and other meteorological
conditions often hamper the retrieval of methane column data, and so does the
presence of aerosols in the atmosphere [22]. Consequently, sparse data over
regions of interest must be averaged over months or years to obtain suitable
coverage for emission budgets [19, 20, 21]. In contrast, TROPOMI observations
of nitrogen dioxide are not subject to the same difficulties posed by cloud
masks [23]. Large excesses of nitrogen dioxide are often linked to regions of
rapid urban expansion and industrial activity [24].

Oil and gas production is often associated with nitrogen dioxide emissions as
well. For example, lit methane flare stacks emit combusted natural gas which
contains nitrogen dioxide [25, 26]; it is also known that the combustion
efficiency in flares is not equal to 100%, leading to the co-emission of
methane into the atmosphere [21, 27]. Similarly, gas-fuelled compressors emit
nitrogen dioxide during operation but they may also emit methane leaking
through their seals. A final example are the pumps and storage tanks which are
an integral part of oil distribution networks: pumps can emit nitrogen dioxide
as a result of fuel combustion, whereas methane can leak from thief hatches in
storage tanks. Indeed, a recent study [27] has shown that there is a
correlation between the nitrogen dioxide and methane concentrations measured
by TROPOMI over the Permian Basin.

In this work, we compensate the missing direct methane data from TROPOMI by
using nitrogen dioxide as a proxy or predictor of methane column density with
the proportionality inferred from sample locations where confident
measurements exist for both species [27]. We develop a Bayesian model to infer
missing methane column data based on co-located column values for nitrogen
dioxide, which expands the spatial coverage of methane observations from
TROPOMI. We use the Permian Basin as a case study; as the most productive
basin in the United States, it produces more than 18,000 million cubic feet of
natural gas and nearly 5 million barrels of oil per day [28]. We fit our model
to TROPOMI observations, test its predictive ability and investigate how the
inclusion of such estimates expands the daily spatial coverage of methane
data. We examine the implications of including inferred values in emission
budgets by calculating the above-background mass of methane observed in the
Permian Basin over the course of the year 2019, and find the inclusion of
inferred methane values results in kilotonnes more of excess methane observed
per day. The spatial augmentation of TROPOMI methane observations will support
inverse methods for estimating methane emission rates on shorter timescales
than currently possible, which will be invaluable as policymakers begin to
require recent and up-to-date methane emission estimates for industrial
regions.

## Results

### Model evaluation

We fit our model to a year’s observations over the Permian Basin from 2019,
using co-located observations contained within a marked study region shown in
Fig. 1. TROPOMI provides one observation per day per location. Our model is a
linear hierarchical model where the relationship between observations of
nitrogen dioxide and methane on a given day $t$ are related to one another by
a linear form with slope parameter $\beta_{t}$ and intercept parameter
$\alpha_{t}$, plus another daily parameter to account for intrinsic scatter
around the linear trend. The daily parameter $\alpha_{t}$ has units of ppbv
and roughly corresponds to the abundance of methane in the study region on day
$t$ that is not associated with lit, burning flare stacks. The daily parameter
$\beta_{t}$ has units of ppbv / (mmol m-2) as we scale all observations of
nitrogen dioxide columns into mmol m-2.

We specify a multivariate normal hierarchical prior on values of $\alpha_{t}$
and $\beta_{t}$ such that

(1) $\displaystyle\begin{pmatrix}\alpha_{t}\\\ \beta_{t}\end{pmatrix}$
$\displaystyle\sim\mathcal{N}\left(\begin{pmatrix}\mu_{\alpha}\\\
\mu_{\beta}\end{pmatrix},\;\begin{pmatrix}\Sigma_{11}&\Sigma_{12}\\\
\Sigma_{21}&\Sigma_{22}\end{pmatrix}\right)$

where $\Sigma_{11}=\sigma_{\alpha}^{2}$, $\Sigma_{22}=\sigma_{\beta}^{2}$ and
$\Sigma_{12}=\Sigma_{21}=\rho\;\sigma_{\alpha}\;\sigma_{\beta}$. The
hyperparameter $\rho$ controls the degree of correlation between $\alpha_{t}$
and $\beta_{t}$, with $\rho\in[-1,\;1]$. The hyperparameter $\mu_{\alpha}$
roughly corresponds to the average background level of methane in the study
region across the year 2019, and the hyperparameter $\sigma_{\alpha}$ is the
deviation that controls the spread of $\alpha_{t}$ around $\mu_{\alpha}$, with
$\sigma_{\alpha}>0$, and similarly for $\mu_{\beta}$ and $\sigma_{\beta}$.

We fit two models to two subsets of data. First, we fit a “data-rich” model to
observations from days that produce numerous highly correlated co-located
observations of methane and nitrogen dioxide in the study region. Afterwards,
we fit a “data-poor” model to observations from all other days where there are
at least two co-located observations of methane and nitrogen dioxide in the
study region. The difference between the two models is that in the data-rich
model, we stipulate a uniform prior independently for each of the
hyperparameters in order for the model to learn their posteriors predominantly
from the data, and in the data-poor model we provide the information learned
from the data-rich model to the hyperparameters as prior information. We do
this so that when the data-poor model is being fit either to sparse or less-
correlated observations we are making full use of the information that can be
obtained from data-rich days. When fitting all models, we monitor the
effective sample size (ESS) and $\hat{R}$ of all parameters in addition to the
energy Bayesian fraction of missing information (E-BFMI) to ensure convergence
and efficient sampling [29, 30]. The data-rich model was fit with an ESS of
212.9 and the data-poor model was fit with an ESS of 860.0.

We examine the joint posterior distribution of the hyperparameters of the
data-rich model, and represent median values of $\mu_{\alpha}$, $\mu_{\beta}$,
$\rho$, $\sigma_{\alpha}$ and $\sigma_{\beta}$ in panel b of Fig. 2. We find
that the correlation $\rho$ between $\alpha_{t}$ and $\beta_{t}$ on data-rich
days for the year 2019 has a median value of -0.18 with a 68% central credible
interval of (-0.30, -0.05). We would expect a value of $\rho<0$ due to the
positive correlation between the amount of methane and nitrogen introduced to
the atmosphere from flare stack combustion chemistry [31, 32, 33, 34], i.e.,
positive correlation in the data implies negative correlation between the
slope and intercept. However, it is important to note that this result was
obtained without stipulating any prior constraint that $\rho$ must be
negative. Multivariate normal hierarchical prior specification and uniform
hyperpriors for hyperparameters in the data-rich model allows us to learn the
degree to which the daily parameters are correlated, and leverage this
information as a prior of appropriate strength in the data-poor model.

Other model hyperparameters that we discuss at this stage are $\mu_{\alpha}$
and $\sigma_{\alpha}$. These two model parameters have physical meaning, where
$\mu_{\alpha}$ can be thought of as the ‘mean’ background of $\mathrm{CH}_{4}$
for the year 2019, and $\sigma_{\alpha}$ as the extent to which the daily
intercept $\alpha$ varies around $\mu_{\alpha}$ from day to day. We find
$\mu_{\alpha}$ to have a median value of 1861.32 ppbv with a 68% central
credible interval of (1859.52, 1863.15), and $\sigma_{\alpha}$ to have a
median value of 14.44 ppbv with a 68% central credible interval of (13.25,
15.73). For the year 2019, the globally-averaged monthly mean value of marine
surface methane ranges from 1858 to 1875 ppbv [2], which reinforces the
physical interpretability of our model.

### Predictive skill

To assess the predictive ability of our models, we performed hold-out testing
by refitting our models to a random selection of 80% of observations in the
study region on each day, using the remaining 20% of TROPOMI methane
observations for comparison against model predictions at those locations. Both
the data-rich and the data-poor model fit to the 80% datasets with
satisfactory E-BFMI, ESS and $\hat{R}$. After generating predictions from the
fitted model and comparing to co-located withheld observations across all
days, we calculate values of reduced chi-squared $\chi_{\nu}^{2}$ for each day
(i.e., chi-squared per degree of freedom), the resulting distribution of which
is shown in panel b of Fig. 3. We also calculate residuals between model
predictions and withheld observations across all days, which have a mean of
0.15 ppbv and standard deviation 12.09 ppbv, correlated with Pearson $R=0.82$
(shown in panel a of Fig. 3). For comparison, first results from TROPOMI were
initially tested against co-located GOSAT observations with a mean difference
of 13.6 ppbv and standard deviation 19.6 ppbv, correlated with Pearson
$R=0.95$ [16]. Panel a of Fig. 3 also demonstrates the tendency of the model
to underestimate methane from high values of observed nitrogen dioxide and
overestimate methane at low values of nitrogen dioxide. This is a result of
the regression dilution caused by uncertain TROPOMI observations of nitrogen
dioxide [35, 36], seen in Panel a of Fig. 2.

###  $\boldsymbol{\beta}$ \- flare stack correlation

After re-fitting the data-rich model to 100% of available observations we
examine posterior estimates of daily parameters as time series. We plot median
values of the estimate for the slope parameter $\beta_{t}$ in Fig. 4 along
with a time series of active flare stacks in the study region on data-rich
days, identified from VIIRS Nightfire [37]. We do not find any correlation
between estimated $\beta_{t}$ and the total number of identified flare stacks
inside the study region on a given day. If the entire extent of the study
region were dominated by identical lit flare stacks burning identical natural
gas, then in atmospherically still conditions we would expect that the
estimated value of $\beta_{t}$ would be close to the methane-nitrogen dioxide
ratio expected from the balance of combustion reactions [31, 33, 32, 34].
Conversely, we would expect an estimated value of $\beta_{t}$ of zero in a
region void of sources of anthropogenic emission. The fact that we do not see
a positive correlation between $\beta_{t}$ and lit flare stack count could
suggest that there are sources of methane in or around the study region that
are not contributing nitrogen dioxide, such as well head completions, leaks,
compressor stations or unlit flare stacks. Some flare stacks may be operating
in an abnormal manner [38, 39, 40], and recent work suggests that super-
emitting individual flares may be accountable for the majority of methane
emissions in the Permian Basin [13]. We also find no observed seasonality in
the time series of $\beta_{t}$.

### Enhancement of spatial coverage

We augment daily observations over the study region by including model
predictions at locations where direct observations are not available from
TROPOMI. Predictions are only made from co-located TROPOMI observations of
nitrogen dioxide with an associated quality assurance value greater than or
equal to 0.75. An example of this procedure is shown in Fig. 5. Prior to
augmentation using predictions from the model, we calculate that the mean
value of daily pixel coverage within the study region for the year 2019 is
$16\%\pm 13\%$. After adding in model predictions at appropriate pixel
locations, the mean value of daily pixel coverage rises to $88\%\pm 18\%$. A
time series demonstrating this rise in spatial coverage is shown in panel a of
Fig. 6. Increasing the spatial coverage over regions of interest like the
Permain Basin can allow for the aggregation of data on shorter timescales than
previously used in performing methane emission estimates, which will be key
for accurately monitoring greenhouse gasses.

### Observed mass of methane hotspots

By subtracting a nominal background level of methane from TROPOMI observations
and model predictions, we can calculate the above-background mass of methane
contained above the study region on a given day. We use the global monthly
marine mean surface value of methane as a far-field reference background [2].

We convert all methane columns within the study region on each day to above-
background masses and integrate over the spatial extent to calculate the
above-background mass. This involves converting all values of methane (which
either returned from TROPOMI or predicted from the model as dry-air column
averaged mixing ratios) into column densities, for which we need a grid of
dry-air column densities at all pixel locations. We calculate a grid of dry
air column density to convert both TROPOMI observations and model predictions
to masses in a consistent manner [41], since dry air column densities are only
provided with the TROPOMI Level 2 CH4 Data Product at the location of observed
pixels [22]. Values of dry air column density from our grid are correlated
with values returned with the TROPOMI Level 2 $\mathrm{CH}_{4}$ Data Product
with a Pearson R of 0.97, and have a mean residual of -1,242.48 kg m-2 with
standard deviation 1,836.32 kg m-2. We use the root mean square deviation of
2,217.16 kg m-2 for error propagation in the calculation of mass of methane
contained with a pixel in the study region.

We find that over the course of the year 2019 the average daily above-
background mass is $0.89\pm 1.50$ kilotonnes. We re-calculate this quantity
including predictions from the model and find a new daily mean of above-
background mass of $4.29\pm 5.08$ kilotonnes. These two quantities are plotted
in panel c of Fig. 6. The choice of a reference background from the far field
does allow for negative values of above-background methane levels, but the
point of the calculation is to demonstrate the affect of including
predictions, and so the choice of background isn’t strictly important.
Previous work has demonstrated that time-integrated TROPOMI observations of
methane can be used to create top-down estimates of methane emission rates
[21]. While we do not currently carry out any inverse analysis, we demonstrate
that the inclusion of predictions of methane from co-located observations of
nitrogen dioxide increases the above-background value of methane loading in
the study region by kilotonnes per day on average, which could indicate that
the inclusion of predictions reveals methane hotspots that were previously
unobserved. By increasing the spatial coverage of observations, we can in
future work attempt to increase the temporal resolution of methane emission
rates in the Permian.

## Discussion

The remote sensing of atmospheric methane from space is crucial for mitigating
anthropogenic climate change. Using the Permian Basin as a case study, we show
that observations of nitrogen dioxide can be used as predictors for methane at
locations for which meteorological factors prevent TROPOMI from making direct
observations of methane from space. We validate our Bayesian model using a
variety of metrics and examine its predictive ability against withheld
observations. The lack of correlation between the daily slope parameter of our
model and the number of active flare stacks in the study region on each day
suggests that there are sources contributing methane without associated
nitrogen dioxide, which could be indicative of the presence of unlit or
malfunctioning flares. When using predicted values of methane observations to
augment daily observations over the Permian, we find that spatial coverage is
increased to nearly full coverage on most days, and the observed mass of
potential methane hotspots is increased on average by kilotonnes per day.
Using a similar algorithm to the work shown here could aid in the inversion of
satellite observations and estimation of methane emission rates on timescales
much shorter than previously achieved.

## Methods

For each day in 2019 we retrieved TROPOMI observations of methane and nitrogen
dioxide over the study region as shown in Fig. 1. Before conducting any
analysis, we reduce the data by ignoring pixels that do not pass the
recommended quality assurance (qa) value. For the TROPOMI Level 2
$\mathrm{CH}_{4}$ Data Product, the threshold qa factor for recommended usage
is 0.5 or greater, and for the TROPOMI Level 2 $\mathrm{NO}_{2}$ Data Product,
the threshold qa factor for recommended usage is 0.75 or greater. We then
classified individual days as being either rich or poor in remaining data,
with data-rich days being those with $N\geq 100$ co-located observations of
methane and nitrogen dioxide in the study region correlated with Pearson’s
$R\geq 0.4$, and data-poor days being all other days with $N\geq 2$. Co-
located observations are determined by linearly interpolating the grid of
$\mathrm{CH}_{4}^{\mathrm{obs}}$ onto the grid of
$\mathrm{NO}_{2}^{\mathrm{obs}}$. We leave the TROPOMI Level 2
$\mathrm{CH}_{4}$ Data Product in units of ppbv, but we scale the TROPOMI
Level 2 $\mathrm{NO}_{2}$ Data Product from mol m-2 to mmol m-2.

### Data-rich model

We next developed a fully Bayesian linear hierarchical model to fit solely to
observations from data-rich days. Data-rich days are indexed via
$t=1,2,...,D$, and co-located observations of methane and nitrogen dioxide
within the study region on day $t$ are indexed via $i=1,2,...,N$.

We relate observed values of methane and nitrogen dioxide to their true latent
values via

(2) $\displaystyle\mathrm{NO}_{2,i}^{\mathrm{obs}}$
$\displaystyle\sim\mathcal{N}\left(\mathrm{NO}_{2,i}^{\mathrm{true}},\,\sigma_{\mathrm{N},i}^{2}\right)$
(3) $\displaystyle\mathrm{CH}_{4,i}^{\mathrm{obs}}$
$\displaystyle\sim\mathcal{N}\left(\mathrm{CH}_{4,i}^{\mathrm{true}},\,\sigma_{\mathrm{C},i}^{2}\right)$

where $\sigma_{\mathrm{N}i}$ and $\sigma_{\mathrm{C}i}$ are the TROPOMI-
provided measurement standard deviations on $\mathrm{NO}_{2i}^{\mathrm{obs}}$
and $\mathrm{CH}_{4i}^{\mathrm{obs}}$ respectively. We also relate latent
values of methane to latent values of nitrogen dioxide via

(4)
$\mathrm{CH}_{4,i}^{\mathrm{true}}\sim\mathcal{N}\left(\alpha_{t}+\beta_{t}\,\mathrm{NO}_{2,i}^{\mathrm{true}},\,\gamma_{t}^{2}\right)$

where we have now introduced the daily model parameters $\alpha_{t}$,
$\beta_{t}$ and $\gamma_{t}$. On a given day $t$, $\alpha_{t}$ and $\beta_{t}$
are respectively the y-intercept and slope of the line of best fit relating
methane to nitrogen dioxide, while $\gamma_{t}$ is the standard deviation of
the scatter around the mean relation.

We stipulate an improper flat prior on latent values of nitrogen dioxide
observations in order to combine equations (2), (3) and (4) into the single
model equation

(5)
$\mathrm{CH}_{4,i}^{\mathrm{obs}}\sim\mathcal{N}\left(\alpha_{t}+\beta_{t}\,\mathrm{NO}_{2,i}^{\mathrm{obs}},\,\beta_{t}^{2}\,\sigma_{\mathrm{N},i}^{2}+\gamma_{t}^{2}+\sigma_{\mathrm{C},i}^{2}\right).$

Writing this equation in this way is desirable because it relates the observed
pixel value of methane to the observed pixel value of nitrogen dioxide
entirely in terms of the associated pixel precisions and the by-day model
parameters $\alpha_{t}$, $\beta_{t}$ and $\gamma_{t}$.

We add hyperparameters to our model by including a multivariate prior
distribution for $\alpha_{t}$ and $\beta_{t}$, shown in equation (1). We
include this multivariate prior to allow for the possibilty of learning the
extent to which $\alpha_{t}$ and $\beta_{t}$ are correlated. Although we
believe it is likely that $\alpha_{t}$ and $\beta_{t}$ are negatively
correlated, we do not encode this belief in the prior. We instead assume a
uniform flat prior on the domain $\left(-\infty,\;\infty\right)$ on each of
$\mu_{\alpha}$, $\mu_{\beta}$, $\Sigma_{12}$ and $\Sigma_{21}$, and assume a
uniform flat prior on the domain $\left(0,\;\infty\right)$ on both
$\Sigma_{11}$ and $\Sigma_{22}$. This lets us learn the extent to which
$\alpha_{t}$ and $\beta_{t}$ are correlated entirely from the posterior
inference. Equations (5) and (1) are the likelihood and prior of our data-rich
model respectively. We write our model in Stan [42] via the interface
CmdStanPy [43] and sample the posterior of our model using Stan’s default
Markov Chain Monte Carlo algorithm NUTS (the No U-Turn Sampler, which is a
variant of Hamiltonian Monte Carlo). When fitting our model we specify the
algorithm to draw samples from the posterior using four separate Markov
chains, each with 500 burn-in iterations with a further 1000 retained draws,
combining for a total of 4000 draws from the posterior for each of our model
parameters. Specifying 1000 post burn-in draws per chain easily allows for a
sufficient ESS.

### Data-poor model

We developed a separate Bayesian model to fit to days that we identified as
being poor in data. Data-poor days were classified as those that were not
data-rich and have $N\geq 2$ co-located observations of methane and nitrogen
dioxide in the study region. The point of developing a second model to fit to
data-poor days after we have already fit a model to data-rich days is to
incorporate information learned from the data-rich model into the data-poor
model as an informed prior.

As in the data-rich model, the likelihood of the data-poor model is given by
equation (5), and the multivariate prior on values of $\alpha_{t}$ and
$\beta_{t}$ is given by equation (1). However, when fitting the hierarchical
model to all data-rich days simultaneously, we’ve learned information about
what sort of values the model hyperparameters “should” take. To account for
this information we’ve learned, we add a prior on the hyperparameters with

(6) $\displaystyle\begin{pmatrix}\mu_{\alpha}\\\ \mu_{\beta}\\\ \Sigma_{11}\\\
\Sigma_{12}\\\ \Sigma_{22}\\\ \end{pmatrix}$
$\displaystyle\sim\mathcal{N}\left(\theta,\,\Upsilon\right).$

We have $\theta$ as the mean vector of the 4000 posterior draws of
($\mu_{\alpha}$, $\mu_{\beta}$, $\Sigma_{11}$, $\Sigma_{12}$, $\Sigma_{22}$)
from fitting the data-rich model to data-rich days. $\Upsilon$ is the 5x5
covariance matrix of $\mu_{\alpha}$, $\mu_{\beta}$, $\Sigma_{11}$,
$\Sigma_{12}$ and $\Sigma_{22}$, constructed using their 4000 posterior draws
from the data-rich model. We fit the data-poor model to all data-poor days
again using Stan and the NUTS MCMC algorithm.

### Reparameterisation

When fitting the data-rich and data-poor models to actual observations, the
MCMC sampler was confronted with an abundance of divergent transitions. To
remove these transitions, we reparameterise our models into effectively
equivalent mathematical representations that present a simpler posterior
geometry [44, 42]. We do so making use of the Cholesky decompositions of our
covariance matrices [45].

In order to decrease fitting time of the data-rich model, we replace the per-
pixel precisions returned with the TROPOMI data products with daily averages,
defined by

(7) $\displaystyle\sigma_{\mathrm{C}t}$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\sigma_{\mathrm{C}i,t}$ (8)
$\displaystyle\sigma_{\mathrm{N}t}$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\sigma_{\mathrm{N}i,t}.$

This yields a new likelihood equation given by

(9)
$\mathrm{CH}_{4i}^{\mathrm{obs}}\sim\mathcal{N}\left(\alpha_{t}+\beta_{t}\,\mathrm{NO}_{2i}^{\mathrm{obs}},\,\beta_{t}^{2}\,\sigma_{\mathrm{N}t}^{2}+\gamma_{t}^{2}+\sigma_{\mathrm{C}t}^{2}\right).$

In order to determine how this affects the predictive accuracy of the data-
rich model, we fit a data-rich model using per-pixel precisions and a data-
rich model using daily averages of precision to data-rich days for the month
of January 2019, and estimate the difference between the two models’ expected
log pointwise predictive density (ELPD) using the widely available information
criterion (WAIC) [46, 47]. We estimate that the difference in ELPD of the two
models in this case is $6.92\pm 7.09$. As the difference in ELPD is within a
standard error of zero we continue using daily averages of TROPOMI pixel
precision in order to decrease the fitting time of our data-rich and data-poor
models.

### Predictions

We can predict $\mathrm{CH}_{4}^{\mathrm{pred}}$ and precision
$\sigma_{\mathrm{C}}^{\mathrm{pred}}$ from a fitted model using an observation
of nitrogen dioxide $\mathrm{NO}_{2}$ and its associated precision
$\sigma_{\mathrm{N}}$ as input. To obtain predictions, we first sample
$K=1000$ potential values of $\mathrm{CH}_{4}$ from the model values via
$\mathrm{CH}_{4,k}\sim\mathcal{N}\left(\alpha_{t,k}+\beta_{t,k}\,\mathrm{NO}_{2},\,\beta_{t,k}^{2}\,\sigma_{\mathrm{N}}^{2}+\gamma_{t,k}^{2}\right)$
where the subscript $k$ denotes a random draw with replacement of a set of
parameter values from one of the 4000 sets in the posterior chain. We then
have
$\mathrm{CH}_{4}^{\mathrm{pred}}=\frac{1}{K}\sum_{k=1}^{K}\mathrm{CH}_{4,k}$
and
$\sigma_{\mathrm{C}}^{\mathrm{pred}}=\sqrt{\frac{\sum_{k=1}^{K}\left(\mathrm{CH}_{4,k}-\mathrm{CH}_{4}^{\mathrm{pred}}\right)^{2}}{K}}$.

### Dropout testing

We perform dropout testing in order test the predictive ability of our model.
For each data-rich day $t$ we ignore a random selection of 20% of the co-
located observations of methane and nitrogen dioxide and fit the model to the
remaining 80% of the data. We use Stan to encode our model and sample the
posterior using NUTS, with four independent chains each with 500 burn-in draws
and 1000 retained sampled draws for a total of 4000 draws from the posterior.

After the data-rich model has been fit to the retained 80% of the data, we can
compare observed values of methane from the held-out subset of data to
predictions from the model and summarise how well the model has fit each day
using a reduced chi-squared statistic. The reduced chi-squared statistic on
day $t$ is calculated via

(10) $\displaystyle\chi_{\nu,t}^{2}$
$\displaystyle=\frac{\chi^{2}_{t}}{\nu_{t}}$ (11) $\displaystyle\chi^{2}_{t}$
$\displaystyle=\sum_{i=1}^{N}\frac{\left(\mathrm{CH}_{4,i}^{\mathrm{obs}}-\mathrm{CH}_{4,i}^{\mathrm{pred}}\right)^{2}}{\sigma_{\mathrm{C},i}^{2}+\sigma_{\mathrm{C},i}^{\mathrm{pred^{2}}}}$
(12) $\displaystyle\nu_{t}$ $\displaystyle=N-8,$

where again $N$ is the number of co-located observations of methane and
nitrogen dioxide in the study region on day $t$. After examining the resulting
distribution of calculated values of $\chi^{2}_{\nu,t}$, we fit the model to
the entire set of observations in the year 2019, without withholding any
subset of the data. This fitted model is the model that we use for predicting
values of methane for our final results.

### VIIRS

We next retrieved VIIRS Nightfire observations of lit methane flare stacks for
each day in 2019. Nightfire datasets provide the location coordinates,
estimated temperature and source size of identified flares [37]. We reduce the
daily datasets down to flare stacks identified within the study region
regardless of estimated temperature. Daily tallies of lit flare stacks were
collected for eventual comparison against time series of daily fitted model
parameters.

### ERA5

The TROPOMI $\mathrm{CH}_{4}$ Level 2 Data Product provides the dry-air column
average mixing ratio of methane in units of ppbv, defined as the ratio between
the column density of methane to the column density of dry air at each pixel
location [22]. To convert the dry-air column average mixing ratio of methane
back to a column density, we simply multiply by the value of the dry air
column density at the pixel location, which we calculate ourselves from ERA5
data. Although values of dry air column density are supplied at TROPOMI pixel
locations where a retrieval was successfully performed, we calculate our own
grid from ERA5 data in order to avoid any discrepancy between values of dry
air column density at the locations of “original” observed methane pixel
values and the locations of model predictions.

We began by retrieving hourly spatial grids of surface pressure
$P_{\mathrm{S}}$ and vertical column of water vapour
$VC_{\mathrm{H}_{2}\mathrm{O}}$ that cover the entire extent of the study
region for each day in 2019 [41]. $P_{S}$ is given in [Pa] and
$VC_{\mathrm{H}_{2}\mathrm{O}}$ given in [kg m-2]. From the value
$P_{\mathrm{S}}$ we can calculate the total column density of air
$VC_{\mathrm{air}}$ via $VC_{\mathrm{air}}=\frac{P_{\mathrm{S}}}{\mathrm{g}}$
where $\mathrm{g}=9.8$ m s-2 is gravity. We then calculate the total vertical
column of dry air $VC_{\mathrm{dry}\;\mathrm{air}}$ via
$VC_{\mathrm{dry}\;\mathrm{air}}=VC_{\mathrm{air}}-VC_{\mathrm{H}_{2}\mathrm{O}}$.

When calculating $VC_{\mathrm{dry}\;\mathrm{air}}$ at the location of either
an observation or model prediction of methane, we interpolate the hourly
spatial grids of ERA5 data to the grid of TROPOMI methane observations, and
then linearly interpolate in time between the two adjacent hourly grids
according to the pixel scanline [22].

### Calculation of methane loading

We were interested in observing how the application of our predictive
algorithm altered the amount of observed methane over the Permian Basin for
the year 2019. By virtue of observing more spatial area, the algorithm would
result in an increase of total observed mass of methane To counteract this
effect, we subtract a nominal background level of methane from each pixel and
then describe how the above-background mass of methane in the study region
changes after the application of our predictive algorithm. We linearly
interpolate the global monthly marine mean surface value of methane as a far-
field reference background, provided by the Global Monitoring Laboratory (GML)
at the National Oceanic and Atmospheric Association (NOAA) [2].

The number of moles of methane above the GML/NOAA background contained in a
pixel is calculated via
$\mathrm{mol}\;\mathrm{CH}_{4}=\left(\mathrm{CH}_{4}-\mathrm{B}\right)\times
10^{9}\times VC_{\mathrm{dry}\;\mathrm{air}}\times A$ where $A$ is the pixel
area in square metres, B is the mean global surface value of methane in ppbv,
and $\mathrm{CH}_{4}$ is either the TROPOMI-observed or model-predicted pixel
value of $\mathrm{CH}_{4}$ in ppbv. We then convert from moles to tonnes by
multiplying by the molar mass of methane and scaling into tonnes. The
precision on the moles of above-background methane contained in a pixel is
calculated via
$\sigma_{\mathrm{mol}\;\mathrm{CH}_{4}}=\mathrm{mol}\;\mathrm{CH}_{4}\sqrt{\left(\frac{\sqrt{\sigma_{\mathrm{C}}^{2}+\sigma_{\mathrm{B}}^{2}}}{\mathrm{CH}_{4}-\mathrm{B}}\right)^{2}+\left(\frac{\mathrm{RMSE}_{VC}}{VC_{\mathrm{dry\;air}}}\right)^{2}}$
where $\sigma_{\mathrm{C}}$ is the standard deviation uncertainty on the value
of $\mathrm{CH}_{4}$ and $\sigma_{\mathrm{B}}$ is the standard deviation
uncertainty on the value of B. We take $\mathrm{RMSE}_{VC}$ to be the root
mean squared error on the residuals between our ERA5-calculated value of
$VC_{\mathrm{dry\;air}}$ and the value returned with the TROPOMI Level 2
$\mathrm{CH}_{4}$ Data Product at all pixel locations possible on all data-
rich days in the year 2019.

We calculate the total methane mass above background on a given day in the
study region via
$\mathrm{total\;mol\;CH}_{4}=\sum_{i=1}^{N}\mathrm{mol\;CH}_{4i}$ where $N$ is
the number of pixels in the study region where we have either
$\mathrm{CH}_{4}^{\mathrm{obs}}$ or $\mathrm{CH}_{4}^{\mathrm{pred}}$. In
order to obtain a conservative estimate, we only include predicted values of
methane $\mathrm{CH}_{4}^{\mathrm{pred}}$ in the summation that are predicted
from values of $\mathrm{NO_{2}}$ with $\mathrm{NO}_{2}/\sigma_{N}>2$. The
precision of total mol CH4 is given by
$\sigma_{\mathrm{total\;mol\;CH}_{4}}=\sqrt{\sum_{i=1}^{N}\left(\sigma_{\mathrm{mol\;CH}_{4},i}\right)^{2}}$.

## Data Availability

The data that support this study are available from the Copernicus Open Access
Hub (for TROPOMI $\mathrm{CH}_{4}$ and $\mathrm{NO}_{2}$ Level 2 Products),
the Copernicus Climate Change Service (for ERA5 reanalysis data), the Earth
Observation Group at the Colorado School of Mines (for VIIRS Nightfire) and
the Global Monitoring Laboratory at the National Oceanic and Atmospheric
Association (for mean global marine surface methane).

## Code Availability

All code used in this study to analyse data and generate figures can be found
at the corresponding author’s Github repository titled TROPOMI_LHM (“Linear
Hierarchical Model”).

## References

  * [1] G. Myhre et al.  “Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change”, 2013
  * [2] Ed Dlugokencky  GML/NOAA  URL: https://gml.noaa.gov/ccgg/trends_ch4/
  * [3] P.. Arias et al.  “Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change”, 2021
  * [4] Pieter Tans and Ralph Keeling Keeling  GML/NOAA  URL: gml.noaa.gov/ccgg/trends/
  * [5] Paul Balcombe, Jamie F. Speirs, Nigel P. Brandon and Adam D. Hawkes  “Methane emissions: choosing the right climate metric and time horizon”  In _Environmental Science: Processes & Impacts_ 20.10, 2018, pp. 1323–1339  DOI: 10.1039/C8EM00414E
  * [6] D.. Jacob et al.  “Satellite observations of atmospheric methane and their value for quantifying methane emissions”  In _Atmospheric Chemistry and Physics_ 16.22, 2016, pp. 14371–14396  DOI: 10.5194/acp-16-14371-2016
  * [7] J. Palmer  “The pollution detectives”  In _Nature_ 577.7791, 2020, pp. 464–466  DOI: https://doi.org/10.1038/d41586-020-00110-8
  * [8] United Nations Environment Programme (UNEP)  “Mineral Methane Initiative OGMP2.0 Framework”, 2020  URL: https://www.ccacoalition.org/en/resources/oil-and-gas-methane-partnership-ogmp-20-framework
  * [9] H. Bovensmann et al.  “SCIAMACHY: Mission Objectives and Measurement Modes”  In _Journal of the Atmospheric Sciences_ 56.2, 1999, pp. 127–150  DOI: 10.1175/1520-0469(1999)056¡0127:SMOAMM¿2.0.CO;2
  * [10] Akihiko Kuze, Hiroshi Suto, Masakatsu Nakajima and Takashi Hamazaki  “Thermal and near infrared sensor for carbon observation Fourier-transform spectrometer on the Greenhouse Gases Observing Satellite for greenhouse gases monitoring”  In _Applied Optics_ 48.35, 2009, pp. 6716–6733  DOI: 10.1364/AO.48.006716
  * [11] R. Glumb, G. Davis and C. Lietzke  “The TANSO-FTS-2 instrument for the GOSAT-2 greenhouse gas monitoring mission”  In _2014 IEEE Geoscience and Remote Sensing Symposium_ , pp. 1238–1240  DOI: 10.1109/IGARSS.2014.6946656
  * [12] Y.. Liu et al.  “The Advanced Hyperspectral Imager: Aboard China’s GaoFen-5 Satellite”  In _IEEE Geoscience and Remote Sensing Magazine_ 7.4, 2019, pp. 23–32  DOI: 10.1109/MGRS.2019.2927687
  * [13] Itziar Irakulis-Loitxate et al.  “Satellite-based survey of extreme methane emissions in the Permian basin”  In _Science Advances_ 7.27, pp. eabf4507  DOI: 10.1126/sciadv.abf4507
  * [14] J.. Veefkind et al.  “TROPOMI on the ESA Sentinel-5 Precursor: A GMES mission for global observations of the atmospheric composition for climate, air quality and ozone layer applications”  In _Remote Sensing of Environment_ 120, 2012, pp. 70–83  DOI: https://doi.org/10.1016/j.rse.2011.09.027
  * [15] André Butz et al.  “TROPOMI aboard Sentinel-5 Precursor: Prospective performance of CH4 retrievals for aerosol and cirrus loaded atmospheres”  In _Rem. Sens. Env., Sentinel-5 special issue_ 120, 2011  DOI: 10.1016/j.rse.2011.05.030
  * [16] Haili Hu et al.  “Toward Global Mapping of Methane With TROPOMI: First Results and Intersatellite Comparison to GOSAT”  In _Geophysical Research Letters_ 45.8, 2018, pp. 3682–3689  DOI: https://doi.org/10.1002/2018GL077259
  * [17] Daniel Clery  “Europe to lead in monitoring carbon from space”  In _Science_ 366, 2019, pp. 1176–1177  DOI: 10.1126/science.366.6470.1176
  * [18] Ramón A. Alvarez et al.  “Assessment of methane emissions from the U.S. oil and gas supply chain”  In _Science_ 361, 2018, pp. 186–188  DOI: 10.1126/science.aar7204
  * [19] Oliver Schneising et al.  “Remote sensing of fugitive methane emissions from oil and gas production in North American tight geologic formations”  In _Earth’s Future_ 2.10, 2014, pp. 548–558  DOI: https://doi.org/10.1002/2014EF000265
  * [20] A.. Turner et al.  “A large increase in U.S. methane emissions over the past decade inferred from satellite data and surface observations”  In _Geophysical Research Letters_ 43.5, 2016, pp. 2218–2224  DOI: https://doi.org/10.1002/2016GL067987
  * [21] Yuzhong Zhang et al.  “Quantifying methane emissions from the largest oil-producing basin in the United States from space”  In _Science Advances_ 6.17, 2020, pp. eaaz5120  DOI: 10.1126/sciadv.aaz5120
  * [22] O. Hasekamp et al.  “Algorithm Theoretical Baseline Document for Sentinel-5 Precursor Methane Retrieval”, 2021
  * [23] J…. Geffen, H.. Eskes, K.. Boersma and J.. Veefkind  “TROPOMI ATBD of the total and tropospheric NO2 data products”, 2021
  * [24] B.. Duncan et al.  “A space-based, high-resolution view of notable changes in urban NOx pollution around the world (2005-2014)”  In _Journal of Geophysical Research-Atmospheres_ 121.2, 2016, pp. 976–996  DOI: 10.1002/2015jd024121
  * [25] Christopher D. Elvidge et al.  “The potential role of natural gas flaring in meeting greenhouse gas mitigation targets”  In _Energy Strategy Reviews_ 20, 2018, pp. 156–162  DOI: https://doi.org/10.1016/j.esr.2017.12.012
  * [26] Christopher D. Elvidge et al.  “Methods for Global Survey of Natural Gas Flaring from Visible Infrared Imaging Radiometer Suite Data”  In _Energies_ 9.1, 2016, pp. 14  URL: https://www.mdpi.com/1996-1073/9/1/14
  * [27] Joost A. Gouw et al.  “Daily Satellite Observations of Methane from Oil and Gas Production Regions in the United States”  In _Scientific Reports_ 10.1, 2020, pp. 1379  DOI: 10.1038/s41598-020-57678-4
  * [28] “Permian Region Drilling Productivity Report”, 2021
  * [29] Michael Betancourt  “Diagnosing Suboptimal Cotangent Disintegrations in Hamiltonian Monte Carlo”, 2016
  * [30] Aki Vehtari et al.  “Rank-Normalization, Folding, and Localization: An Improved $\widehat{R}$ for Assessing Convergence of MCMC (with Discussion)”  In _Bayesian Analysis_ 16.2, 2021, pp. 667–71852  URL: https://doi.org/10.1214/20-BA1221
  * [31] J.. Sonibare and F.. Akeredolu  “A theoretical prediction of non-methane gaseous emissions from natural gas combustion”  In _Energy Policy_ 32.14, 2004, pp. 1653–1665  DOI: https://doi.org/10.1016/j.enpol.2004.02.008
  * [32] K. Deetz and B. Vogel  “Development of a new gas-flaring emission dataset for southern West Africa”  In _Geosci. Model Dev._ 10.4, 2017, pp. 1607–1620  DOI: 10.5194/gmd-10-1607-2017
  * [33] O. Ismail and G. Umukoro  “Modelling combustion reactions for gas flaring and its resulting emissions”  In _Journal of King Saud University - Engineering Sciences_ 28.2, 2016, pp. 130–140  DOI: https://doi.org/10.1016/j.jksues.2014.02.003
  * [34] G. Umukoro and O. Ismail  “Modelling emissions from natural gas flaring”  In _Journal of King Saud University - Engineering Sciences_ 29.2, 2017, pp. 178–182  DOI: https://doi.org/10.1016/j.jksues.2015.08.001
  * [35] Jennifer A Hutcheon, Arnaud Chiolero and James A Hanley  “Random measurement error and regression dilution bias”  In _BMJ_ 340, 2010, pp. c2289  DOI: 10.1136/bmj.c2289
  * [36] Stefano Andreon and Merilee Hurn  “Measurement errors and scaling relations in astrophysics: a review”  In _Statistical Analysis and Data Mining: The ASA Data Science Journal_ 9, 2013, pp. 15  DOI: 10.1002/sam.11173
  * [37] Christopher D. Elvidge, Mikhail Zhizhin, Feng-Chi Hsu and Kimberly E. Baugh  “VIIRS Nightfire: Satellite Pyrometry at Night”  In _Remote Sensing_ 5.9, 2013, pp. 4423–4449  URL: https://www.mdpi.com/2072-4292/5/9/4423
  * [38] C.. Rella et al.  “Measuring Emissions from Oil and Natural Gas Well Pads Using the Mobile Flux Plane Technique”  In _Environmental Science & Technology_ 49.7, 2015, pp. 4742–4748  DOI: 10.1021/acs.est.5b00099
  * [39] Mark Omara et al.  “Methane Emissions from Conventional and Unconventional Natural Gas Production Sites in the Marcellus Shale Basin”  In _Environmental Science & Technology_ 50.4, 2016, pp. 2099–2107  DOI: 10.1021/acs.est.5b05503
  * [40] Anna M. Robertson et al.  “Variation in Methane Emission Rates from Well Pads in Four Oil and Gas Basins with Contrasting Production Volumes and Compositions”  In _Environmental Science & Technology_ 51.15, 2017, pp. 8832–8840  DOI: 10.1021/acs.est.7b00571
  * [41] Hans Hersbach et al.  “The ERA5 global reanalysis”  In _Quarterly Journal of the Royal Meteorological Society_ 146.730, 2020, pp. 1999–2049  DOI: https://doi.org/10.1002/qj.3803
  * [42] Stan Development Team  “Stan Modeling Language Users Guide and Reference Manual”, 2021  URL: https://mc-stan.org
  * [43] Stan Development Team  “CmdStanPy”, 2021  URL: https://pypi.org/project/cmdstanpy
  * [44] Omiros Papaspiliopoulos, Gareth O. Roberts and Martin Sköld  “A General Framework for the Parametrization of Hierarchical Models”  In _Statistical Science_ 22.1 Institute of Mathematical Statistics, 2007, pp. 59–73  DOI: 10.1214/088342307000000014
  * [45] M. Betancourt and Mark Girolami  “Hamiltonian Monte Carlo for Hierarchical Models”, 2013  DOI: 10.1201/b18502-5
  * [46] Sumio Watanabe  “Asymptotic Equivalence of Bayes Cross Validation and Widely Applicable Information Criterion in Singular Learning Theory”  In _J. Mach. Learn. Res._ 11, 2010, pp. 3571–3594
  * [47] Aki Vehtari, Andrew Gelman and Jonah Gabry  “Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC”  In _Statistics and Computing_ 27.5, 2017, pp. 1413–1432  DOI: 10.1007/s11222-016-9696-4

## Acknowledgements

C.R. acknowledges financial support from Shell Research Ltd. through the
Cambridge Centre for Doctoral Training in Data Intensive Science.

## Author Contributions

C.R. retrieved all data, wrote all code and conducted the analysis of this
work. O.S. supervised the project and guided the analysis. K.M., P.J. and M.J.
guided the construction of the Bayesian model. R.I. and B.H. contributed
discussions of the data, satellites and industrial context.

## Competing Interests

We report no competing interests.

(a)  Figure 1. TROPOMI observations of methane (a) and nitrogen dioxide (b)
over the Permian Basin on January 31st, 2019. Only methane and nitrogen
dioxide pixels which pass the recommended quality assurance threshold are
shown [22] and all other pixels are masked. The study region is shown by the
red dashed rectangle in b, superimposed over county lines which are shown in
grey solid lines in both panels. Also shown with green triangles in b are the
locations of lit methane flare stacks identified by VIIRS Nightfire [37],
demonstrating their co-location with atmospheric overabundances of methane and
nitrogen dioxide.

(a)  Figure 2. A graphical representation of how $\alpha_{t}$ and $\beta_{t}$
are estimated (a), with a demonstration of results accross all data-rich days
in 2019 (b). a Co-located TROPOMI observations of methane column average dry
air mixing ratio plotted against nitrogen dioxide column densities within the
study region on January 31st, 2019, shown in red. Error bars on the
observations (shown in blue) are the single-sigma precisions provided in the
TROPOMI Level 2 Data Products. Superimposed over the scatterplot of
observations we plot a random 1,000 selections of $(\alpha_{t},\;\beta_{t})$
pairs from the 4,000 sampled draws from the posterior of the fitted data-rich
model, shown in lime green. In the bottom right corner we show median
estimated posterior values of $\alpha_{t}$ [ppbv] and $\beta_{t}$ [ppbv /
(mmol m-2)] with $\pm$ indicating the extent of the 68% central credible
interval. b Median estimated posterior values of $\alpha_{t}$ and $\beta_{t}$
(shown in red) across all data-rich days with their estimated 68% central
credible intervals (shown in blue). Shown underneath the plotted values of
$\alpha_{t}$ and $\beta_{t}$ is a red ellipse constructed from the median over
random samples drawn from the joint posterior chain of $\mu_{\alpha}$,
$\mu_{\beta}$, $\rho$, $\sigma_{\alpha}$ and $\sigma_{\beta}$. This ellipse
indicates the information that is eventually supplied to the data-poor model
as prior information for $\alpha_{t}$ and $\beta_{t}$, although in the data-
poor model we also account for the spread of posterior values of the
hyperparameters instead of providing median values without any uncertainty.

(a)  Figure 3. Results of the dropout testing. a Predictions of methane
$\mathrm{CH}_{4}^{\mathrm{pred}}$ plotted against actual TROPOMI observations
$\mathrm{CH}_{4}^{\mathrm{obs}}$. Predictions come from a model that was
fitted without using the observed values shown on the x-axis. We plot in red
the line $y=x$ to demonstrate the effect that uncertainty on nitrogen dioxide
observation has on our model; regression dilution results in slight
underestimation of methane at high values of nitrogen dioxide and
overestimation at low values of nitrogen dioxide, though there is still good
agreement overall. b A distribution of reduced chi-squared values for data-
rich days in the year 2019.

(a)  Figure 4. A demonstration that the number of total lit methane flares in
the study region identified by VIIRS does not appear to correlate with
estimated values of $\beta_{t}$. b Median estimated posterior values of
$\beta_{t}$ on data-rich days as a time series with 68% central credible
regions as error bars. Also shown in this time series are the numbers of
active lit methane flares identified by VIIRS Nightfire on data-rich dates,
shown with red x’s. In panel a we plot the same values and error bars on
$\beta_{t}$ as a function of the identified number of lit methane flares on
that day.

(a)  Figure 5. A comparison between original (a) and augmented (b) TROPOMI
observations of methane columns. a The same TROPOMI observation of methane
columns as in panel a in Fig. 1. b The same observations from a, augmented
with predictions from the fitted model at all locations where we have a
TROPOMI observation of nitrogen dioxide that passes the recommended quality
assurance thresholds.

(a)  Figure 6. (Continued on the following page.)

An examination of the effects of including model predictions of methane
columns in addition to original TROPOMI observations. a-c In each panel we
plot two time series. Plotted in blue are quantities calculated purely from
TROPOMI observations that pass the recommended quality assurance threshold.
Plotted in red are the same quantities but calculated with the inclusion of
predictions of methane from the fitted model at pixel locations in the study
region where direct TROPOMI observations of nitrogen dioxide are available but
no direct observations of methane are available. For both time series in each
panel, full circles indicate data-rich days and open circles indicate data-
poor days. a Percentage of usable pixels in the study region, demonstrating
that the application of the predictive algorithm augments spatial coverage to
nearly 100% of the study region on most days. b Median observed pixel value in
the study region, demonstrating that the inclusion of predictions does not
skew the median pixel value in the study region to higher or lower values away
from the original median observed pixel value. c Total observed above-
background mass of methane over the study region (in reference to the NOAA
background). We calculate uncertainty on the quantity plotted in c, but error
bars would so narrow as to not be visible when plotted on this scale. Panel c
demonstrates that the inclusion of predictions could potentially account for
extra kilotonnes of excess methane over the Permian Basin that would nominally
be unobserved.