Technical note: Isolating methane emissions from animal feeding operations in an interfering location

Agricultural emissions, including those from concentrated animal feeding operations (CAFOs) for beef and dairy cattle, make up a large portion of the United States' total greenhouse gas (GHG) emissions. However, many CAFOs reside in areas where methane (CH₄) from oil and natural gas (ONG) complicates the quantification of CAFO emissions. Traditional approaches to quantify emissions in such regions often relied on inventory subtraction of other known sources. We compare the results of two approaches to attribute the CAFO CH₄ emission rate from the total CH₄ emission rate derived from an aircraft mass balance technique. These methods make use of the mixing ratio data of CH₄, ethane (C₂H₆), and ammonia (NH₃) that were collected simultaneously in-flight downwind of CAFOs in northeastern Colorado. The first approach, the subtraction method (SM), is similar to inventory subtraction, except the amount to be removed is derived from the observed C₂H₆ to CH₄ ratio rather than an inventory estimate. The results from this approach showed high uncertainty, primarily due to how error propagates through subtraction. Alternatively, multivariate regression (MVR) can be used to estimate CAFO CH₄ emissions using the NH₃ emission rate and an NH₃ to CH₄ ratio. These results showed significantly less uncertainty. We identified criteria to determine the best attribution method; these criteria can support attribution in other regions. The final emission estimates for the CAFOs presented here were 13 ± 3 g of CH₄ per head per hour and 13 ± 2 g of NH₃ per head per hour. These estimates are higher than the inventory of the US Environmental Protection Agency (EPA) and previous studies highlighting the need for more measurements of CH₄ and NH₃ emission rates.

1 Introduction

Livestock produce large amounts of greenhouse gases (GHGs) and reactive nitrogen species, including methane (CH₄) and ammonia (NH₃), through enteric fermentation and waste generation. Ruminant animals (e.g., cattle, buffalo, sheep, goats, and camels) constitute a significant source of CH₄ as their digestive systems break down coarse plant material through microbial fermentation in their rumen stomach (large frontal stomach) and subsequently release the produced gas (CH₄). From 1990 to 2019, CH₄ emissions from enteric fermentation grew 8.4 % in the United States, making agriculture the largest source of US anthropogenic CH₄ emissions in 2020 (EPA, 2022). Waste and manure management are also significant emission sources of CH₄ and NH₃. Together, enteric fermentation and manure management account for more than 30 % of US anthropogenic CH₄ emissions (Maasakkers et al., 2016; EPA, 2022); CH₄ is a crucial GHG due to its high global warming potential (Myhre et al., 2013; Moumen et al., 2016; Smith et al., 2021).

Large uncertainty remains around the magnitude of livestock emissions. Beef and dairy cattle are confined to feedlots where they are fed and kept in tight areas to process the animals efficiently. These feedlots are known as concentrated animal feeding operations (CAFOs). The compactness of CAFOs has been shown to create significant emissions (Golston et al., 2020; Hacker et al., 2016; Eilerman et al., 2016; Staebler et al., 2009). However, observations indicated large variability in these CH₄ emissions, which creates uncertainty in cumulative estimates of agricultural emissions (Golston et al., 2020). There are many different factors determining the amount of CH₄ released from enteric fermentation and manure management, and practices may vary from farm to farm (EPA, 2022; Maasakkers et al., 2016). However, interfering sources of CH₄ and NH₃, such as oil and natural gas (ONG), waste pools, and landfills, may also complicate measurements from individual CAFOs. Improving the methodology for isolating CAFO signals from other interfering sources will allow more accurate measurements and provide new information to constrain greenhouse gas emissions from the agriculture sector.

When CH₄ concentration data are used alone for the quantification of emissions in regions with multiple CH₄ sources, there is not enough information to distinguish the contributing sources. However, previous studies have attributed CH₄ emissions in complicated regions using a variety of methods. A few examples include attribution by subtracting inventory data (Caulton et al., 2014; Peischl et al., 2015, 2018), collecting ground-based isotope data to attribute CH₄ signals (Townsend-Small et al., 2016), using C₂H₆ as a tracer to subtract or attribute the ONG fraction (Mielke-Maday et al., 2019), and using multivariate regression (MVR) with independent gas tracers to attribute sectors (Kille et al., 2019; Pollack et al., 2022).

There are many considerations when determining which attribution method is most appropriate, and often attribution is done based on the data available rather than the ideal methodology. Attribution using the subtraction of a value determined by inventory or calculated ratio is best done for larger regions due to the variability of emissions from individual sources. For ONG, individual sites can release different amounts of CH₄ and ethane (C₂H₆) compared to other wells and compressors (Yacovitch et al., 2014; Zimmerle et al., 2022). When quantifying emissions in a small region (as in this work), an estimate of emissions from each type of ONG source or a local C₂H₆ : CH₄ ratio near the CAFO is required to accurately separate the CH₄ emissions into contributions from the CAFO and nearby ONG activities.

Another concern with using a tracer like C₂H₆ in isolation to estimate the contribution of ONG is that in complicated regions, individual CH₄ signals become mixed, making it possible that the ratio observed by the aircraft is not representative of the original ratio at the ground. This is the theory behind tracer release, for example, where a tracer gas is released at a known rate near a source of interest and used to back out the source emission rate (Roscioli et al., 2015). The added gas does not have to be introduced exactly at the source emission point (which may be unknown), provided sampling occurs far enough downwind where the species are well mixed. Typically, this means that the tracer gas must be released within 100 m of the source and measured > 500 m downwind (Roscioli et al., 2015). For airborne data, there are many situations where we would expect signals to be mixed, complicating the use of airborne ratio analysis. Townsend-Small et al. (2016) circumvented this by combining ground-based isotope ratios (which show distinct ratios for particular sources) with aircraft data (which showed only one ratio). However, because they used a single isotope ratio, there is large uncertainty in their results stemming from the single isotope ratio that must be attributed to multiple contribution signals (one equation, multiple unknowns). This is similar to the difficulties that Smith et al. (2015) encountered using C₂H₆ as a tracer in a complicated region with multiple C₂H₆ : CH₄ ratios.

On the other hand, an approach like MVR, which makes use of multiple tracer gases, requires sufficient data and is subject to its own sensitivities (Kille et al., 2019). MVR is best used when there are multiple tracer gases that can be treated as independent variables and one dependent variable (i.e., CH₄). The more independent variables that are included, the more data are needed to produce statistically significant results. These independent variables can be thought of as independent source terms. Although MVR is generally an effective method, it is important to note that it may not be appropriate in all situations. Specifically, using MVR with too many independent variables in a region where there are only a few sources may result in misleading or inaccurate results. However, using MVR with a limited number of independent variables that correspond to known sources can yield reliable results for those specific sources. Even in cases where there is a possibility of additional sources, fewer independent variables will produce more robust results and other sources can be treated in the extra term or background.

Here, we demonstrate a methodology to isolate and quantify emission rates for CAFOs in the northeastern Colorado Front Range (NCFR) using airborne measurements of CH₄, C₂H₆, and NH₃. We investigated two methods for CAFO CH₄ emission isolation: (1) a subtraction method (SM) using the C₂H₆ : CH₄ ratio and (2) an MVR method using CH₄, C₂H₆, and NH₃. This study focuses on the NCFR, where there is a high density of large CAFOs. Figure 1 shows a map of the NCFR with CAFOs for beef cattle, dairy cattle, chickens, sheep, and swine in the area. The area is dominated by CAFOs for beef cattle and dairy cattle. However, the NCFR contains a large mixture of CH₄ emissions due to the high production of ONG. In this region, prior estimates indicate that natural gas accounts for 38.5 % of the state-wide CH₄ emissions, while agriculture accounts for 22.3 % of the state-wide CH₄ emissions (Arnold et al., 2014). The NCFR is home to the Denver–Julesburg Basin (DJB), with over 52 000 ONG wells, and it has an abundance of compressors and processing plants, many of which are in close proximity to CAFOs (Higley and Cox, 2007).

Figure 1 Map of northeastern Colorado. CAFOs are colored by animal type and sized by maximum animal units. Animal units are equivalent to 1 live beef cattle such that 1 head of beef cattle = 0.7 dairy cattle = 2.5 swine = 10 sheep = 100 poultry (CDPHE, 2017). Note that the maximum animal units represent the maximum animal capacity of a given facility and not necessarily the actual number of animals present at that facility at the time of sampling. CAFO data as of 2017 are registered with the Colorado Department of Public Health and Environment (CDPHE, 2017).

2 Materials and methods

2.1 Data collection

This study was conducted in the NCFR near Greeley, CO; the analysis includes data collected over Weld, Morgan, Logan, Larimer, and Washington counties. There are many CAFOs within these five counties with an area-wide maximum capacity of > 1 million heads of cattle (United States Department of Agriculture [USDA], 2018). This study used the University of Wyoming King Air (UWKA), which is a relatively small research aircraft capable of flying at low altitudes (100 m a.g.l.) and slow flight speeds (95 m s⁻¹). The UWKA is a national aircraft research facility owned and operated by the University of Wyoming. Flights departed from and returned to the Laramie Airport in Laramie, Wyoming (KLAR). Data outside of the NCFR were removed to ensure regional data only. This was done by screening out everything north and west of 41^∘ latitude, −105.25^∘ longitude. Three flights were performed in November 2019, departing around 12:00 MST (Mountain Standard Time) and lasting 2–4 h.

The UWKA was instrumented to measure CH₄, NH₃, and C₂H₆ mixing ratios. A Picarro G2401-m flight-ready analyzer measured CH₄, carbon monoxide (CO), carbon dioxide (CO₂), and water vapor (H₂O) at 0.25 Hz through infrared cavity ring-down spectroscopy (Crosson, 2008). This Picarro and other Picarro models were tested previously and found to be stable and suitable for airborne field measurements (Richardson et al., 2012). Two separate Aerodyne commercial quantum–cascade tunable infrared laser direct absorption spectrometers (QC-TILDAS) measured NH₃ and C₂H₆. These instruments are described in detail in Yacovitch et al. (2014), Pollack et al. (2019), and Pollack et al. (2022). Measurements of NH₃ are collected at 10 Hz and averaged to 1 Hz for reporting and this analysis; C₂H₆ measurements are collected and reported at 1 Hz. All chemical data were adjusted for time lag between instruments and further averaged to 0.25 Hz for emission calculations. Instruments were calibrated on the ground before and after flights. Instrument zeros were routinely measured in flight by overblowing the instrument inlets with a bottled source of synthetic “zero” air. Other in situ measurements from the UWKA standard instrument package included pressure, temperature, three-dimensional winds, GPS position, aircraft altitude, and heading. The wind speed is measured at 25 Hz then averaged to 0.25 Hz. The precision on the reported wind speed is 0.14 m s⁻¹ with an expected wind direction precision of ∼ 5^∘ (Strauss et al., 2015).

2.2 Flight patterns

Flights were designed to identify the best flight patterns to simultaneously quantify CH₄ and NH₃ emission fluxes, CH₄ to NH₃ ratios, and NH₃ deposition downwind of CAFOs. Prior to each flight, forecast meteorology was used to identify the ideal CAFOs to sample based on prevailing wind direction, isolation from other CAFO plumes, and other logistical constraints (e.g., proximity to urban areas, towers, and airports). Once airborne, the pilot conducted a vertical profile to characterize the mixed boundary layer (MBL). Selected CAFOs were located by coordinates and, once close enough, by sight. Once a selected CAFO was identified by sight, the UWKA pilot would perform a visual safety inspection of the area and then fly up or down in a spiral pattern centered on the selected CAFO. The aircraft proceeded to circle the target CAFO at a low altitude to confirm in situ enhancements of CH₄ and NH₃ mixing ratios and the direction of the outflow plume. Flight altitudes ranged from 0.1 to 3 km a.g.l., and near-CAFO flight altitudes depended on safety constraints.

Ammonia (NH₃) deposition calculations are not the focus of this work; however, observations of deposition require multiple downwind observations. NH₃ : CH₄ ratios should be calculated near the source, as documented in Pollack et al. (2022), which used the same dataset as this work. We investigated both spiral (which can be completed quickly) and horizontal (which require much more flying time) transect flight patterns, shown in Fig. 2. The spiral patterns can be executed in quick succession, provide information as to the upwind background, and provide multiple downwind distances for analysis of NH₃ deposition. However, the spiral transects were found to be undesirable for quantifying emissions because the aircraft could not get far enough outside the plume to characterize background conditions. For example, the further-downwind spiral transects had high enhancements of CH₄ and NH₃ for the full width of the spiral, indicating that the UWKA did not leave the plume. The horizontal transects, while ideal for emission calculation, did not provide much information about the evolution of downwind NH₃ deposition. In order to sample the full plume efficiently and provide multiple downwind observations, racetrack patterns or boxes were later identified as a preferred approach for future sampling.

Figure 2 Observed CAFO flight path color-filled with CH₄ (ppbv). Animal operations are indicated by the different colored circles as in Fig. 1. Green dots represent ONG wells; data of ONG are as of 2015 (Colorado Department of Natural Resources Oil & Gas Conservation Commission, 2016).

2.3 Emission calculations

CAFOs suited for emission calculations were identified based on the following requirements: (1) enhancements of CH₄ and NH₃ above background conditions, (2) flight path includes multiple transects downwind at different altitudes within the MBL, and (3) enough data near the CAFO, but outside the CAFO plume, to characterize background mixing ratios. Only one CAFO sampled on 13 November 2019 satisfied the stated requirements to be suitable for emission quantification and is used in the remainder of this work. The flight on 13 November 2019 (denoted as F2 from hereon) occurred during a period with strong and steady winds (average wind speed of 8.4 m s⁻¹) from the north–northeast (average wind direction 32^∘). The MBL was well mixed with a top at 1200 ± 150 m a.g.l. (Fig. S2 in the Supplement). The flight track was carefully planned to target a specific feedlot. The aircraft performed spiral transects at distances ranging from 4–14 km downwind of the feedlot, and stacked horizontal transects were conducted 12 km downwind. The horizontal transects were extended to account for high NH₃ and CH₄ values observed during the flight. However, it was found that this extended area included more than one feedlot, as determined by NOAA HYSPLIT back trajectories performed on the two peaks of the plume (Fig. S1). The analysis showed that the entire width of the plume contained emissions from multiple feedlots. The combined maximum capacity of all the feedlots within the plume was estimated to be 173 800 maximum heads based on the 2017 data from the Colorado Department of Public Health and Environment (Colorado Department of Natural Resources Oil & Gas Conservation Commission, 2016). The sampling region during F2 is surrounded by many ONG wells and is located 16 km southeast of Greeley, CO, an urban area with a population > 100 000.

We used a mass balance approach to calculate emissions. Aircraft mass balance has been used to quantify emissions from a variety of source types (examples include Cambaliza et al., 2014; Caulton et al., 2014; Karion et al., 2015; Peischl et al., 2015, 2018). Briefly, the mass flow rate of a species through a crosswind plane downwind of the source is approximated by the integration of enhancement above a background concentration over the width and height of the plume. The emissions are derived using Eq. (1) shown here:

$$\begin{array}{}\text{(1)}& M\left(u\right)=\underset{\mathrm{0}}{\overset{Z_{\mathrm{MBL}}}{\int }}\underset{-x}{\overset{x}{\int }}(C_u-C_{\mathrm{b}})\times U_{⟂}\times \mathit{\rho }\left(z\right)\mathrm{d}x\mathrm{d}z.\end{array}$$

In Eq. (1), M represents the molar flux (mol s⁻¹) of a gas downwind of the source. To find the enhancement, the local background concentration, C_b (ppbv), is subtracted from the measured concentration, C_u. The ideal gas law is used to calculate the air density (ρ) at every data point, using the universal gas constant at ∼ 8.31 J mol⁻¹ K⁻¹. The values ±x are the horizontal limits of the plume width from the center point, and Z_MBL is the top of the MBL. The full plume's width was identified by consistent NH₃ values above 5 ppbv. The plume's height is constrained to be from ground level to the top of the MBL. A vertical profile near the source was used to identify the MBL level using H₂O vapor and calculated potential temperature (Fig. S2). We used a constant value to represent the regional background. The background region was determined from the edges of the horizontal transects for the observed CAFO and used to calculate average values for CH₄, NH₃, and C₂H₆.

We transformed all observations surrounding the target CAFO sampled during F2 (Fig. 2) onto a polar coordinate system (r,θ) using the center of the target CAFO as the origin, following the process described in Nathan et al. (2015). The location of the data point on the polar coordinate system (θ) is perpendicular to the theoretical flux surface. Thus, U_⟂ (m s⁻¹ ) is the corrected perpendicular wind found by taking the cosine of the location of the data point on the polar coordinate system (θ), subtracted by wind direction (ϕ) and multiplied by the wind magnitude (V):

$$\begin{array}{}\text{(2)}& U_{⟂}=\cos (\mathit{\theta }-\mathit{\varphi })\times V.\end{array}$$

2.4 Uncertainty analysis

We conducted an uncertainty analysis using a Monte Carlo approach, creating a pseudo-distribution of the data using observed means and standard deviations and recalculating the emission rates. Due to the nature of the plume, not all variables may be represented using a pseudo-distribution. We define the final emission uncertainty as the magnitude of the change in the emission from a combined function of the pseudo-distributions of the following five parameters: background value, perpendicular wind speed, density, MBL depth, and instrumental uncertainty. The Monte Carlo recalculations were first done for individual parameters, then as a combined pseudo-distribution of perpendicular winds, density, background values, and MBL depth to calculate the final uncertainty. The uncertainty for the attribution methods followed the same approach but is addressed separately.

For the uncertainty analysis, we used a Gaussian distribution based on the mean and standard deviation of the original background value to select a new background value randomly. We found the density to have a pattern of increased values to one side of the horizontal transects. Therefore, to account for the pattern, we created a Gaussian distribution using the standard deviation found in each transect and added that onto a moving mean of 1 min. The pseudo-distribution of the MBL height was formed using a uniform randomly selected value of 2600 ± 150 m a.m.s.l. Throughout the day, the MBL height changed by an average of ∼ 150 m. The perpendicular winds consist of two variables: wind speed and wind direction. For each separate transect, we created a pseudo-Gaussian distribution of wind speed and wind direction based on the transects' mean and standard deviations for each variable. These Gaussian distributions of wind speed and direction were recalculated through Eq. (2) to produce the new perpendicular wind speeds. The emission recalculations were done 1000 times using all 4 parameters for CH₄, C₂H₆, and NH₃. The uncertainty of final emission estimates were calculated by using the 95 % confidence interval (CI) from the Monte Carlo approach divided by the average emission.

When considering the effect of instrumental uncertainty, we only examine the factors that actually impact the measurements used in the mass balance equation: the delta value (C_u-C_b). There are two factors that impact this value: (1) the intercept (bias) and (2) the slope (calibration factor) of the instrument. The bias, which can be affected by drift in the zero reading, is quantified for NH₃ and C₂H₆. To further clarify why these are the only instrumental uncertainty factors we evaluate and why total uncertainty does not affect our results, we consider the following hypothetical situation: two otherwise identical instruments are used to quantify the same enhancement. One instrument reads a background value of 1×10⁶ ± 0.1 and peak value of 1.000001×10⁶. The second reads a background value of 1±0.1 and an enhancement of 2. The resulting delta values have the same absolute value and uncertainty because of how error propagates through subtraction regardless of which instrument was accurate. We assume that the variability in the bias is dwarfed by real variability in the background, or if the bias is actually large, it similarly affects the background reading. Thus, we expect our analysis of background variability to be the appropriate metric to account for bias.

The uncertainty of the calibration factor is also possible to analyze and include. As far as we are aware, this source of uncertainty does not appear to be routinely reported or incorporated into uncertainty analysis in the mass balance literature. We assumed that the applied calibration factor can vary randomly and applied a randomly picked factor to the delta value to create a pseudo-distribution of possible enhancements. Note that the accuracy of the sensor does not matter. If the uncertainty of the delta value is calculated in the fashion described, there would be equivalent uncertainty of an accurate or inaccurate reading provided the variability of the background and calibration factor was the same. The variability of the calibration factor was <0.1 % for CH₄ and 9 % for C₂H₆. For NH₃, Pollack et al. (2019) found the variability to be 2 %.

The effects of instrumental precision, which affect both background values and enhancements, are neglected because we average a large area to calculate the background and bin average 5 s of data for the enhancements. Averaging has the effect of decreasing the random error. The 1 Hz in-flight precision is already low at 1 ppbv for CH₄ (Picarro, 2023), 200 pptv for C₂H₆, and 60 pptv for NH₃ (Pollack et al., 2019). Comparatively, the observed variability in the background values dwarfs the error from precision at 33 ppbv, 3 ppbv, and 4 pbbv for CH₄, C₂H₆, and NH₃, respectively. Thus, we would expect the error introduced from precision to be negligible.

2.5 CH₄ attribution

The two methods described here were applied to either all the data within the MBL during the flight (abbreviated F2) or only the downwind transects used for emission calculation (abbreviated Transect). This provides a total of four scenarios that will be analyzed. Ratios are reported in percentage (ratio ×100 %) across all methods.

2.5.1 Subtraction method

The first attribution approach, referred to as the subtraction method (SM), removes the CH₄ emissions related to ONG and attributes the excess CH₄ emission to the CAFOs. Rather than using inventory estimates for interfering sources, we calculated observation-based C₂H₆ : CH₄ ratios. The F2 ratio (11 ± 0.02 %) is from Pollack et al. (2022). Briefly, NH₃ values >5 ppbv were used to screen out data points associated with CAFOs. The remaining data points were assumed to be associated with ONG sources, which is consistent with observations in this region that CAFOs and ONG are the dominant sources (Kille et al., 2019). The slope was calculated using least squares orthogonal distance regressions (ODRs) and is equivalent to the C₂H₆ : CH₄ enhancement ratio. Such a ratio can also be calculated from the CAFO transects as there was a region of CH₄ and C₂H₆ signals that did not include NH₃ interference located on the northwest end of the transect. The resulting ratio for the Transect data following the identical calculation is 14.7 ± 0.6 %. The ONG CH₄ emissions were removed through Eq. (3):

$$\begin{array}{}\text{(3)}& M_{{{\mathrm{CH}}_{\mathrm{4}}}_{\mathrm{ag}}}=M_{{\mathrm{CH}}_{\mathrm{4}}}-{\mathrm{CH}}_{\mathrm{4}}:{\mathrm{C}}_{\mathrm{2}}{\mathrm{H}}_{\mathrm{6}}\times M_{{\mathrm{C}}_{\mathrm{2}}{\mathrm{H}}_{\mathrm{6}}}.\end{array}$$

The inverse of the C₂H₆ : CH₄ ratio is used with the emission rate of C₂H₆ (mol s⁻¹) to retrieve the portion of CH₄ associated with ONG; that value is subtracted from the total CH₄ emission rate to result in the CH₄ emission associated with the CAFOs. Finally, the molar emission rate (mol s⁻¹) is multiplied by the molar mass to return the emission estimate in grams per second (g s⁻¹).

2.5.2 Multivariate regression

An alternative to the SM is to directly calculate the CAFO CH₄ emissions from an NH₃ : CH₄ ratio and the NH₃ emission rate using MVR. Because of the widespread ONG activity in this region, it is not always possible to have clear regions with which to calculate a ratio using a traditional regression approach (as in the SM). Indeed, there is also concern that the ratios calculated in the SM may not be accurate due to the influence of diffuse CAFO or agricultural signals. Figure 3 presents C₂H₆ and NH₃ mixing ratios plotted against CH₄. In both plots, there are elevated regions of the other species (i.e., regions of elevated C₂H₆ in the NH₃ versus CH₄ plot). Instead, MVR using CH₄, C₂H₆, and NH₃ data can be used to calculate NH₃ : CH₄ and C₂H₆ : CH₄ ratios, as described in Pollack et al. (2022). Briefly, C₂H₆ and NH₃ are assumed to be independent tracers (associated with ONG and CAFOs, respectively), and CH₄ is the dependent variable, as shown in Eq. (4).

$$\begin{array}{}\text{(4)}& {\mathrm{CH}}_{\mathrm{4}}=a+b\times {\mathrm{NH}}_{\mathrm{3}}+c\times {\mathrm{C}}_{\mathrm{2}}{\mathrm{H}}_{\mathrm{6}}\end{array}$$

In this equation, a is the background CH₄ mixing ratio, b is the inverse effective NH₃ : CH₄ ratio, and c is the inverse SM C₂H₆ : CH₄ ratio. Unlike Kille et al. (2019) and Pollack et al. (2022), we did not subtract a background mixing ratio from the observed CH₄, NH₃, or C₂H₆ mixing ratios, thus the a variable actually represents the local background, and we can compare its value to the observed background that we identified at the edges of the transect. Kille et al. (2019) performed sensitivity analyses on their MVR results. Following the guidance of Kille et al. (2019), we only use the MVR results when all three variables are positive, R²>0.5, and all variables are statistically different from 0. We also tested scenarios with background-subtracted CH₄, NH₃, and C₂H₆ to compare with the approaches used in Kille et al. (2019) and Pollack et al. (2022). We performed this analysis on the entire F2 dataset versus the Transect only data. The results of these sensitivity analyses are reported in Tables S1–S2 in the Supplement. The choice of background made no difference to the ratios and only affected value a. This is consistent with Kille et al. (2019) and Pollack et al. (2022). Generally, scenarios where NH₃ and C₂H₆ were not background-subtracted produced a values more consistent with the observed CH₄ background and were used for the remainder of the analysis. Slight differences in the MVR results were observed from Pollack et al. (2022). This is attributed to the differences in the area used for MVR between these studies. Pollack et al. (2022) isolated specific source regions for MVR analyses while this study uses all of the data in the study region.

Figure 3 Scatter plots of (a–b) C₂H₆ versus CH₄ colored by the mixing ratio of NH₃ and (c–d) NH₃ versus CH₄ colored by the C₂H₆ mixing ratio sampled during the horizontal transects on 13 November 2019.

3 Results and discussions

3.1 CAFO emissions

Figure 4 shows curtain plots of the measured CH₄, NH₃, and C₂H₆ mixing ratios from the horizontal transects. The curtain plots display the data from the flight paths in color-filled boxes representing the area used for the calculation in Eq. (1). In general, elevated NH₃ coincides with elevated CH₄, indicating the presence of an agricultural plume. However, a high C₂H₆ enhancement appears to be embedded in the agricultural plume, verifying mixed sources of CH₄ in this region. To isolate the agricultural plume and minimize the influence of other sources of CH₄, we created a mask around the NH₃ signal and screened out points that were <5 ppbv of NH₃ (Fig. 4e). Further, the plume was limited to regions with contiguous NH₃ values above 5 ppbv. This removed some isolated areas of NH₃ enhancement at the lowest altitude. Because the lower transects showed a consistent drop in NH₃, we extrapolated that boundary upward to the highest transect as indicated by the vertical dashed line on the left side of Fig. 4. This removes some signal that is primarily downwind of Greeley, CO. The horizontal dotted line represents the boundary layer, and points above this layer were excluded. This threshold was identified by Pollack et al. (2022) and was calculated as 2 times the maximum NH₃ mixing ratio observed outside of the plumes of individual CAFOs and within the MBL. The resulting mask was used to limit the area integrated for emission estimates according to Eq. (1). The calculated emission rates are 3330, 609, and 542 g s⁻¹ for CH₄, NH₃, and C₂H₆, respectively.

Figure 4 Vertical curtain plots of horizontal transects for (a) CH₄, (b) NH₃, and (c) C₂H₆ ∼12 km downwind of the observed CAFO. For reference, the surface is located at ∼ 1400 m.s.l. Vertical curtain plots of (d) CH₄, (e) NH₃, and (f) C₂H₆, where NH₃>5 ppbv downwind of the observed CAFO. Points above the boundary layer were removed for calculation. Additionally, points on the left side of the vertical dashed line, above 5 NH₃ ppbv, were removed due to their likely origin from Greeley, CO.

We present an in-depth look at the uncertainty contributed by each parameter for CH₄ emissions only. Perpendicular winds, density, and MBL depth are consistent in the calculations of NH₃ and C₂H₆ emissions; thus, we only discuss total uncertainty for these. Total uncertainty distribution plots are shown in Fig. S3.

Variability in density had little impact (<1 %) on the final CH₄ emission rate. The winds during the plume had a standard deviation of ±1.3, ±1.9, and ±1.8 m s⁻¹ for lowest to highest transect altitude, respectively, with small changes in the wind direction with increased altitude. The wind speed was recalculated into the perpendicular wind and showed minimal change to the final CH₄ emissions by ∼ 4 %. The location of the background in this study had a relatively consistent CH₄ mixing ratio (mean = 1933 ppbv, standard deviation = 1 ppbv). The background is similar to the regional background (1990 ppbv) identified by Pollack et al. (2022). Background variation affects the final CH₄ emissions by 5 %. NH₃ had higher variations in background values leading to increased uncertainty overall. MBL height was the largest driver of uncertainty and was associated with an 8 % change in CH₄ emissions. The significant uncertainty due to MBL depth is expected as changes in MBL depth would result in the interpolation of a different area without a response in concentration measurements. Thus, if the MBL depth were to increase, there would likely be a corresponding decrease in concentrations. The absence of data between the lowest altitude transects and the surface may affect the accuracy of the results, but the associated uncertainty cannot be quantified because we do not have data at the surface. Instrumental uncertainty had little impact, changing the emission rates by <1 % for all: CH₄ (0.3 %), NH₃ (0.6 %), and C₂H₆ (0.5 %).

Cambaliza et al. (2014) carried out a detailed analysis on uncertainty from aircraft mass balance calculations. They found that MBL depth and background mixing ratios may have uncertainties up to 19 % and 21 %, respectively, on the overall emission rate estimates. Other studies have confirmed that uncertainty associated with the MBL depth makes a large contribution to overall uncertainty for these types of calculations (Karion et al., 2013; Peischl et al., 2015). The winds in the mass balance are assumed to be consistent from the release of the emissions to the location of the measurement. Therefore, the natural variability of winds may contribute to the uncertainty more with less consistent winds (including direction and speed) and may have a large effect on the final uncertainty similar to MBL and background value (Karion et al., 2013; Peischl et al., 2015). During F2, winds were consistent and made a smaller contribution to the overall uncertainty on this particular day than noted by these other studies.

The total uncertainty for CH₄ emissions calculated from horizontal transects was ±10 %. For the C₂H₆ and NH₃ emission estimates, the total uncertainty was ±14 % and ±17 %, respectively. Other factors that may influence the accuracy of the emissions include smaller-scale variations in the mixing ratios and regions of the plume that the flights did not sample. As it is not possible to sample the entire MBL from top to bottom, the vertical spacing between the horizontal transects may result in errors. Errors associated with interpolation were not explored here but are expected to be small (Cambaliza et al., 2014). Further, times when the bulk of the mixing ratio is near the ground are seen to have the most uncertainty (Gordon et al., 2015). The elevated mixing ratio values of CH₄, NH₃, and C₂H₆ were seen throughout the MBL. The uncertainty estimates are specific to this flight. Different meteorological conditions can lead to varying uncertainty estimates. These uncertainties are specific to each flight and can be influenced by factors such as smaller-scale variations in mixing ratios, unexplored errors associated with interpolation, and the proximity of the bulk mixing ratio to the ground. In particular, MBL uncertainty or growth may significantly affect the magnitude of the uncertainty. The uncertainty analysis methodology presented here can be applied to any future mass balance emission estimates.

3.2 CAFO CH₄ attribution: comparison of methods

Results from the SM and MVR CAFO CH₄ attribution and associated uncertainty are presented in Table 1. Reported uncertainties in Table 1 represent results from combining the ratio uncertainty to the rest of the pseudo-distributions as described in Sect. 2.4. Both the magnitude of the attributed CAFO CH₄ emission and its uncertainty vary between approaches. These methods represent typical approaches that can be undertaken to isolate the emissions for a given facility. In order to identify the optimal method, we identify four criteria: (1) the relationships should predict the total CH₄ well; in other words, traditional goodness-of-fit (GOF) values should be optimal (the average residual of the fit should be near 0 and R² values should be high); (2) the relationships for NH₃ : CH₄ and C₂H₆ : CH₄ should be consistent with observations; (3) the method should err on the side of being conservative, meaning it should be more likely to under-attribute than over-attribute the CAFO emissions; and (4) uncertainty of the result should be low. As currently presented, these criteria are qualitative because each analysis is unique. We will discuss the implications and further refinements of these criteria later on. These criteria are presented in order of their importance to the final recommendation. This approach ensures that low uncertainty is not the primary deciding factor.

Table 1 CH₄ Attribution sensitivity.

^a SE = standard error. ^b Pollack et al. (2022). ^c Criteria are defined in Sect. 3.2.

To investigate criterion 1, we calculated total fits, residuals, and R² values for all four methods and assessed how well the NH₃ : CH₄ and C₂H₆ : CH₄ ratios represented the Transect data. The central assumption is that regardless of the attribution method, all methods should predict the observed CH₄ values well as evaluated by traditional metrics of goodness of fit (e.g., R², residuals). MVR directly produces both ratios such that predicted CH₄ can be calculated from the NH₃ and C₂H₆ time series. To construct a CH₄ prediction for the SM, which does not require calculation of the NH₃ : CH₄ ratios, we inferred the effective NH₃ : CH₄ ratio from the ratio of the calculated NH₃ emissions to the CH₄ attributed to the CAFO by the SM. The resulting ratio is then used in Eq. (4) to produce an effective MVR prediction. Table S3 shows the full range of variables for each scenario. GOF statistics are calculated for the Transect data only.

The results of this comparison are shown in Fig. 5, with the prediction limited to the Transect data, because these are the data that are actually used to calculate the CH₄ emission rate. The Transect MVR approach was the method with the best GOF statistics including an average residual closest to 0 and the highest R² value; it is also the only approach that passes criterion 1. In some ways, this is unsurprising because fits are calculated by minimizing the sum of the residual squared, and this is the only method that directly fits the Transect data (Skoog et al., 2004). Still, this analysis does provide the means to compare the other methods and can be useful when site-specific MVR is not possible. The background values were consistent in all methods and consistent with observations. The three other methods all overpredicted CH₄ and had average negative residuals and lower R² values. The overprediction suggests that these methods produce an incorrect relationship between CH₄ and one (or both) of the tracer gases, but it is not possible to identify which relationship is incorrect. The method that best simulates the data after the Transect MVR is the F2 MVR.

Figure 5 (a) Residual CH₄ (observation–prediction) versus predicted CH₄ (ppbv). The ratios corresponding to the transect MVR results are shown in red, the F2 MVR results are shown in green, the transect SM results are shown in blue, and the F2 SM results are shown in cyan. The black line in both plots represents the zero line. (b) Violin plots of residuals for the four different methods: the black lines are the calculated mean residuals, the red lines are the median residuals, and the closed black dots are the sum of the residuals squared (normalized by the minimum value ×10). The violin plots are ordered by increasing absolute average residual.

The effective NH₃ : CH₄ ratios produced for this method can also be compared directly in support of criterion 2. Criterion 2 is difficult to evaluate because the actual ratios would ideally be known. We can compare these to literature values, compare values between scenarios, and look at the underlying assumptions of the calculations to evaluate this criterion. The C₂H₆ : CH₄ ratios calculated for the SM and MVR varied by 40 %, which substantially affects how much CH₄ is attributed to ONG. However, the range of ONG ratios is broadly consistent with other observations in this region (Kille et al., 2017, 2019; Yacovitch et al., 2017; Peischl et al., 2018).

As a plume moves downwind from its source, it is expected that the plume will disperse over the MBL. NH₃ has a short lifetime (∼ hours) and may be removed through dry deposition and transition into the particle phase, even at relatively short distances downwind from sources (Staebler et al., 2009; Miller et al., 2015). Thus, there may be lower concentrations downwind of gaseous NH₃ than expected from dilution alone. The range of reported NH₃ : CH₄ ratios in the literature is quite variable and includes observations from 1 %–200 % (Eilerman et al., 2016; Golston et al., 2020). We used Pollack et al. (2022) as our primary reference for comparison. Comparison to this study revealed that the transects located ∼ 12 km downwind of the observed target CAFO do not have the same NH₃ : CH₄ ratio as transects collected closer to the CAFO. There is also a wide variety of ratios reported in Pollack et al. (2022) ranging from 80 % to 270 % across different facilities which illustrate why the MVR results from the full F2 data may not be appropriate to represent an individual site. In general, the SM produced ratios that were lower than MVR and observations near the source (Pollack et al., 2022). The SM Transect ratio is particularly low (45 %), about half the value of any other ratio produced from the other attribution methods or results from Pollack et al. (2022). This indicates that this approach actually fails criterion 2.

On a theoretical basis, the MVR ratios calculated from the F2 data include locations near other CAFOs and locations near ONG sources. Therefore, the F2 MVR results may not accurately represent the area around this particular plume. The SM using F2 data has similar challenges. The SM assumes that C₂H₆ : CH₄ is constant throughout the area and that there are no other sources of CH₄ other than the CAFO and ONG. However, we know there may be other trace amounts of CH₄ from vehicles, waste areas, and other small sources. Generally, we would expect data nearest the source to produce the most accurate ratios and pass criterion 2. This would generally exclude the F2 results and for this reason, we assume that these results fail criterion 2.

For criterion 3, the two methods are different in approach and it is quite easy to identify which approach is more conservative. The SM attributes anything outside of ONG to the CAFOs plume; thus, there is no unattributed CH₄ for these methods. The SM is, therefore, expected to be less conservative and an upper limit of the CAFOs CH₄ emissions due to the possibility of other sources in the region.

In the MVR approach, the R² value of the fit gives an estimate how much of the CH₄ signal is explained by the chosen tracers. The MVR fits had R² values of 0.72 for F2 and 0.74 for Transects, which leaves an amount of CH₄ that is not well correlated to NH₃ or C₂H₆. In this scenario, this excess CH₄ is left unattributed in the MVR approach (56 % for F2 and 24 % for Transect). A visual representation of the unattributed portion of the CH₄ signal for the transect MVR is presented in Fig. S4. The unattributed CH₄ is broadly distributed and cannot be attributed to any particular source. The MVR approach avoids over-attributing the CH₄ signal to the CAFOs.

In general, we expect the Transect MVR approach, which uses the data closest to the source, to provide the best estimate of these CAFO emissions in accordance with the first three criteria. This approach produces local ratios on par with the literature, has the closest average residual to 0, and is conservative in its approach to attribution. Transect MVR has slightly greater uncertainty than F2 MVR and both have lower uncertainty than SM because they do not require subtraction. The impact of subtraction on uncertainty is evident in Table 1 and is the basis for including criterion 4. While the ONG ratios calculated either by SM or MVR using the F2 data differ by ∼ 40 %, the uncertainty of the results using MVR attribution is a factor of 3 lower. This is despite the fact that the SM ratio uncertainty is actually a factor of 5 lower than the MVR ratio. If we had used a comparable error on the C₂H₆ : CH₄ ratio to MVR with SM, the error would approach 100 %. The primary reason for the difference in uncertainty between the SM and MVR approaches relates to how error is propagated through subtraction versus multiplication or division. For subtraction, absolute errors add in quadrature, while for multiplication or division, relative errors add in quadrature (Skoog et al., 2004). The net effect is that error on the SM will be very high when the absolute amount to be subtracted is large. The ONG ratio is not directly used for attribution in the MVR approach, which is why the CAFO-attributed emissions by the F2 MVR and SM approaches vary by ∼ 30 %. However, the relative uncertainty (standard error or ratio) of the NH₃ : CH₄ MVR ratio is comparable to the relative uncertainty of the MVR ONG ratio.

The subtraction method is very similar to inventory subtraction in principle. Error estimates are often not provided for inventory data, but the error introduced to emission estimates by subtracting inventories may be substantial. Zimmerle et al. (2022) recently showed how ONG sources could vary and produce uncertainties in excess of estimates relying on traditional inventory data. The MVR approach is generally attractive to avoid such errors when interfering signals are high and there are enough data to produce robust MVR results. The plume from the observed CAFOs was broad, allowing for multiple data points within the plume, and produced results with low uncertainty (20 %). For sites with narrower and lower enhancements, local MVR may not be possible. Our criteria would suggest that using MVR on the full flight data may be an appropriate proxy when local ratios with MVR cannot be calculated. However, the results from the MVR Transect and F2 data are statistically different in this case. Additional analysis may need to be done to ensure that the MVR results from larger datasets are applicable to the site of interest. Repeat measurements or ground observations would help to constrain such results.

3.3 Additional case studies

To test the methodology we have developed so far for methane attribution in complicated regions more thoroughly, we have applied our strategy to a case study from a larger dataset of observations in the region. These data were collected as part of the Transport and Transformation of Ammonia (Trans2Am) campaign during the summer of 2021 and 2022, where updated flight plans were performed for CAFOs in the same region as this work. While full results from that campaign are outside the scope of this work, we compare our results to one research flight (research flight 13, denoted as RF13 in the remaining text) in the same vicinity as the facilities measured during the original F2. The updated flight plans, shown in Fig. 6, included multiple downwind transects to allow computation of emissions at consecutive distances over the course of ∼ 1.5 h. We calculated CH₄ and NH₃ emission rates and uncertainties identically to the previously described methodology. We would expect the CH₄ emission rate to be conserved at different downwind distances. The NH₃ emission may be affected by deposition to the surface and thus may remain the same or decrease downwind from a facility.

Figure 6 Path of Trans2Am RF13, colored by CH₄ (ppbv). Black numbers represent the corresponding transect number. Animal operations are indicated by the different colored circles as in Fig. 1. Green dots represent ONG wells; data of ONG are as of 2015 (Colorado Department of Natural Resources Oil & Gas Conservation Commission, 2016).

Table 2 shows results from RF13, which was collected on 23 August 2021 near Greeley, CO. The target CAFO was the same target feedlot in F2. The top of the MBL was located at 2000 ± 200 m a.g.l. An average temperature of 24 ^∘C and southwesterly winds with an average wind speed of 4.1 m s⁻¹ were observed during the flight, and a total of five downwind transects were performed. Background values and uncertainties for CH₄ were obtained from the MVR fit and varied for each transect, ranging from 1886 to 1927 ppbv. The background NH₃ value (3.023 ppbv) was calculated by constraining observations to those associated with the target CAFO (Fig. 6) and identifying the minimum NH₃ concentration. Enhanced concentrations of CH₄ and NH₃ were observed during all transects, spanning about 17 km downwind of the target facility. CAFO CH₄ was attributed using the transect MVR approach for all five transects described previously.

Table 2 CH₄ Attribution for Trans2Am RF13.

^* SE = standard error.

Attributed CH₄ ranged between 189 ± 55 to 533 ±162 g s⁻¹ with a total average of 365 ± 89 g s⁻¹. The average emission rates for this site convert to 12 ± 3 g of CH₄ per head per hour and 13 ± 3 g of NH₃ per head per hour. There is a decreasing trend in the NH₃ : CH₄ ratio observed as the transect distance increases, indicating that NH₃ depletion is occurring, likely by dry deposition. The per head values use the combined maximum number of cattle from all facilities within the plume (109 500 cattle), according to the CAFO Permit Database 2021 (CDPHE, 2022). The total average uncertainty range (95 % CI) for this site was calculated from the distribution of the sum of all five transect Monte Carlo simulations. The total average relative uncertainty for this case (24 %) is higher than the uncertainty for the previous case (20 %) due to different environmental conditions. While there is considerable variability in the CH₄ emissions from this site, we note that the NH₃ emissions show similar variability. There is a factor of 2.7 between the minimum and maximum CH₄ emission observations and a factor of 2.4 between the minimum and maximum NH₃ emission observations. Since the NH₃ emission rate is not affected by any attribution method, it is assumed that the variability in CH₄ is primarily caused by on-site or environmental conditions and not the MVR attribution method. Comparisons between the MVR and SM attribution methods for this facility show (Table S4) that the SM method produces more variable results (factor of 3.9 min–max range) and overpredicts the total average methane by about 49 %, while the total average relative uncertainty associated with the SM method for this site is almost double the uncertainty associated with the MVR method. These results are in agreement with our previous results. The average per head CH₄ and NH₃ results for this site are not statistically different from the November 2019 results despite the different wind direction, background values, etc.

3.4 Comparison to literature

A few recent studies have reported CH₄ emission rates from CAFOs. The emission factors calculated from the Transect MVR method in this work are 13 ± 3 g of CH₄ per head per hour and 13 ± 2 g of NH₃ per head per hour in November 2019 (F2) and 12 ± 3 g of CH₄ per head per hour, and 13 ± 3 g of NH₃ per head per hour August 2021 (RF13). The calculation of emission rates as a value per head is based on the maximum amount of cattle allowed in the CAFOs, although the actual number of cattle in each CAFO at the time of the measurement is unknown. Our maximum number of cattle identified for the F2 data consisted of mostly beef cattle at 161 500 heads, and only 12 300 heads were dairy cattle (total 173 800 heads). The Environmental Protection Agency (EPA) estimated CH₄ emissions for beef cattle is ∼7.2 g per head per hour and previous studies range from 7–9 g per head per hour (EPA Annex, 2021; Golston et al., 2020; Hacker et al., 2016). Dairy cattle have high rates of CH₄ emissions due to higher exertion on the animal with rates ranging from 14–39 g per head per hour and the EPA national average is 18 g per head per hour (EPA Annex, 2021; Leytem et al., 2011; Bjorneberg et al., 2009; Griffith et al., 2008). Golston et al. (2020) conducted a ground-based study in the NCFR and found the highest emission rates in dairy cattle at 39.32 ± 3.92 g of CH₄ per head per hour. In contrast, their reported emissions from beef cattle were much lower and similar to other studies (9.48 ± 0.93 g of CH₄ per head per hour). However, they noted that there was a significant difference in emission rates of repeat observation for an individual CAFO.

The difference between the CH₄ emissions per head of cattle in this work compared to prior works could be due to the management of the CAFOs. Enteric fermentation can be altered through the composition of the food given to the animals. For instance, the type and maturity of diets provided to the animals may modify the nutrients and digestibility of the food (Archimède et al., 2011). CH₄ emissions from cattle are also known to depend on exertion on the animals including exercise and stress (EPA, 2022). In this study, the food source and feeding schedule were unknown. We also note that one CAFO was extremely large with 98 000–100 000 reported maximum cattle heads during the time of study and one of the largest in Colorado. Most emission estimates are based on considerably smaller CAFOs.

The NH₃ emissions per head reported here are also slightly higher than previous estimates, which ranged from 2–12 g per head per hour (Staebler et al., 2009; Golston et al., 2020; Hacker et al., 2016). The NH₃ emission per head calculation in our work does not require removal of interfering signals like the CH₄ emissions. The NH₃ emissions are highly affected by meteorological conditions including seasonality and time of day. Waste that is stored outside may exacerbate NH₃ emissions in warmer temperatures (Montes et al., 2013). There is also evidence of a diurnal pattern of NH₃ emissions from CAFOs with a peak during noon local time, near the time of sampling in this work (Shonkwiler and Ham, 2017). NH₃ : CH₄ ratios follow a similar pattern with higher ratios during the midday and lower ratios at night (Eilerman et al., 2016).

4 Conclusions

We demonstrate an approach to isolate and quantify CH₄ plumes of agricultural CAFO sources with interfering sources by using the SM and MVR. The SM uses NH₃ as a tracer to identify a CAFO plume before using a C₂H₆ : CH₄ ratio to remove any ONG CH₄ emissions interfering with the CAFO source. It is an optimal method when there is little ONG influence to keep the error introduced from subtraction small. The MVR method uses C₂H₆ and NH₃ as tracers and provides lower uncertainty of emission estimates. This approach is appropriate when tracer data are available and there is enough signal to produce statistically significant relationships and thus site-specific ratios. The criteria to identify the best approach may be useful for isolating and attributing emissions from specific sources in other regions. Overall, our best estimates of emissions from the observed CAFOs are 13 (±2) g per head per hour for CH₄ and 13 (±2) g per head per hour for NH₃. These findings are significantly higher than EPA inventory estimates and previous studies highlighting the need for more observations and estimates of CAFO CH₄ and NH₃ emission rates.