Avoiding methane emission rate underestimates when using the divergence method

Methane is a powerful greenhouse gas, and a primary target for mitigating climate change in the short-term future due to its relatively short atmospheric lifetime and greater ability to trap heat in Earth's atmosphere compared to carbon dioxide. Top-down observations of atmospheric methane are possible via drone and aircraft surveys as well as satellites such as the TROPOspheric Monitoring Instrument (TROPOMI). Recent work has begun to apply the divergence method to produce regional methane emission rate estimates. Here we show that when the divergence method is applied to spatially incomplete observations of methane, it can result in negatively biased time-averaged regional emission rates. We show that this effect can be counteracted by adopting a procedure in which daily advective fluxes of methane are time-averaged before the divergence method is applied. Using such a procedure with TROPOMI methane observations, we calculate yearly Permian emission rates of 3.1, 2.4 and 2.7 million tonnes per year for the years 2019 through 2021. We also show that highly-resolved plumes of methane can have negatively biased estimated emission rates by the divergence method due to the presence of turbulent diffusion in the plume, but this is unlikely to affect regional methane emission budgets constructed from TROPOMI observations of methane. The results from this work are expected to provide useful guidance for future implementations of the divergence method for emission rate estimation from satellite data -- be it for methane or other gaseous species in the atmosphere.


Introduction
Methane is a powerful greenhouse gas, with a far greater warming potential (84 times greater on a 20-year timescale) and shorter atmospheric lifetime (9 years instead of centuries) than carbon dioxide [1,2].These attributes make methane an attractive target for mitigating the short-term effects of climate change, and have been the focus of recent climate summits and global commitments towards emission reductions [3].Nevertheless, in recent years, the rate of increase of atmospheric methane has itself increased.[4,5].Roughly 30% of anthropogenic methane emissions are attributed to the fossil fuel industry [6,7], making increased monitoring and accounting of emissions from this sector an important factor in meeting national commitments towards methane emission reductions.
Satellite observations are a powerful tool for monitoring atmospheric methane abundances [8], with remote sensing of methane from space providing opportunities for repeated and unscheduled monitoring of emissions.The era of greenhouse gas observing satellites began with the SCanning Imaging Absorption SpectroMeter for Atmospheric CHartographY (SCIA-MACHY) [9] in 2003, and subsequent generations of satellites have given rise to instruments with ever increasing capabilities.The TROPOSpheric Monitoring Instrument (TROPOMI) provides daily global coverage of methane observations with an updated 5.5x7 km 2 pixel resolution [10], whilst other instruments such as GHGSat provide intermittent, targeted methane observations down to 50x50 m 2 resolution [11].Many greenhouse gas-observing satellites lack the spatial resolution to calculate asset-level emissions, and so aircraft and drone surveys are used to bridge this gap.These instruments can image and estimate facility-level methane emission rates [12,13], and such facility-level measurements can be used for reporting under the Oil & Gas Methane Partnership 2.0 (OGMP 2.0) framework.This is a multi-stakeholder initiative launched by the United Nations Environment Programme (UNEP) and the Climate and Clean Air Coalition, aimed at improving the accuracy and transparency of methane emissions reporting in the oil and gas sector [14].Recently, "top-down" methane emission estimates calculated from aircraft observations have been found to be in disagreement to "bottom-up" emission estimates reported from industrial activity [15,16], and more work is required to reconcile these standards of reporting methane emissions.
There are a variety of methods for constructing top-down emission estimates from satellite observations of methane.Some analyses use Bayesian methods to estimate regional methane emission rates [17]; Bayesian methods require priors to be specified on the spatial distributions of emission rates, which are sometimes constructed from bottom-up emission estimates.Other methods such as the inverse Gaussian plume method [18,19,20], the source pixel method [8,21], the cross-sectional flux method [22,23], and the integrated mass enhancement (IME) method [24,25] rely directly on column measurements of methane (either remotely sensed or in-situ) and wind field data.These methods are particularly attractive for construction top-down emission budgets as they are entirely divorced from prior estimates that are informed by self-reported production statistics.
Another top-down methodology for emission estimation is the divergence method [26].Like the inverse Guassian plume, source pixel, cross-sectional flux and IME methodologies, the divergence method relies solely on wind field and column observations of the quantity of interest.In contrast to these other top-down methodologies, the divergence method can also be used for regional-scale emission rate estimation [26,27,28].In the divergence method, emission is calculated via the spatial gradients of zonal (i.e., advective) fluxes.Mathematically, the total sources and sinks of emission E are calculated via the emission equation where ⃗ F adv is the advective flux of a quantity of interest (e.g., methane or other pollutants).The divergence method for estimating the spatial distribution of methane emissions is attractive because it is entirely data-driven and does not rely on prior estimates of spatial emission distributions as extensively as Bayesian methods do.Unlike the source pixel method, the divergence utilises wind field and concentration information away from source locations [29].Originally presented as a method for estimating the location and emission rates of sources of nitrogen dioxide, this methodology is now being used to estimate regional-level methane emission rates [27,28].
It is important to make explicit some important assumptions that are currently intrinsic to this methodology.First, plumes of any gas (including methane) propagate through the atmosphere not only by advection, but also molecular and turbulent diffusion [30,31,32]; in atmospheric transport, turbulent diffusion is usually the dominant effect over molecular diffusion in practice.The emission equation E = ∇ • ⃗ F adv that is central to many regional-level methane budget estimates does not take the effect of turbulent diffusion into account.Correcting for turbulent diffusion requires the usage of a modified emission equation where ⃗ F dif is the turbulent diffusive flux of a quantity of interest.Second, regional methane emission rate estimates are often time-averaged.The linear property of the divergence operator means that time-averaged estimated emission rates could be calculated in two different ways: first, by time-averaging daily estimated emission rates, or second, by taking the divergence of time-averaged daily advective fluxes.With perfect observations, these procedures should yield identical results, but satellite observations of methane are often spatially disrupted with missing coverage.To the best of our knowledge, no work has yet been done to examine the consequence of this choice of order of operations when the divergence method is applied to methane observations.
In this work, we derive analytical expressions for ⃗ F dif , and generate synthetic simulated satellite observations of Gaussian plumes to examine under what physical scenarios it becomes important to include ⃗ F dif in the emission equation.We find that when a plume of methane is relatively diffusive, methane emission estimates via the divergence method can become inaccurate if ⃗ F dif is excluded from the divergence calculation (under certain conditions).In this case, the estimated emission rate of the source is underestimated, and the spatial distribution of emissions is incorrectly distributed.We also demonstrate that time-averaged emission estimates calculated from spatially incomplete observations may be inaccurate if care is not taken to use a time-averaged advective flux in the emission equation (as opposed to taking the time-average of daily emission estimates via the divergence method).We compare the results of our synthetic study to a case study of the Permian basin, using TROPOMI observations of methane [33].We find that it is unlikely that a regional methane emission budget calculation of the Permian via the divergence method would be negatively biased due to the effects of any turbulent diffusion.We do find though that the sparse nature of the TROPOMI methane data product can result in negatively biased time-averaged methane emission rate estimates when the divergence method is used, specifically affecting cases where daily emission rate estimates which are time-averaged to produce the time-averaged emission rate estimate.When using the divergence method in conjunction with spatially sparse observations, it is important to take the divergence of time-averaged daily fluxes to obtain a time-averaged emission estimate.

Synthetic case study
We generate simulated satellite observations of ideal steady-state Gaussian plumes [30] resulting from isolated point sources with known emission rates, and use them in an investigative synthetic case study.These plumes are characterised by a known emission rate Q true [kg s −1 ], wind speed w [m s −1 ], wind angle θ relative to the x-axis of the observation grid, and constant of turbulent diffusion K [m 2 s −1 ].Simulated plume observations are generated via Eq. 9, which is derived in S1.Note that when we calculate grid cell values, we numerically spatially integrate Eq. 9 over the area of the grid cell, and so the results of our synthetic study are independent of grid cell resolution.This is important to bear in mind as real observational instruments will have pixel resolutions spanning from meter to kilometer scales, which changes the amount of fine-scale structure that will be visible in observed plumes.Fig. 1 shows some examples of how these parameters affect plume morphology.For our simulated observations, we use the divergence method to estimate the spatial emission field both with and without including the diffusive flux term in the emission equation.We find in our synthetic study that there are a variety of scenarios where the use of the divergence method results in negatively biased estimated emission rates, i.e., when turbulent diffusion is neglected in some cases, or when time-averaged estimated emission rates are calculated without time-averaging daily flux fields when daily observations are spatially incomplete.

Underestimating emission rates due to turbulent diffusion in the plume
Fig. 2 demonstrates the application of the divergence method to a simulated satellite observation of a plume.We begin by using the standardly formulated emission equation where the turbulent diffusive flux term term is not included, i.e., where C is the simulated column density of methane in a grid cell and ⃗ w is the simulated wind vector in a grid cell.We find that when estimating the emission field without including the turbulent diffusion term in the emission equation, emissions are incorrectly spatially distributed in the estimated emission field.Rather than estimating a single point of positive emission, we instead find that "arrowhead" shapes of positive emission are estimated, in conjunction with negative emissions (or sinks) downwind within the plume shadow (see Fig. 2b).To obtain a total estimated emission rate Q est , we spatially integrate the estimated emission field within a circle of fixed radius centered on the source location.Using Eq. 3, we find that the total estimated emission rate of the emission field is underestimated, the extent of which depends on the conditions discussed in the following two paragraphs.
When increasing the opening angle of a plume (which scales as a function of K/w), the total estimated emission rate decreases for a fixed radius of spatial integration over the estimated emission field.This demonstrates that the emission rates of "diffuse" plumes are poorly estimated when using the divergence method, relying solely on advective fluxes.To first order, the underestimation is proportional to K/w.It is important to note in our synthetic study that we measure estimated emission rates as percentages of the true emission rate, and so results are independent of the mass emission rate of the point source.For a real instrument, higher emission rates mean a higher measured signal-to-noise ratio.
The total area for spatial integration of the estimated emission field also influences Q est .In Fig. 3 we alter the radius r of the cirular area of integration for a single simulated observation, and find that increasing the radius of integration increases Q est , i.e., the total estimated emission rate is less negatively biased.Although emissions are incorrectly distributed in the estimated emission field, they are distributed such that integrating the estimated emission field over a larger area improves the total estimated emission rate.Thus, we determine that in our procedure for estimating the emission rate of a plume, r and K/w are two independent parameters which determine the extent to which Q est is underestimated.In Fig. 4 we vary r and K/w over physically realistic values (10 to 400 m 2 s −1 , see Sec. 2.2 and 4.4 for a description of methods) and examine how greatly Q est is underestimated for an ideal Gaussian plume resulting from a point source of emission.We find that in the most extreme cases, Q est can be underestimated by more than 40%.
In practice, different plume measurement methodologies will correspond to different parameter locations on Fig. 4, and have characteristic biases associated with them.We indicate three examples on the right hand edge of Fig. 4 of regions where certain instruments tend to lie.Global coverage satellites such as TROPOMI tend to have very large fields of view as well as lower pixel resolutions [9,10,34].Regional emission budgets performed on the scale of tens of kilometers are unlikely to experience a high level of negatively biased emission estimates due to the lack of a diffusive term in the divergence method (e.g., see Sec. 2.2).More targeted satellites have higher pixel resolutions and are capable of imaging plumes on the scale of tens of meters [11].Although such satellites have fields of view that can still exceed ten kilometers or more, it would still be possible (given the high pixel resolution) to spatially integrate estimated emission fields over small enough areas to experience the outlined negative bias of the uncorrected divergence method, an issue that could arise if attempting to spatially isolate one plume from another adjacent source.Lastly, plume imaging surveys based from aircraft or drones [12,13] would likely have the smallest field of view and would consequently be most likely to experience negative biases in emission estimates if they were to use the divergence method for emission rate estimation.
The significance of diffusion on the accuracy of the divergence method can be seen most directly when a diffusive term is included in the emission equation, i.e., Eq. 12. Making this addition to the method, we show in Fig. 5 that when the expression for diffusive flux (as calculated by Eq. 13) is included in the emission equation, then the estimated emission field is correctly constrained to a point source (allowing some slight deviation due to numerical derivative effects).Q est is found precisely to equal Q true .
A caveat to the efficacy of including diffusive flux in the divergence method to restore accurate emission estimates is that it is in practice difficult to estimate the constant of turbulent diffusion K.In our synthetic study, it is possible to know and choose the precisely correct value of K to calculate ⃗ F dif , but for real data this is much harder to calculate.In Sec.2.2, we estimate values of K from high-resolution observations of real plumes.

Miscalculating emission rates due to missing data
Beirle et.al. (2019) [26] was the first study to showcase the divergence method for estimating emission rates, and used TROPOMI observations of nitrogen dioxide to estimate emission rates from cities and power plants.They pointed out that, due to the linear properties of the divergence operator, it was sensible to time-average the daily fluxes of nitrogen dioxide first and then take the divergence of the time-averaged flux to obtain a time-averaged estimate of nitrogen dioxide emissions.Beirle et. al. (2019) outlined that this was a sensible procedure as the time-averaged nitrogen dioxide advective flux would have a smooth spatial distribution and thus allow for a more accurate calculation of spatial derivatives.TROPOMI observations of methane differ significantly from that of nitrogen dioxide, in that the spatial coverage of the methane data product is much less complete than the nitrogen dioxide data product [35,36].Fig. 6 demonstrates the problem this poses as we look to apply the divergence method to TROPOMI observations of methane; when we randomly mask 10% of the observational data of a single simulated plume, the corresponding emission field estimated via the divergence method is missing a significantly larger amount of spatial coverage, and dramatically underestimates the total emission rate.This is because the numerical methods for calculating spatial derivatives [37,38] (see Eqs. 15 and 16) require eight valid neighboring data points, and so any missing observational data is magnified into even more missing spatial coverage of the estimated emission field.Time-averaging observations prior to numerically calculating spatial derivatives, and choosing to describe average emission rates over time domains rather than at specific time points, hopefully allows a chance at circumnavigating the problem of truncated spatial observations.Thus, in this section we investigate two different methodologies for calculating time-averaged emission estimates using the divergence method.In the first methodology (which we denote by E 1 ), we calculate daily estimates of emission fluxes and then temporally average them [27,28].Mathematically, this is represented via where the subscript t denotes advective fluxes on day t in the time period of interest.In the second methodology (which we denote by E 2 ), we temporally average the daily advective fluxes to create a time-averaged advective flux, and then take the  [28], i.e., calculated using the methods of Schneising et.al. 2020 [22], but using the same data as Veefkind 2023 [28].
divergence of the time-averaged advective flux to obtain a time-averaged estimated emission field [26].Mathematically, this is represented via An explanation of these methodologies is given again briefly in Sec.4.3 and in greater detail in supplementary section S4.In this section we do not correct the estimated emission fields of our simulated plumes for the effects of turbulent diffusion as we did in Sec.2.1.1,and focus only on the differences between the E 1 and E 2 methodologies.
In Fig. 7 we simulate a time-averaged study of 30 steady-state Gaussian plumes.Plume parameters are left unvaried, and each of the 30 repeated simulated observations has a random 30% of its pixels masked.We then display the resulting estimated emission fields obtained via E 1 and E 2 , and find that the emission field of E 2 is spatially complete, whereas E 1 is still missing some spatial coverage.The total integrated emission rate of E 1 is also severely underestimated, but the total integrated emission rate of E 2 retrieves the correct time-averaged emission rate (apart from the slight negative bias due to the presence of diffusion in the plume, which was discussed in the previous section).Figures showing the difference in resulting estimated emission fields under non-static conditions are shown in the supplement in Figs.S2, S3, and S4.
We investigate how the difference in performance of the E 1 and E 2 methodologies varies as a function of amount of daily missing data, and when plume parameters are allowed to vary in time.These results are shown in Fig. 8.We find that as the amount of missing observational data increases, both methodologies underestimate the true average emission rate, but that method E 2 (i.e., averaging daily fluxes of C and taking the divergence once) allows for estimates that are far more robust against missing data.This holds true for both static and time-varying simulated plume observations, although the time-averaged estimated emission rates for the time-varying plumes have more variance than the static plumes.This simulation differs from realistic physical scenarios in that data is randomly masked in a physically uncorrelated manner (which would not be the case with cloud cover), but it nonetheless demonstrates that the way in which time-averaged emission estimates are calculated using the divergence method is not trivial.Additionally, Fig. 8 demonstrates that it is also possible to over-estimate the source emission rate, though this typically only occurs at a critical "turn over" point where the fraction of daily missing data begins to dominate over the number of repeated observations.Past this point, we find that estimated emission rates will only be underestimated as a consequence of spatially incomplete data, but the exact location of this critical value is highly dependant on the number of repeated observations and the spatial distribution of missing data.

Permian basin case study
The Permian basin is the largest oil and gas producing region in the United States, producing nearly 6,000 barrels of oil a day as of January 2023 [40].Due to its prominence and size the Permian is frequently a target of ground-based, airborne, and space-based campaigns monitoring methane emissions [41,42].We grid three years of daily TROPOMI methane observations of the Permian basin (2019-2021) onto a 0.2 x 0.2 latitude-longitude grid using an area-weighted oversampling [43] and calculate yearly emission estimates.
We use the Copernicus Sentinel-5P TROPOMI Level 2 methane data product [33], which provides daily global methane observations as total column-averaged methane mixing ratios X(CH 4 ).We convert X(CH 4 ) [ppbv] to total column densities n Table 2. Comparison of yearly methane emission rate estimates [Tg/year] for the Permian when using the two different time-averaging methodologies.In the first column we show yearly methane emission rate estimates for the Permian via E 2 , when we average over advective fluxes for the year and take the divergence to yield a time-averaged emission estimate.In the next column we show the same yearly emission rate estimate calculated via E 1 , where we calculate daily methane emission estimates and time-average them.In the penultimate column we show the difference of the total estimated emission rate between the two methodologies, when the estimated emission fields are only spatially integrated over the intersection of the two estimated emission fields.This is to examine whether the difference in the estimated emission rate is driven by differences in spatial coverage or not.Supplementary figures 5, 6, and 7 plot these estimated emission fields.Also shown in the last column is the average daily spatial coverage of the Permian basin by our regridded TROPOMI methane observations.[mol m −2 where p is the dry surface pressure in Pa, g is the constant of gravity (9.81 m s −2 ), and m air is the molar mass of dry air (28.96 × 10 −3 kg mol −1 ) [28].We determine dry surface pressure from the ERA5 reanalysis dataset from the European Centre for Medium-Range Weather Forecasts (ECMWF) [44].The total column density n must further be background-reduced before applying the divergence method.Following the methodology of Veefkend et.al. (2023), methane column enhancements ∆n above background can be determined on a daily basis via where c 0 and c 1 are constants determined from a linear fit of n vs p.This methodology assumes a limited-size area of interest, such as the Permian, and further assumes that the tropopause pressure is constant over the area of interest [28].The 2019-2021 average methane enhancement over the Permian basin is shown in Fig. 9.
Using ERA5 wind data from the ECMWF [44], we calculate daily advective fluxes of methane, and then calculate yearly time-averaged methane emission maps of the Permian basin using both the E 1 and E 2 methodology.As a reminder: in the E 1 methodology, we use the divergence method to estimate daily emission fields, which are then time-averaged, and in the E 2 methodology, we first time-average daily advective fluxes methane, and use the divergence method to estimate the emission field from the time-averaged fluxes.As in Sec.2.1.2,for simplicity E 1 and E 2 both exclude the turbulent diffusive flux term, in order to purely examine the differences arising from the choice of how to order the time-averaging operations (though later in this section, we examine the extent to which turbulent diffusion matters for this observing system).The yearly estimated emission fields produced via the two methodologies are shown in the supplement, and the estimated emission fields produced via the two methodologies for the entire time period 2019-2021 is shown in Fig. 10.Our yearly total estimated methane emission rates for the Permian are shown in Tables 1 and 2. We find good agreement between our time averaging methodology E 2 and other Permian emission estimates from previous work, but find that the E 1 methodology (in which daily emission estimates for the Permian are time-averaged) significantly underestimates the time-averaged emission rates when compared to previous estimates in the literature.It should be noted that some of the results of previous work quoted in Table 1 use TROPOMI methane columns calculated via the Weighting Function Modified Differential Optical Absorption Spectroscopy (WFM-DOAS) algorithm version 1.5 [39], an alternative data product which has a greater spatial coverage than the traditional Copernicus Sentinel-5P TROPOMI Level 2 methane data product.Consequently, these results are less likely to be susceptible to the weaknesses we note here in using the divergence method with spatially corrupted data.
The difference in results between our E 1 and E 2 methodologies is likely due to the sparse daily spatial coverage of the TROPOMI methane data product over the Permian basin.For the years 2019-2021, the average daily coverage of our regridded Table 3. Yearly estimates of methane emission from the Permian basin [Tg/year].We present yearly estimates (calculated using both the E 1 and E 2 time-averaging methodologies), and compare results between when using the fourth-order central finite difference and the second-order central finite difference to calculate numerical derivatives.The second-order central finite-difference requires fewer valid neighbors to calculate derivatives, and so could potentially lessen the discrepancy between the results of the E 1 and E 2 methodologies.We find that for the year 2019 (which had the poorest average spatial coverage over the Permian by the TROPOMI methane data product), the difference between the yearly methane emission budgets estimated via E 1 and E 2 is slightly decreased, but the gap between the two is not bridged within error.Also shown in this table is the percentage area coverage of the Permian basin by the estimated emission field produced by the E 1 methodology.As expected, the percentage area coverage is improved by the usage of the second order central finite difference in calculating derivatives, and the greatest improvement is seen in the year 2019.However, this increase in spatial coverage of the estimated emission field is not sufficient alone in bridging the discrepancy of results produced by the E 1 and E 2 methodologies.

Year
4th order CFD 2nd order CFD TROPOMI observations of the Permian basin never exceeds 50% (Table 2).In Table 3, we examine whether the choice of order of central finite difference influences the results obtained when calculating time-averaged emission rate estimates for the Permian.Although the spatial coverage of the estimated emission field produced via E 1 improved by using the second order central finite difference to calculate derivatives instead of the fourth order central finite difference, we did not find any significant change in the total estimated emission rates.
As the total estimated emission rate for the Permian is obtained via spatially integrating the estimated emission field over the extent of the Permian, it follows that a methodology which produces greater spatial coverage (E 2 ) will produce a larger total estimated emission rate than a method with less spatial coverage (E 1 ), assuming that grid cells are generally estimated to have positive emissions.However, in Table 2 we show the difference in total estimated emission rate between E 1 and E 2 for when we only spatially integrate over the intersection of their spatial coverages.In this case, any difference in total estimated emission rate for the overlapping coverage is due to the differing methodologies of E 1 and E 2 , rather than the difference in spatial coverage of the Permian they produce.We find that for the year 2019, the difference in average total estimated emission rate between E 1 and E 2 can potentially be explained entirely by the fact that the estimated emission field produced via E 2 covers a significantly larger geographic area than that produced via E 1 .This is shown in Fig. S4.However, for the years 2020 and 2021, the spatial coverage of the estimated emission fields obtained via E 1 and E 2 are complete, and thus the difference in total estimated emission rates cannot be explained by a difference in spatial coverage.Increasing the spatial coverage of the region of interest by the methane data [39,45] may close the gap in the results obtained between the two methods.This indicates that an analytical determination of the constant of turbulent diffusion K from observations may only be suitable for large, plume-like emitters and would not be applicable to large-scale diffuse emission environments.
We additionally calculate the percent change of our yearly estimated methane rates for the Permian when including the diffusive flux calculation of Eq. 14.Based on a number of real methane plumes measured by GHGSat [46], a typical value of K was determined to be highly variable and in the the range between 10 and 400 m 2 s −1 (see Sec. 4.4 for a description of methods).We find in all three year-averaged emission estimates for the Permian that the total estimated emission rate is increased by less than a millionth of a percent when K = 400 m 2 s −1 , the maximum value for the constant of turbulent diffusion that we consider in this work.Furthermore, the inclusion of the turbulent diffusion flux term (even at this strength) does not change the spatial distribution of estimated emissions, and still shows emission enhancements over the Midland and Delaware basins.This is not unexpected given the results shown in Fig. 4, which suggests that integrating over large areas sufficiently corrects for any negative bias introduced by neglecting turbulent diffusion in the divergence method.However, the plume observations [46] used to calculate the potential range of values of K may not be indicative of typical plume morphology for other geographic locations or meteorological conditions.If satellite observations are not screened for low wind fields, greater emission underestimates may be obtained via the divergence method.Smaller plumes than those shown in this example study may exhibit different behaviour, though they may not be resolvable on TROPOMI-pixel scales.

Discussion
In this work, we examine the conditions under which the divergence method for estimating emission rates may prove to produce negatively biased results.Using a simulation study with synthetic satellite observations of ideal Gaussian plumes, we showed that highly-resolved, diffuse plumes may have negatively biased emission rates when their estimated emission fields are spatially integrated over narrow fields of view.Our simulation study suggests that this affect would only be of concern for observations obtained by high-resolution, narrow field of view methodologies, e.g., drones, or where high-resolution satellite data has clipped the area into which a plume extends.In contrast, our case study of the Permian basin with TROPOMI methane observations does not find that yearly estimated methane emission budgets are impacted by including even a high estimate of turbulent diffusion.In the future, as satellites become more capable of resolving individual plumes, it will become important to correct for turbulent diffusion when estimating the spatial distribution of emissions via the divergence method.In these scenarios where fine-scale plume structure can be spatially resolved, emission rate estimation is likely better performed via the IME or plume rotation methods, as the divergence method (like the inverse Gaussian plume method) is unable to adequately take turbulence into account [29].Conditions where diffusion may be dominating over advection may also be identified by screening for very low wind speeds [28].
In our synthetic study, we do not model the effects of instrumentation error or averaging kernel sensitivities [47]; this is in order to isolate and examine the effects of missing data and the turbulent diffusion flux term.We do give an example supplementary figure (Fig. S1) showing a simulated plume observation and resulting estimated emission field where some added noise is included in the observation.Previous work [29] also shows that the Gaussian plume model does not accurately model the fine-scale turbulent eddy structure for highly resolved plumes; although this means our synthetic plume observations are more idealised than real plumes, it allows us to test the assumptions of the divergence method systematically against a simple test dataset in a robust manner.We do not mean to suggest that the inclusion of the turbulent diffusion flux term in the divergence equation would entirely correct the calculation for real plumes, but that it is likely a necessary first-order correction.
We also examine two possible methodologies for calculating time-averaged emission estimates using the divergence method.Using simulated spatially incomplete plume observations, we find that time-averaging daily emission rate estimates produced via the divergence method will consistently underestimate the true average emission rate.We compare these results to an alternative methodology previously described for nitrogen dioxide observations [2].By this method, daily advective fluxes are time-averaged, and the divergence is taken thereof in order to obtain a time-averaged emission estimate.We find in our synthetic study that this latter methodology yields robust emission estimates even in the face of spatially incomplete observations.We compare these two methodologies by constructing yearly methane emission budget estimates for the Permian basin for the years 2019-2021, using the Copernicus Sentinel-5P TROPOMI Level 2 methane data product [33].We find that these two methodologies do not produce congruent emission rate estimates, and that the latter methodology produces estimates in agreement with previous top-down estimates for methane emission in the Permian basin [17,28].Methods and datasets exist that can augment the spatial coverage of the TROPOMI methane data product [39,45], which in turn would augment the spatial coverage of the advective flux field of methane prior to the application of the divergence method.Spatial smoothing and interpolation could also be used to try and make the spatial coverage of the advective methane flux field more complete.In this work, we do not explore these avenues further, preferring to examine the differences between the E 1 and E 2 methodologies.When using the divergence method for methane emission rate estimation in the Permian, we find areas bordering the Delaware and Midland basins that are estimated to have negative emission rates.We do not expect this to truly be the case.Even in the "perfect" case in our synthetic study when the point source of emission is correctly estimated (see Fig. 5, panel c), we still estimate some grid cells to have negative emission rates.In this case, this is an effect of the discretisation of the emission equations and the manner in which numerical derivatives are calculated over our grids.Therefore, in some areas, especially in the vicinity of large methane sources, the divergence method will return negative emissions.These are local artefacts; the area-integrated emissions remain positive both in our synthetic study and our case study of the Permian basin (and in this latter case, our results are in good agreement with previously estimated methane emission rates).Whilst the TROPOMI methane observations over the Permian can not be considered to be ideal plumes, it may be the case that the regions of negative estimated emission are analogous to those in the synthetic study, as we know that that they are bordering known regions of strong positive emission.Other work also demonstrates that some regions of negative emission estimated via the divergence method in the Permian can be related to changes in orogoraphy, surface albedo, or convergent wind fields [28].Other works demonstrate that time-averaged emission calculations from the divergence method are unlikely to be dominated by convergent wind fields [27], and other formulations of the divergence method calculate advective flux without any contribution from wind field divergence [48,49].One could develop a model that prohibits the estimation of negative methane emissions (other works do so in a Bayesian framework [17,50]), though at this stage this would no longer purely be the "divergence method", which is driven entirely by the data and the principle of the conservation of mass.
We conclude that the divergence method for estimating methane emissions (as described in this work) would be best applied to regional analyses where the affects of turbulent diffusion are unlikely to dominate over the scale of advective methane fluxes.Under these scenarios, we have demonstrated that it is possible to obtain accurate total estimated emission rates due to the large area of spatial integration.For such larger-scale analyses, it is important to be able to accurately determine the spatial distribution of emissions in order to identifying regions of high emissions, which may not always align with reported bottom-up inventories.In contrast, emission quantification for point sources with known locations (e.g., facility-scale emission audits) may be better achieved using alternative methodologies such as IME or the plume rotation method.When planning emission measurement campaigns, it is important to identify whether the primary goal is emission location or quantification, as different methodologies may prove better suited to one goal or the other.Whenever possible, spatially completely methane data products should be used.When using spatially incomplete datasets, it may be the case that taking the divergence of time-averaged advective fluxes of methane will produce more accurate methane emission rate estimates.

Generating simulated satellite observations of plumes
For our synthetic studies we generate simulated top-down observations of ideal Gaussian plumes [30] via the equation where and θ is the wind angle relative to the x-axis (in the anti-clockwise direction).Eq. 9 is derived in S1 .There are a variety of assumptions that are fundamental to the ideal Guassian plume equation [30], but most important of note here is that diffusion is assumed to be dominated by advection, and thus diffusion only takes place perpendicular to the wind vector characterised by w and θ .

Calculating emissions and flux terms
We calculate spatially varying estimated emission fields E using both simulated plume observations and the TROPOMI L2 methane data product.We calculate E via two different emission equations.
The first emission equation (commonly found in literature [26,28]) is where is the advective flux of some column density C. In the case of our synthetic plume observations, C is generated via Eq. 9.For TROPOMI observations of methane, we convert column-averaged mixing ratios to above-background column density enhancements [28].⃗ F adv is then given by where ⃗ w is a spatially varying wind vector with magnitude w and angle θ relative to the x-axis of our grid.For our synthetic studies we specify w and θ ourselves.For our work with TROPOMI observations of the Permian basin, we take ⃗ w to be the ERA5 wind data on multiple pressure levels, temporally averaged daily over a wind history at 1700, 1800 and 1900 hours, and then averaged vertically to an altitude of 500m to account for changes in wind vector through the boundary layer [28].
To examine the extent to which turbulent diffusion influences estimated emission rates via the divergence method, we also calculate E via a second emission equation 9/19 ⃗ F dif is the turbulent diffusive flux of some column density C, and for an ideal Guassian plume is given by where θ is the wind angle relative to the x-axis of our grid.C is again either generated via Eq. 9 or calculated from TROPOMI satellite observations of methane.Eq. 13 is derived under the assumption that diffusion only takes perpendicular to the wind vector ⃗ w.If, however, we choose to ignore this assumption (but still assume that K is constant in space), that we can work directly in the (x, y) grid and state that Eqs. 10, 11, 12, 13, and 14 are derived in S2.
We need to calculate spatial derivatives over a cartesian grid to fully obtain E in Eqs. 10 and 12.For first derivatives, we use the fourth-order central finite difference [37] where V is a spatially varying quantity and d is the grid spacing in coordinate p.This numerical recipe is commonly used for emission estimates via the divergence method [26,28].For second derivatives, we use the fourth order discretization [38]

Calculating time-averaged emission rates
We calculate time-averaged estimated emission fields using two methodologies.In the first (denoted by E 1 ), we calculate daily estimated emission fields and time average them to obtain E 1 .In the second methodology (denoted by E 2 ), we time-average daily fluxes of C, and take the divergence of the time-averaged flux to obtain E 2 .With spatially complete observations of C over an entire time period, the two methods yield identical results, but for TROPOMI observations of methane, data is often spatially masked due to cloud cover and albedo effects.Detailed equations describing these two methodologies is given is S4.

Estimating the constant of turbulent diffusion K
It is in practice difficult to estimate constants of turbulent diffusion.If K is assumed to be constant in space and time, then the standard deviation "width" of a Gaussian plume can be described via where σ is the width of the plume [m], x is the downwind distance in the plume [m], u is the wind speed [m s −1 ] and K is the constant of turbulent diffusion [m 2 s −1 ] [30,51].We take multiple GHGSat scenes of isolated methane plumes [46] and measure values of σ at multiple downwind locations within each plume.We then fit a linear function to σ 2 against x for each plume using the method of least squares.The slope of the fitted function yields 2 K/u, and thus K can be determined as u is known for each scene.We determine using these plumes that K can vary between 10 and 400 m 2 s −1 .As the plume in a was constructed using a point source of positive emission located at the origin, we can see in b that the estimated emission field has an incorrect spatial distribution."Positive" emissions are incorrectly distributed in a horseshoe-shaped pattern around the origin, and "negative" emissions (or sinks) downwind in the shadow of the plume.When the estimated emission field is integrated over a circular region of radius r = 500 m centered on the origin, the total emission rate is slightly underestimated.For a fixed r and true point source emission rate Q true , the total estimated emission rate Q est will be underestimated to a greater extent if K/w is increased.When K/w increases for the plume, the plume becomes more "diffuse", and the emission equation E = ∇ • ⃗ F adv becomes less able to explain the resulting column density as arising from a point source of emission.A contour plot demonstrating the extent to which a plume's emission rate will be underestimated as a function of radius of integration r and K/w, the ratio of the constant of turbulent diffusion to the wind speed.Shown on the right had side of the plot are general regions where certain observing methodologies would lie on this spatial scale.The range of values of K/w (i.e., the x-axis) are chosen to represent a range of physically realistic values, estimated from high-resolution GHGSat scenes containing methane plumes.Global coverage satellites like TROPOMI tend to have larger pixel sizes and large fields of view, and hence are less likely to be affected by the negative bias.Targeted satellites can have narrower fields of view on the order of km scales, but due to their higher pixel resolutions can still image plumes on the scales of hundreds of meters.Consequently, these instruments may be affected by the negative bias of estimated emission rates, depending on how they are used.Surveys conducted via drone or plane may be especially susceptible to this bias if the divergence method is used to estimate emissions.The point source of emission is located at the origin with actual emission rate Q true = 10 kg s −1 .10% of the pixels in a are randomly masked to simulate the effect of poor data or otherwise unavailable observations, as is often the case for satellite observations of methane column density.b The resulting estimated emission field for the plume observation in a, obtained using the divergence method and emission equation E = ∇ • ⃗ F adv .Our algorithm uses the fourth-order central finite difference to calculate numerical gradients [37], which means that every unavailable pixel in a results in a 4x4 cross of pixels in b where we are unable to calculate E. Consequently, any amount of missing data in a quickly becomes a large amount of missing coverage in b.Integrating over the available pixels in b yields a severely underestimated estimated emission rate Q est .Scaling up Q est by the fraction of missing coverage yields a scaled total emission estimate of 16.2 kg s −1 , and does not correct the underestimation.In Fig. S7 in the supplement we demonstrate that this scaling procedure is not sufficient to correct for the problem of missing data, and can actually significantly overestimate the total emission rate (as is the case in this particular example), depending on the fraction of daily missing data and number of repeated observations that go into the time-averaged calculation.In all other estimates of total emission rate in this study, we quote "unscaled" results, i.e., the answer obtained by spatially integrating the estimated emission field over available grid cell values.It is also important to note that individual or time-averaged emission fields estimated via the divergence method in our synthetic study will vary depending on the spatial distribution of missing data, which is a random process.We do not cherry pick particular realisations of spatial distribution for missing data, and examine the variance that is inherent to this random process in Fig. 8.Although the extended three-year time period means that complete spatial coverage is achieved for all the basins, the intermittent missing data on a daily basis means that the methodology of E 1 estimates a significantly lower emission rate than the methodology of E 2 .We find the three-year estimated emission fields E 1 and E 2 to be correlated with r = 0.87 .

Figure 1 .
Figure 1.Simulated column-integrated ideal Gaussian plumes.All the plumes shown have the same emission rate of Q true = 10 kg s −1 and wind angle θ = 45 degrees.a Simulated plume observation with w = 8 m s −1 and K = 40 m 2 s −1 .bSame simulated observation as in a, this time with doubled K. c As in a, this time with doubled K and doubled w.Note that panels a and b are plotted with the same colorscale with the colorbar shown on the left hand side, and that panel c is plotted with the colorbar on the right hand side.The simulated plume in b has a wider opening angle than in a due to the increased ratio of K to w.The plume in c has the same opening angle as the plume in a because K/w = 5 m for both plumes.However, the column density of the plume in c is half that of the column density in a, because according to Eq. 9, column density scales as ∼ (K w) −1/2 .

Figure 2 .
Figure 2. a A simulated column-integrated Gaussian plume, with wind angle θ = 45 degrees and other parameters as indicated on the plot.bThe estimated emission field of the plume using the divergence method and emission equation E = ∇ • ⃗ F adv .As the plume in a was constructed using a point source of positive emission located at the origin, we can see in b that the estimated emission field has an incorrect spatial distribution."Positive" emissions are incorrectly distributed in a horseshoe-shaped pattern around the origin, and "negative" emissions (or sinks) downwind in the shadow of the plume.When the estimated emission field is integrated over a circular region of radius r = 500 m centered on the origin, the total emission rate is slightly underestimated.For a fixed r and true point source emission rate Q true , the total estimated emission rate Q est will be underestimated to a greater extent if K/w is increased.When K/w increases for the plume, the plume becomes more "diffuse", and the emission equation E = ∇ • ⃗ F adv becomes less able to explain the resulting column density as arising from a point source of emission.

Figure 3 .
Figure 3. a A simulated column-integrated Gaussian plume.b The estimated emission field using the divergence method, integrating over a circular area of radius r = 250 m to obtain the total estimated emission rate Q est .c The same estimated emission field as in b, this time integrated over a circular region with r = 500 m to obtain Q est .As r increases, Q est is underestimated to a lesser extent.

Figure 4 .
Figure 4.A contour plot demonstrating the extent to which a plume's emission rate will be underestimated as a function of radius of integration r and K/w, the ratio of the constant of turbulent diffusion to the wind speed.Shown on the right had side of the plot are general regions where certain observing methodologies would lie on this spatial scale.The range of values of K/w (i.e., the x-axis) are chosen to represent a range of physically realistic values, estimated from high-resolution GHGSat scenes containing methane plumes.Global coverage satellites like TROPOMI tend to have larger pixel sizes and large fields of view, and hence are less likely to be affected by the negative bias.Targeted satellites can have narrower fields of view on the order of km scales, but due to their higher pixel resolutions can still image plumes on the scales of hundreds of meters.Consequently, these instruments may be affected by the negative bias of estimated emission rates, depending on how they are used.Surveys conducted via drone or plane may be especially susceptible to this bias if the divergence method is used to estimate emissions.

Figure 5 .
Figure 5.A simulated column-integrated Gaussian plume, and resulting estimated emission fields via the divergence method.a A simulated observation of a plume.b The emission field of the plume estimated via the divergence method, using the emission equation E = ∇ • ⃗ F adv .c The emission field of the plume estimated via the divergence method, using the emission equationE = ∇ • ⃗ F adv − ⃗ F dif .The term ⃗ F dif is calculated according to Eq. 13.The emission field in c is now correctly estimated as a point source, with some small perturbations due to the numerical derivative.The correct total estimated emission rate Q est is now obtained.

Figure 6 .
Figure 6. a A simulated observation of a plume with wind angle θ = 45 degrees, wind speed w = 5 m s −1 , and K = 60 m 2 s −1 .The point source of emission is located at the origin with actual emission rate Q true = 10 kg s −1 .10% of the pixels in a are randomly masked to simulate the effect of poor data or otherwise unavailable observations, as is often the case for satellite observations of methane column density.b The resulting estimated emission field for the plume observation in a, obtained using the divergence method and emission equation E = ∇ • ⃗ F adv .Our algorithm uses the fourth-order central finite difference to calculate numerical gradients[37], which means that every unavailable pixel in a results in a 4x4 cross of pixels in b where we are unable to calculate E. Consequently, any amount of missing data in a quickly becomes a large amount of missing coverage in b.Integrating over the available pixels in b yields a severely underestimated estimated emission rate Q est .Scaling up Q est by the fraction of missing coverage yields a scaled total emission estimate of 16.2 kg s −1 , and does not correct the underestimation.In Fig.S7in the supplement we demonstrate that this scaling procedure is not sufficient to correct for the problem of missing data, and can actually significantly overestimate the total emission rate (as is the case in this particular example), depending on the fraction of daily missing data and number of repeated observations that go into the time-averaged calculation.In all other estimates of total emission rate in this study, we quote "unscaled" results, i.e., the answer obtained by spatially integrating the estimated emission field over available grid cell values.It is also important to note that individual or time-averaged emission fields estimated via the divergence method in our synthetic study will vary depending on the spatial distribution of missing data, which is a random process.We do not cherry pick particular realisations of spatial distribution for missing data, and examine the variance that is inherent to this random process in Fig.8.

Figure 7 .
Figure 7. a A time-averaged observation of 30 static plumes with constant plume parameters.Each individual plume has a random 30% of its pixels masked (similarly to panel a in Fig. 6).b The resultant average estimated emission field, constructed by time-averaging the estimated emission fields of each of the 30 individual plumes.c The average estimated emission field obtained by first time-averaging all the advective fluxes of the 30 individual plumes, and then taking the divergence of the time-averaged fluxes.The estimated emission field in c is much more resilient to the missing data, and provides a much more accurate estimate of the average emission rate than in b.

Figure 8 .
Figure 8. Underestimation of time-averaged estimated emission rate, plotted as a function of the percentage of observational data present for each individual plume that goes into into the time-averaged estimation.In a, each time-averaged calculation uses 30 individual plumes, each one static with constant wind speed, angle, point source location and emission rate, etc.Each of the 30 simulated observations differs only in the spatial distribution of the masked pixels.In b, each time-averaged calculated uses 30 individual plumes, but the wind speed, wind angle, and source location are randomised.Each of the 30 simulated observations then have the same percentage of pixels randomly masked.This scenario represents a higher degree of complexity, where time-invariant assumptions about sources of emission no longer hold.We find that for both scenarios, it is better to time-average advective fluxes and take the divergence once, as opposed to time-averaging individually estimated emission fields (from the perspective of accurate emission flux estimation).

Figure 9 .
Figure 9. Average methane enhancement over the Permian basin for the years 2019, 2020, and 2021 as observed by TROPOMI.a The Delaware basin.b The Central basin.c The Midland basin.d The Ozana Arch basin.e The Val Verda basin.Basin boundaries taken from the Energy Information Administration.

Figure 10 .
Figure 10.Estimated emission fields for the five Permian subbasins for three years of observations, 2019-2021.a Estimated emission field obtained when daily estimated emission fields are averaged.b Estimated emission field when daily advective fluxes are averaged, with the divergence taken thereof to obtain the estimated emission field.Although the extended three-year time period means that complete spatial coverage is achieved for all the basins, the intermittent missing data on a daily basis means that the methodology of E 1 estimates a significantly lower emission rate than the methodology of E 2 .We find the three-year estimated emission fields E 1 and E 2 to be correlated with r = 0.87 .

Table 1 .
Estimates of Permian methane emission rates [Tg/year].We compare our yearly estimates via the E 2 methodology (where daily methane fluxes are averaged) to those in other literature and find good agreement.Uncertainties on our time-averaged emission estimates are calculated via the algebraic propagation of the daily variance of advective flux of methane at each grid cell.This methodology is described in supplementary section S4.Year This work Veefkind et.al. 2023 * Schneising et.al. 2020 || Liu et.al. 2021 Zhang et.al. 2020 + Calculated using TROPOMI methane mixing ratios from the Weighting Function Modified Differential Optical Absorption Spectroscopy (WFM-DOAS) algorithm version 1.5 [39].1σ uncertainty of 25%.+ Date range for this paper was March 2018 -March 2019.|| Value taken from Veefkind et.al. 2023 * This is average percentage coverage of the Permian basin by the Copernicus Sentinel-5P TROPOMI Level 2 methane data product for each day within the given year.This is not the percentage coverage of the Permian basin by the time-averaged estimated emission field for the year. * This is the percentage coverage of the Permian basin of the time-averaged estimated emission field for this year, and not the average daily coverage of the Permian basin for the year by the methane data product.