Methane emissions in the United States, Canada, and Mexico: Evaluation of national methane emission inventories and sectoral trends by inverse analysis of in situ (GLOBALVIEWplus CH4 ObsPack) and satellite (GOSAT) atmospheric observations

Xiao Lu1,2, Daniel J. Jacob2, Haolin Wang1, Joannes D. Maasakkers3, Yuzhong Zhang2,4,5, Tia R. Scarpelli2, Lu Shen2, Zhen Qu2, Melissa P. Sulprizio2, Hannah Nesser2, A. Anthony Bloom6, Shuang Ma6, John R. Worden6, Shaojia Fan1, Robert J. Parker7,8, Hartmut Boesch7,8, Ritesh Gautam9, Deborah Gordon10,11, Michael D. Moran12, Frances Reuland13, Claudia A Octaviano Villasana14 10 1 School of Atmospheric Sciences, Sun Yat-sen University, Zhuhai, Guangdong Province, China 2 Harvard John A. Paulson School of Engineering and Applied Sciences, Harvard University, Cambridge, MA, USA 3 SRON Netherlands Institute for Space Research, Utrecht, the Netherlands 4 School of Engineering, Westlake University, Hangzhou, Zhejiang Province, China 15 5 Institute of Advanced Technology, Westlake Institute for Advanced Study, Hangzhou, Zhejiang Province, China 6 Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA, USA 7 National Centre for Earth Observation, University of Leicester, Leicester, UK 8 Earth Observation Science, Department of Physics and Astronomy, University of Leicester, Leicester, 20 UK 9 Environmental Defense Fund, Washington, DC, USA 10 RMI, New York, NY, USA; 11 Watson Institute for International and Public Affairs, Brown University, Providence, RI, USA 12 Environment and Climate Change Canada, Toronto, ON, Canada 25 13 RMI, Boulder, CO, USA 14 Instituto Nacional de Ecología y Cambio Climático (INECC), Mexico City, Mexico


60
Atmospheric methane (CH4) is the most important anthropogenic greenhouse gas after carbon dioxide (CO2). Natural emissions are mainly from wetlands. Anthropogenic emissions are from many sectors including the oil/gas supply chain, coal mining, livestock, and waste management. Individual countries must report their anthropogenic methane emissions by sector to the United Nations in accordance with the United Nations Framework Convention on Climate Change (UNFCCC, 1992). These national 65 emission inventories are mainly constructed by bottom-up methods as the product of activity data and emission factors, following methodological guidelines from the Intergovernmental Panel on Climate Change (IPCC). The emission factors are highly variable and have large uncertainties, leading to errors in estimating national emissions, their trends, and the contributions of different sectors (Kirschke et al., 2013;Saunois et al., 2020). Top-down methods involving inversion of atmospheric methane observations 70 can usefully diagnose these errors (Houweling et al., 2017) Maasakkers et al. (2016) to enable its evaluation using top-down methods. Results using analysis of atmospheric methane measurements from ground, aircraft, and satellite platforms show larger methane emissions than reported in the GHGI, particularly for the oil/gas industry 80 (Alvarez et al., 2018;Zhang et al., 2020;Lu et al., 2021;Maasakkers et al. 2021;Qu et al., 2021) and for livestock Yu et al., 2021). Atmospheric observations also suggest an increasing trend of US anthropogenic emissions over the past decade (Turner et al., 2015;Sheng et al., 2018a;Lan et al., 2019;Maasakkers et al., 2021), while the GHGI indicates a decrease (EPA, 2021). 85 Anthropogenic methane emissions for Canada are reported yearly by Environment and Climate Change Canada (ECCC) (ECCC, 2020a; as part of the National Inventory Report (NIR). Atmospheric observations again indicate an underestimate of emissions from oil/gas production (Atherton et al., 2017;Johnson et al., 2017;Chan et al., 2020;Baray et al., 2021;Lu et al., 2021;Tyner and Johnson, 2021) but a decrease of these emissions over the past decade Maasakkers et al., 2021) for selected years (INECC and SEMARNAT, 2018). The last communication to the UNFCCC was in 2015 and this inventory was allocated to a 0.1 o ×0.1 o grid by Scarpelli et al. (2020). A recent inverse analysis of satellite data finds oil/gas emissions to be underestimated by a factor of 2 over eastern Mexico .
The above top-down studies except for Baray et al. (2020) and Lu et al. (2021) used either in situ or satellite observations but not both. Satellite observations have better data coverage but are less sensitive to emissions and have larger uncertainties, particularly at high latitudes. In a previous inverse analysis , we showed that in situ and satellite observations provide complementary global information for inverse analyses of methane emissions. That inversion was conducted at 4°×5° resolution, 105 too coarse for specific evaluation of national inventories.
Here we apply extensive in situ observations from surface sites, towers, ships, and aircraft (GLOBALVIEWplus CH4 ObsPack data compilation) together with the Greenhouse Gases Observing Satellite (GOSAT) observations, in an inverse analysis for 2010-2017 to optimize methane emissions and 110 their year-to-year variability at up to 0.5×0.625° resolution for North America. We use as prior estimates the gridded national emission inventories from EPA (US), ECCC (Canada), and INECC (Mexico), so that our results can inform inventory improvement planning at the emission sector level. Following Lu et al. (2021), we use an analytical inversion method that provides closed-form characterization of error statistics and information content of the inverse solution, and also allows us to compare quantitatively the 115 information from the in situ and satellite observations.

Methods
We use methane observations from the GLOBALVIEWplus CH4 ObsPack in situ data (Section 2.1) and/or GOSAT satellite retrievals (Section 2.2) with the GEOS-Chem chemical transport model (Section 2.4) as 120 the forward model, to optimize a state vector of mean methane emissions for individual years (Section 2.3) covering the North American continent at a spatial resolution of up to 0.5°×0.625°. We derive posterior estimates of the state vector and the associated error covariance matrix by analytical solution to the Bayesian optimization problem (Section 2.5). Our base inversion uses GOSAT + in situ observations and our best choices of inversion parameters. We also present results from an ensemble of sensitivity 125 inversions using observation subsets (in situ or GOSAT) and varying inversion parameter assumptions (e.g. different error distributions). We attribute inversion results to different methane emission sectors with the methodology described in Section 2.6. 130 We use the comprehensive database of in situ (surface, tower, shipboard, and aircraft) methane observations over North America for 2010-2017 from the GLOBALVIEWplus CH4 ObsPack v1.0 product compiled by the National Oceanic and Atmospheric Administration (NOAA) Global Monitoring Laboratory (Cooperative Global Atmospheric Data Integration Project, 2019). Following Lu et al. (2021), data from surface and tower sites are sampled only during daytime (10-16 local time) and averaged as 135 daytime mean values on individual days for use in the inversion. For sites with standard deviations larger than 30 ppb, we exclude data points that depart by more than two standard deviations from the mean because such conditions are difficult to simulate with the transport model. For other sites we exclude data points that depart by more than three standard deviations from the mean. We also exclude aircraft measurements higher than 9 km above sea level as these measurements would have weak sensitivity to 140 surface fluxes. The in situ observations thus include 49742 data points from surface sites, 15285 from towers, 56 from ship cruises, and 26620 from aircraft campaigns over North America and adjacent waters (Figure 1a). The number of available in situ observations per year increases from 10830 in 2010 to 13593 in 2017. All 145 these in situ data points are used in the base inversion to optimize methane emissions for individual years. We also conduct sensitivity inversions by only using surface and tower sites with continuous eight-year records for trend analyses.

150
The GOSAT satellite launched in 2009 measures the backscattered solar radiation from a sun-synchronous orbit at around 13:00 local time (Kuze et al., 2016). Methane is retrieved in the 1.65 µm shortwave infrared absorption band. We use the column-averaged dry-air methane mixing ratios from the University of Leicester version 9.0 Proxy XCH4 retrieval (Parker et al., 2020a). Comparison with ground-based methane observations from the Total Carbon Column Observing Network (TCCON) shows that the 155 retrieval has a single-observation precision of 13 ppb and an overall global bias of 9 ppb that is removed from the Proxy XCH4 data (Parker et al., 2020a). Here we use a total of 205875 (25734 per year on average) GOSAT retrievals for 2010-2017 over North America in the inversion, excluding glint data over the oceans and data poleward of 60 o which are not representatively sampled and for which errors are large ( Figure 1b).

Prior emission inventories
We use as prior estimates of anthropogenic methane emissions the gridded versions of the official national inventories for the US (EPA, 2016), Canada (ECCC, 2020), and Mexico (INECC and SEMARNAT, 2018) (Maasakkers et al., 2016;Scarpelli et al., 2020. These emissions are listed in Table 1 for individual 165 countries and the spatial distributions for major sectors are shown in Figure 2. We assume no year-to-year trend in the prior emissions, so that trends from the inversion are solely driven by observations. Prior anthropogenic emissions for the contiguous US (CONUS) are 28.7 Tg a -1 . Anthropogenic US emissions outside CONUS (mostly Alaska, not optimized in the inversion) account for only 0.3 Tg a -1 according to Maasakkers et al. (2016). The latest GHGI report from EPA (2021)  from the nine highest-performance members of the WetCHARTs v1.3.1 inventory ensemble , selected for their fit to the global GOSAT inversion results of Zhang et al. (2021). This choice of prior estimate effectively corrects the large overestimates of wetland emissions for North America previously found in inversions of GOSAT and aircraft data when using the overall mean of the WetCHARTs v1.0 ensemble (Sheng et al., 2018b;Maasakkers et al., 2021). We do not include interannual 180 variability from WetCHARTs because it is highly uncertain and we prefer to have it informed by the observations. Unlike in our global inversion , we do not optimize the relative seasonal https://doi.org/10.5194/acp-2021-671 Preprint. Discussion started: 19 August 2021 c Author(s) 2021. CC BY 4.0 License. variation of wetland emissions and instead have it imposed by the prior estimate (Parker et al., 2020b).
Prior estimates of open fire emissions are the daily values for individual years from the Global Fire Emissions Database (GFED) version 4s (van der Werf et al., 2017). Other small natural emissions 185 (seepages, termites) are as described in Lu et al. (2021).

The GEOS-Chem forward model
We use the nested version of the GEOS-Chem 12.5.0 chemical transport model (http://geos-chem.org, last access: 6 April 2021) (Wecht et al., 2014) as the forward model for the inversion. The model is driven 190 by MERRA-2 reanalysis meteorological fields at their native 0.5° × 0.625° resolution (Gelaro et al., 2017). Methane loss from atmospheric oxidation is as described in Lu et al. (2021) but is inconsequential here because it is negligibly slow compared to the timescale for ventilation of the North American domain. Soil uptake of methane is from the MeMo model v1.0 (Murguia-Flores et al., 2018) but is very small and therefore not optimized in the inversion.

195
The GEOS-Chem model simulation is conducted at 0.5° × 0.625° resolution over the North America domain (130-55°W, 15-65°N) ( Fig.1) for the 2010-2017 period, with initial and dynamic boundary conditions archived every 3 hours from a global 2010-2017 simulation at 4°×5° resolution using methane emissions and sinks previously optimized with GOSAT observations . This means that 200 GOSAT observations over North America are used twice, once for the global inversion (along with other observations worldwide) and once for the North American inversion, but this is inconsequential because the sole purpose of the global optimization is to avoid biases in boundary conditions that would cause spurious corrections to emissions within the inversion domain (Wecht et al., 2014). Lu et al. (2021) show that their optimized simulation is unbiased in comparison to global zonal mean observations for 2010-205 2017 but we still find some residual bias for individual years. We therefore optimize the mean boundary conditions for individual years on each side of the domain (north, south, west, east) as part of the North American inversion.

210
Our state vector x to be optimized in the inversion includes spatially resolved emissions in North America and boundary conditions for each year of 2010-2017. Although we could technically optimize methane emissions for each of the 0.5°×0.625° native model grid elements, the observations do not have sufficient coverage to constrain emissions everywhere at that resolution and doing so would introduce large smoothing errors in the inversion (Wecht et al., 2014). Following Turner and Jacob (2015) and 215 Maasakkers et al. (2021), we use instead a Gaussian mixture model (GMM) to determine the emission patterns that can be constrained effectively by the inversion. This is done by projecting the nativeresolution methane emissions onto 600 Gaussian functions optimized to fit the location, magnitude, and distribution of methane emissions for different sectors as given by the prior estimates. Optimal construction of the GMM aggregates regions with weak or homogeneous emissions while preserving

225
The inversion finds the optimal estimate of by minimizing the Bayesian cost function ( ) (Brasseur and Jacob, 2017): where is the prior estimate of , denotes the prior error covariance matrix, is the observation vector, denotes the observation error covariance matrix, is a regularization factor (see below), 230 and ( ) represents the GEOS-Chem simulation of . The GEOS-Chem forward model ( ) as implemented here is strictly linear (because methane sinks are not optimized), so that the model can expressed as = + , where = / represents the Jacobian matrix and is a constant.
Minimizing the cost function (Eq.1) by solving ( ) = 0 yields closed-form posterior estimates of the state vector �, its error covariance matrix � , and the averaging kernel matrix (Rodgers, 2000;235 Brasseur and Jacob, 2017): where in Eq.2 is the gain matrix, The averaging kernel matrix A in Eq. 4 quantifies the sensitivity of the posterior estimate to changes in the true value, and therefore measures the information content provided by the observing system for correcting the prior estimates and returning the true values as posterior estimates. We refer to the diagonal elements of A as the averaging kernel sensitivities, and to the trace of A as the degrees of freedom for 245 signal (DOFS), representing the number of pieces of independent information on the state vector obtained from the observing system (Rodgers, 2000). Our inversion returns the posterior estimates of mean emissions and averaging kernel sensitivities for each Gaussian, and these can be mapped additively to the 0.5 o ×0.625 o grid using their spatial distributions on the grid.

250
Analytical solution to equation (2), and inference of error statistics and information content from equations (3)-(4), requires explicit construction of the Jacobian matrix . We construct by conducting GEOS-Chem simulations where each element of the state vector is perturbed separately. This is readily done computationally as an embarrassingly parallel problem. Analytical solution has several advantages relative to the more widely used variational (numerical) approach. (1) It identifies the true 255 minimum in the cost function. (2) It provides complete explicit forms of the posterior error covariance and averaging kernel matrices. (3) It enables a range of sensitivity analyses at no significant computational cost modifying the inversion parameters and adding/subtracting observations.
To construct the prior error covariance matrix , we assume a 50% error standard deviation for 260 individual Gaussians in the base inversion (and we test the sensitivity to that assumption as will be described later), with no spatial error covariance so that is diagonal. There is necessarily some spatial covariance in the prior estimates since the Gaussians have spatial overlap, and there is also some spatial covariance in the forward model error contributing to , but these are difficult to quantify. The former would underestimate the information content of the observations while the latter would overestimate it. 265 We effectively correct for this using the regularization parameter γ as described below, and we further rely on our inversion ensemble rather than the posterior error covariance matrix to characterize the error in our posterior solution.
The standard assumption of Gaussian error statistics in the cost function of equation (1) is required to 270 achieve an analytical solution but may lead to unphysical negative emissions (Miller et al., 2014) and fail to capture the heavy tail of the emission distribution Frankenberg et al., 2016;Alvarez et al., 2018). We solve this problem by optimizing for ln ( ) instead of , with the error on ln (x) following a normal Gaussian distribution, i.e., lognormal errors for (Maasakkers et al., 2019). The forward model is then nonlinear, so that the solution must be solved iteratively with a transformed 275 Jacobian matrix ′ = ∂ / ∂ln ( ) at each iteration N. Once the original Jacobian matrix = / for the linear model has been computed, we can derive ′ immediately at any iteration by ∂y / ∂ln (x ) = x ∂y / ∂x , where i and j represent the indices of the observation and state vector elements, respectively. The iterative solution is obtained with the Levenberg-Marquardt method (Rodgers, 2000) for each iteration N: , where ′ = ln ( ) with the initial value ′ from the prior estimate, and = 10 is a coefficient for the iterative approach to the solution (Rodgers, 2000). ′ (with diagonal elements denoted by s ′ ) is the prior error covariance matrix for the inversion in log space, and can be derived from the original prior 285 error covariance matrix (with diagonal elements denoted by s ) following (Maasakkers et al., 2019): We adopt as convergence criterion that the maximum difference between ′ + and ′ elements be The above describes our base inversion. We also conduct sensitivity inversions using different error 300 assumptions. This includes 1) using the quadrature sum of error variances for all sectors contributing to a given Gaussian with a cap of 50% following Maasakkers et al. (2021), resulting in a 43% error on average; 2) to 4) using the normal error distributions (then with the linear Jacobian matrix) with 50%, 95%, and the quadrature sum of errors for individual Gaussians as error variances; 5) assuming an error standard deviation of 5 ppb for boundary conditions.

305
The observation error covariance matrix includes contributions from measurement and forward model errors. We compute it following the residual error method originally described by Heald et al. (2004) and previously used by Lu et al. (2021) observation points on the model grid (Heald et al., 2004). The variance of provides the diagonal terms of . The resulting observation error standard deviations average 13 ppb for GOSAT, 26 ppb for surface sites, 39 ppb for towers, 19 ppb for ships, and 22 ppb for aircraft. The observation error is larger for in situ than for satellite observations, even though the in situ measurements are more precise, because the forward model error is larger for local points than for atmospheric columns (Turner et al., 2015). The 320 observation error for in situ observations is dominated by the forward model error while that for GOSAT is dominated by the measurement error.
We do not have sufficient objective information to quantify the error correlation structure of SO and we therefore assume it to be diagonal. This may underestimate because of correlated transport and source 325 aggregation errors in the forward model, as noted above. We follow Zhang et al. (2018) to introduce a regularization factor for the observation terms in the cost function ( ) (Eq. 1) to avoid either overfits or underfits that would result from underestimates of and , respectively. Lu et al. (2021) showed that the optimal value of this regularization factor can be selected such that the sum of the n prior terms in the posterior estimate of the cost function ( (�) = (� − ) − (� − )) has a value ≈ n, which is 330 the expected value from the Chi-square distribution with n degrees of freedom. Here we find that γ = 1 is best for both the in-situ and GOSAT inversions (i.e., no weighting in the inversion). We also conduct a sensitivity inversion using γ = 0.5 for the GOSAT observation terms as adopted in Maasakkers et al. (2021).
335 Table 2 summarizes the settings of our base inversion (in bold) and the inversion ensemble. The ensemble comprises 33 inversions using the different combinations of settings in the Table. The base inversion including GOSAT + in situ data represents our best estimate, but we will compare it prominently to the GOSAT-only and in-situ-only inversions with the same inversion parameters in order to evaluate the contributions from the different observing platforms for optimizing emissions. We will use the other 340 ensemble members to discuss the sensitivity of inversion results to the choices of observations and inversion parameters, and to define the range of uncertainty in the inversion results.

Sectoral attribution and aggregation of inversion results
The inversion returns the posterior estimates of mean emissions for each of the Gaussians, and we allocate 345 these emissions to the native 0.5 o ×0.625 o model grid by summing the contributions of all Gaussians on the grid. This defines a correction factor f0 to total prior emissions for each 0.5 o × 0.625 o grid cell and including the contributions from all q emission sectors (in our case q = 12, cf. Table 1). For sectoral attribution of this total correction factor we need to derive the correction factors fi to the individual sectors ∈ [1, ] contributing to f0. We use two alternative methods for this purpose. The first method simply 350 takes fi = f0 for all i, thus assuming that all sectors contribute equally to the grid-level correction factor Lu et al., 2021;Zhang et al., 2021) . The second method  accounts for emissions from different sectors having different prior error standard deviation σi and therefore contributing differently to f0. Following Shen et al. (2021), is then given by:

355
where is the fraction of total emissions in the grid cell contributed by sector i, is the prior error standard deviation for total emissions in the grid cell, and η = 2 ∑ 2 =1 2 is a normalization factor. For the prior error standard deviations σi on the 0.5 o ×0.625 o grid we use the scale-dependent adaptation by Maasakkers et al. (2016) of EPA sectoral national error estimates. This results in prior error standard deviations of 43% for rice, 66% for wastewater, 51% for landfills, 38% for livestock, 18% for coal, 30% 360 for gas, and 87% for oil emissions. We further use 70% for wetlands (Bloom et al., 2017) and 100% for all other natural sources. These error estimates are solely used to infer fi values in equation (8), so that more uncertain emissions will contribute more to the correction. We use the second method in our base attribution of posterior estimates to emission sectors but will also use the results from the first method to contribute to error ranges in these sector-attributed posterior estimates. 365 We also need to aggregate posterior emission estimates nationally and by sector for comparison to the national emission inventories. Following Maasakkers et al. (2019), this is done by a transformation from the posterior full-dimension state vector � to the reduced state vector � (national emission for a given sector) with a summation matrix W: . The posterior error covariance and averaging kernel matrices for the reduced state vector are then given by Calisesi et al., 2005). � enables us to determine whether national correction factors for individual sectors are affected by error correlations between sectors. enables us to determine the ability of the observing system to quantify national emissions from a particular sector   Figure 3a shows the gridded posterior correction factors from the base inversion averaged over 2010-2017, i.e., the multiplicative factors applied to the total prior emissions from Figure 2, mapped on the 0.5 o ×0.625 o model grid. Figure 3b shows the corresponding averaging kernel sensitivities, indicating the 385 dependence of the posterior solution on the prior estimate (0 = total dependence, 1 = no dependence). The DOFS = 114 can be placed in the context of the 600 Gaussian state vector elements used to optimize the spatial distribution of emissions. We see that the observations provide considerable information to optimize methane emissions but we also see that a finer resolution for the inversion would not be justified on the continental scale.

390
Figures 3c-f show the results from the GOSAT-only and in-situ-only inversions, enabling us to compare the information contents and consistency of the two data sets. The GOSAT-only inversion yields a DOFS of 68, while the in-situ-only inversion yields a DOFS of 80, even though there are 50% fewer in situ observations than GOSAT observations. This is because the sensitivities of surface observations to 395 emissions are an order of magnitude higher than those of satellite observations (Cusworth et al., 2018). The GOSAT observations have the advantage of broader coverage. Thus we find that the in-situ observations dominate the information content of the base inversion over California, the upper Midwest, and Canada; whereas GOSAT dominates the information content in Mexico (where there are no in-situ observations) and most of the western US. GOSAT and in-situ observations contribute comparably in the We evaluated the ability of the base GOSAT + in situ inversion to fit the two observational data sets by comparing 2010-2017 GEOS-Chem simulations with posterior versus prior emissions and boundary conditions. Results are shown in Figure 5. The posterior simulation reduces the model mean bias (MB) at surface and tower measurements from -11 ppb in the prior simulation to -5 ppb, and also narrows the rootmean-square error (RMSE) from 24 to 14 ppb. For GOSAT the improvement is less apparent from the 425 comparison statistics, because the prior simulation already has a low mean bias MB = -0.5 ppb, and the prior RMSE is only 6.9 ppb (which decreases to 6.5 ppb). However, we see from Figure 5 a significant whitening of the noise with reduction of regional-scale biases.

430
Tables 1a-c summarize our inversion results for national 2010-2017 methane emissions by sector in CONUS, Canada, and Mexico. Our best estimates of total anthropogenic + natural emissions from the base inversion are 46.3 (40.2-48.4) Tg a -1 for CONUS, 16.2 (13.5-17.4) Tg a -1 for Canada, and 6.8 (5.4-6.9) Tg a -1 for Mexico. The ranges given in parentheses are from the 33 inversion ensemble members (Table 2). Averaging kernel sensitivities for these total national emissions (the diagonal elements in , 435 section 2.6) are 0.72 for CONUS, 0.60 for Canada, and 0.40 for Mexico, indicate that the GOSAT + in situ observation system informs 72% of the total methane emissions in CONUS, 60% in Canada, and 40% in Mexico. The lower values for Mexico are due to the lack of in situ observations.
We partition these national totals into different sectors as described in Section 2.6, and use the posterior 440 error covariance matrix (equation (10)) to evaluate the ability of the inversion to separate between sectors. This is shown in Figure 6 as the posterior error correlation matrix, displaying the error correlation coefficients (r) in the inversion results for all sector pairs. Error correlation coefficients are generally lower than 0.2 for CONUS, indicating successful separation, except for small sources (termites, seeps, other anthropogenic). The same holds for Canada except for error correlation between landfills and 445 wastewater treatment, both associated with urban areas. Anthropogenic emissions in Canada are well separated from the large wetland emissions. Error correlations are higher in Mexico, because emissions from different sectors tend to be concentrated in Mexico City and the eastern part of the country (Scarpelli et al., 2020), but even there the error correlation coefficients are generally less than 0.4. Optimization of the oil/gas sector is well separated from the other sectors in all three countries. 450 We find that anthropogenic methane emissions for all three countries are larger in our inversion results than in the national inventories submitted to the UNFCCC. Our best estimate of the mean 2010-2017 anthropogenic methane emission for CONUS is 36.9 (32.5-37.8) Tg a -1 , which is 30% higher than the 28.7 Tg a -1 in the 2016 version of the EPA GHGI used as prior estimate (EPA, 2016), and 42% higher Our best estimate of the mean 2010-2017 anthropogenic methane emission for Canada is 5.3 (3.6-5.7) Tg a -1 , which is 43% higher than the 3.7 Tg a -1 in the ECCC NIR (2020 version) used as prior estimate, and 33% higher than the 4.  used as prior estimates. We find that emissions from all major sectors except coal and wastewater are lower in the national inventories than our inversion results, with the largest underestimates for fugitive emissions from the oil sector. The total CONUS oil and gas emissions in our inversion are 4.6 and 9.9 Tg a -1 , respectively, 109% and 45% higher than the EPA (2016) inventory used here as prior estimate, and 177% and 65% higher than the most recent EPA (EPA, 2021) inventory for the 2010-2017 mean. The EPA 480 inventory reports an uncertainty of -24 to +29% for oil and -15 to +14% for natural gas emissions (EPA, 2021). Our estimates are also higher than those in Maasakkers et al. (2021), which are 3.6 and 8.0 Tg a -1 respectively for oil and gas emissions in 2010-2015. They are consistent with the Alvarez et al. (2018) estimates for total CONUS oil and gas emissions of 13 (11-15) Tg a -1 in 2015 based on field measurements within oil and gas basins, scaled up to derive a national value. 485 We mentioned previously that the lower estimates in Maasakkers et al. (2021) could reflect their use of GOSAT observations only, the difference in time frame, and their high prior estimate for wetlands, but another factor is their assumption of normal distributions for prior emission error standard deviations. We find from our inversion ensemble that assuming a log normal distribution instead (as in our base inversion)  (Table 1a) but still represents our best estimate.
Our inversion increases the oil emissions over Canada by more than a factor of two to 1.8 Tg a -1 compared to the ECCC inventory. The total posterior oil and gas emissions for Canada are 2.9 (1.6-3.3) Tg a -1 . This is in good agreement with a recent inversion study ( We further compared our oil/gas inversion results for CONUS, Canada, and Mexico to the TRACE bottom-up inventory aggregating data from individual assets up to the country level (Climate TRACE, 2021). This inventory uses lifecycle assessment emissions models for production, processing, refining, 510 and shipping (Gordon et al., 2015;Masnadi et al., 2018;Gordon and Reuland, 2021). The TRACE oil/gas total emission estimates for CONUS (9.6 Tg a -1 ), Canada (1.8 Tg a -1 ), and Mexico (0.5 Tg a -1 ) are similar to the prior estimates from EPA, ECCC, and INECC respectively (Table 1) and correspondingly lower than our best posterior estimates of 14.5 Tg a -1 for CONUS, 3.2 Tg a -1 for Canada, and 1.3 Tg a -1 for Mexico. The bottom-up oil and gas modeling in TRACE assesses routine methane emissions from normal 515 operations, assuming normal fugitive emissions. Recent flyover work, however, shows that methane emissions are highly intermittent (Cusworth, et. al., 2021) and this is not well captured in bottom-up estimates. over the oil/gas production regions of Alberta and Saskatchewan, again consistent with previous studies (Johnson et al., 2017;Baray et al., 2018;Chan et al., 2020). For Mexico the upward correction is concentrated in the onshore Sureste Basin which is the largest oil field in the country, but with a downward correction for offshore operations. This is consistent with aircraft and TROPOMI satellite observations, suggesting that methane from offshore oil platforms is piped onshore and inefficiently flared (Zavala-   Figure 9 shows the 2010-2017 time series of total anthropogenic methane emissions from CONUS, Canada, and Mexico, and the contributions from the dominant sectors (oil, gas, coal, livestock, and 550 landfills). We include no trend in the prior estimates so that the trends in Figure 9 are solely driven by the observations. Table 3 gives  (Figure 9). In contrast, the most recent EPA GHGI inventory reports a steady decreasing trend of -0.8 % a -1 in US anthropogenic methane emissions 565 over the 2010-2017 period mostly driven by coal (-5.4% a -1 ) and landfills (-1.6% a -1 ) (EPA, 2021). The decrease for gas is more pronounced in our inversion than in the EPA inventory (-0.4% a -1 ). The EPA inventory shows no significant trend for oil emissions. Figure 10 shows the spatial distributions of the linear regression fits to the 2010-2017 trends for the major 570 anthropogenic sectors, i.e., the equivalent linear trends over the period. We find that the oil increases are mostly driven by major basins in the south-central US including the Permian and Eagle Ford basins. The gas decreases are mostly driven by fields in the western US (Niobrara) and southeastern US (Haynesville).
Livestock emissions show variable regional patterns of increase and decrease that could reflect variations in animal populations. The increase in Iowa is consistent with a previous study of GOSAT trends by Sheng ). The inversion also suggests a decreasing trend in Mexican anthropogenic methane emissions by -3.3 (-3.4 --1.7) % a -1 , but this is mainly driven by a decrease from 2010 to 2011. We find very large relative decreases of oil emissions (-11.6 (-15.0 --3.5) % a -1 ) in particular for offshore Sureste, consistent with increasing utilization of associated gas (Zhang et al., 2019).

2010-2017 wetland methane emissions and trends
Our inversion shows strong ability to optimize wetland emissions over CONUS and Canada (averaging kernel sensitivity of 0.57). Wetland emissions in Mexico are much smaller and are not efficiently optimized by the inversion, as shown in Table 1c. Posterior wetland emissions are 8.4 (6.4-10.6) Tg a -1 for CONUS and 9.9 (7.8-12.0) Tg a -1 for Canada, compared to 7.5 and 12.0 Tg a -1 in the prior estimate 595 from the WetCHARTs v1.3.1 high-performance subset for North America (Ma et al., accepted). There are larger regional upward (southeast US) and downward (Upper Midwest) corrections even with this highperformance subset, as shown in Figure 11, pointing to major gaps in our understanding. Figures 11 and 12 show the 2010-2017 trends of wetland emissions for CONUS and Canada. We find a 600 significant increase of 2.6 (1.7-3.8) % a -1 in wetland methane emissions over CONUS in 2010-2017, in particular after 2014, and this is consistent with but higher than the WetCHARTs trend estimates of 1.3 % a -1 (not used in the inversion). The trends over CONUS are mostly driven by increases in the southeast US (Fig.11b). Fluctuations in emissions for temperate and boreal wetlands are mostly modulated by temperature, snow melt, precipitation, and drought events (Watts et al., 2014). We find a significant 605 correlation of 0.89 between the CONUS wetland emissions and annual precipitation in the CONUS wetland regions, and a strong 2010-2017 increase in the later that may drive the wetland trends (Fig.12c). Wetland emissions over Canada do not show significant trends in the inversion. The 2016 peak is consistent with WetCHARTs and may be explained by high precipitation (Fig.12d). range) Tg a -1 for CONUS, of which 36.9 (32.5-37.8) Tg a -1 is anthropogenic. This anthropogenic emission is 30% higher than the EPA inventory of 28.7 Tg a -1 used as prior estimate (EPA, 2016), and 42% higher than the 2010-2017 mean of 26.0 Tg a -1 in the most recent version of the EPA inventory (EPA, 2021). These upward corrections are largely attributed to the oil (4.6 Tg a -1 ) and gas (9.9 Tg a -1 ) sectors, which are respectively 177% and 65% higher than the EPA (2021) estimates. The upward corrections of the oil 630 and gas sectors are mainly in large basins of the south-central US. The inversion also shows upward corrections of livestock emissions to 10.6 Tg a -1 , 15% higher than the EPA estimate (9.2 Tg a -1 ), and of landfill emissions to 7.2 Tg a -1 , 24% higher than the EPA estimate (5.8 Tg a -1 ).
We estimate a mean 2010-2017 emission for Canada of 16.2 (13.5-17.4) Tg a -1 , of which 5.3 (3.6-5.7) Tg 635 a -1 is anthropogenic. This anthropogenic emission is 43% higher than the 3.7 Tg a -1 in the ECCC (2020) national inventory used as prior estimate. Most of this difference is due to oil emissions which we estimate at 1.8 Tg a -1 , more than twice the ECCC estimate, and mainly from production in Alberta and Saskatchewan.

640
We estimate a mean 2010-2017 emission for Mexico of 6.8 (5.4-6.9) Tg a -1 , of which 6.0 (4.7-6.1) Tg a -1 is anthropogenic. This anthropogenic emission is 20% higher than the 5.0 Tg a -1 in the INECC (2018) national inventory used as prior estimate. Again, most of the underestimate is due the oil sector and specifically to oil production in the Sureste onshore region. Offshore oil emissions are lower than the INECC estimate, suggesting that the associated gas is piped onshore and then vented, perhaps because of 645 inefficient flaring.
We find from the inversion that anthropogenic emissions in CONUS peaked in 2014 and had no net trend over the 2010-2017 period (0.1 (-0.1-0.3) % a -1 ), in contrast with the EPA inventory that reports a steady decreasing trend of -0.8 % a -1 over this period. The net trend in the inversion reflects compensating effects 650 from increases in emissions from the oil and landfill sectors, decreases from the gas sector, and flat emissions from the livestock and coal sectors. We find a decreasing trend in Canadian anthropogenic emissions of -2.3 (-2.5 --1.6) % a -1 over the 2010-2017 period, mainly driven by oil and gas production. We also find a decreasing trend in Mexican anthropogenic methane emissions (-3.3 (-3.4 --1.7) % a -1 ) over the 2010-2017 period, mostly driven by the oil sector and in particular by offshore operations.

655
Wetlands are the main natural source of methane in all three countries. Starting from the high-performance subset of the WetCHARTs inventory ensemble as prior estimate, our inversion yields mean wetland emission estimates for 2010-2017 of 8.4 (6.4-10.6) Tg a -1 for CONUS, 9.9 (7.8-12.0) Tg a -1 for Canada, and 0.6 (0.4-0.6) Tg a -1 for Mexico. Wetland emissions in CONUS show a significant increase of 2.6 (1.7-   O., and Wunch, D.: Estimating global and North American methane emissions with high spatial resolution using GOSAT satellite data, Atmos. Chem. Phys., 15, 7049-7069, http://doi.org/10.5194/acp-15-7049-2015, 2015. Turner, A. J., Jacob, D. J., Benmergui, J., Brandman, J., White, L., and Assessing the capability of different satellite observing configurations to resolve the distribution of methane emissions at kilometer scales, 915 Atmos. Chem. Phys., 18, 8265-8278, http://doi.org/10.5194/acp-18-8265-2018, 2018. Tyner, D. R., andJohnson, M. R          = 50% (log normal) for emissions, and = 10 ppb for boundary conditions; (2) the same as (1) except that γ= 0.5 for GOSAT observations; (3)-(6) the same as (1), except that for emissions uses the other 4 options in the Table; and (7) is the same as (1), except that = 5 ppb for boundary conditions. Similarly, the GOSAT + in situ (long-term) and GOSAT inversions have 7 ensemble members, respectively. The in situ and in situ (long-term) inversion have 6 ensemble members, respectively. This adds up to 33 inversion ensemble members. Sectoral attribution is done by two alternative 1010 methods (see text in Section 2.6), resulting in a total of 66 members. c Including only long-term surface and tower sites with observations for all years of the 2010-2017 record. d Adding the errors from individual sectors in quadrature following Maasakkers et al. (2021).         Figure 6. Posterior error correlation coefficients (r) between sectoral methane emissions in the contiguous US (CONUS), Canada, and Mexico, using the sector-aggregated error covariance matrix as described in Section 2.6. Error correlation coefficients indicate the ability of the inversion to separate emissions between sectors (0 = perfectly, ± 1 = not at all).     inversion so that the trend in the inversion results is solely from the atmospheric observations. The right panels show the annual precipitation over the wetland regions of CONUS and Canada, as determined by weighting precipitation amounts with the WetCHARTs wetland emission fluxes on their native 0.5 o x0.5 o grid. The gridded precipitation data are from the ERA-Interim re-analyses (https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era-interim/, 0.5°×0.5°). 1115