Spatiotemporal Industrial Activity Model for Estimating the Intensity of

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


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


Page 3 of 28

37 38


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


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


Page 5 of 28


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


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


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)



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


Page 9 of 28


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



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.


Page 11 of 28

213


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


Intensity (Unit-less)


Page 12 of 28

Intensity (Unitless)


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.


Page 13 of 28


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



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

Spatiotemporal Industrial Activity Model for Estimating the Intensity of

Recommend Documents