Apportionment and Inventory Optimization of Agriculture and Energy Sector Methane Emissions Using Multi-Month Trace Gas Measurements in Northern Colorado

we

series into contributions from the energy and agriculture sectors.In contrast to the more common static linear regression employed by other studies with hour-to-day length time series, we adopted a dynamic linear model (DLM) approach to capture the expected temporal variations in regression coefficients more accurately (Kille et al., 2019;Pollack et al., 2022;Yacovitch et al., 2014Yacovitch et al., , 2015)).A Bayesian inversion then used the DLM-derived energy and agriculture sector methane observations and an atmospheric transport model to optimize energy and agriculture methane inventory fluxes within the study's 850 km 2 area of sensitivity (for derivation of the sensitivity area, see Section 2.5).
Results from the energy sector are synthesized with previous regional measurements and historical data.Since the first study of WF energy emissions in 2008, barrel of oil equivalent (BOE) energy production within both the study's sensitivity area and the larger WF has increased several-fold while emissions have slightly decreased.This trend is the result of declining mean WF emissions factors (EFs), which, after dropping 2.9 ± 0.4 kg CH 4 / BOE (∼75%) between 2008 and 2017 have stagnated at 0.4 ± 0.2 kg CH 4 /BOE since.As a consequence, further WF production increases may yield increasing methane emissions.In contrast, inferred agricultural methane fluxes were 3.5× greater than inventory estimates.We demonstrate that this discrepancy arises partially from the spatial distribution of livestock which is not captured in the inventory model.Our work highlights that tracer gas measurements can separate emissions from different sectors even in complex emissions environments while avoiding sectoral misallocation, and reinforces the importance of further monitoring to refine inventory models as EFs change.

Materials and Methods
First, we discuss the collection of time series methane and tracer gas data, and subsequent sector apportionment using a DLM.Next, we give a brief description of the atmospheric transport model and sector-resolved emissions inventory used in this work.Finally, we describe the Bayesian inversion approach which generates the optimized posterior emissions inventories.

Observational Data Collection
Methane (CH 4 ), ethane (C 2 H 6 ), and water (H 2 O) concentrations were measured at the Platteville Atmospheric Observatory (PAO, {40.182, −104.725})from 1 November 2021 to 17 January 2022 with an open-path MIR-DCS instrument; ammonia (NH 3 ) was measured with a commercial CRDS.(While two instruments were used in this work, in the future all four gases could be measured using a single DCS instrument with adequate spectral coverage (Herman et al., 2021).) Figure 1c shows the dry air CH 4 , C 2 H 6 , and NH 3 mole fraction time series reported in ppm ([μmol/mol]).Subsequent analysis relies on periods when all three species were measured.The study's estimated sensitivity area (black dashed rectangles, Figures 1a and 1b) encompasses 850 km 2 around PAO and denotes the area within which measurements substantially constrain methane emissions.Further information on the DCS system and experimental setup at PAO are provided in Texts S1 and S2 in Supporting Information S1.

Dynamic Linear Model Tracer Gas Analysis
Energy and agriculture contributions in a methane time series can be extracted using correlations with ethane and ammonia (Kille et al., 2019).Generally this is achieved by fitting the methane data to a linear regression model comprised of energy sector methane ( ), a background term (β 0 ), and a Gaussian noise term (ϵ): This model is appropriate because the majority of methane emissions within the study's sensitivity area are from energy and agriculture.While landfills emit substantial volumes of methane, estimated landfill fluxes within this work's sensitivity area are <1% of predicted contributions from energy and agriculture (Text S3 in Supporting Information S1).
Fluctuations in the β 0 , β 1 , and β 2 regression coefficients are expected; the background methane concentration β 0 varies diurnally as the boundary layer height changes, and the two tracer gas coefficients, β 1 and β 2 , change as emissions from different sources are transported to PAO.Since a static linear regression cannot model all such variations without sub-dividing the ∼2-month time series into arbitrarily smaller segments, we instead perform the tracer gas analysis using a DLM (West & Harrison, 1997).Methane data are modeled with the observation equation, and the system equation, where t is an index representing data time steps.Tracer gas observations, along with a constant unity term which models the intercept, are represented by the regression vector Observations are assumed to be subject to Gaussian noise ν t with a mean of zero and a variance V t (defined here as the variance of the point-wise difference of the methane time series).The state vector θ t = (β 0,t ,β 1,t ,β 2,t ) evolves over time as a function of the θ t-1 state vector and the evolution variance vector W t .Because the variance is difficult to directly estimate and may not be time-invariant, DLMs are often solved using a discount factor δ instead as a proxy for the "memory" of the system over time (West & Harrison, 1997).The discount factor is defined as δ = P t /(W t + P t ), where P t is the prior variance corresponding to a state vector with zero stochastic change (W t = 0).In that limiting case, δ = 1 (irrespective of the actual value of P t ) and the DLM is identical to a static linear regression model.An optimal discount factor can be determined through minimizing the model's mean standard error, but in practice this minimization becomes expensive for large data sets such as ours.Instead, 100 DLM fits were performed over the full times series data with discount factors sampled from a random uniform distribution spanning [0.98, 0.999]; the mean values from the 100 DLM fits are used throughout.(Discount values below 0.98 lead to numerical instability; data where the fractional variance of either β 1 or β 2 was greater than 100% of the fit value are excluded in subsequent analysis.DLM-derived β 0 values were consistent with region-wide background methane concentrations; see Text S4 in Supporting Information S1).

Atmospheric Transport Modeling
Influence footprints in a 6° × 6° domain centered on PAO were calculated with the STILT-R atmospheric transport model and 3-km High Resolution Rapid Refresh (HRRR) meteorological data (Benjamin et al., 2016;Fasoli et al., 2018;Lin, 2003).Each influence footprint H(z r ,T r |z i ,T i ) (units of [ppm m 2 s/μmol CH 4 ]) connects sector-specific emissions throughout the spatial domain, at location z i and time T i , to observed sector-apportioned methane mixing ratios at PAO (z r ) at time T r .Footprints were calculated for each hour in an 8-week period of observations from November and December 2021.Each footprint is the sum of a 48-hr duration back trajectory of 100 particles originating from PAO, calculated at 0.1° × 0.1° resolution and hourly step size with hyper near field effects enabled.

Emissions Inventories
Energy and agriculture emissions are estimated using a 0.1° × 0.1° resolution sector-resolved methane inventory derived from the 2012 US EPA GHGI (Maasakkers et al., 2016).The energy sector, x Energy (units of [μmol/m 2 s]), is the sum of IPCC categories 1B2b (Natural Gas Production + Processing + Transmission + Distribution) and 1B2a (Petroleum); coal methane emissions are not considered (IPCC, 1996).The agriculture inventory, x Agri , is the sum of IPCC categories 4A (Enteric Fermentation) and 4B (Manure Management).

Bayesian Inversion
Sector-resolved methane time series, (y Energy , y Agri ), can be modeled as the product of time-independent methane inventories, (x Energy , x Agri ), and the combined set of time-varying influence footprints, H, plus an error term ϵ, Energy =  Energy +  Agri =  Agri +  Bayesian inverse modeling uses observational constraints (   Obs Energy∕Agri ) to generate maximum a posteriori probability (MAP) estimates,   Posterior Energy∕Agri , using the prior information provided by the inventories,   Prior Energy∕Agri (Cusworth et al., 2020).The observation vector   Obs Energy∕Agri are the hourly mean mixing ratios of energy and agriculture methane averaged from the 2-min time series.Following other studies, analysis is restricted to measurements within the hours of 11:00-16:00 local time when the simulation's boundary layer is well mixed and better captured by the meteorological models (Bianco et al., 2022;Fasoli et al., 2018;Kunik et al., 2019;McKain et al., 2015;Sargent et al., 2018).This restriction yielded a total of 238 valid data points for each observation vector.The H matrix contains the corresponding STILT footprint for each valid hour, where each footprint is restricted to a 5.8° × 5.8° domain centered on PAO at 0.1° resolution for a total of 3,422 state vector elements; footprints were flattened and stacked to yield the final H matrix with shape (238 × 3,422).Variances for the diagonal prior and observational error covariance matrices were estimated using a restricted maximum likelihood (RML) approach (Michalak, 2004;Michalak et al., 2005).Analysis of the averaging kernel sensitivity matrix indicates the posterior inventory is constrained by observations in an 850 km 2 area centered around PAO.This sensitivity area is highlighted with a dashed rectangular outline in Figures 1, 3, and 4. Further details on the inverse analysis are provided in Text S5 in Supporting Information S1.

Time-Resolved Sector Apportioned Methane
We first examine the DLM tracer gas results which provide key observational constraints for the Bayesian inversion.Three examples, shown in Figure 2, demonstrate how the DLM analysis captures the evolution of tracer gas coefficients and associated uncertainties as different sources are transported to PAO.During periods with a low tracer gas concentration or little variation in the tracer gas, uncertainty in the respective coefficient increases.Additionally, an increased correlation between methane and one tracer gas reduces the respective coefficient's uncertainty.
DLM-derived coefficients can provide insight into emission source characteristics.Figures 2d and 2e show kernel density estimates of the energy (β 1 ) and agriculture (β 2 ) tracer gas coefficients over the multi-month observation period.The β 1 coefficient has been observed to vary as natural gas is extracted, processed, and transported (Cardoso-Saldaña et al., 2019;Peischl et al., 2013).Ethane and methane mole fractions for natural gas samples collected after 2010 in the WF by the Colorado Oil and natural gas Conservation Commission (COGCC) provide a direct comparison to our estimates for β 1 (Figure 2d) (Colorado Oil and Gas Conservation Commission, 2022).These data are collected from a range of sample locations, including well casings (bradenheads, well tubing, and surface, intermediate, and production casings), produced gas, and separators and water tanks.The β 1 values determined from the PAO data overlap with the lower end of well casing and the higher end of separator and water tank data, while being most consistent with measurements of produced gas.
Similarly, β 2 is expected to vary as emissions from different livestock species can have substantially different ratios of methane and ammonia concentrations (Golston et al., 2020).Other sources of variation could include atmospheric chemistry effects such as deposition and reactivity (primarily for NH 3 ).We compare our β 2 results with two mobile measurement studies in Figure 2e.While extensive sampling of ammonia/methane ratios throughout Colorado are not available, studies in both the San Joaquin Valley of California and Northern Colorado overlap well with β 2 results obtained at PAO, indicating a consistent, if broad, distribution of β 2 values for agriculture across the western United States (Eilerman et al., 2016;Miller et al., 2015).
Significant variations in tracer gas coefficients observed in this analysis emphasizes the difficulty in determining a unique set of energy and agriculture coefficients, even for measurements conducted in a single location.Despite these complexities, the DLM approach successfully generates energy and agriculture sector-apportioned methane time series (   Obs Energy ,  Obs Agri ) which provides observational constraints for inventory optimization.

Methane Inventory Optimization
In order to provide emissions estimates with quantified uncertainties, an ensemble of 100 inverse analyses were performed for each sector, with each inversion using prior and observational error variances drawn from an optimal range determined using an RML approach.Inversion results were evaluated on several metrics, including the posterior's reduced chi-square (χ 2 ) value, and changes in the coefficient of determination (ΔR 2 ) and root mean square error (ΔRMSE) between observed (y Obs ) and predicted (y Prior and y Posterior ) mixing ratios (Kunik et al., 2019;Tarantola, 2005).Optimal inverse results were identified as having a χ 2 close to 1 (Text S5 in Supporting Information S1).Finally, we compare mean fluxes from x Prior and x Posterior within the 850 km 2 sensitivity area identified by the averaging kernel sensitivity matrix.All prior uncertainties quoted for comparison to posterior results are calculated following (Maasakkers et al., 2016).(Mean observed, prior, and posterior diurnal energy and agriculture methane mixing ratios are shown in Text S6 in Supporting Information S1).

Energy Sector
Inverse analysis results for the energy sector are shown in Figure 3. On average, the inversions achieved a χ 2 = 0.98, a RMSE reduction of 12%, and a 13% increase in R 2 .Posterior mean methane fluxes within the study's sensitivity area (78 ± 33 nmol CH 4 m −2 s −1 ) agree within uncertainty with the prior (100 ± 53 nmol CH 4 m −2 s −1 ), although posterior emissions were slightly reduced north-east of PAO (Figure 3c).A constant error in the DCS measurement of ethane cannot account for the comparable mean prior and posterior fluxes (Text S7 in Supporting Information S1).From these results, we next calculate an EF for the study's sensitivity area.EFs quantify the amount of methane emitted per volume of energy produced by a process, and are useful for generating and improving emissions inventories.Using the posterior energy flux and oil and gas production volumes within our study's 850 km 2 sensitivity area, we estimate a mean EF of 0.4 ± 0.2 kg CH 4 emitted per BOE produced (Skinner et al., 2022).For context, we also calculate historical EFs for the WF using data from multiple airplane mass balance, flask sample, and satellite inversion studies spanning 2008 to 2021 (Alvarez et al., 2018;Cusworth et al., 2022;Fried & Dickerson, 2023;Lu et al., 2023;Peischl et al., 2018;Pétron et al., 2012Pétron et al., , 2014;;Riddick et al., 2022;Shen et al., 2022).Finally, we estimate inventory EFs for the WF and our study's sensitivity area (Maasakkers et al., 2016).Combined, this synthesis (Figure 3d, data provided in Text S8 in Supporting Information S1) demonstrates a clear decline in mean WF EFs from 2010 to 2017, with post-2017 EFs remaining steady.Our study, though sensitive to a smaller central region of the overall field, yields EFs consistent with the larger WF.
From this analysis, we estimate how total energy emissions have changed over time in the WF by calculating the product of each EF and the Field-wide energy production (Figure 3e).While we note that applying this work's EF to the entire Field is an extrapolation, trends in both new well installations and energy production are closely mirrored in our study's sensitivity area and the larger WF (see Text S9 in Supporting Information S1).A declining trend in total energy sector methane emissions since 2010 emerges from this ensemble of independent measurements, although current emissions are comparable to 2008 levels.State and federal air quality regulations likely contributed to this decline by encouraging the adoption of emissions reduction practices and technologies.Changing well infrastructure may have also contributed (Text S9 in Supporting Information S1).Determining the relative importance of regulations and infrastructure requires further analysis.

Agriculture Sector
Results of the agriculture optimization are shown in Figure 4.The mean posterior had a χ 2 = 1.04, and a RMSE reduction of 22% and a R 2 increase of 41% compared to the prior.Within the sensitivity area, fluxes around PAO increased from a prior mean of 14 ± 16 nmol CH 4 m −2 s −1 to a posterior mean of 49 ± 22 nmol CH 4 m −2 s −1 (Figures 4a and 4b).This 3.5× ± 2.4× increase is surprising given that the total permitted livestock population around PAO has remained roughly constant since 2012 (National Agricultural Statistics Service, n.d.).While a threefold error in livestock EFs is possible, we instead investigated whether a spatial misallocation of emissions could explain the enhanced posterior flux in the sensitivity area.A comparison of the prior (Figure 4a) to registered concentrated animal feeding operation sites (CAFOs, Figure 1b) demonstrates that fluxes are not localized around CAFOs.This is a result of methodology: the agriculture inventory was generated by probabilistically distributing county-level livestock headcounts throughout each county using multiple livestock occurrence probability maps (Maasakkers et al., 2016).For some livestock, such as beef cattle which graze in pastures for parts of the year, this is a logical approach; however, poultry and dairy cattle are often on CAFOs throughout the animal's lifespan.
To determine if localizing inventory emissions at CAFOs improved agreement with observations, county-level emissions were redistributed to CAFO sites (Text S10 in Supporting Information S1) proportionate to the fraction of total county-level animal equivalent emissions units located at each CAFO (Golston et al., 2020).Total county level emissions were unchanged, reflecting our assumption that agricultural emissions have remained constant.The redistributed prior had a mean flux of 30 ± 34 nmol CH 4 m −2 s −1 within the sensitivity area, which increased in   Posterior Redist Agri to 43 ± 29 nmol CH 4 m −2 s −1 , comparable to the original   Posterior Agri results.The redistributed posterior had a similar χ 2 value (0.96) but a smaller RMSE reduction of 11% and a smaller R 2 increase of 25%, consistent with the smaller differences (Figure 4f) between the redistributed prior (Figure 4d) and posterior (Figure 4e), compared to those observed with the original inventory (Figure 4c).

Conclusions
We constrain energy and agriculture methane emissions in a ∼850 km 2 area in Northern Colorado by analyzing measurements of methane, ethane, and ammonia with a DLM and Bayesian inversion.Comparison with the 2012 gridded EPA inventory within the study area showed a small decrease in energy sector methane emissions which was consistent with a decrease in energy EFs from 2010 to 2017 observed by other studies.State and federal regulations and changing energy infrastructure likely contributed to this decline.While current energy emissions are lower than 2010 levels, they are comparable to 2008 observations, which indicates that further reductions are necessary to meet Colorado's greenhouse gas emissions reduction targets of 50% by 2030 and 90% by 2050 relative to 2005 levels.A significant increase in posterior agricultural methane emissions helped identify issues in the spatial distribution of agricultural fluxes.Redistributing emissions to CAFO sites improved agreement between the redistributed prior and posterior, although posterior agriculture emissions in the sensitivity area remain ∼50% higher than the redistributed prior.Refining the spatial distribution of emissions inventories is critical for regional scale studies using aircraft or satellite observations where multiple tracer gas observations are not present (Cusworth et al., 2021;Peischl et al., 2018).While conclusions from our single-sensor study can be improved with a distributed sensor network, it is noteworthy this approach can refine sector-resolved methane emission across areas comparable to the footprints of many methane observing satellites (Cusworth et al., 2021;Ware et al., 2019).by NIST through the Innovations in Measurement Science and Technology Maturation Acceleration Programs, and by the National Aeronautics and Space Administration Earth Science Technology Office's Instrument Incubator Program (Grant NNG20OB05A).The authors thank NOAA (Eric Williams) for access to the Platteville Atmospheric Observatory, Andy Neuman for loan of the ammonia CRDS instrument, and Dan Zimmerle, Dan Bon, and Chad DeVolin for CAFO information.The authors thank Newton Nguyen, Israel Lopez-Coto, Zachary Grey, and Subhomoy Ghosh for helpful suggestions and comments, as well as NIST reviewers Christopher Flood and Kristina Chang.Finally, the authors thank both anonymous reviewers for their constructive and knowledgeable critiques of this work.Official contribution of the National Institute of Standards and Technology; not subject to copyright in the United States.

Figure 1 .
Figure 1.Energy and agriculture sector methane sources are intermingled around the Platteville Atmospheric Observatory (PAO) measurement site.The study's sensitivity area is outlined in the black dashed rectangle centered on PAO, while county borders are outlined in black.(a) Thousands of wellheads (shown as a density map) extract oil and gas from the Wattenberg Field (WF, red outline).Locations of other down-stream components of the extraction process are not shown.This work's sensitivity area covers ∼19% of the WF.(b) Major agricultural developments called concentrated animal feeding operations (CAFOs, color coded by livestock and scaled to relative expected emissions magnitude), are widely distributed and spatially overlapped with energy infrastructure.(c) The full multi-month methane, ethane, ammonia, (expressed as dry mixing ratios) and water time series recorded at PAO.

Figure 2 .
Figure 2. Three methane plumes (a-c) illustrate how the DLM apportions methane into contributions from the energy and agriculture sectors.The tracer gas coefficients (dashed lines, left axis) are shown with uncertainties in gray shading.In addition, the top panel shows both the full modeled methane concentration (solid line, right axis) and the measured methane concentration (red dots, right axis).The second and third rows show the ethane and ammonia measurements (blue and green dots, right axis).Panel (d) compares the DLM-derived β 1 coefficients from the full time series with β 1 coefficients calculated from COGCC sample data.Panel (e) compares β 2 coefficients at PAO with other studies performed in Northern Colorado (Eilerman et al., 2016) and California's San Joaquin Valley(Miller et al., 2015).

Figure 3 .
Figure 3.The spatial distribution of energy sector methane emissions are optimized with a Bayesian inversion using the energy sector methane time series observed at PAO. (a-b) Prior Prior Energy and posterior Posterior Energy energy sector methane flux maps remain largely similar in both distribution and magnitude of emissions within the sensitivity region of the study area (black dashed rectangle).(c) The Posterior-Prior difference indicates a slight reduction in emissions north-east of PAO.(d) Comparison of this work's estimated mean emissions factor to other studies of the larger WF.A notable decrease in emissions factors is apparent from 2010 to 2017.(e) Estimated WF methane emissions are calculated using the EFs from (d) and WF energy production volumes (black line).Despite increasing production from the region, total emissions have slightly declined in the past decade.

Figure 4 .
Figure 4. Comparison of the original (top row) and re-distributed (bottom row) agriculture sector methane inventory and posterior results supports localizing emissions at CAFO sites.(a-b) Posterior Posterior Agri agriculture methane within the study area's sensitivity region (black dashed rectangle) are more localized around PAO than in the prior Prior Agri .(c) Difference between prior and posterior emissions are significant, with a several-fold increase in emissions to the north-west.The change in agriculture methane emissions between the redistributed prior (d) and posterior (e) are smaller than with the original inventory prior, as illustrated by the difference map (f).