Inter‐Annual Variability in Atmospheric Transport Complicates Estimation of US Methane Emissions Trends

US natural gas production increased by ∼43% between 2005 and 2015, but there is disagreement among existing studies on whether this growth led to increased methane emissions. We evaluate the likely contributions of atmospheric transport to an upward trend in atmospheric methane enhancements during 2007–2015, defined as the contribution of North American emissions to atmospheric observations across the US. We find that interannual variability (IAV) in transport yields an apparent upward trend in enhancements across much of the US during this time and can explain disagreements among existing studies over emissions trends. We further find that enhancements at satellite and in situ monitoring sites are 19% higher during El Niño than La Niña, possibly because air masses spend more time over North America on average during some years. The results show that accurate modeling of IAV in transport is a key prerequisite to quantifying emissions trends.

examine trends in atmospheric observations from a site in Oklahoma and from the Greenhouse Gases Observing Satellite (GOSAT). They estimate that US emissions increased by 2.5%-4.7% per annum between years 2010 and 2014, depending on the observations analyzed. Sheng et al. (2018), also using GOSAT, report a similar upward emissions trend of 2.5 ± 1.4% per annum between years 2010-2016. By contrast, Bruhwiler et al. (2017) point out that global inverse modeling studies show no upward emissions trend for the US and provide several possible explanations for this discrepancy, including the impacts of sampling variability, the estimated methane boundary condition, atmospheric transport, and/or the limited sensitivity of column-averaged methane observations to surface emissions. Additional studies focused on the US reach similar conclusions; Lan et al. (2019) report a trend in US emissions of 0.7 ± 0.3% per annum (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) using in situ aircraft observation sites, Maasakkers et al. (2021)  The purpose of this work is to help explain the disparate trends reported by recent studies that use atmospheric methane observations and examine the role of variability in atmospheric transport. In a previous study, Bruhwiler et al. (2017) simulate methane timeseries at two aircraft sites in the northeastern US and show substantial inter-annual variations in methane mixing ratios, even without a trend in emissions. The results suggest that atmospheric transport may have a large imprint on IAV in atmospheric methane observations. Indeed, many existing studies emphasize the salient role of transport in determining the atmospheric distribution of GHGs and of numerous other air pollutants-from hourly to annual time scales and site-level to global spatial scales (e.g., Barnes et al., 2016;Denning et al., 1999;Keppel-Aleks et al., 2011;Li, 2018;Li et al., 2010;Liu et al., 2015;Lu et al., 2019;Pal et al., 2020;Samaddar et al., 2021;Sweeney et al., 2015;Torres et al., 2019). In the current manuscript, we compare how emissions versus atmospheric transport likely impact AME-for both in situ and satellite observations across the US. We specifically explore whether this IAV in transport is likely to average out or cancel out from one region to the next. We then evaluate what meteorological factors may be driving this methane variability. Lastly, we explore connections between methane variability and large-scale climate patterns like the El Niño Southern Oscillation (ENSO) and North Atlantic Oscillation (NAO).

Atmospheric Modeling
We model atmospheric methane enhancements (AME) between years 2007 and 2015 at 8 tower measurement sites in the continental US that are part of the National Oceanic and Atmospheric (NOAA) Global Monitoring Laboratory (GML) Cooperative Air Sampling Network (Figure 1a; Andrews et al., 2014). Tall tower observations in the US greatly expanded in 2007, and the 8 tower sites included in this study have observations available during most or all years of the study period. Sites include Argyle, Maine (AMT); Erie, Colorado (BAO); Park Falls, Wisconsin (LEF); Billings, Oklahoma (SGP); Sutro Tower, San Francisco, California (STR); West Branch, Iowa (WBI), Walnut Grove, California (WGC), and Moody, Texas (WKT) (Andrews et al., 2014). We further model AME at 80,914 GOSAT sounding locations across the continental US (CONUS) between years 2009 and 2015. GOSAT sounding locations are specifically taken from the UoL Proxy XCH 4 Retrieval Version 9 (Parker et al., 2020), which is used in several recent studies of methane emissions (Maasakkers et al., 2021;Sheng et al., 2018).
We model AME at these locations using simulations from the Stochastic Time-Inverted Lagrangian Transport model (STILT) (e.g., Lin et al., 2003). STILT is a regional particle trajectory model; it tracks a large set of tracer particles (500 in this study), and the dispersion of those particles in the atmosphere is used to generate an influence footprint (in the units of ppb per unit of emissions). We model AME at each location and time by multiplying each of these footprints by a methane emission estimate (described below). Because STILT is regional in nature, it only estimates the methane increment or enhancement at the observation sites due to emissions in the model domain, in this case North America. A methane boundary condition or background can be added to this enhancement to estimate total methane mixing ratios at the observation sites, though variability in the background is not the focus of the present study.
The STILT simulations used here were specifically generated as part of the NOAA CarbonTracker-Lagrange project (e.g., Hu et al., 2019) and are driven by meteorology from the Weather Research and Forecast (WRF) 10.1029/2022GL100366 3 of 10 model (Skamarock et al., 2008). To date, WRF-STILT has been used for atmospheric transport in numerous existing regional methane and greenhouse gas modeling studies (e.g., Hegarty et al., 2013;Hu et al., 2019;Miller et al., 2013Miller et al., , 2014Nehrkorn et al., 2010). Each STILT simulation is run 10 days back in time, and footprints have a spatial resolution of 1° latitude by 1° longitude.
We further use several methane emissions estimates in the STILT simulations. Specifically, we use the US EPA gridded inventory across CONUS (available only for year 2012 at the time of writing; Maasakkers et al., 2016) and the Emission Database for Global Atmospheric Research (EDGAR) gridded methane emissions version 5 (Crippa et al., 2019) for anthropogenic emissions outside CONUS. We additionally use wetland methane emissions calculated using the model in Pickett-Heaps et al. (2011) (and as used in Miller et al., 2014    . Panel (b) further displays trends in modeled atmospheric methane enhancements (AME) (with uncertainties) for different scenarios at the in-situ observation sites (years 2008-2015). Lastly, panel (c) displays the difference in estimated trends in AME at Greenhouse Gases Observing Satellite observations sites (years 2009-2015) between modeling scenario 3 (interannual variability (IAV) in transport and no trend in emissions) and scenario 4 (no IAV in transport or emissions).

Modeling Scenarios and Trend Fitting
We analyze four modeling scenarios: one with trends in emissions and IAV in atmospheric transport (scenario 1), one with trends in emissions and without IAV in atmospheric transport (scenario 2), one without trends in emissions and with IAV in atmospheric transport (scenario 3), and one without trends in either emissions or IAV in atmospheric transport (scenario 4).
The emissions scenarios are generated based on the methane emissions estimates described in Section 2.1. For the scenario with no emissions trend, we use the monthly US EPA inventory estimate, monthly wetland emissions, and daily QFED emissions for year 2012 in all years of the study. For the scenarios with an emissions trend, we scale EPA oil and gas emissions in each state relative to monthly state-wide dry natural gas production data (US EIA, 2022) from years 2007-2015. This scaling results in an emissions trend of 2.3% per annum from the oil and gas sectors, changing from 8.3 Tg yr −1 2007, to 10 Tg yr −1 in 2015. By comparison, total US anthropogenic and natural emissions in all simulations are 47.9 Tg CH 4 for year 2012. Note that we do not add a trend to other methane source types because we are primarily interested in how a plausible trend in oil and gas emissions would manifest at the atmospheric observation sites, all else being constant. Some studies argue that US methane emissions trends are likely being driven by the oil and gas sectors (e.g., Sheng et al., 2018;Turner et al., 2016), and we therefore create a hypothetical emissions scenario that focuses on that sector.
We further generate meteorology scenarios that include IAV in atmospheric transport and scenarios that do not. For the former scenarios, we run WRF-STILT using standard protocols (Section 2.1). For the latter scenarios, we average footprints from different years to remove IAV in transport. Specifically, at each in-situ monitoring site, we average the footprints from each month of the year across all years of modeling simulations (see Text S2 in Supporting Information S1). In other words, we average the WRF-STILT footprints from all Januarys (across 2007-2015), across all Februarys, etc. This approach preserves seasonal variability in the footprints but removes IAV. For the GOSAT observations, we group the observations into 4° latitude by 4° longitude grid boxes across the United States. Within each box, we average the footprints from each month as described above.
We subsequently fit trend lines to the model estimates of AME for each scenario. We specifically fit trend lines using the procedures outlined in Lan et al. (2019) for in situ observations and Sheng et al. (2018) for GOSAT observations. We use line-fitting procedures from these studies to ensure that the results presented here are directly comparable to existing research (Text S2 in Supporting Information S1).

Atmospheric Transport Can Explain Conflicting Estimates of US Emissions Trends
Atmospheric transport yields an apparent upward trend in AME at all tall tower observation sites during the study period. Figure 1b displays the results of the four modeling scenarios at these sites, and the individual bars in the plot display the trend in AME estimated using a linear regression (Section 2). We find an upward trend in AME at all sites during the study period (2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015), irrespective of whether we include a trend in emissions (e.g., scenarios S1 and S3). Even without an emissions trend, the model outputs often display an upward trend between 2% and 5% per annum, ranging from 0.2% per annum (at Argyle) to 4.9% per annum (at Erie) (scenario 3). By contrast, when we remove IAV in meteorology, the upward trend in AME largely disappears, even when we include a trend in emissions. Specifically, the differences in estimated trends between scenarios 1 and 3 is between 0.1% and 1.1% per annum, comparable in magnitude to the modeled emissions trend (0.5% of total US emissions per annum).
These numbers are comparable in magnitude to the US methane emissions trend estimated by several atmospheric studies (e.g., Sheng et al., 2018;Turner et al., 2016). These studies attribute trends in observed atmospheric mixing ratios to emissions, while our results suggest that IAV in atmospheric transport can yield comparable numbers. Note that we estimate small negative trends in a few scenarios. In most of these cases, the standard error bars encompass zero. In two other instances (S2 at STR and WGC), the negative trend estimate occurs at sites that have a very large seasonal cycle in AME and have sustained data gaps; the combination makes trend estimation at these sites prone to error (Text S3 in Supporting Information S1).
More broadly, atmospheric transport has a large impact on inter-annual variability (IAV) in AME. For example, we calculate the maximum and minimum values in annually-averaged AME at each observation site in the tall tower network. For these calculations, we use anthropogenic emissions that do not contain any trend (scenario 3 described above), such that IAV in AME does not reflect variability in emissions. We find, on average, that IAV in AME at the observation sites is equal to 40% of the total average AME from North America (e.g., Figures S1 and S2 in Supporting Information S1). At some sites, particularly sites that are close to large agricultural or oil and gas emissions sources, this IAV is as high as 59% of average AME (e.g., at Eerie, Colorado, site BAO). In a previous study, Bruhwiler et al. (2017) argue that atmospheric transport can yield a trend of up to 1.5 ppb/yr in total column methane at two regular NOAA aircraft sites in the northeastern US. By contrast, we evaluate the impacts of transport across the entire US for both in situ and GOSAT satellite observations and find IAV as large as 25 ppb or over half the average annual AME at sites near large methane sources (e.g., BAO, SGP, WGC, and WKT; Figure S2 in Supporting Information S1).
We find similar results for simulated GOSAT methane enhancements. For example, Figure 1c shows the difference between estimated trends in scenarios 3 and 4 (simulations with and without IAV in transport), and most locations show a large positive difference in estimated trends (over 3% per annum), indicating that variability in atmospheric transport (in scenario 3) is likely driving much of that trend. Note that a small number of grid boxes yield unrealistic trend estimates (Figure 1c, e.g., coastal northern California and northern Vermont). These grid boxes contain a limited number of observations that are not evenly distributed across seasons and years, making trend estimation challenging.
Note that the apparent trends described above are specific to the study period (years 2008-2015). AME during other time periods exhibits similar IAV. However, the apparent trend in AME is downward as often as it is upward. A handful of in situ monitoring sites in the continental US began operations before 2007, and Figure S4 in Supporting Information S1 displays observational timeseries from three sites with long data records (AMT, LEF, and WKT). Even over two decades, the IAV in these timeseries is larger than any apparent long-term trends, implying that it would be challenging to separate emissions from atmospheric transport without using an atmospheric transport model.
We further conduct two sensitivity tests for in-situ observation sites in oil and gas producing regions (SGP and WKT)-one test explores the impact of the meteorological product used in STILT and one explores the impact of observation sampling time and frequency (Text S6 in Supporting Information S1). In simulations using both meteorology products, the impact of a trend in emissions is small relative to IAV in atmospheric transport, though the models do not always agree on the exact magnitude of AME in specific months. In the second test, we find that variations in sampling time have little impact on AME at one site (WKT) but do impact the results at another site (SGP); hence, we cannot rule out the role of observation sampling frequency and time on estimates of AME.

Links Between Methane Enhancements and Transport Variability
We hypothesize that IAV in AME could be caused by two broad mechanisms. First, IAV in wind direction could advect methane from source regions to the tower or GOSAT observation locations more frequently in some years than in other years, leading to IAV in AME. Second, regional-scale IAV in wind speeds or vertical mixing could alter the impact of emissions on downwind observations, irrespective of wind direction. For example, faster surface wind speeds and more vigorous vertical mixing could cause methane to be ventilated out of source regions more quickly. Similarly, these phenomena could also change the amount of time air masses spend over the North America. In other words, it could change the average number of days between when particles enter the North America model domain and when they reach an observation site.
We find evidence for the second explanation in the STILT model simulations. As evidence of this relationship, we see a close correlation between the overall magnitude of the STILT footprints and IAV in AME at both in situ and GOSAT observation sites (Figures 2a and 2b). The STILT model releases a set of imaginary particles from the observation site, and those particles travel backward in time following estimated wind fields. These particles indicate where air masses traveled before reaching the observation site. The magnitude of the STILT footprint indicates how long those particles spent in the surface mixed layer, and how concentrated or diluted the mixed layer was. The longer particles spend in the surface mixed layer, the larger the footprint. Figures 2a and 2b indicate that in years with higher AME, the average footprint magnitude is larger. In other words, the STILT particles have a longer residence time near the surface and have more intensive interactions with the surface during years with higher AME.
We also explore IAV in the time air masses spend over North America. Specifically, we calculate the average number of days STILT particles travel over North America before reaching the observation sites (Figures 2c  and 2d). At the end of the study period (years 2014 and 2015) when AME is larger at almost all sites, STILT particles spend more time in North America. By contrast, at the beginning of the study period when AME is generally smaller (years 2008 and 2009), STILT particles spend less time in North America. Remarkably, nearly all in situ sites show the same general pattern, indicating that IAV in atmospheric transport does not average out or cancel out from location to location across the US.
Unlike the footprint, Figures 2c and 2d do not indicate whether the STILT particles interacted with the surface where emissions sources are located or whether the particles traveled high in the atmosphere far above emissions sources. With that said, this panel yields parallel results as the footprint analysis.

Links Between Methane Enhancements and Large-Scale Climate Patterns
There appears to be a relationship between variability in AME and larger-scale climate patterns like the El Niño-Southern Oscillation (ENSO) and NAO. The first half of the study period is characterized by two La Niña events and persistent negative NAO anomalies (NWS, 2022a), while the end of study period features a strong El Niño event with consistent positive NAO anomalies (NCEP, 2022) ( Figure 3). Modeled AME is 19% larger (mean difference) at both tower and GOSAT observation sites during El Niño than during La Niña (Figures 3a-3c). Bruhwiler et al. (2017) also hypothesize that ENSO could lead to variability in methane mixing ratios-by changing the direction of air inflow into the North American continent and impacting the methane background or boundary condition. In theory, this variability could be accounted for by setting an appropriate background estimate. By contrast, we see variability in AME during different ENSO phases that is more pervasive and cannot be fully accounted for by setting appropriate background values.
This connection with ENSO may also be linked to variability in the amount of time air masses spend over North America (e.g., Figure 2d). Specifically, we find that STILT particles spend more time over North America before reaching the observation sites during El Niño years than during La Niña years-an average of 1.2 days or 15% longer at the tower sites. Relatedly, several studies report that wind speeds across the continental US generally decrease during El Niño and increase during La Niña (Enloe et al., 2004;Harper et al., 2007;St. George & Wolfe, 2009;Yu et al., 2015), except in California where the opposite appears be true (Berg et al., 2013;Mohammadi & Goudarzi, 2018). El Niño is also associated with easterly wind anomalies in the central and Eastern US (Ning & Bradley, 2015). Yu et al. (2015) further find that the primary mode of IAV in 80-m wind speeds across the United States is not only correlated ENSO but also with the NAO. Decreased wind speeds during El Niño and during positive NAO anomalies are consistent with decreased ventilation of local methane emissions and a longer residence times over the continental US. Unfortunately, few (if any) studies evaluate the impact of ENSO or NAO on other parameters that govern tracer transport across the US, parameters like planetary boundary layer heights.
Note that the primary time period of this study (2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) only covers a limited number of ENSO cycles-two El Niño events and three La Niña events (Figure 3). With that said, we see evidence for these relationships across longer time periods, though these relationships are more apparent in the southern US (e.g., at WKT) than in the northern US (e.g., at AMT and LEF). Figure S4 in Supporting Information S1 shows methane observed at WKT, AMT, and LEF minus two different methane boundary or background estimates. These sites have observations extending back to the 1990s or early 2000s. Values at WKT are 15-18ppb higher during El Niño than La Niña in this extended time series, depending on the boundary condition used (27%-73% of the overall AME). By contrast, values at AMT and LEF are roughly similar between ENSO cycles (within 2ppb or 15% of AME). These results are perhaps unsurprising; WKT is proximal to large emissions and may therefore be more sensitive to variations in ventilation of air from the region. LEF and AMT, by contrast, are remote sites, and enhancements in these locations depend on advection of air from distant source regions. Furthermore, meteorology at WKT may be more impacted by ENSO than AMT and LEF due to its southerly location.

Conclusions
Our results show that IAV in AME largely reflects atmospheric transport variability and that climate patterns like ENSO may be associated with large variations in AME, at least during the study period. Relatedly, we find that IAV in atmospheric transport does not necessarily average out or cancel out from one observation location to another: even observation sites at very different locations across the US show similar inter-annual variations in AME and in the amount of time STILT particles spend over North America. This variability is in addition to plausible transport-related IAV in the methane background (e.g., Bruhwiler et al., 2017).
We further argue that this transport variability can explain the disagreement over methane emissions trends among existing studies. Specifically, existing studies fall into two categories: studies that directly interpret trends in atmospheric observations (e.g., Lan et al., 2019;Sheng et al., 2018;Turner et al., 2016) and studies that estimate . Modeled atmospheric methane enhancements (AME) is higher during both El Niño and positive North Atlantic Oscillation cycles (a). Note that the right side of panel (a) displays the range for Greenhouse Gases Observing Satellite observations that have been averaged into 4° latitude by 4° longitude grid boxes, as in Figure 1c. Additionally, panels (b) and (c) show de-seasonalized modeled and observed AME at two tower sites in oil and gas production regions. emissions using inverse modeling, which accounts for variability in transport (e.g., Benmergui et al., 2015;Lu et al., 2022;Maasakkers et al., 2021). The former studies often report an upward emissions trend during a similar time period as the present study (2.5%-4.7% per annum) while the latter studies find little upward emissions trend (e.g., 0.1%-0.7% per annum).
These findings pose inherent challenges for the detection of greenhouse gas emissions trends, especially given the limited time span of many existing in situ and satellite observation records. With that said, the use of inverse modeling and reanalysis products that are true to observed IAV in meteorology provides the best path forward for quantifying emissions trends.