The Pennsylvania State University The Graduate School College of Earth and Mineral Sciences BENZENE OBSERVATIONS AND SOURCE APPOINTMENT IN A REGION OF OIL AND NATURAL GAS DEVELOPMENT A Dissertation in Meteorology by Hannah Selene Halliday © 2016 Hannah Selene Halliday Submitted in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy December 2016 The dissertation of Hannah Selene Halliday was reviewed and approved* by the following: William H. Brune Distinguished Professor of Meteorology Co-Chair of Committee, Dissertation Co-Adviser Anne M. Thompson Adjunct Professor of Meteorology Co-Chair of Committee, Dissertation Co-Adviser Special Member George S. Young Professor of Meteorology and GeoEnvironmental Engineering Guido Cervone Associate Professor of Geography Johannes Verlinde Professor of Meteorology Associate Head, Graduate Program in Meteorology *Signatures are on file in the Graduate School. iii Abstract Benzene is a primarily anthropogenic volatile organic compound (VOC) with a small number of well characterized sources. Atmopspheric benzene affects human health and welfare, and low level exposure (< 0.5 ppbv) has been connected to measureable increases in cancer rates. Benzene measurements have been increasing in the region of oil and natural gas (O&NG) development located to the north of Denver. High time resolution measurements of VOCs were collected using a proton-transfer-reaction quadrupole mass spectrometry (PTR-QMS) instrument at the Platteville Atmospheric Observatory (PAO) in Colorado to investigate how O&NG development impacts air quality within the Wattenburg Gas Field (WGF) in the Denver-Julesburg Basin. The measurements were carried out in July and August 2014 as part of NASA’s DISCOVER-AQ field campaign. The PTR-QMS data were supported by pressurized whole air canister samples and airborne vertical and horizontal surveys of VOCs. Unexpectedly high benzene mixing ratios were observed at PAO at ground level (mean benzene = 0.53 ppbv, maximum benzene = 29.3 ppbv), primarily at night (mean nighttime benzene = 0.73 ppbv). These high benzene levels were associated with southwesterly winds. The airborne measurements indicate that benzene originated from within the WGF, and typical source signatures detected in the canister samples implicate emissions from O&NG activities rather than urban vehicular emissions as primary benzene source. This conclusion is backed by a regional toluene-to-benzene ratio analysis which associated southerly flow with vehicular emissions from the Denver area. Weak benzene-to-CO correlations confirmed that traffic emissions were not responsible for the observed high benzene levels. Previous measurements at the Boulder Atmospheric Observatory (BAO) and our data obtained at PAO allow us to locate the source of benzene enhancements between the two atmospheric observatories. Fugitive emissions of benzene from O&NG operations in the Platteville area are discussed as the most likely causes of enhanced benzene levels at PAO. A limited information source attribution with the PAO dataset was completed using the EPA’s positive matrix factorization (PMF) source receptor model. Six VOCs from the PTR-QMS measurement were used along with CO and NO for a total of eight chemical species. Six sources were identified in the PMF analysis: a primarily CO source, an aged vehicle emissions source, a diesel/compressed natural gas emissions source, a fugitive emissions source, and two sources that have the characteristics of a mix of fresh vehicle emissions and condensate fugitive emissions. 70% of the benzene measured at PAO on the PTR-QMS is attributed to fugitive emissions, primarily located to the SW of PAO. Comparing the PMF source attribution to source calculations done with a source array configured from the literature returns a contradictory result, with the expected sources indicting that aged vehicle emissions are the primary benzene source. However, analysis of the contradictory result indicates that the toluene to benzene ratio measured for PAO is much lower than the literature values, suggesting that the O&NG source emissions have a lower ratio of toluene to benzene than anticipated based on studies of other regions. Finally, we propose and investigate an alternative form of the source receptor model using a constrained optimization. Poor results of the proposed method are described with tests on a synthetic testing dataset, and further testing with the observation data from PAO indicate that the proposed method is not able to converge the best global solution to the system. iv Table of Contents List of Figures .............................................................................................................................................. vi List of Tables .............................................................................................................................................. xii Preface ........................................................................................................................................................ xv Acknowledgments ...................................................................................................................................... xvi Chapter 1 Introduction ................................................................................................................................. 1 Chapter 2 Benzene Analysis in Colorado .................................................................................................... 5 2.1 Introduction ......................................................................................................................................... 5 2.2 Methods............................................................................................................................................... 5 2.2.1 Site Description ............................................................................................................................ 5 2.2.2 PTR-QMS benzene measurements .............................................................................................. 6 2.2.3 TOGA VOC Collection ............................................................................................................... 7 2.2.4 Canister Whole Air Sample VOC Collection .............................................................................. 7 2.2.5 Ozonesonde Boundary Layer Height Calculations ...................................................................... 7 2.3 Results ................................................................................................................................................. 8 2.3.1 Spatial Distribution of Benzene ................................................................................................... 8 2.3.2 Benzene Observations at PAO ..................................................................................................... 8 2.3.3 Vertical Distribution of Benzene from Aircraft Measurements ................................................. 11 2.3.4 Polar Frequency Analysis of Benzene ....................................................................................... 13 2.3.5 Hydrocarbon Tracer Correlations from PAO VOC Canister Samples ...................................... 13 2.3.6 Toluene-to-Benzene Ratio ......................................................................................................... 17 2.3.7 Polar Frequency Analysis of CO and Benzene-to-CO Correlation ............................................ 18 2.4 Discussion ......................................................................................................................................... 20 2.5 Conclusions ....................................................................................................................................... 22 Chapter 3 PMF Source Identification ........................................................................................................ 23 3.1 Introduction ....................................................................................................................................... 23 3.1.1 Previous Results ......................................................................................................................... 24 3.2 Source Receptor Models and Methods ............................................................................................. 26 3.2.1 Source Receptor Models in Air Quality ..................................................................................... 26 3.2.2 PMF Model Description and Methods ....................................................................................... 28 3.2.3 Source Receptor models with VOCs, limited information ........................................................ 30 3.3 Expected Sources and VOC Tracer Selection ................................................................................... 32 3.4 PMF Analysis.................................................................................................................................... 34 3.4.1 Analysis with PMF: Full dataset results..................................................................................... 34 v 3.4.2 PMF Source Profile Identification ............................................................................................. 40 3.4.3 Benzene Analysis from the PMF Results................................................................................... 48 3.5 Discussion: The challenges of using source receptor models on limited data .................................. 54 3.6 Conclusions ....................................................................................................................................... 56 Chapter 4 A Proposed Alternative Source Receptor Model Methodology and Analysis of the Expected Sources ........................................................................................................................................................ 58 4.1 Introduction and Motivation ............................................................................................................. 58 4.2 Proposed Statistical Source Identification Method ........................................................................... 59 4.2.1 Proposed Constrained Optimization Methods ........................................................................... 59 4.2.2 Comparison of PMF and Proposed Method ............................................................................... 63 4.3 Analysis of the Expected Sources ..................................................................................................... 65 4.3.1 Analysis of the Expected Sources using a Simplified Linear Algebra Solution ........................ 65 4.3.2 Discussion: Why are the CMB Type Source Results Inverted from the PMF Source Results? 71 4.4 Conclusions ....................................................................................................................................... 73 Chapter 5 Conclusions and Summary ........................................................................................................ 74 References ................................................................................................................................................... 76 Appendix A Supplementary Material for Chapter 2 .................................................................................. 86 Appendix B Supplementary Material for Chapter 3 .................................................................................. 91 Appendix C Supplementary Material for Chapter 4 ................................................................................ 101 vi List of Figures Figure 1-1. Overview map of the study region for DISCOVER-AQ 2014. The six instrumented ground sites are marked in red and labeled with the site name. The measurements used in this dissertation were primarily collected at the PAO site, which is located within the WGF. The boundary of the WGF is shown in black. The locations of gas wells and waste pits are marked on the plot. Data for the O&NG infrastructure was supplied by the Colorado Oil and Gas Conservation Commission (cogcc.state.co.us). .. 3 Figure 2-1. Map of the 2014 DISCOVER-AQ study area. The urban areas are shown in grey (data courtesy of the United States Census Bureau, http://www.census.gov/geo/maps-data/data/tiger.html). The boundary of the WGF is shown in black along with the gas wells (brown points) (data courtesy of the Colorado Oil and Gas Conservation Commission, http://cogcc.state.co.us/). The DISCOVER-AQ ground sites are plotted and colored using the mean benzene volume mixing ratio measured during the aircraft spirals over each site. Benzene statistics were calculated using data from the bottom 1 km AGL for each site. ................................................................................................................................................................ 6 Figure 2-2. PTR-QMS benzene data collected at PAO during the DISCOVER-AQ 2014 campaign. (a) Ground level benzene time series for the entire campaign. (b) Box and whisker benzene mixing ratios organized by day, midnight to midnight MDT. Outlier points are not shown. Open points are the daily means, and these mean values are calculated from all data.The boxes indicate the locations of the 25th, 50th and 7th quartiles. The whiskers exptend to 1.5 times the interquartile range. (c) Diurnal cycle of benzene plotted in log scale. Outlier points are shown in black. Open points are the hourly mean values. All dates and times are in local time (MDT). ............................................................................................... 9 Figure 2-3. Benzene relationship to PBLH. Purple triangles are PBLH calculated from ozonesonde data, using 15 minute mean benzene data starting at balloon launch time. Gold points are mean benzene associated with PBLH derived from MPL data (binned by 50 m, collected between 0700 and 1900 MDT). Boxes and whiskers for benzene associated with the MPL derived PBLH data are shown in grey. Outlier points have been excluded for clarity. All PBLH heights are calculated AGL. ......................................... 10 Figure 2-4. Vertical distribution of benzene observed over PAO as a function of time of the day (MDT). The aircraft data are binned into 50 m boxes and plotted in box-and-whisker format. Mixing ratios are cut off at 0.60 ppbv for clarity. A total of 42 spirals have been used to construct this plot, and the number of each profile is labeled within each facet (n = number of spirals). .............................................................. 11 Figure 2-5. Polar frequency plots of benzene as measured at PAO during the DISCOVER-AQ 2014 campaign. Wind observations are binned by 10 degrees wind direction and 0.5 ms-1 wind speed. Wind speed and direction bins are plotted on a polar plot and colored by mean benzene mixing ratio for each bin. .............................................................................................................................................................. 12 Figure 2-6. Regression plots of various hydrocarbons as measured in the canister samples during the DISCOVER-AQ 2014 campaign. Propane is an O&NG tracer and acetylene is a tracer for vehicular emissions. i-pentane, n-pentane, n-butane, ethane and benzene, respectively, show a stronger correlation with propane than with acetylene. ............................................................................................................... 14 vii Figure 2-7. Regression plot of i-pentane vs. n-pentane as measured as measured in the canister samples during the DISCOVER-AQ 2014 campaign and during previous studies in NE Colorado. (aBaker et al., 2009; bGilman et al., 2013; cSwarthout et al. 2013; dThompson et al., 2014; eThis work, 2016) ............... 15 Figure 2-8. Toluene-to-benzene ratios as measured from the NCAR C-130 during the FRAPPÉ campaign at all altitudes. Data points are colored by the toluene-to-benzene ratio, with ratios < 2 colored with in blue tones, and ratios  2 presented in warm colors. All observations are sized by the measured benzene mixing ratio. High benzene levels with toluene-to-benzene  2 are only observed over the Denver urban area and close to Greeley. PAO is marked with a cyan square. .................................................................. 16 Figure 2-9. Regression plot of benzene vs. toluene as observed by PTR-QMS at PAO. (a) Regression plot of benzene vs. toluene for the entire DISCOVER-AQ 2014 campaign. The points are color coded by the toluene-to-benzene ratio (grey < 2, orange  2; the latter being associated with fresh vehicular emissions). (b) Polar frequency plot of the toluene-to-benzene ratio showing that high ratios were associated with southerly winds. .......................................................................................................................................... 17 Figure 2-10. Polar frequency plots of CO measured at PAO during the DISCOVER-AQ 2014 campaign. The wind observations are binned by 10 degrees wind direction and 0.5 ms-1 wind speed. Enhanced CO has a directional dependence associated with SW flow over the site, similar to the results for benzene (Figure 2-5). ................................................................................................................................................ 19 Figure 2-11. Regression plots of benzene vs. CO, faceted by wind direction. Calm winds had a measured wind speed less than 1 ms-1 for the minute average. The CO axis scaling changes based on wind direction. .................................................................................................................................................................... 20 Figure 3-1. An cartoon overview of the basic SRM method. Two sources (red and blue) are measured at the receptor location as atmospheric concentrations. Each source outputs VOC 1, VOC 2, and VOC 3, but in different ratios, which are shown in the known Source Type Array. The known concentrations and the known source arrays are used to calculate the unknown source contributions at each step, referred to as the weights. ................................................................................................................................................. 28 Figure 3-2. A cartoon schematic of the steps taken in the PMF SRM method. The PMF method is seeded with a random source array (S) based on the measurements, solved the linear system for the weight contributions using a non-negative solution, and uses the ME-2 to optimize the randomized starting source arrays. .............................................................................................................................................. 29 Figure 3-3. Observational time series of the species used in the source identification at PAO. All data was collected at the surface at PAO and reported at 1-minute resolution. All mixing ratios are reported in parts per billion by volume (ppbv). ..................................................................................................................... 35 Figure 3-4. Source contribution time series for the six sources identified in the p6 results from the PMF analysis. All source contributions are shown on the same scale to show relative importance throughout the source identification time period................................................................................................................. 41 Figure 3-5. The diurnal cycle of the source contributions from the p6 PMF results displayed in with a log scale on the y axis. The x axis shows hours of day in MDT (local time). All data is shown as a box and viii whisker plot, with the box showing 25th, 50th, and 75th quartiles, the whiskers extend to the highest and lowest points within 1.5 times the inter-quartile range. Any points beyond the whisker ranges are considered outlier data and shown as points. Each source is shown on its own scale. ............................... 42 Figure 3-6. Mean polar frequency analysis of the source contributions determined with the p6 PMF analysis. The wind speeds and directions are binned by 0.5 m/s and 10 degrees. The mean value of each source contribution is shown as a color map. The plot shows the mean contributions for A. source 1, B. source 2, C. source 3, D. source 4, E. source 5, and D. source 6. ............................................................... 43 Figure 3-7. A time series comparison of the source weight contributions from source 2 and 6 from the p6 PMF analysis. The top panel compares source 2 and 6. The bottom panel compares source 3 and the difference between source 2 and 6. Units can be interpreted as ppbv. ...................................................... 45 Figure 3-8. A polar frequency comparison of the maximum source contributions from A. source 2 and B.6 from the p6 PMF analysis. B. The maximum polar frequency analysis of the difference between source 2 and source 6, considered to the fugitive emission contribution to source 2. D. The maximum polar frequency plot for source 3, identified as fugitive emissions. ........................................................... 46 Figure 3-9. The time series of the benzene mixing ratios by source type, from the source array solution to the p6 PMF results. Each source type is shown with its own scale. The contributions to benzene by each source type are colored by source, and faceted by type for clarity. The labels for each source indicate the assigned source, and the percentage of each source that is made up of benzene. ....................................... 47 Figure 3-10. The benzene mixing ratio measurements separated by source for the full source separation dataset. All data is shown as a box and whisker plot, with the box showing 25th, 50th, and 75th quartiles, the whiskers extend to the highest and lowest points within 1.5 times the inter-quartile range. Any points beyond the whisker ranges are considered outlier data and shown as points. The plots are shown on a log scale to show the variability between sources at low mixing ratios. .......................................................... 48 Figure 3-11. The benzene observations separated by campaign section: early, middle, and late. A. The observed benzene time series, colored by campaign period. B. The daily (midnight to midnight) benzene statistics, colored by campaign period. ....................................................................................................... 49 Figure 3-12. The benzene mixing ratio measurements by source, colored by collection period for the full source separation dataset. The plots are shown in log scale for better comparison of the low values. All data is shown as a box and whisker plot, with the box showing 25th, 50th, and 75th quartiles, the whiskers extend to the highest and lowest points within 1.5 times the inter-quartile range. Any points beyond the whisker ranges are considered outlier data and shown as points. The plots are shown on a log scale to show the variability between sources at low mixing ratios. The box plots are colored by collection period. .................................................................................................................................................................... 50 Figure 3-13. Percentage contribution to the total calculated benzene by source. The 0% and 100% contributions are shown as dotted lines for reference. Each source is colored and faceted for clarity, and shown on its own scale. .............................................................................................................................. 51 ix Figure 3-14. Percentage contribution to the total calculated benzene by source, with source contributions less than 0% removed from source 3 and source contributions great than 100% removed from source 4. The 0% and 100% contributions are shown as dotted lines for reference. Each source is colored and faceted for clarity, and shown on its own scale. ......................................................................................... 52 Figure 3-15. The total calculated benzene correlated against total calculated CO, unmodified (blue), and modified to remove the primarily CO source (source 1) (red). The plots are faceted by wind direction; each direction encompasses 90 degrees. The calm winds are any observations with a wind speed less than 1 m/s for the minute average. Note that the x-axis change scale with wind direction, with the largest values of CO observed from the SW. The correlation values for each wind facet are listed, and color coded by type of CO. .................................................................................................................................. 54 Figure 4-1. A cartoon schematic of the constrained optimization algorithm. ............................................ 62 Figure 4-2. Observational time series of the species used in the source identification at PAO. All data was collected at the surface at PAO and reported at 1 minute resolution. All mixing ratios are reported in parts per billion by volume (ppbv). ..................................................................................................................... 64 Figure 4-3. The diurnal profiles of the expected source contributions at PAO, plotted by hour of day (local time, MDT). The labels are abbreviated to maintain plot clarity: vehicle emissions (VE), aged vehicle emissions (AVE), fugitive emissions (FE), condensate emissions (CE), compressed natural gas (CNG), and carbon monoxide (CO). ........................................................................................................... 67 Figure 4-4. Time series of the calculated benzene CMB type solution with the expected sources, colored by source contribution. A. The total time series of all calculated benzene mixing ratios. B. The time series from -1 ppbv to 5 ppbv calculated benzene, to more clearly show the source dynamics on the low end of the scale where the sources other than aged vehicle emissions contribute to the calculated benzene mixing ratios. ........................................................................................................................................................... 69 Supplementary Figure 2-1. Regression plot of 2-min average PTR-QMS benzene (with the mean benzene window being centered on the canister opening time) vs. benzene in the canister samples. The grey dashed line is the 1-to-1 line. The black line is the linear regression best fit line, which has a slope of 1.003, and an R2 of 0.9205. The points are colored by their hour of collection, in local time (MDT). ...... 86 Supplementary Figure 2-2. PBLH diurnal profile. The box and whisker plots are created from MPL observations between 22 July 2014 and 11 August 2014. The open points show the mean PBLH derived from MPL data for each hour. The red points show the mean PBLH derived from ozonesonde data for each hour. The daytime-type boundary layer starts developing between 900 and 1000 MDT. The MPL PLBH data is based on aerosol retrievals, which remain lofted into the atmosphere after the collapse of the PBLH at night, so we consider the post-sunset data to be unreliable for PBLH measurements. .......... 87 Supplementary Figure 2-3. Vertical benzene profiles and ground benzene measurements collected at PAO during the DISCOVER-AQ 2014 study. Each point is plotted based on the time stamp and GPS altitude data from the P3-B spirals. The altitudes are reported above ground level (AGL). Each point is colored by measured benzene mixing ratio. The color plot is based on the quartiles of the benzene data. .................. 88 x Supplementary Figure 2-4. Benzene mean profile over PAO from the aircraft spirals over the site. The points are the mean mixing ratios binned by 50 m bins. The error bars show the standard deviation of the mean for each vertical bin. .......................................................................................................................... 89 Supplementary Figure 2-5. Waste pit locations and statuses provided by COGCC. The cyan points show the locations of the PAO and BAO sampling sites for the campaign. The black outline shows the political boundary of the WGF as defined by the COGCC. The red triangles show the locations of the active waste pit locations as reported from the COGCC data. The green square shows the location of the Platte Valley Gas Plant. .................................................................................................................................................... 90 Supplementary Figure 3-1. Comparison of the observations vs calculated mixing ratios for the p5 PMF model run. Each species is shown on its own scale. Blue is observed data, red is calculated mixing ratios from the p6 PMF solution. .......................................................................................................................... 94 Supplementary Figure 3-2. Comparison of the observations vs calculated mixing ratios for the p6 PMF model run. Each species is shown on its own scale. Blue is observed data, red is calculated mixing ratios from the p6 PMF solution. .......................................................................................................................... 95 Supplementary Figure 3-3. The diurnal cycles of the source profiles for the p6 PMF results. All data is shown as a box and whisker plot, with the box showing 25th, 50th, and 75th quartiles, the whiskers extend to the highest and lowest points within 1.5 times the inter-quartile range. Any points beyond the whisker ranges are considered outlier data and shown as points. Each source is shown on its own scale. The log scale version of this plot is shown in Supplementary Figure 3-3. .............................................................. 96 Supplementary Figure 3-4. Median polar frequency analysis of the p6 source contributions. The wind speeds and directions are binned by 0.5 m/s and 10 degrees. The median value of each source contribution is shown as a color map. The plot shows the median contributions for A. source 1, B. source 2, C. source 3, D. source 4, E. source 5, and D. source 6. .............................................................................................. 97 Supplementary Figure 3-5. Maximum polar frequency analysis of the p6 source contributions. The wind speeds and directions are binned by 0.5 m/s and 10 degrees. The maximum value of each source contribution is shown as a color map. The plot shows the maximum contributions for A. source 1, B. source 2, C. source 3, D. source 4, E. source 5, and D. source 6. ............................................................... 98 Supplementary Figure 3-6. Correlation between the source contributions from Source 2 and Source 3 in the p6 PMF results. The coefficient of determination is 0.6625, and the slope is 0.87 source 3 to source 2. .................................................................................................................................................................... 99 Supplementary Figure 3-7. Correlation between the source contributions from Source 2 and Source 6 in the p6 PMF results. The coefficient of determination is 0. 0.5787, and the slope is 0.5866 source 6 to source 2. .................................................................................................................................................... 100 Supplementary Figure 4-1. The time series source weight contributions calculated from the expected sources from Table 1 using the CMB type solution. The labels are abbreviated to maintain plot clarity: vehicle emissions (VE), aged vehicle emissions (AVE), fugitive emissions (FE), condensate emissions (CE), and compressed natural gas (CNG). ................................................................................................ 102