Enhanced monitoring of atmospheric methane from space over the Permian basin with hierarchical Bayesian inference

TROPOMI provides one observation of nitrogen dioxide and methane 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 β_t and intercept parameter α_t , plus another daily parameter to account for intrinsic scatter around the linear trend. We fit our model to a year's observations of methane and nitrogen dioxide over the Permian basin from 2019, using co-located observations contained within a marked study region shown in figure 1. The daily parameter α_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 nitrogen dioxide. The daily parameter β_t has units of ppbv (mmol m⁻²)⁻¹ as we scale all observations of nitrogen dioxide columns into mmol m⁻².

We specify a multivariate normal hierarchical prior on values of α_t and β_t such that

$$\begin{align} \begin{pmatrix} \alpha_t \\ \beta_t \end{pmatrix} &\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) \end{align} \tag{ 1 }$$

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 α_t and β_t , with $\rho \in [-1,\;1]$. The hyperparameter µ_α roughly corresponds to the average background level of methane in the study region across the year 2019, and the hyperparameter σ_α is the deviation that controls the spread of α_t around µ_α , with $\sigma_\alpha \gt 0$, and similarly for µ_β and σ_β .

We fit two models, each to a specific subset 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. Note that the latter model can be fit to data very well and is referred to as the 'data-poor' model because it is fit to data from days that have relatively few co-located TROPOMI observations of methane and nitrogen dioxide. The difference in construction 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. The effect of requiring that data-rich days be those with both large amounts of co-located observations of nitrogen dioxide and methane as well as a high degree of correlation is that model predictions of methane from the fitted data-poor model may be more dominated by the prior when observations are extremely sparse, but we find that this is typically not the case. 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 [30, 31]. 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, indicative of sufficient mixing in both cases.

We examine the joint posterior distribution of the hyperparameters of the data-rich model, and represent median values of µ_α , µ_β , $\rho$, σ_α and σ_β in panel (b) of figure 2. We find that the correlation $\rho$ between α_t and β_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 [32–35], 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 model parameters α_t and β_t are correlated, and leverage this information as a prior of appropriate strength in the data-poor model.

Zoom In Zoom Out Reset image size

Other model hyperparameters that we discuss at this stage are µ_α and σ_α . These two model parameters have physical meaning, where µ_α can be thought of as the 'mean' regional background of methane in our study region for the year 2019, and σ_α as the extent to which the daily intercept α varies around µ_α from day to day. We find µ_α to have a median value of 1861.32 ppbv with a 68% central credible interval of (1859.52, 1863.15), and σ_α to have a median value of 14.44 ppbv with a 68% central credible interval of (13.25, 15.73).

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 figure 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 figure 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]. Recent work has derived the observation standard deviation over the Permian to be 11 ppbv [21]. In an extreme case where all 11 ppbv is independent of our model error, the standard error of our model predictions would be roughly $\sqrt{12^2 + 11^2} \approx 16$ ppbv. Panel (a) of figure 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 regression dilution caused by the relatively high uncertainty of the TROPOMI observations of nitrogen dioxide [36, 37], seen in panel (a) of figure 2. Correcting for regression dilution is not necessary in predictive modelling scenarios. Validation studies have shown that the TROPOMI nitrogen dioxide data product can be biased low by as much as 50% over highly polluted regions when compared to ground-based observations [38]. Since we fit our models using nitrogen dioxide TROPOMI observations that are not bias corrected, the resulting predicted methane abundances will not be degraded when non-bias-corrected TROPOMI nitrogen dioxide observations are used as input for the methane predictions.

Zoom In Zoom Out Reset image size

We also assess the predictive ability of our models against the weighting function modified differential optical absorption spectroscopy (WFMD) CH₄ data product [39]. This data product has much better spatial coverage than the operational Sentinal-5P CH₄ data product that we fit our models to. We find that predictions of methane from our fitted model correlate nearly as well with co-located observations of the Permian basin from the WFMD CH₄ data product as the 'original' TROPOMI methane observations (figure 4). TROPOMI observations of methane exhibit a standard deviation of 15.5 ppbv around their line of best fit with co-located observations from the WFMD CH₄ data product red with a correlation of Pearson R = 0.73, whilst predictions from our model exhibit a standard deviation of 15.8 ppbv around the same line red with a Pearson R = 0.58.

Zoom In Zoom Out Reset image size

After re-fitting the data-rich model to 100% of available observations we examine posterior estimates of daily model parameters as time series. We plot median values of the slope parameter β_t in figure 5 along with a time series of active flare stacks in the study region on data-rich days, identified from VIIRS Nightfire [40]. We do not find any significant correlation between estimated β_t and the total number of identified flare stacks inside the study region on a given day, which is not unexpected as recent work suggests that super-emitting individual flares may be accountable for the majority of methane emissions in the Permian basin [13]. However, we do find that higher values of β_t tend to occur in the summer months. We investigate this further by examining the seasonality of nitrogen dioxide in the study region, shown in of figure 6(c). Nitrogen dioxide has a shorter atmospheric lifetime in the warmer summer months, which could explain the higher inferred values of β_t shown in figure 5. However, in figure 6 we also investigate the correlation between the average value of nitrogen dioxide in the study region and number of active flares, and find that nitrogen dioxide correlates with number of flares to a significant level in the warmer summer months. Future work is needed to determine the origin of this correlation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 figure 7. 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\%\pm13\%$. After adding in model predictions at appropriate pixel locations, the mean value of daily pixel coverage rises to $88\%\pm18\%$. A time series demonstrating this rise in spatial coverage is shown in figure 10(a). We find that the spatial coverage of our augmented TROPOMI observations at least matches and usually exceeds that of the WFMD CH₄ data product (figure 8). Increasing the spatial coverage over regions of interest like the Permian 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 gas emissions in near-real time.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The cities of Odessa and Midland are contained within the study region and each have populations that exceed 100 000 people. Cities are known to emit large quantities of nitrogen dioxide [25], and may emit co-located methane at a rate that differs from the rural surrounding oil and gas producing regions. We here examine the extent to which including these cities in the study region may affect predictions of methane in non urban areas. To do this, we define a new sub-region within the study region that contains the two cities. Figure 9 shows comparative distributions of various TROPOMI observations and model predictions, segregated according to whether or not they are located in the urban sub-region containing the cities. Using a Kolmogorov–Smirnov test, we determine these pairs of distributions to be significantly different from one another [41]. We calculate the fractional difference between the mean value of observed nitrogen dioxide columns over the urban sub-region and the mean value of observed nitrogen dioxide columns over the surrounding rural areas and find that they differ by 31.28%. This indicates that nitrogen dioxide emissions tend to be higher on average over the city sub-region than over the surrounding rural regions, and that as a result it may be the case that methane emissions are slightly over-predicted by our model over the city sub-region. However, the difference between the mean values of predicted and observed methane over the urban sub-region is less than one percent. The same is true for the fractional difference between the mean values of predictions and observations of methane over the rural surroundings. Although it is beyond the scope of this work to fully investigate this result, it may be the case that our model under-predicts methane emissions in the surrounding rural areas of the study region where methane sources like livestock are not directly correlated with nitrogen dioxide emissions. We note that the majority of locations of data in our sample correspond to rural areas, and section 2.2 demonstrated that in general our model predictions agree well with the data. In future work, this methodology could be applied to smaller study areas where infrastructure remains more homogenous, though this would come at a cost of less data.

Zoom In Zoom Out Reset image size

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 [42], since dry air column densities are only provided with the TROPOMI Level 2 CH₄ 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 −1242 kg m⁻² with standard deviation 1836 kg m⁻². We use the root mean square deviation of 2217 kg m⁻² 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 of methane is $0.9\pm1.5$ kilotonnes. We re-calculate this quantity including predictions from the model and find a new daily mean of above-background methane mass of $4.3\pm5.1$ kilotonnes. These two quantities are plotted in panel (c) of figure 10. 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 is not 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.

Zoom In Zoom Out Reset image size

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,\ldots,D$, and co-located observations of methane and nitrogen dioxide within the study region on day t are indexed via $i = 1,2,\ldots,N$.

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

$$\begin{align} \mathrm{NO}_{2,i}^{\mathrm{obs}} &\sim \mathcal{N}\left(\mathrm{NO}_{2,i}^{\mathrm{true}},\,\sigma_{\mathrm{N},i}^2\right) \end{align} \tag{ 2 }$$

$$\begin{align} \mathrm{CH}_{4,i}^{\mathrm{obs}} &\sim \mathcal{N}\left(\mathrm{CH}_{4,i}^{\mathrm{true}},\,\sigma_{\mathrm{C},i}^2\right) \end{align} \tag{ 3 }$$

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

$$\begin{equation} \mathrm{CH}_{4,i}^{\mathrm{true}} \sim \mathcal{N}\left(\alpha_t + \beta_t\,\mathrm{NO}_{2,i}^{\mathrm{true}},\,\gamma_t^2\right) \end{equation} \tag{ 4 }$$

where we have now introduced the daily model parameters α_t , β_t and γ_t . On a given day t, α_t and β_t are respectively the y-intercept and slope of the line of best fit relating methane to nitrogen dioxide, while γ_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)–(4) into the single model equation

$$\begin{align} \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). \end{align} \tag{ 5 }$$

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 α_t , β_t and γ_t .

We add hyperparameters to our model by including a multivariate prior distribution for α_t and β_t , shown in equation (1). We include this multivariate prior to allow for the possibility of learning the extent to which α_t and β_t are correlated. Although we believe it is likely that α_t and β_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 µ_α , µ_β , Σ₁₂ and Σ₂₁, and assume a uniform flat prior on the domain $\left(0,\;\infty\right)$ on both Σ₁₁ and Σ₂₂. This lets us learn the extent to which α_t and β_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 [43] via the interface CmdStanPy [44] 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.

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 \geqslant 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 α_t and β_t is given by equation (1). However, when fitting the hierarchical model to all data-rich days simultaneously, we have learned information about what sort of values the model hyperparameters 'should' take. To account for this information we have learned, we add a prior on the hyperparameters with

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

We have θ as the mean vector of the 4000 posterior draws of (µ_α , µ_β , Σ₁₁, Σ₁₂, Σ₂₂) from fitting the data-rich model to data-rich days. Υ is the $5 \times 5$ covariance matrix of µ_α , µ_β , Σ₁₁, Σ₁₂ and Σ₂₂, 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.

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 [43, 45]. We do so making use of the Cholesky decompositions of our covariance matrices [46].

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

$$\begin{align} \sigma_{\mathrm{C}t} & = \frac{1}{N}\sum_{i = 1}^{N} \sigma_{\mathrm{C}i,t} \end{align} \tag{ 7 }$$

$$\begin{align} \sigma_{\mathrm{N}t} & = \frac{1}{N}\sum_{i = 1}^{N} \sigma_{\mathrm{N}i,t}. \end{align} \tag{ 8 }$$

This yields a new likelihood equation given by

$$\begin{align} \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). \end{align} \tag{ 9 }$$

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 [47, 48]. 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.

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}}$. If predictions CH$_4^{\mathrm{pred}}$ were to be used in an inversion analysis, $\sigma_\mathrm{C}^{\mathrm{pred}}$ would be the error associated to CH$_4^{\mathrm{pred}}$.

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

$$\begin{align} \chi_{\nu,t}^2 & = \frac{\chi^2_t}{\nu_t} \end{align} \tag{ 10 }$$

$$\begin{align} \chi^2_t & = \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}}}} \end{align} \tag{ 11 }$$

$$\begin{align} \nu_t & = N - 8, \end{align} \tag{ 12 }$$

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.

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 [40]. 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.

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 [42]. P_S is given in (Pa) and $VC_{\mathrm{H}_2\mathrm{O}}$ given in (kg m⁻²). 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⁻² 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].

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}\gt2$. The precision of total mol CH₄ is given by $\sigma_{\mathrm{total\;mol\;CH}_4} = \sqrt{\sum_{i = 1}^{N}\left(\sigma_{\mathrm{mol\;CH}_{4},i}\right)^2}\label{eq:total_mol_precision}$.