Subscriber access provided by Warwick University Library
Article
A Spatiotemporal Industrial Activity Model for Estimating the Intensity of Oil and Gas Operations in Colorado William B. Allshouse, John L. Adgate, Benjamin D. Blair, and Lisa M. McKenzie Environ. Sci. Technol., Just Accepted Manuscript • DOI: 10.1021/acs.est.7b02084 • Publication Date (Web): 17 Jul 2017 Downloaded from http://pubs.acs.org on August 3, 2017
Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.
Environmental Science & Technology is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.
Page 1 of 28
Environmental Science & Technology
1
A Spatiotemporal Industrial Activity Model for Estimating the Intensity of Oil and Gas
2
Operations in Colorado
3
William B. Allshousea, John L. Adgatea, Benjamin D. Blaira, Lisa M. McKenziea*
4 5
a
6
Aurora, Colorado, USA
Department of Environmental and Occupational Health, Colorado School of Public Health,
7
8
9
10
11
12
13
14
15
*
16
Health, Colorado School of Public Health, 13001 E 17th Pl, Campus Box B119, Aurora,
17
Colorado 80045 USA. Telephone: 303.724.5557. E-mail:
[email protected] Corresponding Author: Lisa McKenzie, Department of Environmental and Occupational
1 ACS Paragon Plus Environment
Environmental Science & Technology
Page 2 of 28
18 19 20 21
Abstract Oil and gas (O&G) production in the United States has increased in the last 15 years and
22
operations, which are trending towards large multi-well pads, release hazardous air pollutants.
23
Health studies have relied on proximity to O&G wells as an exposure metric, typically using an
24
inverse distance weighting (IDW) approach. Since O&G emissions are dependent on multiple
25
factors, a dynamic model is needed to describe the variability in air pollution emissions over
26
space and time. We used information on Colorado O&G activities, production volumes, and air
27
pollutant emission rates from two Colorado basins to create a spatiotemporal industrial activity
28
model to develop an intensity-adjusted IDW well count metric. The Spearman correlation
29
coefficient between this metric and measured pollutant concentrations was 0.74. We applied our
30
model to households in Greeley, Colorado, which is in the middle of the densely developed
31
Denver-Julesburg basin. Our intensity-adjusted IDW increased the unadjusted IDW dynamic
32
range by a factor of 19 and distinguishes high intensity events, such as hydraulic
33
fracturing/flowback, from lower intensity events, such as production at single well pads. As the
34
frequency of multi-well pads increases, it will become increasingly important to characterize the
35
range of intensities at O&G sites when conducting epidemiological studies.
36
2 ACS Paragon Plus Environment
Page 3 of 28
37 38
Environmental Science & Technology
Introduction Widespread implementation of technological advancements for oil and gas (O&G)
39
extraction, such as horizontal drilling and hydraulic fracturing, have increased the spatial extent
40
of economically viable O&G reserves, leading to more development near where people live.1
41
The state of Colorado, which currently ranks as the 7th largest oil and 6th largest gas producing
42
state, has at least 378,000 people living within 1 mile of an active O&G well.2-4 The increasing
43
intersection of O&G development and homes has led to concern about the possible health
44
consequences of residing near O&G operations. One focus of attention is that a wide range of
45
hazardous air pollutants (HAPs), including diesel exhaust, benzene, ethylbenzene, toluene, and
46
xylenes (BTEX), particulate matter, and polycyclic aromatic hydrocarbons (PAHs), are emitted
47
during each stage of O&G well development.5-7 Exposure to HAPs associated with O&G
48
development may affect the health of people living near the sites.8 In addition, people living near
49
O&G development experience odors, noise, and night-time light, which are associated with
50
increased psychosocial stress, sleep disturbance, and cardiovascular health.9-16
51
In the absence of spatially and temporally resolved air monitoring at multiple locations,
52
comprehensive emission inventories, and detailed meteorological data necessary for robust
53
modeling approaches; proximity models have been used to retrospectively approximate exposure
54
to HAPs and other stressors emitted from or associated with O&G wells.17-21 The first proximity
55
models assumed that exposure to HAP emissions and other stressors were a function of the
56
number of O&G wells within a specific buffer around a home, weighted by distance.17, 18, 20
57
These first studies used an inverse distance weighted (IDW) count to estimate the density of
58
O&G wells within 10 miles of residential address as an individual-level exposure metric. These
59
studies found statistically significant associations between IDW well counts and congenital heart 3 ACS Paragon Plus Environment
Environmental Science & Technology
Page 4 of 28
60
defects, acute lymphocytic leukemia, and low birth weights.17, 18, 20 These preliminary studies did
61
not account for the differences in emissions among the phases of well development, levels of
62
activity among well pads, or the impact of individual well pad infrastructure, such as valves and
63
tanks.
64
More recently, the IDW approach was expanded to include all wells regardless of
65
distance from residence and separate IDW metrics for phases of well development (preparation,
66
drilling, hydraulic fracturing, and production). One study in Pennsylvania used a sum of the
67
normalized, phase-specific activity metrics and found an association between proximity to
68
natural gas activity and both pre-term birth and high-risk pregnancy.19 Due to high collinearity
69
among phase-specific metrics, another Pennsylvania study modelled phase-specific metrics
70
separately and found an association between each phase of natural gas activity and asthma
71
exacerbation.21 These expanded IDW models are able to estimate an individual’s exposure
72
relative to other individuals within each well development phase, but only include one phase at a
73
time and assume equivalent intensity of activity for all phases.
74
However, levels of activities and the magnitude of HAP emissions differ by phase of
75
development and emission controls.8, 22, 23 Drilling and hydraulic fracturing involve numerous
76
truck trips and often use diesel engines to power the process. Flowback after hydraulic fracturing
77
has been associated with some of the largest HAP emissions during the well development life
78
cycle.8, 22, 23 Shortly after the well is completed, it goes into production, which brings
79
hydrocarbons from the target formation to the surface, where they are either stored in tanks or
80
sent into a pipeline. A well pad with a large number of tanks rather than pipelines is likely to
81
have higher air pollutant levels due to tank emissions.24 In addition, the United States
4 ACS Paragon Plus Environment
Page 5 of 28
Environmental Science & Technology
82
Environmental Protection Agency (EPA) estimates that O&G wells with emission controls (e.g.,
83
green completions) can reduce O&G emissions at numerous infrastructure points by up to 95%.25
84
The IDW model approach would be improved by accounting for both temporal inter and
85
intra-well variability in intensity among well development phases, as well as the number of O&G
86
wells and types of equipment on the pad, and the use of emissions controls. Due to the variety of
87
HAPs that are produced at a well pad, as well as other stressors, selecting either one or multiple
88
pollutants to model may not adequately approximate exposure to all stressors associated with
89
O&G development. However, estimating the relative intensity at each well pad according to the
90
type of activity (i.e., phase of well development), number of wells, and volume of production
91
provides an estimate of cumulative exposure to multiple stressors and adds a temporal
92
component to O&G exposure assessment that has not been fully developed to date. This
93
temporal component adds intra-well variability across phases and within the production phase,
94
making for a more dynamic cumulative exposure model.
95
In this manuscript, we outline development and implementation of an industrial activity
96
model for retrospectively estimating the temporal and spatial intensity of O&G development
97
using data on activities and production volumes from the Colorado Oil and Gas Information
98
System (COGIS). This approach builds upon previously used IDW models to better delineate
99
exposure categories for health studies, especially during the short-lived, critical exposure
100
windows that are theorized to exist for many developmental and childhood health effects.26-28
101
102
Materials and Methods
103
Data Sources 5 ACS Paragon Plus Environment
Environmental Science & Technology
Page 6 of 28
104
We built a relational database of production volumes, well pad activities, and associated
105
metadata from information in the COGIS to build a spatiotemporal model of O&G activity. The
106
Colorado Oil and Gas Conservation Commission (COGCC) maintains the COGIS, which
107
contains a wide variety of data about upstream O&G extraction in the state.29 From the COGIS,
108
we obtained the location, dates of activities, and the amount of O&G produced by month for
109
each well in Colorado. Dates of activities include the spud date (when drilling begins), treatment
110
dates (when hydraulic fracturing occurs), whether green completion techniques were used to
111
control flowback emissions, first production date, shut in dates (temporarily plugging a well),
112
and abandonment date (permanently plugging a well). Based on a unique well pad identifier, we
113
then determined which wells were co-located on multi-well pads. For pads with more than 2
114
wells, we obtained the numbers of tanks on each well pad from the 2012 inspection report
115
available in the COGIS database and verified this information using the Google Earth satellite
116
image of the site that was closest in time to the inspection report. If a 2012 inspection report was
117
not available, we used the closest report to the year 2012 and the Google Earth image that was
118
the closest match to the report date. For pads with 1 or 2 wells, we assumed that equipment
119
would not differ significantly at the site based on observations of these sites using Google Earth
120
and data contained in inspection reports within the COGIS.
121
Emissions data from two studies, one conducted in the Piceance Basin in Garfield
122
County, Colorado between 2013 and 2015 and the other in the Denver-Julesburg Basin (DJB)
123
along Colorado’s Front Range area between 2014 and 2016, were used to calibrate the relative
124
phase specific intensity factors in our model.22, 23 These studies conducted experiments around
125
several well pads undergoing a particular activity to determine emission rates for specific air
126
pollutants during that activity. Activities included 5 well pads being drilled, 11 that were 6 ACS Paragon Plus Environment
Page 7 of 28
Environmental Science & Technology
127
undergoing hydraulic fracturing/flowback, and 11 in production. Specific pollutant emission
128
rates were determined by releasing a tracer of acetylene at the well pad, collecting a downwind
129
ambient air sample, identifying the tracer in the downwind sample, and using the ratio of a
130
specific background-corrected pollutant concentration, also measured in the downwind sample,
131
to the tracer at the downwind location to derive the emission rate of the pollutant at the source.
132
We obtained the addresses of all residential properties in Weld County, Colorado from
133
the DataQuick Information Systems database30 (as of 2012). The Google Maps Geocoding
134
Application Programming Interface (API) was used to geocode the addresses of these properties
135
and we kept those whose geocode could be returned with “Rooftop” accuracy and was within 2
136
miles of the Greeley, Colorado city limit (N=32,737). Greeley sits in the middle of the DJB and
137
has notable O&G operations in and around city limits.2 We modelled the intensity of O&G
138
activities within 10 miles of this population (Eq. 1) and compared the modelled O&G activity to
139
the IDW exposure metric that does not adjust for intensity.
140
Model Development and Parameters
141
We used the information from COGIS to inform the construction of our spatiotemporal
142
model (Eq. 1) that estimates monthly intensity of O&G activity at each well pad, based on the
143
specific events (pad preparation, drilling, hydraulic fracturing and flowback, and production),
144
production volume, and O&G infrastructure. The time-varying intensity from a well pad can be
145
estimated as:
146
∑ = ∑ , + , + , + + +
(1)
7 ACS Paragon Plus Environment
Environmental Science & Technology
Page 8 of 28
147
, where ZL(tm) is the intensity score Z for well pad L during month tm. This intensity score is a
148
summation of activities at each of the nw number of wells on the pad over the nd number of days
149
in the month. The variable IP,j indicates whether the pad is being prepared on day j, and we
150
assume that this occurs 30 days prior to the first well spud date on the pad.19, 21 The variable ID,ij
151
indicates whether well i is being drilled on day j, with the assumption that drilling begins on the
152
spud date and whose duration is dependent on the year and the total well depth. Based on
153
available information, we assume that it took 48 hours to drill 1000 feet in 2008 and 20 hours to
154
drill 1000 feet in 2013.31, 32 We use these two time points to estimate the time needed to drill a
155
well in a particular year, based on its measured depth. For wells that we did not have a measured
156
depth, we used the average number of days needed to drill a well for that year based on the wells
157
that did have a measured depth. The variable IF,ij indicates whether well i is undergoing hydraulic
158
fracturing/flowback on day j. Hydraulic fracturing and flowback typically last 2-4 days and 7-12
159
days, respectively.22 Based on the data available in the COGIS, we cannot distinguish between
160
these two activities and assume that both fracturing and flowback begins on the (hydraulic
161
fracturing) treatment date from COGIS, and that the process lasts 14 days. The parameter αi
162
represents whether green completion was used during the hydraulic fracturing/flowback phase.
163
This parameter is a function of control device efficiency, rule effectiveness, and capture
164
efficiency. Using the Colorado Department of Public Health and Environment’s parameters for
165
the 2008 state implementation plan of 95% control device efficiency, 83% for rule effectiveness
166
and 75% capture efficiency for tank emissions, we calculate that αi=0.41 for green completions
167
and αi=1 for uncontrolled completions.33 The variables Oi and Gi represent the percentile of
168
O&G produced by well i, respectively, compared to oil (or gas) production of all active wells by
8 ACS Paragon Plus Environment
Page 9 of 28
Environmental Science & Technology
169
month. The variable TL is the number of tanks on the well pad as determined by using inspection
170
reports and Google Earth.
171
Model Intensity Factors
172
The intensity vector β was based on the sum of BTEX, propane, ethane, and n-hexane
173
emission rates from the previously described O&G phase-specific emissions experiments.22, 23
174
We derived our intensity factors using respective medians of the log-transformed sample
175
emission rates for these eight O&G pollutants for each activity phase (drilling, hydraulic
176
fracturing/flowback, and production). Monthly oil (β3) and gas (β4) production intensity were
177
each set to 5, and other phase intensities were calculated relative to the production intensity using
178
the emission rate data. Given that a well with an average O&G production would have a
179
production intensity of 5, we then used the ratio of the median log-transformed drilling emission
180
rate to that of the production emission rate to obtain our drilling intensity factor and used the
181
same method for deriving our hydraulic fracturing/flowback intensity. These calculations yielded
182
a relative drilling intensity of β1=25 and hydraulic fracturing/flowback relative intensity of 60 for
183
controlled flowback. Since all flowback emissions experiments were conducted at sites using
184
green completion techniques, we used the green completion emission reduction factor (αi), to
185
obtain a relative intensity β2 =146 for uncontrolled flowback. See Supporting Information for
186
more detail on deriving intensity factors.
187
The tank intensity (β5) was estimated using data on total uncontrolled volatile organic
188
compounds (VOC) emissions and the corresponding number of tanks subject to an EPA consent
189
decree against Noble Energy.33, 34 The total uncontrolled VOC emissions were divided by the
190
number of tanks in the battery and multiplied by αi to obtain the VOC emissions per tank. Since
9 ACS Paragon Plus Environment
Environmental Science & Technology
Page 10 of 28
191
individual VOCs were not available for tanks, the median of the log-transformed total VOC
192
emissions per tank were compared to the log-transformed total VOC emissions from the phase-
193
specific emissions experiments, yielding β5 = 0.2.
194
Model Validation
195
To validate our model, we utilized results from two studies that collected ambient air
196
samples in areas of dense O&G development and analyzed the samples for non-methane
197
hydrocarbons. The first was a Garfield County study that collected 24-hour average ambient air
198
samples within 350 feet in each cardinal direction of four well pads each undergoing drilling and
199
completions in 2008.35 This study also collected samples at eight “background” sites, which
200
given the high density of O&G development in that area can be considered representative of the
201
production sites in their vicinity. We supplemented the Garfield County data with results from
202
ambient air grab samples collected at nine locations across the DJB in 2014.36 We included
203
locations where samples were collected on multiple days and were not next to gas plants,
204
injection wells, or within 150 feet of O&G well equipment. The results for the sum of the eight
205
O&G pollutants collected at each site were averaged over the number of samples collected. We
206
used the model described in Eq. 1 to calculate an intensity-adjusted IDW (IA-IDW) at the 25
207
locations where ambient air samples were collected. The IA-IDW well count for each sampling
208
location was determined by identifying all O&G well pads within 10 miles of the sample;
209
dividing the corresponding intensity factors (ZL) by the square of the distance between the
210
sample and pad; and summing over the pads. The Spearman correlation was used to compare the
211
log-transformed IA-IDW to the log-transformed average sum of measured O&G pollutant
212
concentrations.
10 ACS Paragon Plus Environment
Page 11 of 28
213
Environmental Science & Technology
Model Implementation
214
The contribution of each model parameter was evaluated by setting the corresponding β
215
parameter to 0 to evaluate the effect of the variable’s absence on the distribution of the monthly
216
well pad intensity scores for 2004-2011. We modelled the intensity of O&G activities for all
217
O&G wells within 10 miles of residences in Greeley and calculated an IA-IDW for each
218
residence. We then compared the IA-IDW to an unadjusted IDW exposure metric. The
219
unadjusted IDW calculation is the same as the IA-IDW, but uses an intensity factor of 1 for each
220
active O&G well (i.e. a well pad with 10 wells would have ZL =10).The relative change in
221
variance was used for comparison as our goal is to account for differences in intensities.
222
223
Results Table 1 summarizes each parameter’s contribution to the intensity model. The hydraulic
224 225
fracturing/flowback term had the largest effect on the variability of model outputs due to the
226
magnitude of its intensity and its relatively short 14 day duration. The production stage variables
227
had the biggest impact on the central tendency of the distribution due to both the intensity and
228
long duration, which lasts years as opposed to months for the other stages. Drilling had a smaller
229
impact on the model due to its short duration, less frequent re-occurrence, and lower intensity of
230
activities compared to hydraulic fracturing/flowback.
231
232
Model
Mean
Median
Variance (σ2) of
Range of 11
ACS Paragon Plus Environment
Environmental Science & Technology
Intensity (Unit-less)
Intensity (Unit-less)
Page 12 of 28
Intensity (Unitless)
Intensity (Unit-less)
Full Model
6.33
4.20
225
0.003-1440
Pad Preparation (β0=0)
6.34
4.20
225
0.003-1440
Drilling (β1=0)
6.19
4.18
213
0.003-1440
Hydraulic Fracturing/Flowback (β2=0)
5.25
4.18
49.7
0.003-259
Oil Production (β3=0)
4.40
2.36
201
0.003-1440
Gas Production (β4=0)
4.90
2.84
253
0.032-1413
Tanks (β5=0)
6.28
4.19
223
0.003-1440
233 234 Table 1. Each parameter’s contribution to the industrial activity model outputs by setting that factor (β) to 0. 235
The distribution of intensity scores for active well pads is shown in Figure 1. This
236
distribution has two main components: 1) smaller producing well pads whose intensity scores
237
range from 0.003 to 10 (-5.8 to 2.3 when log-transformed); and 2) pads with wells being drilled,
238
hydraulically fractured/flowback, or large multi-well production pads whose intensity scores
239
range from 10 to 1440 (2.3 to 7.3 when log-transformed). The full model produced an overall
240
mean intensity of 6.33 with a variance of 225, and the mean monthly intensity generally
241
increased over time, corresponding to the increase in well starts and related activity.
12 ACS Paragon Plus Environment
Page 13 of 28
Environmental Science & Technology
242 243
Figure 1. The distribution of the log-transformed monthly modeled well pad intensity for O&G well pads in
244
Colorado from January 2004-December 2011. Intensity scores below 2.3 (10 on a linear scale) generally reflect
245
smaller pads in the production phase whereas scores above 2.3 represent large multi-well production pads or those
246
pads undergoing drilling or hydraulic fracturing/flowback.
247
Figure 2 shows an example of how the relative intensity at a well pad can change over
248
time using our model. The intensity profile reflects the transformation of a single well pad into a
249
large 21-well production pad. The initial intensity is based on the relatively small volume of
250
O&G produced at the single well site. There are 20 additional wells drilled over the period of 5
251
months and then these wells are all hydraulically fractured in the same month, leading to a spike
252
in intensity. Following the hydraulic fracturing/flowback period, the 21 wells go into production
253
and reflect an intensity of the increased production volume and additional equipment on site.
254
13 ACS Paragon Plus Environment
Environmental Science & Technology
Page 14 of 28
255 256
Figure 2. An example of the change in temporal intensity at a well pad in Garfield County. This well pad was a
257
single well production pad and transformed into a large 21-well production pad. Over the period of 5 months, these
258
20 additional wells were drilled and then they were hydraulically fractured in the same month, where the intensity
259
spiked.
260
As depicted in Figure 3, the Spearman correlation coefficient between the log-
261
transformed IA-IDW and log-transformed measured average sum of pollutant concentrations was
262
0.74 (p