High-resolution US methane emissions inferred from an inversion of 2019 TROPOMI satellite data: contributions from individual states, urban areas, and landfills

We quantify 2019 annual mean methane emissions in the contiguous US (CONUS) at 0.25° × 0.3125° resolution by inverse analysis of atmospheric methane columns measured by the Tropospheric Monitoring Instrument (TROPOMI). A gridded version of the US Environmental Protection Agency (EPA) Greenhouse Gas Emissions Inventory (GHGI) serves as the basis for the prior estimate for the inversion. We optimize emissions and quantify observing system information content for an eight-member inversion ensemble through analytical minimization of a Bayesian cost function. We achieve high resolution with a reduced-rank characterization of the observing system that optimally preserves information content. Our optimal (posterior) estimate of anthropogenic emissions in CONUS is 30.9 (30.0–31.8) Tg a⁻¹, where the values in parentheses give the spread of the ensemble. This is a 13 % increase from the 2023 GHGI estimate for CONUS in 2019. We find emissions for livestock of 10.4 (10.0–10.7) Tg a⁻¹, for oil and gas of 10.4 (10.1–10.7) Tg a⁻¹, for coal of 1.5 (1.2–1.9) Tg a⁻¹, for landfills of 6.9 (6.4–7.5) Tg a⁻¹, for wastewater of 0.6 (0.5–0.7), and for other anthropogenic sources of 1.1 (1.0–1.2) Tg a⁻¹. The largest increase relative to the GHGI occurs for landfills (51 %), with smaller increases for oil and gas (12 %) and livestock (11 %). These three sectors are responsible for 89 % of posterior anthropogenic emissions in CONUS. The largest decrease (28 %) is for coal. We exploit the high resolution of our inversion to quantify emissions from 70 individual landfills, where we find emissions are on median 77 % larger than the values reported to the EPA's Greenhouse Gas Reporting Program (GHGRP), a key data source for the GHGI. We attribute this underestimate to overestimated recovery efficiencies at landfill gas facilities and to under-accounting of site-specific operational changes and leaks. We also quantify emissions for the 48 individual states in CONUS, which we compare to the GHGI's new state-level inventories and to independent state-produced inventories. Our posterior emissions are on average 27 % larger than the GHGI in the largest 10 methane-producing states, with the biggest upward adjustments in states with large oil and gas emissions, including Texas, New Mexico, Louisiana, and Oklahoma. We also calculate emissions for 95 geographically diverse urban areas in CONUS. Emissions for these urban areas total 6.0 (5.4–6.7) Tg a⁻¹ and are on average 39 (27–52) % larger than a gridded version of the 2023 GHGI, which we attribute to underestimated landfill and gas distribution emissions.

1 Introduction

All projected pathways that prevent global warming above 1.5 °C require methane emission reductions (IPCC, 2022). The Global Methane Pledge, launched at a 2021 meeting of the United Nations Framework Convention on Climate Change (UNFCCC), aims to achieve a 30 % global reduction in methane emissions from 2020 to 2030 (About the Global Methane Pledge, 2023). Consistent with that goal, the US government is targeting methane emission decreases from oil and gas, livestock, and landfills (The White House, 2021). The UNFCCC requires member parties to report their anthropogenic methane emissions including sectoral contributions from oil and gas, coal, livestock, rice, landfills, and wastewater. The bottom-up approaches used to generate these emission inventories use information on sectoral activity levels and emission factors, but considerable uncertainty can exist in these values. Top-down evaluations of bottom-up inventories use observations of atmospheric methane to infer emissions, often through inverse analyses using a chemical transport model. These top-down emission estimates are most useful for the evaluation of emission inventories and emission mitigation efforts if they achieve high spatial resolution consistent with the information content of the observation–model system. Here, we use column methane observations from the Tropospheric Monitoring Instrument (TROPOMI) aboard the Sentinel-5 Precursor satellite in a reduced-rank analytical inversion to infer methane emissions and the associated information content at 0.25° × 0.3125° (≈ 25 km × 25 km) resolution over the contiguous US (CONUS) for 2019, allowing for detailed analysis of sectoral, state, and urban emissions.

Satellite observations of atmospheric methane column concentrations inferred from measurement of backscattered sunlight in the shortwave infrared have been used extensively in inverse analyses of methane emissions (Streets et al., 2013; Jacob et al., 2022). Previous satellite instruments were limited by large pixel sizes (SCIAMACHY, 2003–2012) or sparse observations (GOSAT, 2009–present). TROPOMI provides daily, global observations of atmospheric methane columns at 5.5 km × 7 km nadir pixel resolution (Hu et al., 2018) with a ∼ 3 % success rate limited by cloud cover, optically dark surfaces, and heterogeneous terrain (Hasekamp et al., 2021). Inversions of TROPOMI data allow for high-resolution quantification of methane emissions but require an understanding of the information content of the observations.

Inverse analyses optimize methane emissions (the state vector) by fitting observations to simulated concentrations from a chemical transport model (CTM) that serves as the inversion forward model. The optimization is typically done by minimizing a Bayesian cost function regularized by a prior emission estimate given by a bottom-up inventory. When a linear relationship exists between emissions and concentrations, as in the case of methane, the optimal (posterior) solution and the associated error covariances and information content can be found analytically (Brasseur and Jacob, 2017). However, this requires the computationally expensive but embarrassingly parallel construction of the Jacobian matrix that represents the relationship between emissions and concentrations in the CTM. This matrix is typically constructed by conducting a CTM perturbation simulation for each optimized emission element, limiting either the spatial resolution of the optimized emissions or the size of the inversion domain. Nesser et al. (2021) demonstrated an alternative method that approximates the Jacobian matrix by perturbing emission patterns that are optimally informed by both the prior emissions and the observations. This approach optimally exploits the information content of the observations, quantifying emissions at the highest resolution possible where the satellite–model observing system provides a constraint and defaulting to the prior estimate elsewhere.

Many inverse studies that have quantified US methane emissions using surface, aircraft, or satellite observations have found large discrepancies with the US Environmental Protection Agency's (EPA) Greenhouse Gas Emissions Inventory (GHGI), which is the bottom-up emission estimate reported by the US to the UNFCCC (EPA, 2023b). Wecht et al. (2014a) found livestock emissions that were 40 % larger than the GHGI for the summer of 2004. Miller et al. (2013) inferred emissions that were 50 % larger than the GHGI for 2007 and 2008, which they attributed to underestimated oil, gas, and livestock emissions. Turner et al. (2015) found similar results for 2009 to 2011. Maasakkers et al. (2021) inferred oil and gas emissions that were 35 % and 22 % higher than the GHGI, respectively, for 2010 to 2015. Lu et al. (2022) found mean 2010–2017 anthropogenic emissions that were 42 % larger than the GHGI, which they attributed largely to oil and gas emissions.

Higher-resolution regional studies have targeted specific aspects of US methane emissions, including contributions from different sectors, states, and urban areas. Karion et al. (2015) found oil and gas emissions in the Barnett Shale in eastern Texas that were consistent with the GHGI when scaled by the region's relative contribution to national gas production but larger than reported by most basin facilities to the EPA's Greenhouse Gas Reporting Program (GHGRP). A series of studies inferred much higher emissions in the Permian Basin than implied by a spatially allocated (gridded) version of the GHGI (Zhang et al., 2020; Schneising et al., 2020; Liu et al., 2021; Chen et al., 2022; Varon et al., 2023). Chen et al. (2018) and Yu et al. (2021) found underestimated livestock emissions in the gridded GHGI in the Upper Midwest. Jeong et al. (2016) inferred Californian emissions that were 20 % to 80 % larger than a state inventory from the California Air Resources Board (CARB). Plant et al. (2019) found methane emissions from six East Coast urban areas in 2012 to be more than 2 times larger than the gridded GHGI.

Here, we use the reduced-rank method of Nesser et al. (2021) in an analytical inversion of 2019 TROPOMI observations to quantify annual mean emissions at 0.25° × 0.3125° resolution over North America using national emission inventories reported by the US, Mexico, and Canada to the UNFCCC as prior estimates. The reduced-rank approach decreases the computational cost by an order of magnitude compared with conventional methods while maximizing the information content from TROPOMI. We focus our analysis on CONUS, paying particular attention to emissions from individual landfills, states, and urban areas. We compare our results to the 2023 GHGI, including estimates for individual states (EPA, 2023b, c). Our inversion provides the first observational evaluation for these state inventories. We also compare our results to inventories prepared by individual states and cities.

2 Data and methods

We conduct an ensemble of inversions of 2019 TROPOMI methane observations over the North American domain shown in Fig. 1 (9.75–60° N, 130–60° W) using the nested GEOS-Chem CTM at 0.25° × 0.3125° resolution as the forward model. The m= 2 919 358 TROPOMI observations are fit to simulated GEOS-Chem concentrations to optimize annual mean methane emissions for 2019 at the 0.25° × 0.3125° GEOS-Chem resolution. This corresponds to n= 23 691 emission grid cells with prior methane emissions larger than 0.1 Mg km⁻² a⁻¹, accounting for over 99 % of North American methane emissions. In a subset of the ensemble, we optimize boundary conditions for the nested GEOS-Chem simulation for each of the four cardinal directions (north, south, east, and west). Methane chemical and soil sinks are not optimized because they are relatively small compared with emissions.

Figure 1 Bottom-up methane emission inventories used as prior estimates for the inversion. Panels show annual mean methane emissions for different sectors. Anthropogenic sectors are given by the gridded versions of the national inventories of Canada (ECCC), the US (EPA GHGI), and Mexico (INECC) reported to the UNFCCC (Maasakkers et al., 2016; Scarpelli et al., 2020, 2021). US oil and gas emissions are updated as described in Sect. 2.2. Wetland emissions are given by the high-performance subset of the WetCHARTs version 1.3.1 wetlands inventory ensemble (Ma et al., 2021), excluding two ensemble members as described in Sect. 2.2. Emissions are shown on the 0.25° × 0.3125° GEOS-Chem grid used for the inversion.

2.1 Reduced-rank analytical inversion

The inversion uses m observed concentrations arranged in a vector y to optimize n gridded emissions arranged in the state vector x by minimizing a Bayesian cost function J assuming normal errors and regularized by the prior emission estimate x_A (Rodgers, 2000):

$$\begin{array}{}\text{(1)}& \begin{array}{rl}J\left(\mathit{x}\right)=& {\left(\mathit{x}-{\mathit{x}}_{\mathrm{A}}\right)}^{\mathrm{T}}{\mathbf{S}}_{\mathrm{A}}^{-\mathrm{1}}\left(\mathit{x}-{\mathit{x}}_{\mathrm{A}}\right)\\ & +\mathit{\gamma }{\left(\mathit{y}-\mathbf{K}\mathit{x}\right)}^{\mathrm{T}}{\mathbf{S}}_{\mathrm{O}}^{-\mathrm{1}}\left(\mathit{y}-\mathbf{K}\mathit{x}\right).\end{array}\end{array}$$

The prior and observing system error covariance matrices S_A and S_O, respectively, are assumed to be diagonal in the absence of better information. The regularization factor γ corrects for the absence of covariance in S_O (Chevallier, 2007). We generate an eight-member inversion ensemble using a range of prior error variances and γ values to capture the inversion's sensitivity to uncertainty in these parameters (Sect. 2.7). The reduced-rank Jacobian matrix $\mathbf{K}=\partial \mathit{y}/\partial \mathit{x}$ represents the sensitivity of concentrations to emissions in the CTM. We construct a rank-k Jacobian matrix for the 0.25° × 0.3125° GEOS-Chem grid by perturbing in the CTM the k emission patterns that best capture the prior emissions and the information content of the TROPOMI observations (Sect. 2.6).

Analytical minimization of the cost function following Rodgers (2000) yields the optimal (posterior) state vector estimate $\widehat{\mathit{x}}$, error covariance matrix $\widehat{\mathbf{S}}$, and information content given by the averaging kernel matrix $\mathbf{A}=\partial \widehat{\mathit{x}}/\partial \mathit{x}=\mathbf{I}-\widehat{\mathbf{S}}{\mathbf{S}}_{\mathrm{A}}^{-\mathrm{1}}$, which describes the sensitivity of the posterior estimate to the true state vector. However, this solution requires inverting the cost function Hessian, which produces numerical instabilities due to the rank reduction of the Jacobian matrix. Here, we use a reduced-rank approximation of the posterior solution following Bousserez and Henze (2018) to solve the inversion on an orthonormal basis that optimally spans the information content of the satellite–forward-model observing system. The basis is given by the eigendecomposition of the prior-preconditioned Hessian of the cost function:

$$\begin{array}{}\text{(2)}& {\widehat{\mathbf{H}}}_p={\mathbf{S}}_{\mathrm{A}}^{\mathrm{1}/\mathrm{2}}{\mathbf{K}}^{\mathrm{T}}{\mathbf{S}}_{\mathrm{O}}^{-\mathrm{1}}{\mathbf{KS}}_{\mathrm{A}}^{\mathrm{1}/\mathrm{2}}=\mathbf{V}\mathbf{\Lambda }{\mathbf{V}}^{\mathrm{T}},\end{array}$$

where the columns of V are the eigenvectors and Λ is a diagonal matrix with entries equal to the eigenvalues. The calculation of ${\widehat{\mathbf{H}}}_p$ requires substantial memory for large m and n, for which we use Dask, a Python parallelization package (Dask Development Team, 2016). The reduced-rank posterior approximation is then generated using the largest k eigenvalues Λ_k and the associated eigenvectors V_k (Bousserez and Henze, 2018):

$$\begin{array}{}\text{(3)}& {\displaystyle }{\mathbf{A}}_{\mathbf{K}}=\mathit{\gamma }{\mathbf{S}}_{\mathrm{A}}^{\mathrm{1}/\mathrm{2}}{\mathbf{V}}_k{\mathbf{\Lambda }}_k{\left({\mathbf{I}}_k+\mathit{\gamma }{\mathbf{\Lambda }}_k\right)}^{-\mathrm{1}}{\mathbf{V}}_k^{\mathrm{T}}{\mathbf{S}}_{\mathrm{A}}^{-\mathrm{1}/\mathrm{2}},\text{(4)}& {\displaystyle }{\widehat{\mathbf{S}}}_{\mathbf{K}}=\left({\mathbf{I}}_n-{\mathbf{A}}_{\mathbf{K}}\right){\mathbf{S}}_{\mathrm{A}},\text{and}\text{(5)}& {\displaystyle }{\widehat{\mathit{x}}}_{\mathrm{FR}}={\mathit{x}}_{\mathrm{A}}+\mathit{\gamma }{\widehat{\mathbf{S}}}_{\mathbf{K}}{\mathbf{K}}^{\mathrm{T}}{\mathbf{S}}_{\mathrm{O}}^{-\mathrm{1}}\left(\mathit{y}-\mathbf{F}\left({\mathit{x}}_{\mathrm{A}}\right)\right).\end{array}$$

Here, ${\widehat{\mathit{x}}}_{\mathrm{FR}}$ approximates the full-rank (FR) posterior $\widehat{\mathit{x}}$ by minimizing the difference between the two, and ${\widehat{\mathbf{S}}}_{\mathbf{K}}$ and A_K are the optimal posterior error covariance and averaging kernel matrices, respectively, for an inversion solved with a reduced-rank forward model. We set k to match the rank of the reduced-rank Jacobian matrix, which is chosen to maximize information content within the available computational resources (Sect. 2.6). The diagonal elements of A_K are often referred to as averaging kernel sensitivities and are a measure of the dependence of the optimized emissions on the prior estimate. Their sum (trace of A_K) gives the degrees of freedom for signal (DOFS) that represent the number of pieces of information independently quantified by the observing system (Rodgers, 2000). The reduced-rank inversion and Jacobian matrix do not attempt to optimize emissions in areas with low information content, so we default to the prior estimate for grid cells with averaging kernel sensitivities less than 0.05 (Nesser et al., 2021).

2.2 Prior estimates and errors

Figure 1 shows the annual-average prior emission estimates for different sectors. Anthropogenic emissions are given by the spatially disaggregated (gridded) versions of the 2016 EPA GHGI for the US for 2012 (Maasakkers et al., 2016), the Instituto Nacional de Ecología y Cambio Climático (INECC) inventory for Mexico for 2015 (Scarpelli et al., 2020), and the Environment and Climate Change Canada (ECCC) inventory for Canada for 2018 (Scarpelli et al., 2021). We update the distribution and magnitude of GHGI oil and gas emissions to the 2020 GHGI for 2018 following Shen et al. (2022) and use the Environmental Defense Fund's inventory for the Permian Basin for 2019 (Zhang et al., 2020), where GHGI estimates are known to be too low (Zhang et al., 2020; Schneising et al., 2020; Liu et al., 2021; Chen et al., 2022; Varon et al., 2023). We treat oil and gas as a single sector in our analysis due to significant source co-location and uncertainty in the partitioning of oil and gas wells. The magnitude of GHGI livestock, landfill, and wastewater emissions changed by less than 10 % from 2012 to 2019, while coal emissions decreased by 26 %. The distribution of these sources is unlikely to have changed significantly. Anthropogenic emissions for Central America and the Caribbean islands are from the EDGAR v4.3.2 global emission inventory for 2012 (Janssens-Maenhout et al., 2019). Anthropogenic emissions are assumed to be aseasonal except for manure management and rice cultivation, for which we apply monthly scaling factors as described by Maasakkers et al. (2016) and Zhang et al. (2018), respectively.

Prior monthly emissions for wetlands are given by the high-performance subset of the WetCHARTs ensemble version 1.3.1, which includes the nine ensemble members that best match global GOSAT inversion results (Ma et al., 2021). In an inversion of GOSAT data over North America, Lu et al. (2022) found that this high-performance subset overestimated wetland methane emissions, particularly at high latitudes. We remove the two members (WetCHARTs models 1923 and 2913; Bloom et al., 2017) that are most responsible for this overestimate from the ensemble. Other natural methane emission sources are minor and include open fires, termites, and geological seeps, for which we follow the emissions described in Lu et al. (2022). Methane losses from oxidation and soil uptake are prescribed as in Maasakkers et al. (2019) and are not optimized in the inversion.

We assume uniform relative error standard deviations for the prior emissions of between 50 % and 100 % for the different members of our inversion ensemble, with no error covariance between grid cells. Previous inversions that optimized methane emissions over North America assumed prior error standard deviations up to 50 %. We inflate errors up to 100 % in our ensemble to account for increased errors at high resolution (Maasakkers et al., 2016). Errors for each ensemble member are chosen as described in Sect. 2.7.

2.3 Forward model

We use the nested version of the GEOS-Chem CTM 12.7.1 (https://doi.org/10.5281/zenodo.3676008, Developers of GEOS-Chem, 2020) at 0.25° × 0.3125° resolution over North America as the forward model for the inversion. Earlier versions of the methane simulation have been described by Wecht et al. (2014a) and Turner et al. (2015). The model is driven by GEOS-FP meteorological fields from the NASA Global Modeling and Assimilation Office (Lucchesi, 2017). Methane sinks from OH, Cl, soil uptake, and stratospheric oxidation are as described in Maasakkers et al. (2019). Initial conditions for 1 January 2019 and 3-hourly boundary conditions for the year are specified by methane concentration fields from a global GEOS-Chem simulation at 2° × 2.5° resolution using optimized emissions from a global inversion of TROPOMI observations (Qu et al., 2021).

2.4 TROPOMI observations

TROPOMI has provided daily, global observations of dry column methane mixing ratios at 7 km × 7 km nadir pixel resolution since May 2018 and at 5.5 km × 7 km nadir pixel resolution since August 2019 (Lorente et al., 2021a). TROPOMI measures backscattered solar radiation in the 2.3 µm methane absorption band from a Sun-synchronous orbit with a local overpass time of 13:30 (Veefkind et al., 2012). Methane concentrations are inferred from a full-physics retrieval with a ∼ 3 % success rate limited by cloud cover, variable topography, low or heterogeneous albedo, and high aerosol loading (Hasekamp et al., 2021). We use retrieval v14 as described by Lorente et al. (2021a), which has a −3.4 ± 5.6 ppb bias relative to the Total Carbon Column Observing Network (TCCON). We use only retrievals with quality assessment flag equal to 1.

Previous analyses of TROPOMI data have identified surface artifacts (Barré et al., 2021) and spatially variable biases relative to the more accurate but sparser GOSAT data (Jacob et al., 2022). We filter the data to remove snow- and ice-covered scenes using blended albedo, an empirical parameter developed by Wunch et al. (2011) and suggested for the TROPOMI data by Lorente et al. (2021a). We remove scenes with a blended albedo greater than 0.75 in non-summer seasons. We also remove scenes with albedo in the shortwave infrared of less than 0.05 following de Gouw et al. (2020), which account for most of the remaining unphysical TROPOMI observations (methane mixing ratio less than 1700 ppb), and scenes north of 50° N in winter.

Figure 2 shows the final m= 2 919 358 observations used for the inversion on the GEOS-Chem 0.25° × 0.3125° grid. The data are dense and seasonally consistent across high-emission regions of CONUS (Fig. S1 in the Supplement). The filters preserve 69 % of the high-quality retrievals of TROPOMI v14 and increase the GOSAT–TROPOMI correlation in all seasons, with the largest increases in winter and spring (Fig. S2). Seasonal regional biases decrease by between 7 % and 21 % and are always within the 1 standard deviation range of both the TROPOMI and GOSAT data. Comparison to a GEOS-Chem simulation driven by the prior emissions (Fig. S3) shows a mean aseasonal (GEOS-Chem–TROPOMI) bias of ξ= 9.1 ppb over North America that we attribute to errors in the boundary conditions. This bias can also be fit as a linear function of degrees latitude θ as $\mathit{\xi }=-\mathrm{5.40}+\mathrm{0.39}\mathit{\theta }$. We correct the bias in our inversion ensemble members by removing either the continental mean bias or the latitude-dependent correction from the GEOS-Chem concentrations.

Figure 2 TROPOMI methane observations in 2019. Panel (a) shows the annual average column dry methane mixing ratios for 2019 averaged on the 0.25° × 0.3125° GEOS-Chem grid. Panel (b) shows the number of observations for the year on the same grid. The observations have been filtered as described in Sect. 2.4.

2.5 Observing system errors

The observing system error covariance matrix S_O includes contributions from forward-model, instrument, and representation errors (Brasseur and Jacob, 2017). Forward-model errors include contributions from transport and from random temporal variability unresolved by the prior emissions estimate. We calculate the total observing system error variances using the residual error method (Heald et al., 2004). This method assumes that the mean difference between the TROPOMI observations and the prior GEOS-Chem simulation, calculated here on a seasonal 2° × 2° grid, is caused by errors in emissions that will be corrected by the inversion. The standard deviation of the residual errors after subtracting the mean gridded errors then defines the standard deviation of the observing system errors. We set a minimum error standard deviation of 10 ppb, which applies to 32 % of observations. We find a mean observing system error standard deviation of 11.5 ppb, with the largest errors in winter and at high latitudes. The resulting error variances are the diagonal elements of S_O. Off-diagonal terms are assumed to be zero in the absence of better information, which we account for with the regularization factor γ (Chevallier, 2007). We describe the choice of γ in Sect. 2.7.

2.6 Jacobian matrix

Constructing the Jacobian matrix K for our inversion would normally require conducting a 1-year perturbation simulation for each of the n= 23 691 grid cells optimized. This is computationally intractable. We construct the Jacobian matrix at substantially decreased computational cost using the reduced-rank method introduced by Nesser et al. (2021), which takes advantage of the heterogeneous information content of the TROPOMI observations. This method updates an initial, low-cost estimate of the Jacobian matrix by perturbing the patterns that best explain the information content of the observing system, constructing a reduced-rank Jacobian matrix that optimally preserves information content.

We construct the initial, low-cost estimate of the Jacobian matrix K⁽⁰⁾ using the mass-balance approach described by Nesser et al. (2021). We assume that a perturbation of methane emissions Δx_j in grid cell j (with units kg m⁻² s⁻¹) produces column mixing ratio enhancements Δy_i over grid cell i according to

$$\begin{array}{}\text{(6)}& \mathrm{\Delta }{y}_i={\mathit{\alpha }}_{ij}{\displaystyle \frac{M_{\mathrm{air}}}{M_{{\mathrm{CH}}_{\mathrm{4}}}}}{\displaystyle \frac{Lg}{Up}}\mathrm{\Delta }{x}_j,\end{array}$$

where ${\mathit{\alpha }}_{ij}\in [\mathrm{0},\mathrm{1}]$ is a dimensionless coefficient providing a crude representation of turbulent diffusion; M_air and $M_{{\mathrm{CH}}_{\mathrm{4}}}$ are the molecular weights of dry air and methane, respectively; L is a ventilation length scale equal to the square root of the grid cell area; g is gravitational acceleration; U is the wind speed, taken here as 5 km h⁻¹; and p is the surface pressure, taken here as 1000 hPa. The use of α_ij produces off-diagonal structure in K⁽⁰⁾, which Nesser et al. (2021) found to be necessary for an effective first estimate. We apply a simple isotropic turbulent diffusion scheme in which the influence of emissions spreads linearly to concentric rings of grid cells. This is represented as ${\mathit{\alpha }}_{ij}=\left(\mathrm{8}-∥i-j∥\right)/\mathrm{36}c$, where $∥i-j∥=\left\{\mathrm{0},\mathrm{1},\mathrm{\dots },\mathrm{7}\right\}$ gives the distance in latitude or longitude grid cell index between i and j, 36 is the sum of $∥i-j+\mathrm{1}∥$ values, and c gives the number of grid cells in the corresponding concentric ring. For $∥i-j∥\ge \mathrm{8}$, α_ij=0.

We use K⁽⁰⁾ together with the error covariance matrices S_A and S_O to calculate the initial patterns of information content that are perturbed in the forward model. We calculate the prior preconditioned Hessian (Eq. 2) using K⁽⁰⁾ and perform its eigendecomposition. The resulting matrix of eigenvectors V⁽⁰⁾ is related to the patterns of information content via ${\mathbf{S}}_{\mathrm{A}}^{\mathrm{1}/\mathrm{2}}{\mathbf{V}}^{\left(\mathrm{0}\right)}$, which is equivalent to the eigenvector matrix of the averaging kernel matrix calculated with K⁽⁰⁾ (Bousserez and Henze, 2018). We perturb the k₁= 434 eigenvectors that capture 50 % of the DOFS generated with K⁽⁰⁾. We then apply an optimal operator that restores the original state dimension and minimizes information content loss to yield an updated reduced-rank Jacobian matrix estimate K⁽¹⁾. We recompute the eigenvectors, perturb the k₂= 1952 eigenvectors that explain 80 % of the initial DOFS, and construct the final reduced-rank Jacobian matrix K⁽²⁾. This iterative update scheme optimizes the information content of the posterior solution while reducing the computational cost by an order of magnitude (Nesser et al., 2021).

2.7 Inversion ensemble

The posterior error covariance matrix that results from Bayesian optimization (Eq. 4) does not account for errors in inversion parameters including the prior and observing system error covariance matrices (Houweling et al., 2014). The analytical solution readily allows for the creation of an ensemble of inversions that reflects the sensitivity of the results to the chosen setup, including parameters. Table 1 summarizes our quality-controlled ensemble of inversions. We conduct inversions that do or do not optimize the boundary conditions and apply either a latitudinal or mean bias correction to the prior (model–observation) difference as driven by boundary condition biases. For each inversion, we choose the relative prior error (50 %, 75 %, or 100 %) and regularization factor (between 0.175 and 0.5) so that the prior term of the cost function evaluated at the posterior solution $J_{\mathrm{A}}\left(\widehat{\mathit{x}}\right)={\left(\widehat{\mathit{x}}-{\mathit{x}}_{\mathrm{A}}\right)}^{\mathrm{T}}{\mathbf{S}}_{\mathrm{A}}^{-\mathrm{1}}\left(\widehat{\mathit{x}}-{\mathit{x}}_{\mathrm{A}}\right)$ averages to 1 across all grid cells optimized by the reduced-rank inversion as expected from the chi-square distribution, which J_A(x) follows by definition (Lu et al., 2021). This yields an ensemble of eight quality-controlled inversions with indistinguishable validity. All inversions have few grid cells with negative emissions, most of which are of the same order of magnitude as the soil sink. Unless otherwise noted, our results give the mean posterior emissions for the ensemble, with uncertainty ranges given by the ensemble range.

Table 1 The eight members of the inversion ensemble.

¹ We conduct inversions that either do or do not optimize the boundary conditions. In inversions with optimized boundary conditions, we include four boundary condition elements, corresponding to the northern, eastern, southern, and western borders of the North American domain, in the inversion state vector.
² We also conduct inversions that apply either a latitudinal or mean bias correction to the prior (model–observation) difference. The latitudinal correction fits the bias with a first-order polynomial. In inversions with a mean bias correction, we remove the mean prior (model–observation) difference as driven by boundary condition biases.
³ We balance the prior and observing system errors to avoid overfitting the emissions to the observations. The regularization factor γ is applied to the inverse observing system error covariance matrix ${\mathbf{S}}_{\mathrm{O}}^{-\mathrm{1}}$ so that values less than 1 increase the observing system errors. We choose the value of the regularization factor and the prior error standard deviation for a given inversion so that the prior term of the posterior cost function is approximately 1, as required by chi-square statistics (Sect. 2.7).

2.8 Source attribution

The high resolution of the inversion facilitates the attribution of the posterior emission estimates to individual source sectors or regions, including states and urban areas. We aggregate the native resolution emission and error estimates to the corresponding p sectors, states, or urban areas using a summation matrix $\mathbf{W}\in {\mathbb{R}}^{p\times n}$. The rows of W are given by the relative contribution of each grid cell to each source category. For sectoral attribution, the rows are given by the relative, area-normalized contribution of each grid cell to a given sector in the prior emission estimate. For state attribution, the rows are given by the fraction of each grid cell within a given state. For urban area attribution, the rows have binary values depending on whether the grid cell overlaps with a given urban area. If the grid cell contains multiple urban areas, the fractional contribution of the grid cell to a given urban area is used instead. The reduced-dimension posterior estimate ${\widehat{\mathit{x}}}_{\mathrm{FR},\mathrm{red}}$, posterior covariance matrix ${\widehat{\mathbf{S}}}_{\mathbf{K},\mathrm{red}}$, and averaging kernel matrix A_K,red are then given by

$$\begin{array}{}\text{(7)}& {\displaystyle }{\widehat{\mathit{x}}}_{\mathrm{FR},\mathrm{red}}=\mathbf{W}{\widehat{\mathit{x}}}_{\mathrm{FR}},\text{(8)}& {\displaystyle }{\widehat{\mathbf{S}}}_{\mathbf{K},\mathrm{red}}=\mathbf{W}{\widehat{\mathbf{S}}}_{\mathbf{K}}{\mathbf{W}}^{\mathrm{T}},\text{and}\text{(9)}& {\displaystyle }{\mathbf{A}}_{\mathbf{K},\mathrm{red}}={\mathbf{WA}}_{\mathbf{K}}{\mathbf{W}}^{\ast },\end{array}$$

where ${\mathbf{W}}^{\ast }={\mathbf{W}}^{\mathrm{T}}{\left({\mathbf{WW}}^{\mathrm{T}}\right)}^{-\mathrm{1}}$ is the Moore–Penrose pseudoinverse (Calisesi et al., 2005). In the case of disaggregating our emission estimates to individual landfills, we scale the posterior estimate in the corresponding grid cell by the fraction of emissions attributed to landfills in the prior estimate. These approaches to source aggregation and disaggregation assume that the prior fractional sectoral contributions are correct in each grid cell and that emission sources are evenly distributed in grid cells that cross state lines. The uncertainty of this method is not reflected in the reported error bounds, but the high resolution of our emission estimates decreases the influence of these assumptions relative to coarser-resolution estimates. Newly developed methods use prior and posterior error covariances to improve upon these assumptions (Cusworth et al., 2021).

3 Results and discussion

Figure 3 shows the ensemble mean posterior scale factors relative to the annual average prior emission estimate as described in Sect. 2.2 (left) and the corresponding averaging kernel sensitivities (right). Grid cells unoptimized by the inversion (mean averaging kernel sensitivity less than 0.05) are left blank. We find 772 (421–1279) DOFS for the domain, where the values in parentheses here and elsewhere are the quality-controlled, eight-member inversion ensemble minimum and maximum, respectively. This represents a large increase in information content relative to past inversions over North America: Lu et al. (2022) found 114 DOFS in a joint inversion of data from GOSAT and the National Oceanic and Atmospheric Administration's (NOAA) GLOBALVIEWplus ObsPack in situ data, while Shen et al. (2022) found 201 DOFS in an inversion of TROPOMI observations over 14 oil and gas basins. This increase reflects both the improved coverage from TROPOMI and the benefit of achieving 0.25° × 0.3125° resolution on the continental scale. Of these DOFS, 641 (350–1058) are found for CONUS, 86 (53–134) for Mexico, and 37 (15–69) for Canada. The high information content for CONUS reflects both the large emissions (Fig. 1) and the high density and consistency of TROPOMI observations (Fig. 2). As a result, we focus our discussion on CONUS. We isolate anthropogenic emissions by removing contributions from wetlands and other natural sources following Sect. 2.8. We compare our posterior emissions to the 2023 EPA GHGI inventory for 2019, including to the emission estimates for individual states (EPA, 2023b, c). We remove emissions from Hawaii and Alaska from the GHGI total using these state estimates, which account for less than 0.5 % of the national total.

Figure 3 Optimization of methane emissions for 2019 by inversion of TROPOMI observations. Panel (a) shows the scale factors relative to the prior estimate for the inversion given by gridded versions of the national anthropogenic emissions inventories for the US (EPA GHGI), Mexico (INECC), and Canada (ECCC), with US oil and gas emissions updated as described in Sect. 2.2, and by WetCHARTs wetland emissions (top left panel of Fig. 1). Panel (b) shows the observing system information content as measured by the averaging kernel sensitivities (the diagonal elements of the averaging kernel matrix). Values of 1 indicate that TROPOMI quantifies emissions independently of the prior estimate, while values of 0 indicate that emissions are not optimized by the inversion. The sum of the averaging kernel sensitivities gives the degrees of freedom for signal (DOFS), shown inset, which defines the number of independent pieces of information quantified by the observing system. Grid cells with averaging kernel sensitivities less than 0.05 are left blank.

We evaluate the inversion results by comparing simulated observations from GEOS-Chem driven by either the prior or the mean posterior emissions to TROPOMI observations and to independent in situ surface and tower observations from NOAA's GLOBALVIEWplus CH₄ ObsPack v3.0 database (Schuldt et al., 2019). We follow Lu et al. (2021) and use only daytime ObsPack observations with outliers excluded. We use monthly average ObsPack observations over CONUS to increase consistency with the annual temporal resolution of our inversion and its distribution of information content. Compared to TROPOMI, both the prior and posterior GEOS-Chem simulations produce similar coefficients of determination (R²) and root-mean-square errors (RMSEs). Compared with ObsPack, the posterior simulation improves upon the prior simulation, increasing the R² value from 0.55 to 0.65 and decreasing the RMSE from 80 to 73 ppb, similar to previous inversions of satellite data (Lu et al., 2021). The broad agreement of both simulations with observations reflects the high quality of the prior emission estimate in North America (Maasakkers et al., 2019).

We also compare the TROPOMI v14 data used here to the most recent data (v19), which have improved bias corrections and performance compared with GOSAT in North America (Balasus et al., 2023). We define the grid cells containing observations sensitive to the optimized emissions by calculating the row-wise sum of the Jacobian matrix weighted by the prior and observing system error standard deviations, limited to emission grid cells with averaging kernel sensitivities greater than 0.05. Of the 95 % of observation grid cells most sensitive to the optimized emissions, only 14 % have an average mean-bias-corrected (v14–v19) difference greater than 5 ppb, while less than 2 % have a difference greater than 10 ppb. We similarly define the observation grid cells that influence the optimized emission grid cells using the Jacobian matrix columns. We find that the average mean-bias-corrected (v14–v19) difference for the observational grid cells influenced by optimized grid cells is −0.05 ppb with a standard deviation of 0.1 ppb, indicating that there is little bias in the observations that influence any single grid cell. We finally find no correlation (R²= 0.03) between our posterior scaling factors and the mean (v14–v19) difference, suggesting that biases in the v14 data do not influence our posterior emissions.

3.1 CONUS sectoral emissions

We find posterior anthropogenic methane emissions of 30.9 (30.0–31.8) Tg a⁻¹ for CONUS in 2019, a 13 % increase from the GHGI estimate of 27.3 (25.1–30.6) Tg a⁻¹, where the values in parentheses represent the GHGI 95 % confidence interval (EPA, 2023b). This estimate excludes Alaska and Hawaii, which likely represent a small (∼ 1 %) contribution to the national anthropogenic total (Miller et al., 2016; Konan and Chan, 2010). Lu et al. (2022) found larger anthropogenic emissions of 36.2 (32.1–37.6) Tg a⁻¹ over the same domain for 2017 by optimizing emissions and trends in a joint inversion of GOSAT and in situ observations for 2010 to 2017. Worden et al. (2022) found lower anthropogenic emissions of 27.6 (24.3–30.9) Tg a⁻¹ over the US for 2019 by regridding global inversions of GOSAT data that optimized emissions at 2° × 2.5° resolution using uncertainties for the prior and posterior estimates. Deng et al. (2022) reviewed an ensemble of global inversions and found median US posterior anthropogenic emissions for 2019 of 26.5 (20.8–38.7) Tg a⁻¹ with GOSAT data and 31.9 (23.9–43.1) Tg a⁻¹ with in situ data.

We allocate our national total to individual emission sectors using the attribution method described in Sect. 2.8. From the off-diagonal structure of ${\widehat{\mathbf{S}}}_{\mathbf{K},\mathrm{red}}$ (Eq. 8), we find very low posterior error correlation between all anthropogenic and biogenic sectors (mean error correlation coefficients less than 0.2), indicating that we can accurately separate sectoral emissions. Figure 4 and Table 2 summarize the results compared to the GHGI. Livestock, oil and gas, and landfills account for 89 % of posterior anthropogenic emissions, and all of the aforementioned sectors show increases relative to the GHGI. We find a significant decrease from the GHGI only for coal. For these four sectors, we find sectoral averaging kernel sensitivities of between 0.47 and 0.91, larger than the values found by Lu et al. (2022) from GOSAT and in situ data, indicating that TROPOMI constrains most of the emissions from these sources. We find a small but significant increase in wetland emissions that is consistent with the large range found by Lu et al. (2022). However, the reduced-rank observing system only optimizes about half of wetland emissions, with most of the inferred increase limited to the south eastern coast, including South Carolina, Georgia, and eastern Florida.

Figure 4 Sectoral methane emissions in the contiguous United States (CONUS) for 2019. The 2023 EPA GHGI emissions for 2019 (top bars) and posterior estimates given by inversion of TROPOMI data for 2019 (bottom bars) are shown for different sectors. For wetland emissions, we show the WetCHARTs estimate (top bar). The shading corresponds to emissions in grid cells that are optimized by the inversion (grid cells with averaging kernel sensitivities greater than 0.05), while the white represents emissions not optimized by the inversion so that the posterior defaults to the prior estimate. Error bars on the GHGI emissions represent the GHGI 95 % confidence intervals. Error bars on the posterior emissions are given by the spread of the eight-member inversion ensemble. Also shown are independent sectoral emission estimates from previous inversions.

Table 2 The 2019 methane emissions for the contiguous United States (CONUS).

¹ Inventory estimates of sectoral methane emissions. Anthropogenic emissions are given by the EPA 2023 ​​​​​​​GHGI for 2019, with error ranges inferred from the sum in quadrature of bottom-up subsector errors given as 95 % confidence intervals. Wetland emissions are from a subset of the high-performance WetCHARTs ensemble version 1.3.1; see Sect. 2.2 for details.
² Optimized emissions from the inversion of TROPOMI data, with the range from the eight members of the inversion ensemble shown in parentheses.
³ The sensitivity of the posterior emissions to the observing system is given by the diagonal elements of the sectoral averaging kernel matrix calculated as described in Sect. 2.8. The values in parentheses give the range of the inversion ensemble. Values range from 0 (no sensitivity) to 1 (full sensitivity).

Landfill emissions show the largest relative and absolute increase from the GHGI for 2019. We find posterior emissions of 6.9 (6.4–7.5) Tg a⁻¹, a 51 % increase relative to the GHGI estimate of 4.6 (3.7–5.8) Tg a⁻¹. Lu et al. (2022) found similar posterior landfill emissions of 7.5 (5.9–7.7) Tg a⁻¹ for 2017. We attribute the GHGI underestimate to two components of the GHGRP landfill inventory methodologies that produce key inputs for the GHGI, which we discuss in detail in Sect. 3.2. First, for landfills with gas recovery systems, the GHGRP assumes overly high collection efficiencies. Second, the GHGRP does not account for site-specific operations that may produce anomalous emissions.

Coal mining emissions of 1.5 (1.2–1.9) Tg a⁻¹ exhibit the largest decrease in sectoral emissions relative to the GHGI estimate of 2.1 (1.9–2.5) Tg a⁻¹. Lu et al. (2022) found much larger posterior emissions of 2.9 (2.3–3.4) Tg a⁻¹ for 2017, and Worden et al. (2022) found similar values of 2.8 ± 0.4 Tg a⁻¹ for 2019. Compared with these studies, we achieve a stronger constraint on coal emissions as measured by averaging kernel sensitivities, reflecting the increased coverage from TROPOMI compared with GOSAT. Our lower estimate better reflects the 30 % decrease in CONUS coal production from 2012 to 2019 (EIA, 2021), which is also shown in the 30 % decrease in GHGI coal emissions over the same period (EPA, 2023b). As expected, emissions correlate with underground coal mining: Appalachia generated 56 % of US coal from underground mines in 2019 and 64 % of posterior emissions from coal, while the Illinois Basin yielded 30 % of US underground coal and 20 % of posterior emissions (EIA, 2021).

Livestock emissions show broad agreement with the GHGI, with posterior emissions of 10.4 (10.0–10.7) Tg a⁻¹ representing an 11 % increase from the GHGI estimate of 9.4 (8.5–10.7) Tg a⁻¹. Lu et al. (2022) found similar mean posterior livestock emissions of 10.4 (8.8–11.6) Tg a⁻¹ over CONUS for 2017, and Worden et al. (2022) found similar values of 9.9 ± 0.4 Tg a⁻¹ for 2019. Yu et al. (2021) conducted a seasonal inversion of aircraft observations over the northern central US and southern central Canada for 2017–2018 and found mean posterior livestock emissions of 5.5 (5.1–6.2) Tg a⁻¹, which agrees with our livestock estimate of 5.4 (5.2–5.6) Tg a⁻¹ over the same region. Despite agreement with total GHGI livestock estimates, we find a significant increase in manure management emissions from 2.3 (1.9–2.8) Tg a⁻¹ to 3.1 (2.9–3.2) Tg a⁻¹, which would almost entirely explain the observed discrepancy between the mean GHGI and posterior emissions. The increase in manure management emissions is concentrated over the Californian Central Valley, northern Iowa, and Sampson and Duplin counties in North Carolina. California is home to more dairy cattle than any other state, Iowa is the largest pork-producing state, and Sampson and Duplin counties are the two largest pork-producing counties in CONUS (USDA, 2019). We find no correlation between our inferred increase and dairy cattle or hog populations, which could reflect variability in manure management practices.

Posterior oil and gas emissions are 10.4 (10.1–10.7) Tg a⁻¹, a 12 % increase from the GHGI estimate of 9.3 (8.1–10.6) Tg a⁻¹. Lu et al. (2022) found much larger posterior emissions of 4.8 (3.1–4.9) Tg a⁻¹ for oil and 8.9 (8.0–9.8) Tg a⁻¹ for gas in 2017, and Lu et al. (2023) used the same inversion framework to find even larger total oil and gas emissions of 15.6 (12.8–17.1) Tg a⁻¹ for 2019 driven by increased emissions in the Anadarko, Marcellus, Barnett, and Haynesville shales. Although we find good agreement on average with the basin-level emissions from Lu et al. (2023), we find much smaller emissions in the Anadarko and Marcellus shales, as shown in Fig. S4. This difference likely results in part from the use of lognormal prior errors in Lu et al. (2023). Compared with Lu et al. (2022, 2023), Worden et al. (2022) found smaller 2019 emissions in the US for oil of 2.4 ± 0.3 Tg a⁻¹ and for gas of 7.9 ± 0.9 Tg a⁻¹, and Shen et al. (2022) found oil and gas emissions of 12.6 ± 2.1 Tg a⁻¹ from an inversion of TROPOMI data over 14 North American basins extrapolated to the national scale for May 2018 to February 2020. Both of these emission estimates are within the uncertainty range of our posterior estimate. We also find consistent basin-level results with Shen et al. (2022), as shown in Fig. S4. Emissions for all posterior basins but one are within 0.25 Tg a⁻¹ of Shen et al. (2022) and all but six are within 0.10 Tg a⁻¹. In particular, we find agreement within error bars in the Haynesville, Barnett, and Anadarko shales. Of the basins where posterior emissions exceed the 0.5 Tg a⁻¹ threshold defined by Shen et al. (2022) for successful quantification of basin emissions by TROPOMI, we find significant differences only in the Permian Basin, where we find smaller emissions of 2.8 (2.8–2.9) Tg a⁻¹. Our Permian estimate is consistent within error bars with Lu et al. (2023) and with other recent studies when basin extent differences are accounted for (Zhang et al., 2020; Schneising et al., 2020; Liu et al., 2021; Varon et al., 2023; McNorton et al., 2022; Veefkind et al., 2023).

3.2 Landfill emissions

We consider in more detail the 51 % increase in our posterior landfill emissions relative to the GHGI. GHGI landfill estimates scale up the total emissions reported to the GHGRP to account for non-reporting landfills (EPA, 2023b). The GHGRP reporting requirements applied to 1297 landfills emitting more than 1 Gg a⁻¹ across the US in 2019 (EPA GHGRP, 2019), over 500 of which had gas recovery systems (EPA LMOP, 2019). The GHGRP requires that landfills use two methods to report emissions. Facilities without gas collection use two approaches that rely on landfill attributes and a first-order decay model based on landfilled mass so that emissions peak the year after waste disposal. However, a survey of 128 California landfills with gas recovery systems found that methane was produced at relatively constant rates over time (Spokas et al., 2015). Landfills with gas collection use one of these methods with recovered methane removed from the modeled emissions in addition to a back-calculation approach that estimates emissions as a function of recovered methane given an estimated collection efficiency based on cover and operation methods. A default efficiency of 0.75 is assumed if cover information is unavailable (EPA, 2023b). Both the model and back-calculation methods have high uncertainties and have not been field validated (NAS, 2018).

We compare our posterior landfill emissions to individual GHGRP facilities that reported more than 2.5 Gg a⁻¹ methane in 2019. Of these 531 landfills, we limit our analysis to the 87 0.25° × 0.3125° grid cells where TROPOMI provides an averaging kernel sensitivity of at least 0.20 and where landfills explain at least 50 % of prior emissions so that we are confident of our ability to separate landfill emissions from other sources. We exclude 33 facilities in grid cells containing multiple landfills because we are unable to separate the individual contributions to total emissions. Figure 5 shows the posterior emissions and corrections to the GHGRP for the remaining 70 facilities, Table 3 shows GHGRP and posterior information for the top 10 methane-producing landfills as ranked by posterior emissions, and Table S3 shows GHGRP and posterior information for all 70 facilities.

Figure 5 Methane emissions for 2019 from 70 individual landfills that report methane emissions of 2.5 Gg a⁻¹ or more to the EPA's Greenhouse Gas Reporting Program (GHGRP) for 2019 and for which our TROPOMI inversion provides site-specific information. The left panel shows the location of the landfills, with insets for parts of California (left) and Illinois and Indiana (right). Posterior emissions for each landfill are shown by the size of the marker. The colors show differences (Δ) between the posterior and GHGRP emissions for 2019, with red colors indicating posterior emissions larger than the reported value. Facilities that collect landfill gas are shown as circles, whereas others are shown as diamonds. The numbers (1 to 10) identify the top 10 methane-producing landfills listed in Table 3, and the letters (a to i) identify the nine validation sites listed in the right panel and outlined in gold. Validation sites are landfills with independent estimates from aircraft campaigns, as listed in the legend. Cambaliza et al. (2015) based their estimates on data from 2011, Duren et al. (2019) on data from 2016 to 2018, Smith (2021) on data from 2019 to 2021, and Catena et al. (2022) on data from November 2021. The right panel shows GHGRP (top bars) and posterior (bottom bars) emissions for the validation sites, along with values reported from the aircraft campaigns. The sites are as follows: (a) South Side Landfill, (b) West Miramar Sanitary Landfill, (c) Seneca Meadows Landfill, (d) Kiefer Landfill, (e) Puente Hills Landfill, (f) Frank R. Bowerman Landfill, (g) Altamont Landfill, (h) Newby Island Landfill, and (i) Keller Canyon Landfill.

Table 3 Top 10 methane-producing landfills in CONUS for 2019.

¹ The top 10 landfills with the largest posterior methane emissions from the TROPOMI inversion for 2019. Numbers correspond to the labels in Fig. 5.
² Emissions reported by individual landfills to the EPA GHGRP for 2019 in gigagrams per year.
³ Posterior emissions from the inversion of TROPOMI observations in gigagrams per year. Posterior emissions are allocated to individual facilities as described in Sects. 2.8 and 3.2. Values in parentheses represent the range from the eight-member inversion ensemble.
⁴ For facilities that capture landfill gas, the recovery efficiency is calculated from emissions and recovered methane reported by individual landfills to the EPA LMOP. Facilities that do not capture landfill gas are listed as NA.
⁵ The posterior recovery efficiency is calculated from posterior emissions and the recovered methane reported by individual landfills to the EPA LMOP. Facilities that do not capture landfill gas are listed as NA.

We validate our posterior landfill results by comparison to aircraft-derived estimates for nine facilities, as shown in Fig. 5. Cambaliza et al. (2015), Smith (2021), and Catena et al. (2022) used mass-balance approaches to estimate emissions using observations from 2011, 2019 to 2021, and November 2021, respectively. Duren et al. (2019) used the integrated methane enhancement method with data from 2016 to 2018. We find agreement within error bounds at the Seneca Meadows Landfill in New York (landfill c in Fig. 5; Catena et al., 2022) and at the Kiefer (d), Frank R. Bowerman (f), Altamont (g), Newby Island (h), and Keller Canyon (i) landfills in California (Smith, 2021; Duren et al., 2019). We find much larger emissions than previous studies at the South Side Landfill (a) in Indiana (Cambaliza et al., 2015) and at the West Miramar Sanitary (b) and Puente Hills (e) landfills in California (Smith, 2021; Duren et al., 2019). The discrepancy at the South Side Landfill could reflect changed emissions since 2011, including the construction of a large landfill gas recovery facility beginning in June 2019 (EPA LMOP, 2019). Methane concentrations of 8662 ppm were recorded at a leak at the West Miramar Sanitary Landfill in November 2019 (San Diego Air Pollution Control District, 2019), suggesting that estimates from other years may not be representative of 2019 emissions. The Puente Hills Landfill closed in 2013 but was previously one of the largest landfills in CONUS (EPA GHGRP, 2019). Our landfill attribution approach, which relies on a prior estimate from 2012, may therefore misallocate emissions to the Puente Hills Landfill instead of to co-located oil and gas operations.

We find mean facility emissions of 13 Gg a⁻¹ compared with the GHGRP mean of 7.2 Gg a⁻¹ for the 70 landfills considered here, with a median 77 % increase in reported emissions. The largest increases occur for facilities that capture landfill gas, for which we find a median 204 % increase from the reported values. As reflected in Table 3, we find no correlation (R²= 0.00) between GHGRP emissions and our posterior estimates; this does not improve when we consider only facilities that do or do not capture landfill gas. This implies that the bottom-up approaches used for emissions estimation have little predictability.

For the 38 facilities that recover gas, we use captured methane emissions reported to the EPA Landfill Methane Outreach Program (LMOP) in 2019 as well as posterior and GHGRP emissions to calculate a posterior and reported recovery efficiency, respectively. The average posterior recovery efficiency of 0.50 (0.33–0.54) is much smaller than the GHGRP mean of 0.61, and both are much smaller than the 0.75 default (EPA, 2023b). Across the six landfill gas facilities at the top 10 methane-producing landfills, we find a mean posterior recovery efficiency of 0.33, half the GHGRP value of 0.65. Indeed, four of the six facilities report methane emission and recovery values consistent with efficiencies larger than the 0.75 default. We find a similar but still lower efficiency at the Seminole Road MSW Landfill (landfill 8) and a marginally higher recovery efficiency only at Sampson County Disposal, LLC (11). We find a low correlation (R²= 0.17) between the efficiencies that does not depend on facility size but improves slightly for facilities constructed within the last decade (R²= 0.31).

We consider in detail the 33 facilities for which posterior emissions show a significant 50 % difference from the GHGRP. We find larger emissions for 28 of these facilities, with the largest discrepancies occurring in 9 of the top 10 methane-producing landfills. Three of these nine facilities have experienced significant operational changes in the last decade. The South Shelby (landfill 2 in Fig. 5) and South Side (3) landfills constructed large landfill gas facilities in 2019 (EPA LMOP, 2019; Russell, 2019), suggesting that emissions from gas infrastructure development may be large. The City of Dothan Sanitary Landfill (6) has been full since 2014, when it stopped accepting most trash (Wise, 2019). Reported emissions peaked at 7.4 Gg a⁻¹ that year (EPA GHGRP, 2019), a value almost 5 times smaller than our posterior emissions, suggesting that the first-order decay model is inadequate to reproduce methane emissions over time. We also find a record of air quality and landfill standard violations at these 34 facilities. At the West Miramar Sanitary Landfill (10, b), a leak emitting 8662 ppm methane was recorded in November 2019 (San Diego Air Pollution Control District, 2019). The Sussex County Landfill in Virginia was fined USD 99 000 in 2016 for failing to address cracks in the landfill cover (Vera, 2016). Lastly, the Newby Island Landfill (h), received 30 violation notices from 2014 to 2020, including for gas collection system shutdowns (Bay Area Air Quality Management District, 2022).

There are five facilities for which our posterior emissions are significantly smaller by 50 % than the 2019 GHGRP. Three report large decreases in estimated methane emissions from 2019 to 2020 that result from changed methodology (EPA GHGRP, 2019). The updated estimates are consistent with our posterior emissions within error estimates in two cases and within 30 % of our posterior emissions in the third case.

3.3 State emissions

The EPA recently began disaggregating the GHGI by state. The EPA uses the same methods to calculate state emissions as in the national inventory so that the total emissions are the same in both estimates. State estimates are developed without reference to greenhouse gas inventories prepared by state governments, which may result in discrepancies in sectoral or total values due to different methods or accounting (EPA, 2023c). In addition to the GHGI state estimates, the EPA provides references to the independent inventories of 24 states and Washington, DC (EPA, 2023a). Of these, we find that eight produce a methane emission estimate separate from their inventory of total CO₂-equivalent greenhouse gases.

We partition our anthropogenic gridded posterior emission estimates, excluding offshore emissions, to each of the 48 states in CONUS (Sect. 2.8) and compare the results to the GHGI state estimates and to inventories prepared by state governments. Figure 6 shows the results for the 29 states responsible for 90 % of posterior CONUS anthropogenic emissions excluding offshore emissions and ordered by posterior emissions, and Table S1 shows the full results for all 48 CONUS states. TROPOMI provides a strong constraint at this resolution, with most state averaging kernel sensitivities larger than 0.5. Our state emissions are on average 7 % larger than the GHGI estimates and 27 % larger in the top 10 methane-emitting states, which produce 55 % of CONUS posterior emissions. Oil and gas emissions on average generate 37 % of posterior emissions and 48 % of the observed increase relative to the GHGI in these 10 states. In Texas, Oklahoma, New Mexico, and Louisiana, the oil and gas sector explains more than 60 % of posterior emissions, with emissions concentrated in the Permian Basin, the Haynesville Shale, and the Anadarko Shale. Livestock and landfills also play a significant role in these states. Emissions in California and Iowa are dominated by the livestock sector, with much of the observed increase relative to the GHGI attributed to manure management emissions (Sect. 3.1). Landfills account for 41 % of posterior emissions in Illinois and 62 % in Florida. Indeed, three of the 10 largest landfills as reported to the GHGRP in 2019 are in Florida (EPA GHGRP, 2019). Consistent with our sectoral analysis, the largest posterior emission decreases relative to the GHGI are found in coal-producing states, including West Virginia and Pennsylvania. While we find a large decrease compared with the GHGI in Pennsylvania, we cannot confidently attribute the difference to a specific sector due to co-location of oil, gas, and coal facilities at the resolution of our inversion.

Figure 6 Anthropogenic methane emissions in 2019 for the 29 states responsible for 90 % of US anthropogenic posterior emissions. The bottom panel shows 2023 EPA GHGI state estimates for 2019 (left bar) and our posterior estimates from the inversion of TROPOMI data (right bar) divided by sector. States are listed from largest to smallest posterior emissions. The information content from the TROPOMI data, as defined by the reduced-form averaging kernel sensitivities (the diagonal elements of the reduced-form averaging kernel matrix; Sect. 2.8), is shown in the top panel. Values of 1 indicate full sensitivity to TROPOMI, whereas values of 0 indicate no sensitivity. The error bars give the spread from the eight-member inversion ensemble. Also shown are emissions estimates from independent state inventories referenced by EPA (2023c).

We consider in more detail Texas and California, which are responsible for 21 % and 7 % of posterior CONUS anthropogenic emissions, respectively. Our posterior estimate for Texas is 6.3 (6.1–6.5) Tg a⁻¹, a 58 % increase from the GHGI estimate of 4.0 Tg a⁻¹. This increase is attributed almost entirely to the oil and gas sector, which accounts for 69 % of posterior emissions compared to 57 % in the GHGI. The Permian Basin alone explains almost 40 % of Texas' posterior emissions. In California, we find posterior emissions of 2.1 (2.0–2.1) Tg a⁻¹, 53 % of which occur in the San Joaquin Valley Air Basin. Our posterior emissions increase 27 % from the GHGI estimate of 1.6 Tg a⁻¹ and 34 % from an independent estimate produced by CARB of 1.5 Tg a⁻¹ (CARB, 2023). Our posterior estimate is smaller than but consistent within error bars with a value of 2.4 ± 0.5 Tg a⁻¹ found by an inversion of in situ observations in California from June 2013 to May 2014 (Jeong et al., 2016). We find general good agreement with the sectoral partitioning in the GHGI, the CARB inventory, and Jeong et al. (2016). Livestock explain 54 % of emissions in our posterior estimate, 47 % in the GHGI, 55 % in the CARB inventory, and 54 % in Jeong et al. (2016), while landfills explain 25 %, 22 %, 21 %, and 19 % of emissions, respectively. We find slightly smaller relative contributions from the oil and gas sector, which accounts for 11 % of emissions in our posterior estimate compared with 14 %, 17 %, and 18 % in the GHGI, the CARB inventory, and Jeong et al. (2016), respectively. This partitioning differs from that found in an inversion of the 2010 CalNex aircraft campaign observations, where 30 % of emissions were attributed to livestock, 38 % to landfills, and 22 % to oil and gas based on the sectoral distribution of the EDGAR v4.2 methane emission inventory (Wecht et al., 2014b).

We also compare our posterior emissions to independent state greenhouse gas inventories from Pennsylvania, Louisiana, Iowa, and Colorado referenced by EPA (2023a), where we have a strong constraint from the inversion (state averaging kernel sensitivity greater than 0.5). Our posterior agrees with the Pennsylvania estimate (Pennsylvania DEP, 2022), but we find a source shift from fossil fuels (from 76 % in the inventory to 63 % in our work) to landfills (from 3 % in the inventory to 16 % in our work). We find that Louisiana's state inventory (Dismukes, 2021) is too low due to underestimated oil and gas emissions, while Iowa's (Iowa DNR, 2020) is too low due to underestimated livestock emissions, particularly from manure management (Sect. 3.1). Colorado's state inventory (Taylor, 2021) is 65 % larger than our posterior estimate due to oil and gas emissions that are more than twice as large.

3.4 Urban area emissions

Urban areas are home to 81 % of the US population (U.S. Census Bureau, 2010) and are major sources of greenhouse gas emissions, including methane (Gurney et al., 2015; Hopkins et al., 2016). As urban populations grow (Seto et al., 2012), these emissions are likely to increase. Cities are well positioned to address methane emissions through waste-reduction initiatives, leak-detection programs, and strategic contracts with landfill operators and gas utilities. Regulation by air pollution control districts can also aid urban emission reduction efforts (Hopkins et al., 2016). C40, a performance-based coalition of over 100 mayors dedicated to climate change mitigation, recommends that cities target a 50 % reduction in methane emissions by 2030 (C40, 2022b). Numerous cities, including New York City, Los Angeles, and Philadelphia, are working toward these reductions through zero-waste programs (C40, 2022a). The US Methane Emissions Reduction Action Plan intends to work with local governments to set up methane monitoring systems to identify and publicize information about municipal gas distribution leaks. The plan also challenges members of the US Climate Mayors to prioritize pipeline abandonment or replacement (The White House, 2021).

We calculate posterior emissions for 95 urban areas across CONUS with 2010 populations over 1 million and averaging kernel sensitivities from our inversion greater than 0.2, providing the first comprehensive national analysis of urban methane emissions. Quantification of urban emissions depends significantly on the definition of city extent due to the presence of large emitters such as landfills on the urban periphery (e.g., Balashov et al., 2020; Plant et al., 2022). We follow Plant et al. (2022) and others in using the US Census Topographically Integrated Geographic Encoding and Referencing system (TIGER)/Line urban areas to standardize the definition across CONUS (U.S. Census Bureau, 2017). These urban areas are responsible for almost a quarter of the GHGI emissions spatially allocated using the gridded inventory from Maasakkers et al. (2016). The gridded inventory does not include post-meter emissions introduced in later versions of the GHGI, which we distribute by population for this analysis. In an average city, the gridded GHGI emissions originate from landfills (40 %), gas distribution (9 %, including 4 % from post-meter emissions), wastewater (6 %), and other sources that are not specific to urban areas such as livestock and oil and gas production and transmission (45 %).

Anthropogenic posterior emissions in these 95 urban areas total 6.0 (5.4–6.7) Tg a⁻¹, 38 (24–54) % larger than the gridded GHGI value of 4.3 Tg a⁻¹. Individual urban area emissions, listed in Table S2, increase by an average of 39 (27–52) %. These increases are much larger than the 13 % increase that we find in total CONUS anthropogenic emissions relative to the GHGI. We are unable to attribute the increased emissions to individual sectors due to source co-location within urban areas at the 0.25° × 0.3125° resolution of our inversion. However, given that landfills account for 40 % of gridded GHGI emissions in an average urban area and that their emissions increase 51 % relative to the GHGI, it is likely that they are responsible for a large fraction of the observed discrepancy. It is also likely that gas emissions, which represent less than 20 % of gridded GHGI emissions in an average urban area but explain between 32 % and 100 % of methane emissions in many cities based on field measurements of methane : ethane ratios (Plant et al., 2019; Floerchinger et al., 2021; Sargent et al., 2021), are significantly underestimated. Finally, recent studies have shown large underestimates of methane emissions from wastewater treatment in the GHGI (Moore et al., 2023; Song et al., 2023) and over urban areas (de Foy et al., 2023), but increasing wastewater emissions accordingly only accounts for 2 % of our observed discrepancy. City-specific variability prevents further attribution of urban emissions. Indeed, we find no correlation between the posterior emission increase and urban area population, population change from 2000 to 2010, population density, or surface area.

Figure 7 Anthropogenic methane emissions for the largest 10 methane-producing urban areas in the contiguous United States (CONUS) for 2019 as identified by the inversion of TROPOMI data. Urban area extents are given by the US Census Bureau TIGER/Line files (U.S. Census Bureau, 2010). The top bars show prior anthropogenic sectoral emissions from the 2023 EPA GHGI for 2019 spatially allocated following Maasakkers et al. (2016) with post-meter emissions allocated by population. The bottom bar shows posterior emissions from the TROPOMI inversion for 2019. We do not resolve posterior sectoral emissions estimates due to source co-location within urban areas at the scale of the inversion. Total emissions (a), per capita emissions (b), and averaging kernel sensitivities (c) are shown for each urban area. Error bars represent the spread of the eight-member inversion ensemble. Also shown are independent urban emissions estimates.

Figure 7 shows results for the top 10 methane-producing urban areas as ranked by posterior emissions from landfills, gas distribution, and wastewater. These 10 regions explain 35 (34–36) % of anthropogenic posterior emissions across the 95 urban areas considered here. We find a mean increase relative to gridded GHGI emissions of 58 (37–84) %. We also compare our posterior emissions to municipal inventories from New York City and Philadelphia, the only available bottom-up urban methane emission estimates. Our emissions are more than twice as large as these inventories, but this likely results from our consideration of broader urban areas.

Figure 7 also compares our results to 12 top-down studies published since 2015. Most of these focused on New York City or Los Angeles. Almost all the studies used larger definitions of urban area extent, with only Pitt et al. (2022) and Plant et al. (2022) using the US Census designation. Most used aircraft or tower observations to infer emissions by inverting a CTM (Cui et al., 2015; Jeong et al., 2016; Cusworth et al., 2020; Pitt et al., 2022; Yadav et al., 2019, 2023). Kuwayama et al. (2019) used a mass-balance approach, while others used observed methane : CO₂ or methane : CO ratios combined with bottom-up inventories of these gases (Wong et al., 2015; Wunch et al., 2016; Plant et al., 2019). Plant et al. (2022) used the same approach with methane : CO ratios from TROPOMI.

We find generally lower but statistically consistent emissions compared with these studies. Our smaller estimates likely result from our restrictive definition of urban area extent. The only study that used aircraft data to estimate methane emissions within a US Census urban area found 314 ± 96 Gg a⁻¹ in New York City (Pitt et al., 2022), which is very similar to our estimate of 309 (241–417) Gg a⁻¹. Plant et al. (2022) used US Census urban areas but relied on TROPOMI methane : CO ratios. They found slightly larger emissions in Atlanta and Philadelphia and much larger emissions in New York City, but their error bars spanned ranges almost twice as large as the derived emissions, limiting the utility of the comparison. Plant et al. (2019) found larger emissions in New York City and Philadelphia but used larger definitions of urban areas and produced similarly wide error ranges.

We find much lower emissions than these studies only in Los Angeles, a difference that decreases but remains significant when we use the same extent as these studies. We attribute much of the discrepancy to decreasing emissions over time. Methane emissions from the Puente Hills Landfill, previously one of the largest landfills in CONUS, decreased following its closure in 2013 (Yadav et al., 2019). This change is not fully reflected in the estimates of Cui et al. (2015), Wong et al. (2015), or Wunch et al. (2016). Yadav et al. (2023) found that Los Angeles emissions decreased an additional 7 % from January 2015 to May 2020. However, their posterior estimate of 251 ± 5 Gg a⁻¹ for 2019 is still larger than our value of 179 (171–193) Gg a⁻¹.

4 Conclusions

We used TROPOMI atmospheric methane column observations for 2019 to optimize methane emissions at 0.25° × 0.3125° resolution over North America with a focus on the contiguous US (CONUS). The high resolution of our inversion allowed us to quantify emissions from individual landfills, states, and urban areas. We compared our results to the 2023 EPA Greenhouse Gas Emissions Inventory (GHGI), including state-level estimates, for 2019; to emissions reported by individual landfills to the EPA Greenhouse Gas Reporting Program (GHGRP); and to other estimates from states and cities. We find large upward corrections to the GHGI at all scales, which may present a challenge for US climate policies and goals, many of which target significant reductions in methane emissions.

We optimized methane emissions using an analytical inversion of TROPOMI methane observations with the GEOS-Chem chemical transport model run at 0.25° × 0.3125° resolution. The inverse solution, or posterior emission estimate, was obtained through a reduced-rank approximation of the analytical minimum of a Bayesian cost function regularized by a prior emission estimate from a gridded version of the GHGI. The analytical solution characterizes the error and information content of the posterior emissions and supported the generation of an eight-member inversion ensemble. We constructed the Jacobian matrix required for the high-resolution, continent-scale analytical solution by iterative approximation using the emission patterns best informed by the prior emission estimate and the observations. This approach decreases the computational cost of our inversion by an order of magnitude compared with conventional analytical methods while optimally preserving its information content.

We find posterior anthropogenic methane emissions of 30.9 (30.0–31.8) Tg a⁻¹ in CONUS, where the range is given by the inversion ensemble. This is a 13 % increase from the GHGI estimate of 27.3 (25.1–30.6) Tg a⁻¹, where the range is given by the 95 % confidence interval. Emissions for landfills, oil and gas, and livestock explain 89 % of posterior CONUS emissions, and each of these sectors' emissions increase by at least 10 % relative to the GHGI. We find a significant decrease compared with the GHGI only for coal emissions. These increases present a challenge to goals set by the US government to decrease methane emissions from oil and gas, livestock, and landfills.

Most of the total increase from the GHGI to the posterior emissions is attributed to a 51 % increase in landfill emissions. We compare our optimized emissions for 70 individual landfills to those reported to the GHGRP and find a median 77 % increase in emissions relative to reported values. We attribute the underestimated GHGI and GHGRP landfill emissions to standard inventory methods that (1) assume overly high recovery efficiencies at facilities that collect landfill gas and (2) inadequately account for anomalous operating events such as gas leaks or the construction of new landfill gas facilities.

We took advantage of the high resolution of our inversion to quantify emissions for each of the 48 states in CONUS and compare them to the newly available state emission inventories published with the GHGI. We find a 7 % average increase, with a 27 % average increase in the top 10 methane-emitting states. Much of the discrepancy in these 10 states is attributed to increased oil and gas emissions, although livestock and landfills also play significant roles. Texas and California, the two largest methane-producing states, emit 21 % and 7 % of total CONUS anthropogenic emissions in our posterior estimate, respectively. Emissions in Texas increase by 58 % relative to the GHGI almost entirely due to the oil and gas sector. Operations in the Permian Basin alone explain almost 40 % of all posterior emissions in the state. In California, we find a 27 % increase from the GHGI and a 34 % increase from an independent inventory prepared by the California Air Resources Board (CARB). Our sectoral partitioning for California is consistent with both inventories, including 54 % of emissions from livestock, 25 % from landfills, and 11 % from oil and gas.

We also provide a first national analysis of urban methane emissions by calculating emissions for 95 urban areas across CONUS. We find total emissions of 6.0 (5.4–6.7) Tg a⁻¹ across these urban areas, representing a fifth of posterior anthropogenic emissions in CONUS and a 38 (24–54) % increase from the gridded 2023 GHGI value of 4.3 Tg a⁻¹. Urban emissions increase on average by 39 (27–52) % compared with the GHGI. We attribute the observed discrepancy to underestimated landfill and gas emissions. Our urban emission estimates are in general consistent with previous top-down studies except for Los Angeles, which may be attributable in part to decreasing emissions between study periods.