In this study, we evaluate the scale of methane emissions from the U.S. NG-based ammonia fertilizer industry. Since ammonia fertilizer plants are dispersed all over the U.S., a mobile sensing approach was adopted to be able to quickly cover a representative number of plants. Using data collected with a mobile sensor, a Bayesian approach was used to characterize methane emissions from the plants surveyed. These estimates were up-scaled to infer methane emissions from the U.S. ammonia fertilizer industry.

The ammonia fertilizer industry is a major industrial consumer of NG, which is used as feedstock and energy source to produce ammonia and some upgraded products (e.g. urea) ( Elvers et al., 1989 ). In 2015, the U.S. gross production capacity of ammonia (including ammonia as a feedstock for upgraded products) was >12,000 Gigagram per year (Gg NH 3 /yr) ( United States Geological Survey, 2017 ). Due to stable low NG prices resulting from the U.S. shale gas boom, the U.S. ammonia production capacity is expected to increase in the next several years ( United States Geological Survey, 2017 ). Consequently, a large amount of NG is and will continue to be consumed in the production of ammonia fertilizer. According to the Manufacturing Energy Consumption Survey (MECS) conducted by the Energy Information Administration (EIA) ( Energy Information Administration, 2017 ), 383 trillion British thermal units (BTU) of NG (approximately 7,700 Gg of NG) was consumed by the U.S. ammonia fertilizer industry in 2014. Ammonia fertilizer plants are large GHG emissions sources mainly because CO 2 is a major byproduct of ammonia production. However, there is less understanding of methane emissions despite the large amount of NG consumption by this industry.

Material and Methods

Field experiments Two sampling campaigns were conducted, the first from Jun. 16th to 19th, 2015, and the second from Sep. 12th to 29th, 2016. Three and eight ammonia fertilizer plants were sampled, respectively, during the two campaigns. Two of the three plants sampled in the first campaign were re-visited in the second campaign. As shown in Figure 1, the plants surveyed were mostly located in the mid-western U.S. surrounded by relatively flat terrain, resulting in few obstacles impacting plume transport. The downwind measurements were conducted without notifying plant operators with intent of encountering routine plant operation. Water vapor from cooling towers and effluent flares were observed when sampling all the plants, suggesting that they were in production. Useful data were collected along downwind roads at six of the nine plants sampled, all during the second field campaign in 2016. These six plants represent >25% of the total of 23 NG-based ammonia fertilizer plants in the U.S. that were operating in 2015–2016. No useful data were obtained at the other three plants (two were visited during both campaigns) due to poor meteorological and road conditions, and hence were not included in the analysis. The data quality control protocol ruled out some data from the six plants surveyed (Supplemental materials Text S1). Key characteristics and sampling conditions related to the six plants that were successfully surveyed are summarized in Table 1. Facility location

(Latitude, Longitude) Date

[MM/DD/YY] Local time

[HH:MM] N

[-] x s

[m] ū

[m/s] θ m

[deg] Verdigris, OK (36.2338, –95.7190) 09/12/16 10:42–14:48 52 1730 2.90 164 Enid, OK (36.3794, –97.7624) 09/14/16 10:24–14:33 45 350 2.80 8.7 09/15/16 10:02–12:38 48 1230 3.84 184 Dodge City, KS (37.7777, –99.9299) 09/20/16 16:16–18:58 8 2170 5.97 183 09/21/16 10:19–13:49 15 1880 5.41 196 09/22/16 10:17–12:01 14 1800 5.67 195 Beatrice, NE (40.3202, –96.8401) 09/26/16 09:47–12:04 24 590 2.50 291 Creston, IA (41.1170, –94.3551) 09/26/16 17:41–18:41 16 650 2.56 305 09/27/16 10:38–13:08 37 590 3.51 304 09/27/16 14:38–15:38 12 730 3.73 291 Fort Dodge, IA (42.4992, –94.0171) 09/28/16 13:08–15:38 34 1400 3.43 355 09/29/16 10:52–12:22 5 1440 2.77 352 09/29/16 14:52–16:52 9 1410 3.12 347 The mobile sensing approach utilizes a Google Street View (GSV) car equipped with a fast-response and high-precision methane analyzer to repeatedly measure methane mixing ratios (c, in ppm) along public roads downwind of the ammonia fertilizer plants. Here, we only offer a brief description of some components of GSV car that are relevant to this field campaign, and more detail on the GSV car can be found elsewhere (von Fischer et al., 2017). A Picarro cavity ringdown analyzer (Model G2301 from Picarro Inc., Santa Clara, CA, USA) configured to sample at 1 Hz was used in the first sampling campaign, and a fast methane/ethane gas analyzer (Los Gatos Research, Inc., San Jose, CA, USA) sampling at 2 Hz was used in the second sampling campaign. We did not re-calibrate the Los Gatos Research (LGR) fast methane/ethane gas analyzer during deployment. After initial factory calibration, the instruments largely operate on a self-calibrating basis, accounting for variation in performance of the analytical cell using a wavelength that was not absorbed by any gases. Otherwise, the analysis is derived from first-principles (cell dimensions, pressure, temperature, and methane’s per-molar absorptivity) (Gupta, 2012). Together, these concentration measures yield methane values that are accurate to better than 60 ppb (Robert Provencal, LGR, personal communication, 04/16/2019), which is relatively small considering that our methane data vary from ~2.0 ppm (ambient level) to ~30 ppm (peak mole fraction). Though ethane concentration were also measured by the gas analyzer, it was not used for analysis since the concentrations were low and the noise was high, resulting a low signal to noise (SNR). The inlet of the methane analyzer was located at the front bumper of the GSV car (~30 cm above-ground) to avoid measuring exhaust of the car. Real-time location of the GSV car was determined by a roof-mounted GPS unit (Model A100 from Hemispheres GNSS, Scottsdale, AZ, USA) with an acquisition frequency of 1 Hz. A portable 3-D sonic anemometer (Model 81000 from R.M. Young Co., Traverse City, MI, USA) was installed near each visited fertilizer plants in a relatively flat and open location to measure local meteorological conditions (Table 1). The height of the anemometer was ~1.5 m above the ground. During the sampling, the GSV car first circled around the plant to be surveyed to map the presence of elevated methane mixing ratios, and to determine whether the pattern of emissions was clear enough such that attribution of any emissions could be made based on wind direction (similar to triangulation). An example is shown in Figure 2, a well-developed methane plume was observed downwind of the ammonia fertilizer plant, while ambient background methane mixing ratios were found when sampling upwind of the plant. After identifying the facility-introduced plume, the GSV car would pass through the plume multiple times to make repeated samplings. Local emissions, such as small pipeline leaks, can be identified as small spikes in the data (Figure 2). In contrast, emissions from the ammonia plant (sampled approximately 500 to 2,000 m downwind as shown in Table 1) are characterized by a broad plume (Figure 2). In the data post-processing, elevated mixing ratios due to local emissions were excluded from the analysis, and we focused on the broad plumes. A low-pass moving average filter was established based on the low frequency drift of the internal sensor pressure, which is likely caused by partial clogging of the sensor inlet due to dust. The filter was later used to mask and adjust the raw mixing ratio data to remove the adverse effects caused by this issue. The background methane mixing ratio (c b ) was estimated as the 5th percentile of the ranked time series of raw methane mixing ratios (c r ) (Albertson et al., 2016; Brantley et al., 2014), and the above-ambient mixing ratios c = c r – c b . The estimated c b is very close to the ambient methane mixing ratios measured upwind of the plant, suggesting that the determination of c b is robust.

Emission characterization using a Bayesian approach Using data collected during multiple sensor passes, a Bayesian approach was adopted to characterize methane emission from the ammonia fertilizer plants. This approach is fully described in (Albertson et al., 2016); here only a brief introduction is included. Following Bayes’s rule, the posterior probability density function (PDF) of the emission rate Q [kg/h] is (Albertson et al., 2016; Yee, 2008; Yee, 2012): P ( Q | c y , I ) = P ( Q | I ) P ( c y | Q , I ) P ( c y | I ) , M1 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ P\left({Q|{c_y},I} \right) = \frac{{P\left({Q|I} \right)P\left({{c_y}|Q,I} \right)}}{{P\left({{c_y}|I} \right)}}, \] \end{document} (1) where c y (in ppm × m) is the cross-plume integrated methane mixing ratio. Practically, c y can be estimated as c y = ∑cΔx, where Δx [m] is the distance between the geo-referenced mixing ratio data. I is the ancillary information, including source information and the prevailing meteorological conditions. P(Q|I) is the prior PDF, which represents the distribution of Q prior to the observation of c y . P(c y |Q, I) is the likelihood function, which is the probability of observing c y given Q and I. P(c y |I) is the evidence term that simply ensures that P(Q|c y , I) integrates to unity. Before making any measurement, the prior (P(Q|I)) is unknown. Thus, a uniform distribution is adopted as a non-informative prior (Yee, 2007). After the first sensor pass (with a valid measurement of c y ), equation (1) is updated recursively such that P(Q|I) is replaced by the posterior PDF (P(Q|c y , I)) derived from the previous sensor pass (Albertson et al., 2016). P ( Q | I ) = { 1 / ( Q max − Q min ) , j = 1 P ( Q | c y , I ) j − 1 , j > 1 , M2 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ P\left({Q|I} \right) = \left\{{\begin{array}{*{20}{c}} {1/\left({{Q_{max}} - {Q_{min}}} \right),j = 1}\\ {\,\,\,\,\,\,\,P{{\left({Q|{c_y},I} \right)}_{j - 1}},j > 1} \end{array}} \right., \] \end{document} (2) where j is a counter for successive sensor passes, Q max and Q min are the upper and lower bound of Q. All the information provided by the measurement c y about the unknown Q is contained in the likelihood function. Following (Yee, 2008; Yee, 2012) and (Albertson et al., 2016), the likelihood function is formulated as: P ( c y | Q , I ) = 1 2 π σ e exp ( − 1 2 ( c y − c y M ( Q ) σ e ) 2 ) , M3 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ P\left({{c_y}|Q,I} \right) = \frac{1}{{\sqrt {2\pi} {\sigma _e}}}{\rm{exp}}\left({ - \frac{1}{2}{{\left({\frac{{{c_y} - c_y^M\left(Q \right)}}{{{\sigma _e}}}} \right)}^2}} \right), \] \end{document} (3) where c y M ( Q ) M4 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ c_y^M(Q) \] \end{document} is the modeled c y as a function of the candidate emission rate Q. σ e [-] is the “error term” of the likelihood function, which is a measure of the uncertainty when comparing the modeled c y M ( Q ) M5 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ c_y^M(Q) \] \end{document} against the measurement c y . The parameterization of σ e is detailed in the Supplemental materials Text S2. When the vehicle’s path is perpendicular to the wind direction, the effect of plume lateral dispersion vanishes and the modeled c y M ( Q ) M6 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ c_y^M(Q) \] \end{document} can be simplified as c y M ( Q ) = Q U ¯ D z M7 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ c_y^M(Q) = {\textstyle{Q \over {\bar U}}}{D_z} \] \end{document} , where Ū is the plume advection speed, and D z accounts for the plume vertical dispersion (Albertson et al., 2016). A common form of D z = A z ¯ exp [ − ( B z z ¯ ) s ] M8 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {D_z} = {\textstyle{A \over {\bar z}}}\,{\rm{exp}}\,\,[ - {({\textstyle{{Bz} \over {\bar z}}})^S}] \] \end{document} , where z [m] is the height of the sensor inlet, z̄ [m] is the mean plume height, and A [-], B [-], and s [-] are all empirical parameters of atmospheric stability and z̄ (Gryning et al., 1987). In the case of a non-perpendicular sensor pass, c y M ( Q ) M9 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ c_y^M(Q) \] \end{document} can be calculated by numerically integrating along the vehicle’s path (Albertson et al., 2016). By updating the prior term P(Q|I) with the posterior (P(Q|c y , I)) derived from the previous sensor pass, P(Q|c y , I) is calculated recursively to incorporate data collected after each pass. After the final sensor pass (the Nth pass), the expected emission rate (Q e , in kg/h) and the associated uncertainty (σ Q , in kg/h) can be quantified as the expectation and standard deviation of the final posterior PDF P(Q|c y , I). This approach considers the part of the facility that is related to NG-consumption (e.g. chemical reactors and boilers) as a point source, and the emissions are continuous and quasi-steady (i.e. Q is quasi-steady). The point-source assumption is based on the observation that the NG-related equipment is usually clustered within a radius of ~100 m from the center region (as visually inspected from Google Earth), while the source-to-sensor distances are 500 to 2,000 m (Table 1). We considered the point-source assumption as an added uncertainty in the source characterization model, and explicitly included it in the σ e parametrization (Supplemental materials Text S2). The assumption of continuous and quasi-steady emissions comes from the premise that NG emissions originate from the chemical production and associated fuel combustion processes. These processes are often continuous and quasi-steady under normal operational conditions of typical chemical plants in order to maximize production efficiency (Green and Perry, 1973). Based on this assumption, the estimated Q e and σ Q (in the unit of kg/h) can be converted into the unit of Gg CH 4 /yr, considering 340 days per year of effective production (i.e. excluding planned maintenance and normal outrages) as suggested by previous reports and surveys (Paul et al., 1977; The Fertilizer Institute, 2006; United States Geological Survey, 2017).