University of Southern California Dissertations and Theses

Assessing induced seismicity rate increase based on deterministic and probabilistic modeling

(USC Thesis Other)

Assessing induced seismicity rate increase based on deterministic and probabilistic modeling

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content ASSESSING INDUCED SEISMICITY RATE INCREASE BASED ON
DETERMINISTIC AND PROBABILISTIC MODELING
by
Seyed Mehran Hosseini
A Ph.D. Dissertation Presented to the
FACULTY OF THE GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(PETROLEUM ENGINEERING)
May 2018
Copyright 2018 Seyed Mehran Hosseini
Abstract
The substantial increase in earthquake rates in the Central and Eastern United
States starting in 2009, has raised the question of a possible causal relationship
between wastewater disposal and earthquakes. This dissertation focuses on the
study of induced seismicity through probabilistic and deterministic modeling of
involved physical mechanisms and validating theoretically predicted seismicity rate
changes with observations. An accurate description of such rate changes is a key
component in seismic hazard assessment at the proximity of geothermal reservoirs,
carbon sequestration, hydraulic fracturing and wastewater disposal sites. This
dissertation meets the need for a representative model and validated simulator for
injection-induced seismicity.
Wecategorizecurrentinducedseismicitymodelsintotwobroadclassesofdeter-
ministic and probabilistic models. Deterministic models simulate pressure for fault
stability of known faults while probabilistic models focus on the probability of
inducing earthquakes on a random distribution of faults. Both classes of existing
models do not consider uncertainty in faults stresses and possible reservoir bound-
ary conditions due to different geological settings. In this dissertation, in addition
to deterministic induced seismicity modeling, we develop a probabilistic-physics-
based induced seismicity model that captures the effects of uncertainty in faults
stresses, reservoir properties and boundary conditions.
ii
We developed a probabilistic-physics-based model which considers a randomly
distributednetworkoffaultsthatundergofrictionalfailureduetopressurechanges.
We include uncertainty in faults stresses through parameterization of required
pressure increase for failure, based on probability distributions of tectonic stresses
using Mohr-Coulomb failure criterion. Reservoir boundary effects are explored
by including pressure solutions for early-time (i.e., infinite-acting) and late-time
(i.e., pseudo-steady-state) periods. We sample the probability density function
for required pressure-increase for failure through simulated pressure in the domain
and calculate the probability of inducing an earthquake at each point. To obtain
the total number of possible earthquakes, we integrate the calculated probability
over the study domain and multiply it by the fault density. Using frequency-
magnitudedistributionlawsforinducedearthquakes(e.g., Gutenberg-Richterlaw),
we calculate the cumulative number of earthquakes above a specific magnitude.
Using deterministic numerical simulation of pressure, we study several induced
seismicity cases in the U.S. and Europe.
Weshowthatthedifferencesinflowboundaryconditionscandrasticallychange
the seismic response to injection operations. We observe an exponential seismic-
ity rate increase in the case of a bounded reservoir versus a linear increase for
an unbounded reservoir. Our deterministic modeling results suggest that net
injected fluid volume and boundary conditions are the primary key factor, con-
trolling induced seismicity. We conclude that secondary and tertiary oil recovery
may experience less induced seismicity than wastewater disposal, likely as a result
of less net injected fluid volumes.
We observe that a linear relationship between injected volume and cumulative
number of seismic events may not always hold due to injection rate and flow
boundaryeffects. Injectionrateandflowboundaryeffectsmaybemorepronounced
iii
in long-term injections such as wastewater disposal whereas a linear behavior may
be expected in short-term, constant-rate injections such as hydraulic fracturing.
We validate our model against seismicity data from German KTB and Paradox
Valley saltwater disposal project, Colorado. The differences in flow boundary
conditions can drastically change the seismic response to injection operations and
for a bounded reservoir, induced seismicity rate increases exponentially.
Maximum magnitude induced seismicity (MMIS) is one of the critical param-
eters in earthquake hazard estimation. Using the developed induced seismicity
model in this dissertation, we study the relationship between operational injection
parameters and MMIS. Current MMIS models in the literature do not include the
effects of injection rates, permeability, and reservoir boundaries. We developed
a probabilistic-physics-based induced seismicity model for MMIS, which includes
the operational and reservoir parameters. Our model is a physics-probabilistic-
based model which accounts for initial state of stress, reservoir flow and injection
characteristics. We conclude that reservoir permeability, injection rate, and the
boundary conditions can drastically change expected MMIS, therefore some of the
conclusions and models in the literature based on very different cases of wastewa-
ter disposal with very different characteristics may not hold. Instead, we suggest
studying each case on its own due to differences in all the key parameters sur-
rounding MMIS.
We study induced seismicity and microseismicity due to hydraulic fracturing.
Besides the pressure increase due to diffusive processes, stress shadow around the
fracture is capable of changing the effective stress in favor of shear failure in pre-
existingfaults. Wedevelopedasimulatortopropagateahydraulicfracturethrough
a porous medium and capture shear failure among a random distribution of faults
in the domain to study statistics of induced seismicity associated with hydraulic
iv
fracturing. We observed very different stimulated rock volume (SRV) activated
as a result of same hydraulic fracture but with different reservoir geomechanical
conditions. Reservoir deviatoric stress, fracture pressure, existing faults and ini-
tial reservoir pore-pressure are among the key parameters pertaining to SRV and
observed microseismicity. The ideal combination for maximizing SRV and produc-
tive SRV is high fracture pressure, high initial reservoir pressure, high stress differ-
ence between minimum and maximum reservoir stress. One of the key completion
strategies should focus on creating larger SRVs with smaller fracture spacing to
maximize production from the pay zone.
v
Dedication
To my mother, to the memory of my father, and to my brothers, Mehrdad and
Mahmood
vi
Acknowledgment
I would like to express my sincere gratitude to Dr. Fred Aminzadeh for his
guidance, patience, and support throughout my studies at USC. He was always
ready to explore new ideas and helped me a lot with his knowledge and expertise.
I would like to thank Dr. Iraj Ershaghi, the director of the petroleum engi-
neering program at USC for his continuous support, encouragement, and great
technical discussions. His insight helped me in different aspects of this disserta-
tion.
MyspecialthanksgoestoDr. ThomasGoebelforhisgreatdiscussionsthrough-
out my studies at USC. His critical thinking and innovative questions improved
my research.
I am grateful and honored for receiving the Viterbi fellowship which paved my
way to focus on my research. I would like to express my gratitude for financial
support from the Induced Seismicity Consortium (ISC) and Reservoir Monitoring
Consortium (RMC) at USC.
I would like to thank all my Ph.D. committee members including Dr. Charles
Sammis, Dr. Birendra Jha, Dr. Aiichiro Nakano, and Dr. Katherine Shing for
their fruitful discussions.
I would like to thank Dr. Susan Hough, Dr. Nicholas van der Elst, Dr. Morgan
Page, and Dr. Deborah Weiser from USGS for their valuable discussions. I would
vii
liketothankDr. EgillHauksson, Dr. FredericCappa, andDr. JeanPaulAmpuero
from California Institute of Technology for their collaborations and discussions.
I also would like to thank the fellow graduate students and friends in the Mork
FamilyDepartmentofChemicalEngineering &MaterialsScience. Theirfriendship
made my time at USC even more pleasant.
Finally, I would like to thank my family, mother and brothers, especially
Mehrdad, for their support during my Ph.D. studies. This dissertation could not
be possible without them.
viii
Contents
Abstract ii
Dedication vi
Acknowledgment vii
List of Tables xii
List of Figures xiii
1 Introduction 1
1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Scope of Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2 Deterministic Approach to Induced Seismicity Modeling 14
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.2 Wastewater Disposal versus Oil Recovery . . . . . . . . . . . . . . . 16
2.3 The Effect of Permeability Structure; a Case Study from Illinois Basin 17
2.4 Induced Earthquake Numbers and Injected Volume . . . . . . . . . 21
2.5 Induced Seismicity at Tejon Oil Field, California . . . . . . . . . . . 23
2.5.1 Hydrological Model . . . . . . . . . . . . . . . . . . . . . . . 24
2.5.2 Model Parameters . . . . . . . . . . . . . . . . . . . . . . . 25
2.5.3 Initial and Boundary Conditions . . . . . . . . . . . . . . . . 27
2.5.4 Model Parameter Uncertainty . . . . . . . . . . . . . . . . . 28
2.5.5 Modeling Results . . . . . . . . . . . . . . . . . . . . . . . . 28
2.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3 A Probabilistic Approach to Injection-Induced Seismicity Assess-
ment for Different Reservoir Types and Pressure-Diffusion Mod-
els 31
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.2 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.2.1 Fault Distribution and Required Overpressure for Failure . . 34
ix
3.2.2 Overpressure Models . . . . . . . . . . . . . . . . . . . . . . 35
3.2.3 Induced Seismicity Probability and Cumulative Number of
Events . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.3.1 Paradox Valley 1996-2005 . . . . . . . . . . . . . . . . . . . 38
3.3.2 KTB 2004-2005 . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.3.3 Induced Seismicity in Oklahoma . . . . . . . . . . . . . . . . 43
3.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.5 Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4 ProbabilisticBoundsonMaximumMagnitudeInducedSeismicity
(MMIS) 54
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.2 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.2.1 Fault Distribution and Required Overpressure for Failure . . 58
4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.3.1 Permeability . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.3.2 Reservoir Size . . . . . . . . . . . . . . . . . . . . . . . . . . 62
4.3.3 Flow Rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
4.3.4 Net Injected Volume . . . . . . . . . . . . . . . . . . . . . . 64
4.3.5 Periodic Injections . . . . . . . . . . . . . . . . . . . . . . . 64
4.4 Possibility of Universal Maximum Magnitude Induced Seismicity
Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
4.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
5 AProbabilistic, GeomechanicsbasedInterpretationofStimulated
Rock Volume in Hydraulic Fracturing 70
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
5.2 Linear Elastic Fracture Mechanics . . . . . . . . . . . . . . . . . . . 75
5.2.1 Near-Tip Hydraulic Fracture Stress Solution . . . . . . . . . 76
5.2.2 Full-Field Hydraulic Fracture Stress Solution . . . . . . . . . 77
5.3 Shear Failure Criteria . . . . . . . . . . . . . . . . . . . . . . . . . . 80
5.4 Probability of an Event . . . . . . . . . . . . . . . . . . . . . . . . . 82
5.5 Sensitivity Analysis for SRV Shape and Size . . . . . . . . . . . . . 84
5.5.1 Case (1). High deviatoric stress and fracture pressure . . . . 85
5.5.2 Case (2). High deviatoric stress and low fracture pressure . . 87
5.5.3 Case (3). Low deviatoric stress and fracture pressure with
high initial reservoir pressure . . . . . . . . . . . . . . . . . 87
5.5.4 Case (4). Low deviatoric stress and fracture pressure with
low initial reservoir pressure . . . . . . . . . . . . . . . . . . 88
5.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
x
6 Conclusions and Future Research 91
6.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
6.2 Future Research . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
Reference List 96
xi
List of Tables
3.1 Nominal Paradox Valley Parameters . . . . . . . . . . . . . . . . . 47
3.2 Nominal KTB Parameters . . . . . . . . . . . . . . . . . . . . . . . 48
4.1 KTB Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
5.1 Nominal Hydraulic Fracture Simulation Parameters . . . . . . . . . 85
xii
List of Figures
1.1 Cumulative number of earthquakes with M≥ 3 in the Central and
Eastern United States, 1970-2017. The blue part of the curve cor-
responds to the long-term rate of 24 earthquake/year which sharply
increases after 2009 shown in red, with the corresponding average
rate of 330 earthquakes/year. (Inset) Distribution of epicenters in
the study area. Adapted from USGS. . . . . . . . . . . . . . . . . . 2
1.2 USGS map displaying the relative seismicity hazard in the United
States as a result of possible ground acceleration. Adapted from
Petersen et al. (2015). . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3 USGS map displaying the potential of damage from natural and
human-inducedearthquakesin2016. Oklahomashowsalargeincrease
in potential damage as high as California in 2016, which is unprece-
dented in past USGS hazard maps and Oklahoma seismicity history.
Adapted from Petersen et al. (2016). . . . . . . . . . . . . . . . . . 4
xiii
1.4 Location of active class II injection wells in the CEUS are shown
as blue circles and spatiotemporally associated injection wells, are
shown as yellow circles. Pie diagram shows spatiotemporally asso-
ciated injection wells by state. 40% of the associated injection wells
in the CEUS are located in Oklahoma while only 8% of all injection
wells are located in Oklahoma. Adapted from Weingarten et al.
(2015). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.5 locations of M≥ 0.0 earthquakes from 1 January 1973 through 31
December 2014 are shown with white dots for earthquakes that are
not spatiotemporally associated with injection wells and red dots for
the associated earthquakes. Adapted from Weingarten et al. (2015). 6
1.6 OklahomaSeismicityMap(M≥3),from1973toJuly2017. Adapted
from USGS. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.7 (Top) Cumulative number of magnitude 2.5 or greater earthquakes
in Oklahoma from 1997 to 2015. (Bottom) Letf axis shows cumu-
lative injected volume associated with three types of injection. The
right axis shows the magnitude of all Oklahoma earthquakes from
1997 to 2015. Adapted from Walsh and Zoback (2015). . . . . . . . 9
1.8 The spatial distribution of seismicity in the Paradox Valley area
over time: (a) injection tests, 1991-1995; (b) continuous injection,
1996-2000; (c) continuous injection, 2001-2008; and (d) continuous
injection, 2009 to February 2013. Adapted from King et al. (2014). 12
xiv
2.1 PermeabilitystructurefromSAIGUPProject(Matthewsetal. 2008)
(a). Randomuniformdistributionoffaults(b). pore-pressureincrease
due to wastewater disposal (c). pore-pressure increase due to simul-
taneousinjectionandproduction(d). Failedfaultsinredforwastew-
ater disposal (e), and simultaneous injection and production (f). . . 18
2.2 (Left) Geology, permeability and well structure used in the numer-
ical model. (Right) Boxplot of core permeability for several forma-
tions within the Illinois Basin after Zhang et al. (2013). The red
dots show he permeability magnitude used for the base case study. . 19
2.3 (Left)Pore-pressurecontourcorrespondingtoapore-pressureincrease
of 142 psi after ten years of injection with the injection rate of 5455
m
3
/d as a function of different injection layer widths (W). (Right)
Pore-pressure contour corresponding to a pore-pressure increase of
142 psi after ten years of injection with the injection rate of 5455
m
3
/d as a function of different permeability contrasts, (Kr). . . . . 20
xv
2.4 (a) Failed faults are shown in the case of constant injection after
three months, where the required pore-pressure for shear failure is
constant. The colors express the event magnitudes. The warmer
colors correspond to higher magnitude events. (b) Failed faults are
shown for the case of constant injection after six months of injec-
tion where the required pore-pressure for shear failure is constant.
(c) Failed faults are shown for the case of constant injection after
three months of injection where the required pore-pressure for shear
failure has a log-normal distribution. (d) Failed faults are shown
for the case of constant injection after six months of injection where
the required pore-pressure for shear failure has a log-normal distri-
bution. (e) Total number of failed faults plotted versus time for
the case of a constant required pore-pressure for shear failure. (f)
Total number of failed faults plotted versus time for the case of a
log-normally distributed required pore-pressure for shear failure. . . 22
2.5 (a) Seismicity and injection at Tejon Oil Field. (b) Seismicity (gray)
andinjectionrate(blue). (c)Seismicitycross-section. Adaptedfrom
Goebel et al. (2016). . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.6 (a)ModelSetupwithahigh-permeabilityflowpathway. (b)Increase
in Pore-pressure at the hypocenter of the M
W
4.6 Earthquake, 11
km away from the injection point. (c) Seismicity front, signaling a
diffusive process. Adapted from Goebel et al. (2016). . . . . . . . . . 29
xvi
3.1 (A,B) Locations of the induced seismicity cases including Paradox
ValleyinColorado1996-2005, thestateofOklahoma1997-2014, and
German KTB site 2004-2005. (C) Normalized cumulative number
of earthquakes versus normalized cumulative injected volume for
Oklahoma, KTB, and Paradox Valley shows that induced seismic-
ity behavior can not be described by a linear relationship between
injected volume and number of events. (D) Normalized injection
rates versus normalized time for Oklahoma, KTB, and Paradox Val-
ley. OverallinjectionrateincreasesinOklahomawhileinjectionrate
decreasesinParadoxValleyandremainsroughlyconstantattheKTB. 33
3.2 (A) Mohr-Coulomb failure envelope along initial and ultimate Mohr
circles in blue and red respectively. Normal distribution for S
1
and
S
3
are shown. Distribution for required pore-pressure for shear fail-
ure is shown in blue. (B) Number of events versus time in two
cases of bounded versus unbounded reservoir. (C) The probabil-
ity of inducing an earthquake is equal the area beneath its ΔP
p−req
distribution shown for far-wellbore and near-wellbore. (D) Possible
reservoir boundaries. . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.3 (left) Injection rate history from late-1996 to mid-2004 in Paradox
Valley, Colorado, along unbounded model pressure behavior and
calculated pressure data from Ake et al. (2005). (right) Cumulative
number of observed M≥ 0 seismic events from late-1996 to mid-
2004, model predictions for bounded, and unbounded reservoirs in
Paradox Valley, Colorado. . . . . . . . . . . . . . . . . . . . . . . . 40
xvii
3.4 (left) Cumulative injected volume at German KTB from 2004 to
2005 for 335 days, original, and simulated wellhead pressure data.
(right) Cumulative number of seismic events during 2004-2005 injec-
tion experiment at KTB, Germany, along developed bounded and
unbounded model predictions, original Segall and Lu (2015) model,
and bounded Segall and Lu (2015) model. . . . . . . . . . . . . . . 44
3.5 KTB and Paradox Valley parameters (A) Normal distribution of
KTB’s minimum total stress,S
3
. (B) KTB’s maximum total stress,
S
1
. Required pore-pressure for shear failure, ΔP
p−req
, for KTB (C)
and Paradox Valley (D). . . . . . . . . . . . . . . . . . . . . . . . . 49
3.6 Recordedandmodeledvaluesofpressurechangedividedbyflowrate
versus time during phase I (top) and II (bottom) cycles after July
10th, 1997 and July 26th, 1999. Derivative of the same parameters
is also plotted which shows the signature of a radial flow. Adapted
from King et al. (2016). . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.7 Recorded and modeled values of pressure change divided by flow
rate versus time during phase III (top) and IV (bottom) cycles after
January 8th, 2001 and January 14th, 2007. Derivative of the same
parametersisalsoplottedwhichshowsthesignatureofaradialflow.
The orange line shows the emergence of radial flow (Adapted from
King and Block, 2016). . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.8 Schematicdiagram of radial reservoir model and expected seismicity
in the case of bounded and unbounded reservoirs . . . . . . . . . . 53
xviii
4.1 (A) Mohr-Coulomb failure envelope along the initial and ultimate
Mohr circles in blue and red respectively. The normal distribution
forS
1
andS
3
are shown. The distribution for required pore-pressure
for shear failure is shown in blue. (B) The number of events versus
time in two cases of bounded versus unbounded reservoir. (C) The
probability of inducing an earthquake is equal the area beneath its
ΔP
p−req
distribution shown for far-wellbore and near-wellbore. (D)
The possible reservoir boundaries. . . . . . . . . . . . . . . . . . . . 59
4.2 Sensitivity analysis for the effect of permeability onN
total
and
ˆ
M
max
(left) linear log-log plot for the number of events greater thanM
0
=
2 versus the injected volume (right)
ˆ
M
max
plotted versus the loga-
rithm of the injected volume showing a linear behavioir. . . . . . . . 62
4.3 Sensitivity analysis for the effect of reservoir size onN
total
and
ˆ
M
max
(left)Non-linearlog-logplotfornumberofeventsgreaterthanM
0
=
2 versus injected volume (right)
ˆ
M
max
plotted versus log of injected
volume showing a non-linear behavior. . . . . . . . . . . . . . . . . 63
4.4 SensitivityanalysisfortheeffectofinjectionrateonN
total
and
ˆ
M
max
(left) linear log-log plot for number of events greater than M
0
= 2
versus injected volume (right)
ˆ
M
max
plotted versus log of injected
volume showing a linear behavior. . . . . . . . . . . . . . . . . . . . 64
xix
4.5 SensitivityanalysisfortheeffectofperiodicinjectiononN
total
,
ˆ
M
max
andP
p
(left) plot of pressure versus time for three cases of injections
with different flow rates of 159m
3
/day (2) periodic injection of 159
m
3
/day with the intervals of 10 days of injection and 10 days of
no injection and finally (3) reduced rate to half which equals 79.5
m
3
/day. (right) semi-log plot for number of events greater than
M
0
= 2 versus time on the right axis and
ˆ
M
max
plotted versus time
on the left axis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
5.1 Hydraulic fracturing in a horizontal wellbore in a shale layer. . . . . 72
5.2 Top: Interaction of hydraulic fracture and natural network of frac-
tures (from Huang et al. 2014). Bottom: SRV around hydraulic
fracture stages. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
5.3 (A) Mohr-Coulomb failure envelope along initial and ultimate Mohr
circles in blue and red respectively. Normal distribution for S
1
and
S
3
areshown. Distributionforrequiredproximitytofailureisshown
in blue. (B) The probability of inducing an earthquake is equal the
area beneath its ΔPF
req
distribution shown for far-wellbore and
near-wellbore. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
5.4 Fracture geometry and geometric parameters used in the near-tip
solution as well as stress sign convention. . . . . . . . . . . . . . . . 77
5.5 Process zone shown in red and singularity-dominated zone shown in
green, at the crack tip. . . . . . . . . . . . . . . . . . . . . . . . . . 78
xx
5.6 Maximum, S
1
and minimum, S
3
principal stress concentration at
the crack-tip shown on left and right respectively for a fracture
half-length of 3 ft., S
1
of 2000 psi and S
3
of 1000 psi with frac-
ture pressure of 1100 psi. Cooler colors at the tip of the fracture
show a tensile regime while surrounding the face of the fracture,
compressive stress regime in hot colors are dominant. . . . . . . . . 78
5.7 Fracture geometry and geometric parameters used in the full-field
stress solution for hydraulic fracturing. . . . . . . . . . . . . . . . . 79
5.8 A typical shear failure zone at the tip of fracture for a case of S
1
of
8040 psi, S
3
of 1740 psi, fracture pressure of 2240 psi and fracture
half-length of 100 ft. The size of shear failure zone for a constant
failure envelope in a 100 ft fracture is roughly 2 ft wide and 8 ft long. 80
5.9 (A) Mohr-Coulomb failure envelope along initial and ultimate Mohr
circlesinblueandredrespectively. NormaldistributionforS
1−Initial
andS
3−Initial
are shown. Distribution for required pore-pressure for
shear failure is shown in blue. (B) The probability of inducing an
earthquake is equal the area beneath its ΔPF
req
distribution shown
for far-wellbore and near-wellbore. (D) Possible reservoir boundaries. 83
5.10 Case (1). High deviatoric stress and fracture pressure. The change
in proximity to failure is plotted. Warmer colors refer to higher
proximity to failure resulting in SRV while cooler colored areas are
less prone to shear failure and SRV. . . . . . . . . . . . . . . . . . . 86
5.11 Case (2). High deviatoric stress and low fracture pressure. The
change in proximity to failure in plotted. Warmer colors refer to
higher proximity to failure resulting in SRV while cooler colored
areas are less prone to shear failure and SRV. . . . . . . . . . . . . 88
xxi
5.12 Case (3). Low deviatoric stress and fracture pressure with high ini-
tial reservoir pressure. The change in proximity to failure in plotted.
Warmer colors refer to higher proximity to failure resulting in SRV
while cooler colored areas are less prone to shear failure and SRV. . 89
5.13 Case (4). Low deviatoric stress and fracture pressure with a low ini-
tial reservoir pressure. The change in proximity to failure in plotted.
Warmer colors refer to higher proximity to failure resulting in SRV
while cooler colored areas are less prone to shear failure and SRV. . 90
xxii
Chapter 1
Introduction
1.1 Background
Active tectonic regions such as plate boundaries in the western United States
are usually the expected regions for earthquakes. However, earthquakes report-
edly occur throughout the United States with rates significantly below the plate
boundary seismicity rates (Petersen et al., 2015). The presence of faults at the
verge of failure makes it possible for small stress perturbation to trigger earth-
quakes in these regions (Townend and Zoback, 2000; Wesson and Nicholson, 1987;
Report, 2013; Schultz et al., 2018). What makes intraplate seismicity in the U.S.
standout, is the dramatic increase over the past nine years from 2009 (Figure 1.1).
The cumulative number of earthquakes with M≥ 3 in CEUS from 1970 to 2016
are shown in the Figure 1.1. The long-term rate of 24 earthquake/year shown
in blue, sharply increases after 2009 to the average rate of 330 earthquakes/year
shown in red. A 14 fold increase in the rate of earthquakes in CEUS after 2009 is
observed compared to the 1973-2008 period.
States experiencing elevated seismicity levels include Arkansas, California, Col-
orado, Kansas, New Mexico, Ohio, Oklahoma, Texas, and Virginia (Ellsworth,
2013; Goebel et al., 2016; Rubinstein and Mahani, 2015). The probability of miss-
ing M≥ 3 earthquakes in the Central and Eastern United States has been near
zero for decades which rules out the possibility of the recorded increase in seismic-
ity rate being an artifact of seismic network changes (Commission et al., 2012).
1
Figure 1.1: Cumulative number of earthquakes with M≥ 3 in the Central and
Eastern United States, 1970-2017. The blue part of the curve corresponds to the
long-term rate of 24 earthquake/year which sharply increases after 2009 shown in
red, with the corresponding average rate of 330 earthquakes/year. (Inset) Distri-
bution of epicenters in the study area. Adapted from USGS.
According to the U.S. Geological Survey (USGS) “The USGS National Seismic
Hazard Maps display earthquake ground motions for various probability levels
across the United States and are applied in seismic provisions of building codes,
insurance rate structures, risk assessments, and other public policy.” (Figure 1.2).
A greater number of large, damaging earthquakes might result from the increase
in the rate of earthquakes in CEUS which consequently persuaded USGS scientists
to update the 2014 USGS seismic hazard map (Figure 1.2 Petersen et al. (2015))
to the 2016 version of the one-year seismic hazard forecast (Figure 1.3) for CEUS
that includes contributions from both induced and natural earthquakes (Petersen
2
et al., 2016). Comparison between the seismic hazard map and the damage map
in the Figures 1.2 and 1.3, respectively, shows that the potential of damage has
increased drastically in Oklahoma due to the inclusion of induced seismicity.
Figure 1.2: USGS map displaying the relative seismicity hazard in the United
States as a result of possible ground acceleration. Adapted from Petersen et al.
(2015).
Some of the CEUS earthquakes happened in the areas where oil and gas oper-
ations were highly active which motivated scientists to look at the correlation
between oil and gas activities and induced earthquakes (Report, 2013). The U.S.
shale revolution was not possible without massive hydraulic fracturing treatments
which use a massive amount of water. Flow-back water together with produced
water after hydraulic fracturing is saline and re-injected into the deeper high per-
meability layers. Class II wells are designed and regulated under Environmental
3
Figure1.3: USGSmapdisplayingthepotentialofdamagefromnaturalandhuman-
inducedearthquakesin2016. Oklahomashowsalargeincreaseinpotentialdamage
as high as California in 2016, which is unprecedented in past USGS hazard maps
and Oklahoma seismicity history. Adapted from Petersen et al. (2016).
Protection Agency (EPA) to regulate the injection of fluids associated with oil and
gas production. These wells consist of saltwater disposal, enhanced recovery, and
a small number of hydrocarbon storage wells. More than 180,000 Class II wells
are operating in the U.S. and over 2 billion gallons of brine per day were injected.
Weingarten et al. (2015) examined the spatiotemporal correlation between the
earthquakesandinjectionwellsacrosstheCEUS.Theydefinedassociatedinjection
wells as the active injection wells at the time of earthquakes within 15 km radial
distance. They found that 10% of all the active injection wells in CEUS are
associatedwithearthquakeswhichismuchhigherthanthereportedcasesinCEUS.
Associated injection wells are mainly located in Texas and Oklahoma shown in
4
the Figure 1.4 with the corresponding associated and non-associated earthquakes
shown in the Figure 1.5. Out of 187,570 Class II injection wells in CEUS, 56%
are active, and approximately 75% of the active Class II injection wells consist of
enhanced oil recovery wells (EOR) while about 25% are saltwater disposal wells
(Weingarten et al., 2015).
Figure 1.4: Location of active class II injection wells in the CEUS are shown as
blue circles and spatiotemporally associated injection wells, are shown as yellow
circles. Pie diagram shows spatiotemporally associated injection wells by state.
40% of the associated injection wells in the CEUS are located in Oklahoma while
only 8% of all injection wells are located in Oklahoma. Adapted from Weingarten
et al. (2015).
We present a brief review of some of the main reported induced earthquakes
in the U.S. The M
W
4.8 earthquake on 9 August 1967 near Denver, Colorado was
the largest induced earthquake in the U.S. by wastewater disposal (Healy et al.,
5
Figure 1.5: locations of M≥ 0.0 earthquakes from 1 January 1973 through 31
December 2014 are shown with white dots for earthquakes that are not spatiotem-
porally associated with injection wells and red dots for the associated earthquakes.
Adapted from Weingarten et al. (2015).
1968) before the M
W
5.7 Prague, Oklahoma earthquake on November 2011 (Ker-
anen et al., 2013). The M
W
4.0 30 December 2011 earthquake was reported near
Youngstown, Ohio, which might have been induced by wastewater disposal (Kim,
2013). The M
W
4.7 Guy-Greenbrier, Arkansas earthquake on 27 February 2011
was reported to be related to wastewater disposal associated with Fayetteville
shale production (Horton, 2012). The M
W
5.7 6 November 2011 Prague, Okla-
homaearthquakehashappenedintheproximityofseveralhigh-volumewastewater
injection wells (Keranen et al., 2013), TheM
W
4.8 17 May 2012 East Texas earth-
quake near Timpson is the largest recorded earthquake in Eastern Texas and is
known to be associated with an active wastewater injection well (Frohlich et al.,
6
2014). A swarm of earthquake activity has been recorded from November 2013
through January 2014, near Azle and Reno, Texas with two M
W
3.6 events that
were widely felt. Geomechanical studies showed that the swarm activity might be
related to brine production and wastewater disposal near existing faults (Hornbach
et al., 2015). In October 2014, two M
W
4.0 and M
W
4.3 earthquakes happened
near Cushing, Oklahoma close to the largest crude oil facility in the world (McNa-
mara et al., 2015). Based on different studies it is becoming clear that main hazard
of induced seismicity is due to wastewater disposal into deep formations adjacent
to crystalline basement or layers with hydraulic communications with crystalline
basement. The stress change due to such activities in these environments can
transfer to crystalline basement faults which have the possibility of generating
large magnitude earthquakes. A map of earthquakes in Oklahoma from 1973 to
July 2017 is shown in Figure 1.6.
Most of the induced earthquakes in CEUS were located in Oklahoma. The top
portion of the Figure 1.7 shows the cumulative number of M≥ 2 earthquakes in
Oklahoma (Walsh and Zoback, 2015). After 2009 a sharp increase in the number
of earthquakes is observed. The seismicity rate increases each year after 2009 and
the overall behavior of the cumulative number of events versus time is superlin-
ear. Figure 1.7 on the bottom shows the statewide monthly injection volumes for
Oklahoma from approximately 7000 UIC wells. Associated injection rates for three
types of enhanced oil recovery (EOR) wells, saltwater disposal wells (SWD) and
wells of unknown type are shown in the Figure 1.7 at the bottom. The earth-
quakes times and magnitudes are shown with red dots at the bottom part of the
Figure 1.7. The monthly injection rate in Oklahoma almost doubled from 80 mil-
lion barrels/month to 160 million barrels/month. This is mostly due to SWD wells
(Weingarten et al., 2015), which inject into the Arbuckle formation at the vicinity
7
Figure 1.6: Oklahoma Seismicity Map (M≥ 3), from 1973 to July 2017. Adapted
from USGS.
of the crystalline basement (Keranen et al., 2013). Although the injection rate has
almost linearly has doubled from 1997 to 2014, the corresponding seismic activity
abruptlyincreasesin2009withalineartrenduntil2013andseismicityacceleration
after 2013 in a non-linear fashion is observed.
We present three cases on large-scale, long-term, deep wastewater disposal
projects in the U.S. and their associated seismicity. These three cases are widely
discussedintheliteratureandsignificantlycontributedtoourunderstandingabout
induced seismicity. In 1961, a deep injection well was used by U.S. Army at the
RockyMountainArsenal(RMA)northeastofDenver, Coloradotoinjecthazardous
chemical waste (Healy et al., 1968). The injection at the depth of 3.6 km started on
March 1962, and the project was terminated on February 1966, after the residents
8
Figure 1.7: (Top) Cumulative number of magnitude 2.5 or greater earthquakes in
Oklahoma from 1997 to 2015. (Bottom) Letf axis shows cumulative injected vol-
ume associated with three types of injection. The right axis shows the magnitude
of all Oklahoma earthquakes from 1997 to 2015. Adapted from Walsh and Zoback
(2015).
of the northeastern Denver area began to report unprecedented earthquakes. 13
earthquakes with M
b
4 and greater occurred during the operational life of the
project (Ellsworth, 2013). The next year after the termination of the project,
the largest Denver earthquake occurred with M
W
4.8 event on 9 August 1967
approximately 10 km away from the injector (Herrmann et al., 1981). RMA events
highlight the possible pore-pressure transfer through faults to distances far away
from the injection point as well as the possible long delays in the order of months
or years between the injection operations and earthquakes.
9
AftertheRockyMountainArsenalearthquakeobservations, suggestinginduced
earthquakes due to pressure transfer in the nearby faults, USGS scientists began
an experiment in 1969 to test this hypothesis in the Rangely oil field in Northern
Colorado. Chevron Oil Company, the field operator, gave the USGS permission
to regulate pressure in a part of the field which was seismically active. The field
had undergone water flooding since 1957. Using the measurements for fault fric-
tion factor from reservoir rock core samples and the in-situ state of stress, the
critical pore-pressure increase of 25.7 Mpa was calculated to induce earthquakes
(Ellsworth, 2013). Between 1969 to 1973, two cycles of fluid injection and with-
drawal were performed. It was observed that when the fluid pressure in the moni-
toring well exceeded the critical pressure, earthquake activity increased. One day
after the onset of the fluid withdrawal, seismic activity ceased (Ellsworth, 2013).
This experiment showed the importance of pressure and its transmission mecha-
nism to faults in inducing earthquakes.
The Colorado River provides water to about 27 million people and irrigation
water to nearly four million acres of land in the United States. Dolores River, a
tributary of the Colorado River, is a major contributor to the natural salt load in
Colorado River which is picked up from Paradox Valley at the rate of 20000 kilo
tonnes of salt per year (Ake et al., 2005). The Paradox Valley Unit (PVU), is a
salinity control project located in Southwest Colorado, which disposes of brine in
a deep injection well to reduce the amount of highly saline groundwater entering
the Dolores River. The PVU project started large scale injections in 1996 injecting
14,000 to16,000 feet deepinto Precambrianand Paleozoic rocks (Yeck et al., 2015).
Until now more than around eight million cubic meters of fluid have been injected
through this project (Yeck et al., 2015). A seismic network installed in 1985 has
been dedicated to seismicity monitoring in the area. AM
W
3.8 in 2000, aM
W
3.6
10
in 2004, and a M
W
4 earthquake in 2013 are among the largest events that have
happened at the vicinity of the Paradox Valley area (Yeck et al., 2015). Between
1985andJuly1996, only12tectonicearthquakesweredetectedwithin35kmofthe
well. This number jumped to thousands as the injection started (Ellsworth, 2013;
Ake et al., 2005). The Figure 1.8 shows the spatial distribution of seismicity in
the Paradox Valley area over time (King et al., 2014). Through rate and injection
volumecontrol, PVUwasabletoreduceassociatedinducedseismicityearthquakes.
We will discuss the details of PVU operation and history in Chapter 3. PVU is an
excellent example of how rate and volume of wastewater disposal can be managed
to decrease induced seismicity.
1.2 Scope of Work
In Chapter 2, we present a deterministic approach for induced seismicity mod-
eling, which is based on the knowledge of the state of stress, fault friction factor,
and fault trajectory, the amount of pore-pressure increase required for shear fail-
ure is calculated. The deterministic approach is suitable for fault stability analysis
where a known fault is subject to stress changes. We study induced seismicity in
two field cases in the Tejon oil field, California, and the Illinois Basin. Simulating
the pressure increase due to injection, we examine the effect of permeability struc-
ture, operational parameters, and reservoir geometry. We observe that injection
layerwidthandpermeability-contrastrelativetothecrystallinebasementinfluence
the induced seismicity potential. Employing the deterministic approach, we study
the difference between induced seismicity response in wastewater disposal versus
secondary, and tertiary oil recovery. Our results suggest that net injected fluid vol-
ume (i.e., the difference between injected and produced volumes) is one of the key
11
Figure 1.8: The spatial distribution of seismicity in the Paradox Valley area over
time: (a) injection tests, 1991-1995; (b) continuous injection, 1996-2000; (c) con-
tinuous injection, 2001-2008; and (d) continuous injection, 2009 to February 2013.
Adapted from King et al. (2014).
factors controlling induced seismicity potential. We observe that secondary and
tertiary oil recovery, as well as, hydraulic fracturing, have less induced seismicity
potential than wastewater disposal as a result of less net injected fluid volumes.
In Chapter 3, we present a model for the probabilistic assessment of induced
earthquakes that incorporates geomechanical analysis of the state of stress, pore
pressure changes and reservoir boundary conditions. The probabilistic approach is
12
bettersuitedtoaddresstheoverallinducedseismicitybehaviorofafieldandtocap-
ture induced seismicity statistics. The geomechanical model considers a randomly
distributed network of faults that undergo frictional failure due to overpressure-
induced effective stress changes. The model includes the effect of flow properties
and flow boundary conditions on overpressure, as well as the effect of the initial
stress regime on failure. Theoretical predictions are compared to observations of
induced seismicity in Colorado, and the German Continental Deep Drilling Pro-
gram KTB. As expected, differences in flow boundary conditions can drastically
change the seismic response to injection operations. The strongest difference is
observed for a bounded reservoir for which the induced seismicity rate increases
exponentially. The proposed model may provide a pathway to assess the expected
seismic response to fluid injection as a function of operational parameters and
reservoir boundary conditions.
In Chapter 4, we discuss the highlights of the work and discuss future research
plans.
13
Chapter 2
Deterministic Approach to
Induced Seismicity Modeling
2.1 Introduction
The increase in wastewater disposal associated with hydraulic fracturing, and
othersubsurfacefluidinjectionandproductionoperationsaswellastheheightened
public concern regarding possible induced seismicity, motivated us to study the
related issues. In this work, we investigate induced seismicity caused by pressure
changes due to fluid injection. We use a reservoir simulator to model the pressure
change in a realistic case with different injection and production scenarios. This
research shows that the presence of critical faults and the net injected volume
are among the most important risk factors contributing to induced seismicity.
Reservoir characteristics such as geometry, size, and permeability are identified as
the main components that pertain to reservoir pressure response, causing a critical
pressure increase and a subsequent potential seismicity. We apply our method to
the Central Illinois Basin, which is a primary candidate for CO
2
sequestration.
We observe that the pressure changes differ widely, but can easily lead to fault
instability and seismic activity up to 10 km away from the injection well.
The magnitude of pressure increase and the geometry of the affected area
greatly depends upon the permeability structure of the injection formation and
14
surroundings. Further, we define the required pressure buildup for fault instabil-
ity. In addition, we generate a random distribution of faults in the area of interest
where they have a certain distribution for required pore-pressure buildup for fault
slip (criticality). Earthquake magnitudes are sampled from a distribution following
the Gutenberg-Richter law. We introduce a simple homogeneous, isotropic repre-
sentation for the reservoir in order to study the primary parameters that affect
the frequency and magnitude of induced seismic events. Given the absence of a
physics-based method for induced seismicity risk assessment, this chapter proposes
a novel approach for risk assessment. The results can be used to mitigate the fault
reactivation risk.
We investigate the potential for induced seismicity through modeling pore-
pressure change. First, we study simultaneous injection and production, which is
common in water floods, and then injection with no production, which is common
in wastewater disposal (WWD). We show that net fluid volume (i.e., the differ-
ence between injected and produced volumes) is a significant factor in controlling
potential seismicity. Our results suggest that enhanced oil recovery and hydraulic
fracturing have less induced seismicity potential than wastewater disposal as a
result of less net injected fluid volume. This finding is consistent with a National
Research Council report on Induced Seismicity Potential in Energy Technologies
(2013). Second,weexaminetheeffectofpermeabilitystructureandreservoirgeom-
etry in the Central Illinois Basin. We run a suite of models with changing the (a)
injection layer width and (b) injection layer permeability contrast relative to the
crystalline basement. Both of the aforementioned influence the induced seismicity
potential (orders of magnitude change in the extent of the elevated pore-pressure
volume is observed). Finally, we investigate the effect of net injected volume of the
total number of induced earthquakes. We compare two cases of constant injection
15
after three and six months. One of the important aspects of induced seismicity
is the relationship between the number of induced earthquakes and injected fluid
volume. We organize two sets of models where faults are distributed randomly,
and their critical pore-pressure distribution is uniform and log-normal. We dis-
cuss the different flow regimes that the reservoir might experience by presenting
corresponding equations. We observe that independent from the distribution of
required pore-pressure for shear failure with the increase in injected volume, there
will be an increase in both the total number of induced seismic activity and its
extent.
2.2 Wastewater Disposal versus Oil Recovery
Induced seismicity can result from different types of subsurface fluid injection
and production (SFIP) operations and hydraulic fracturing to different degrees.
In this section, we show results from two series of model runs, which explored
pressure changes caused by wastewater disposal versus that of water flooding. We
use a synthetic reservoir model from the Sensitivity Analysis of the Impact of
Geological Uncertainties on Production (SAIGUP) Project (Krogstad et al., 2015).
MATLABReservoirSimulatorToolbox(MRST)isutilizedforpressuresimulations
(Krogstad et al., 2015). Inset (a) of the Figure 2.1 shows the permeability structure
in the model. The synthetic permeability field that is generated includes faults and
natural features typically found in shoreface reservoirs. We uniformly distribute
some faults in the domain as shown in Figure 2.1, Inset (b). We simulate pressure
increase in two scenarios: (1) wastewater disposal with no production and (2)
water flooding where the injection rate equals the production rate. One observes
16
the pressure increase associated with each scenario in Insets (c) and (d) of Figure
2.1, respectively.
As shown, the amount of pressure increase in the case of injection with no
production is much higher than that of simultaneous injection and production. In
themiddleofFigure2.1, theleftsiderepresentsawastewaterdisposalcase, andthe
right side symbolizes a typical case of water flooding. The pressures are very high,
but we only focus on the difference between the two cases. In Insets (e) and (f), red
dots show the failed faults due to injection and the blue dots denote the faults that
have not been failed. In Inset (e), the failed and stable faults are illustrated for
the water flooding case. Due to smaller pressure increase and smaller affected area
with pressure increase, in the case of waterflooding, fewer faults are activated while
more faults are activated due to wastewater disposal in Inset (f). This illustrates
the importance of net injected volume on induced seismicity potential.
Several studies in the literature focus solely on injected volumes, whereas the
produced volumes should be taken into account to have a realistic picture of the
induced seismicity potential. Reservoirs that are already depleted might have an
entirely different seismic response to injection than the reservoirs that are not
depleted. In the next section, we analyze the effect of reservoir permeability struc-
ture through a case study from the Central Illinois Basin.
2.3 The Effect of Permeability Structure; a Case
Study from Illinois Basin
A suite of models were developed to investigate the effect of permeability struc-
ture on injection induced seismicity. We examine a case of fluid disposal in prox-
imity to active faults in the Central Illinois Basin (Zhang et al., 2013). The basal
17
Figure 2.1: Permeability structure from SAIGUP Project (Matthews et al. 2008)
(a). Random uniform distribution of faults (b). pore-pressure increase due to
wastewater disposal (c). pore-pressure increase due to simultaneous injection and
production (d). Failed faults in red for wastewater disposal (e), and simultaneous
injection and production (f).
Mount Simon Sandstone is considered a target for CO
2
Sequestration (Person
et al., 2010). We use the simulation parameters employed by Zhang et al. (2013),
found in Table 1 of their paper. We then add a sensitivity analysis for the effect of
injection layer width (W), and the logarithm of the ratio of permeability contrast
between the injection layer (A1) and the crystalline basement (CB), defined as Kr.
18
The model’s geometry, geology, and layers are shown in Figure 2 on the left.
The name of the layer and the range of the reported permeability in the literature
corresponding to each layer in the Illinois Basin is shown in the boxplot on the
right (after Zhang et al. (2013)). The red dots show the permeability magnitude
used for the base case study.
Figure 2.2: (Left) Geology, permeability and well structure used in the numerical
model. (Right) Boxplot of core permeability for several formations within the Illi-
nois Basin after Zhang et al. (2013). The red dots show he permeability magnitude
used for the base case study.
We investigate the effect of the injection layer width (W) and the permeabil-
ity contrast (Kr) defined as log10(KA1/ KCB) on the spatial-temporal pressure
response to injection. Figure 3 shows the pore-pressure increase front of 142 psi
after ten years of injection with 5455 m
3
/d as a function of different injection
layer widths, W, (left) and permeability contrast, Kr (right). One sees that as the
width of the injection layer increases from 85 m to 850 m, the lateral extent of
19
pressure decreases. Hence, the geometry of the injection layer would be an impor-
tant parameter that can substantially affect the pressure increase at the crystalline
basement where most earthquakes are generated.
In the definition of permeability contrast (Kr), we assume that the permeabil-
ity of the injection layer is constant. We reduce the permeability of the crystalline
basement to get a different Kr. As Kr increases from 1 to 6, we observe that the
pressure increase is more confined in the injection layer due to the low permeability
of the crystalline basement. Therefore, the existence of the low permeability crys-
talline basement ensures that pressure increase is confined to the injection layer,
which might reduce the chances of induced seismicity, as most earthquakes happen
in the crystalline basement.
Figure 2.3: (Left) Pore-pressure contour corresponding to a pore-pressure increase
of 142 psi after ten years of injection with the injection rate of 5455 m
3
/d as
a function of different injection layer widths (W). (Right) Pore-pressure contour
corresponding to a pore-pressure increaseof 142 psi after ten yearsof injection with
the injection rate of 5455 m
3
/d as a function of different permeability contrasts,
(Kr).
20
2.4 Induced Earthquake Numbers and Injected
Volume
One crucial aspect of induced seismicity is the relationship between the number
of induced earthquakes and the injected fluid volume. We run two sets of models
wherein the faults are distributed randomly. Their critical pressure distribution is
uniform and log-normal. At the onset of injection, we have infinite acting behavior.
We start with the flow equations in radial flow in the infinite-acting reservoir.
Infinite acting means that the boundaries are located infinitely distant from the
injection point. Carslaw and Jaeger (1959) solved governing PDEs. We adapt
the formulation presented by Lee et al. (2003). This is based on the Carslaw and
Jaeger (1959) in the field units, presented in Appendix 1.
The results of the two model runs are presented below, and the event magni-
tudes are color coded. Figure 4, Insets (a), and (b), show the failed faults in the
case of constant injection after three and six months, respectively. The required
pore-pressure increase for shear failure is constant in these two cases. The colors
depict the event magnitudes, with the warmer colors corresponding to higher mag-
nitude events. Insets at top right and bottom left show the failed faults in the case
of a constant injection after three and six months, with a log-normally distributed
required pressure increase for shear failure. One perceives that irrespective of the
distribution for required pressure for shear failure with the increase in injected
volume, the total number of induced seismic activity, as well as its extent, both
increase. The induced seismicity pattern in the case of constant pressure required
for shear failure is radial, while for the case of log-normally distributed criticality,
we see in Insets (c), and (d), is distributed. The frequency-magnitude distribution
follows the GutenbergâĂŞRichter law and each event magnitude is sampled from
21
this distribution. Finally, we plot the total number of failed faults versus time in
the case of constant and log-normally distributed criticality, as shown in Figure 4,
Insets (e), and (f), respectively.
Figure 2.4: (a) Failed faults are shown in the case of constant injection after three
months, where the required pore-pressure for shear failure is constant. The colors
express the event magnitudes. The warmer colors correspond to higher magnitude
events. (b) Failed faults are shown for the case of constant injection after six
months of injection where the required pore-pressure for shear failure is constant.
(c) Failed faults are shown for the case of constant injection after three months
of injection where the required pore-pressure for shear failure has a log-normal
distribution. (d) Failed faults are shown for the case of constant injection after
six months of injection where the required pore-pressure for shear failure has a
log-normal distribution. (e) Total number of failed faults plotted versus time for
the case of a constant required pore-pressure for shear failure. (f) Total number of
failed faults plotted versus time for the case of a log-normally distributed required
pore-pressure for shear failure.
22
2.5 Induced Seismicity at Tejon Oil Field, Cali-
fornia
1
In order to assess the effects of subsurface waste water disposal on pressure
change and possible induced seismicity at Tejon oil field, we developed a 3D pres-
sure model for the injection area. Tejon oil field is located at the Southern part
of San Joaquin Basin Province, California (Hosford Scheirer, 2007) which is com-
posed of three main parts, namely, the Western Tejon, the Central Tejon and the
Eastern Tejon (California oil and gas fields, 1998). The history of the Tejon oil
field production dates back to 1935 (Kasline, 1948) and as of now 589 wells still
exist in this field (DOGGR Well Finder, 2015). There are several injection and
production wells at the Western Tejon oil field which inject and produce mainly
from three layers of Transition, Pulv (Reserve) and Valv with average depths of
792, 1402, and 1645 m, respectively (California Oil and Gas Fields, 1998). Due to
proximity of WWD-32 waste water disposal well to the hypocenter of earthquakes
and high injection rate prior to and at the time of earthquakes, we aim to model
pore-pressure change associated with injection in this particular well.
Waste water disposal well, WWD-32 with API. No 030-26630 is a horizontal
well located in Section 32, T11N, R19W which kicks-off from vertical at 551 m
depthandhorizontallyextendsintheTransitionlayerforabout762mtowardseast
at the depth of 813 m (DOGGR Well Finder, 2015). The target layer for injection
is the Transition layer with 30 m thickness. Based on gamma ray, resistivity logs,
and geologic layering interpretation for the reservoir (California oil and gas fields,
1998), we model the reservoir with four distinct layers of combined Alluvium %
1
Parts of this chapter have been published in Goebel et al. (2016).
23
Kern River-Chanac (undifferentiated), Transition, Santa Margarita and Basement
with heights of 800, 30, 300 m and 17.9 km respectively (California oil and gas
fields DOGGR report, 1998). Figure 5, adapted from Goebel et al. (2016), shows
the earthquake locations, magnitudes, injection rates, and seismicity cross-section.
Figure 2.5: (a) Seismicity and injection at Tejon Oil Field. (b) Seismicity (gray)
and injection rate (blue). (c) Seismicity cross-section. Adapted from Goebel et al.
(2016).
2.5.1 Hydrological Model
Thepartialdifferentialequation(PDE)describingtheflowissolvednumerically
using a finite-difference method by ECLIPSE reservoir simulator from Schlum-
berger. ECLIPSE is a reservoir simulator, using finite difference method to solve
fluid flow problems. ECLIPSE 100, or black-oil simulator is a single-phase reser-
voir simulator that is used in this dissertation. In order to solve the PDE, we first
discretize the domain. Then we assign required properties such as porosity, perme-
ability and compressibility to each of the grid-blocks. Next, we assign boundary
24
and initial conditions. We also introduce wells, their location and their injection
schedule. Finally, by assigning fluid rheology and reporting time-steps, the model
setup is completed.
We initially discretize the domain with grid sizes of Dx = 109.75, Dy = 500 and
Dz = 100 meters. At the vicinity of the injection point and the high permeability
pathway, we use a local grid refinement to refine the grid size gradually down to
5.78, 9.55 and 1.58 m in X, Y and Z directions respectively. To eliminate boundary
effects, we locate the injection well 7.8 and 2.5 km away from reservoir boundaries
in X and Y directions at the depth of 0.8 km. For modeling purposes we model
the well as a vertical injector and complete the well at the 30 meters interval of
the transition zone.
2.5.2 Model Parameters
Permeability
Zhang et al. (2013), among many others, show that having a thin low per-
meability layer between two high permeability layers can isolate them from each
other and there would be no pressure transmission between the two layers. In other
words, a low permeability layer can dictate the behavior of the whole domain by
isolating different high permeability layers from each other. Our goal is to model
the reservoir behavior in large scale where there is a lack of detailed description for
every part of the domain. We focus on detailed modeling for the high permeability
Transition layer and high permeability pathways in the reservoir. Encountering a
lowpermeabilitylayer, weassignedthelowpermeabilitytothewholelayersbeyond
that point where the pressure transmission is dependent upon passing through the
low permeability layer. Townend and Zoback (2000) showed that critically stressed
faults are mostly permeable and will experience shear deformation that can result
25
in seismic activity. Following the premise of Shapiro et al. (1997) and others, we
use spatio-temporal distribution of seismic activities to find traces of permeability
structures and their magnitudes.
Due to a high clay content observed in upper Kern River-Chanac and lower
Santa Margarita layers from logs adjacent to the Transition layer, we assign a
typical low permeability of 0.001 mD to these two layers. There might be other
sublayers in these main layers with higher permeability that is assigned to the
whole layer but as we mentioned before, overall model response is dictated by those
low permeability layers adjacent to the Transition layer. Following this logic, for
basement permeability, we assign a constant value of 0.001 mD, for the whole 17.9
km interval.
Porosity
We Compared the porosities, Φ, reported in California oil and gas fields report
(1998) for Transition (0.29), Reserve (0.28), and Valv (0.25 to 0.30) layers respec-
tively. We observe that although these layers are spanning average depths of 792,
1402, and 1645 m, they all have very similar porosities. Therefore, with the lack
of further data for porosity in other layers and our observation regarding the small
change of porosity among the three layers, we assign a uniform porosity of 0.2
to the whole domain. Core and seismic data from the industry can improve this
uniform porosity assumption.
Transition Layer Extent
Transition layer’s 2.1 km
2
horizontal extent is determined from approximate
coverage of Western Tejon reservoir map. This is consistent with reservoir extent
26
for Reserve and Valv layers and well locations in Western Tejon (California oil and
gas fields DOGGR report, 1998).
Rock and Fluid Compressibility
Reservoir rock compressibility, c, is one of the parameters with the least mea-
surements in our study domain. In the absence of measurements for rock com-
pressibility we use the correlation proposed by Hall (1953), who found a correlation
between porosity and rock compressibility for different sandstones and limestones.
BasedonHall’stotalrockcompressibilityplot, byincreasingtherockporositymore
than 13%, rock compressibility curve plateaus to 4× 10
−10
pa
−1
. We use this
magnitude for rock compressibility. Core measurements can enhance the accuracy
of this estimation.
Fluid Properties
In the absence of fluid sample data and since the injected fluid is water, we
assign typical water properties for the injected fluid. The water density is 1000
kg/m
3
at the surface condition. The water compressibility at 30 Mpa (reference
pressure) is 1.1× 10
−9
pa
−1
. The water formation volume factor at the reference
pressure is 1.13 (reservoir m
3
/ standard condition m
3
). The water viscosity at
reference pressure is 0.5 cP.
2.5.3 Initial and Boundary Conditions
Since we study the increase in pressure, not the absolute pressure, the initial
pressure condition in the reservoir may not influence the results. However, to solve
the PDE, the initial pressure for each grid-block is required. In the absence of
initial-pressure data in the study domain, we assumed a pressure gradient of 12000
27
pa/m throughout the domain and pressure of zero at the surface (top boundary).
We impose Neumann boundary conditions to the PDE, where the flow rate at
the domain boundaries are zero (no-flux boundary condition). In order to avoid
artificial effects of boundary domains on the results, we distanced the injection well
from the domain boundaries and chose the size of the domain to be three times the
distance of well and mainshock hypocenter in X direction and two times the depth
of the mainshock in Z direction. The model size is 30 km in horizontal direction,
18 km in vertical direction with a 5 km width.
2.5.4 Model Parameter Uncertainty
The 3D reservoir model size is 30 km in horizontal direction, 18 km in verti-
cal direction with 5 km width. We are relying on public data for geologic layer
interpretations and properties for this large-scale domain. There can always be
improvements to the model by adding 3D, and 4D seismic data, as well as tracer
test data to better understand the geology and permeability structure of the reser-
voir. Typical buildup test, interference and pulse well tests can also be useful,
especially since the field is mature and well developed.
2.5.5 Modeling Results
Model structure, (Figure 6(A)), The amount of pressure increase at the
hypocenter of the M
W
4.6 earthquake as a function of different permeability val-
ues and pathway widths (Figure 6(B)), and the seismicity front (Figure 6(C)) are
shown (Goebel et al., 2016).
28
Figure 2.6: (a) Model Setup with a high-permeability flow pathway. (b) Increase
in Pore-pressure at the hypocenter of the M
W
4.6 Earthquake, 11 km away from
the injection point. (c) Seismicity front, signaling a diffusive process. Adapted
from Goebel et al. (2016).
2.6 Conclusions
Numerical modeling of simultaneous injection and production showed that the
net volume of injected or produced fluid is an key parameter in induced seismicity
occurrence. In the Central Illinois Basin case, we observe that the size of the
reservoir, presence of faults, and permeability contrast relative to the host rock,
strongly influence the pressure changes with distance and time, which may act as
29
a trigger for induced seismicity. These pressure changes fluctuate widely and can
easily lead to fault instability and seismic activity up to a 10 km distance from the
injection well. Finally, by modeling a constant injection scenario with the presence
of randomly distributed faults, we witness that as the injected volume increases,
the extent and total number of induced earthquakes increases. We also observe
that depending upon the criticality distribution for faults, the spatio-temporal
distribution of the seismic cloud would be different. The knowledge obtained from
this chapter can be useful in the selection of wastewater disposal site, injection
schedule, and hazard studies.
30
Chapter 3
A Probabilistic Approach to
Injection-Induced Seismicity
Assessment for Different
Reservoir Types and
Pressure-Diffusion Models
1
3.1 Introduction
It has long been known that injection of fluids into the subsurface might cause
earthquakes(Evans,1966). EarthquakestriggeredatRockyMountainArsenaland
Rangely, Colorado are classic examples of fluid induced seismicity (Healy et al.,
1968; Herrmann et al., 1981). The recent increase in seismicity rate from the
background rate in the Central U.S. starting in 2009, brought much attention to
inducedseismicityassociatedwithwastewaterdisposal(Healy et al.,1968;Frohlich,
2012; Horton, 2012; Ellsworth, 2013; Keranen et al., 2013; Kim, 2013; Zhang et al.,
2013). In California, where the background seismicity rate is high, a new case of
1
Parts of this chapter is under publication at the journal of Geophysical Research Letters.
31
induced seismicity associated with non-geothermal fluid injection was reported
(Goebel et al., 2016) which expands the induced seismicity prone area from U.S.
midcontinent to the western part of the country.
The relationship between injection rate or injected volume and the rate of
inducedseismicityisobservedtobecomplexwithsignificantdeviationfromalinear
relationship. Some examples are shown in Figure 1 where the cumulative number
of earthquakes versus injected volume and injection rate versus time are plotted
for three different datasets: Oklahoma 1997–2014 (Walsh and Zoback, 2015; Wein-
garten et al., 2015), KTB 2004–2005 injection experiment (Shapiro et al., 2006a),
and Paradox Valley 1996–2004 saltwater disposal (Ake et al., 2005; King et al.,
2014). Comparing the number of earthquakes versus injected volume for differ-
ent sites (Figure 1C) suggests that a linear relationship between the cumulative
number of events and injected volume may not always be applicable. Namely, the
Oklahoma and KTB datasets show an acceleration in the seismicity rate which is
difficult to explain by standard models.
Here, we explore plausible explanations for the non-linearity observed in the
relationship between the cumulative number of seismic events and injected vol-
ume. We introduce a probabilistic induced seismicity framework that integrates
the physical processes of pressure diffusion in the reservoir and frictional failure
along faults with the initial state of stress and reservoir boundary conditions.
First, we calculate the probability distribution of pore pressure increase required
to induce frictional failure along faults given the distributions of principal stresses.
Next, we quantify expected seismicity rate changes in a randomly distributed fault
system that undergoes frictional failure due to injection-induced effective stress
changes. We specifically test the consequences of different reservoir boundary con-
ditions and injection settings. We validate our model using the data from the
32
KTB 2004-2005 injection experiment (Shapiro et al., 2006a) and Paradox Valley
1996-2004 saltwater disposal (Ake et al., 2005; King et al., 2014). We propose that
the acceleration in seismicity rate at the KTB and in Oklahoma can partially be
explained by the effects of ‘closed’ (no flow) reservoir boundaries in addition to
changes in injection rates. We conclude that a physics-based, probabilistic frame-
work of induced seismicity can successfully capture the site-specific relationships
among injection parameters, reservoir parameters, and expected seismicity rates.
Figure 3.1: (A,B) Locations of the induced seismicity cases including Paradox
Valley in Colorado 1996-2005, the state of Oklahoma 1997-2014, and German
KTB site 2004-2005. (C) Normalized cumulative number of earthquakes versus
normalized cumulative injected volume for Oklahoma, KTB, and Paradox Valley
showsthatinducedseismicitybehaviorcannotbedescribedbyalinearrelationship
between injected volume and number of events. (D) Normalized injection rates
versus normalized time for Oklahoma, KTB, and Paradox Valley. Overall injection
rate increases in Oklahoma while injection rate decreases in Paradox Valley and
remains roughly constant at the KTB.
33
3.2 Method
We derive a mathematical model to predict the number of seismic events
induced by injection-related overpressure. There are four key elements of the
model: (1) A random spatial distribution of faults in the reservoir, (2) the dis-
tribution of pore pressure change required for failure on these faults (3) the over-
pressure caused by fluid injection and (4) the resulting probability of inducing an
earthquake at each point in the domain, which is used to compute the total number
of possible induced events. Below we describe these elements of the model.
3.2.1 Fault Distribution and Required Overpressure for
Failure
In our model, shear failure is induced by effective stress changes on faults as a
result of changes in reservoir pressure. Based on the Mohr-Coulomb failure criteria
forcriticallystressedfaultsatthevergeoffailure, followingZoback (2010)notation,
the amount of required overpressure for shear failure, ΔP
p−req
, is calculated as:
ΔP
p−req
=
nS
3
−S
1
n− 1
−P
p0
(3.1)
where n = (
q
1 +μ
2
f
+μ
f
)
2
, μ
f
is fault coefficient of friction, P
p0
is initial pore-
pressure, S
1
and S
3
are maximum and minimum total principal stresses, respec-
tively. ΔP
p−req
, is a linear combination of S
1
, S
3
, and P
p0
. Zoback (2010) among
others show that stress distribution in fault networks can be described by a normal
distribution. In other words, the principal stresses, S
3
and S
1
, that drive failure
on individual faults within the network of reservoir-bounded faults are assumed to
follow a normal distribution. In this case, based on probability theory, the linear
34
combination of several normal distributions remains normally distributed. There-
fore, given a normal distribution for maximum and minimum stresses, ΔP
p−req
can
also be described by a normal distribution. In Figure 2 (A), we show the distribu-
tions of S
1
, S
3
, ΔP
p−req
along with the Mohr circles before and after injection.
3.2.2 Overpressure Models
Here, we discuss the expected overpressure ΔP
p
based on models with differ-
ent boundary conditions. The modeling results allow us to distinguish induced
seismicity in a bounded reservoir from that in an unbounded reservoir. We use
a line-source pressure solution with two different boundary conditions to approxi-
mate injection from a single well in a homogeneous, isotropic, cylindrical reservoir.
Schematics of a cylindrical reservoir with spatially random fault distribution, an
injector at the center and corresponding pressure profile is shown in Figure 2 (B).
The details of the pressure model can be found in the supporting information for
unbounded reservoirs with infinite acting flow behavior, and bounded reservoirs
with pseudo-steady flow behavior (Lee et al., 2003). Due to the large spatial scale
ofmostwastewaterdisposalcasesandtheextentofunknowngeologicalparameters
involved, the model introduced in this work does not include crustal heterogeneity
but has the advantage of capturing average large-scale pressure changes in very
different environments. It should be mentioned that more detailed pore-pressure
modeling approaches can be used within the framework introduced here, but they
are outside the scope of this work.
The corresponding number of induced seismic events versus time for two cases
of identical constant injection rate with different bounded and unbounded reser-
voir boundary conditions are illustrated in Figure 2 (B). At early injection times,
the presence of outer reservoir boundaries has no effect on the resolved pressure
35
changes. As injection continues, the reservoir pressure starts to be influenced by
the boundaries and the flow is governed by pseudo-steady flow behavior. Figure
2 (B) shows that for identical injection rates, a bounded reservoir will eventually
experience an acceleration in seismicity while both reservoirs share the same initial
behavior. The observed higher number of induced earthquakes in Figure 2 (B) is
due to higher reservoir pressure increase in the case of a bounded reservoir com-
pared to an unbounded reservoir after the common initial infinite-acting period.
Sealing faults, pinch-outs and low permeability layers such as shale are among pos-
sible reservoir boundaries that strongly inhibit fluid flow (Manzocchi et al., 2008),
while the absence of these elements create unbounded reservoirs where the pressure
behavior is described as infinite-acting (Figure 2, D).
3.2.3 Induced Seismicity Probability and Cumulative
Number of Events
Based on probability theory, one can find the probability of activating pre-
stressed faults due to overpressure, as the integral probability density function
(PDF) of ΔP
p−req
, from zero to ΔP
p
, which is defined as the difference in the
cumulative density functions (CDF) at those points. In this study, we are only con-
cerned with positive values of ΔP
p−req
, which corresponds to pore pressure increase
due to injection, therefore we only focus on the positive part of the ΔP
p−req
dis-
tribution from zero to ΔP
p
. We refer to this integral as G(ΔP
p
), which presents
the probability of activating a fault due to a pore-pressure increase of ΔP
p
and is
shown in Figure 2 (C) as the area beneath the probability curve. Specifically the
integral can be expressed by:
G(ΔP
p
) =
Z
ΔPp
0
f(ΔP
p−req
)dΔP
p
=F (ΔP
p
)−F (0) (3.2)
36
wheref andF are the PDF, and the CDF, respectively. The probability of induc-
ing a seismic event is smaller farther away from the wellbore due to smaller induced
overpressure, ΔP
p
, corresponding to a smaller area beneath the PDF (Figure 2.,
C, left), and larger closer to the wellbore due to the higher ΔP
p
(Figure 2. C.
right). It should be noted that the required pore-pressure change for shear failure
can have any arbitrary probability density function and Equation (2) holds true
for any PDF.
Figure 3.2: (A) Mohr-Coulomb failure envelope along initial and ultimate Mohr
circles in blue and red respectively. Normal distribution for S
1
and S
3
are shown.
Distribution for required pore-pressure for shear failure is shown in blue. (B)
Number of events versus time in two cases of bounded versus unbounded reser-
voir. (C) The probability of inducing an earthquake is equal the area beneath
its ΔP
p−req
distribution shown for far-wellbore and near-wellbore. (D) Possible
reservoir boundaries.
37
So far, we calculated the probability of inducing an earthquake at each point in
the domain. In order to calculate the cumulative number of induced seismic events
as a function of injection, we define ρ
(x,y,z)
as the density of faults in the domain.
MultiplyingG(ΔP
p
)byρ
(x,y,z)
andintegratingoverthevolumeofthedomain, gives
an estimate of the total number of probable seismic events. To get an insight into
the number of events greater than a certain magnitude, we assume a Gutenberg-
Richter-type frequency-magnitude relationship for induced earthquakes with an
exponent, b∼ 1. We define R(M
0
) as N
M(M>M
0
)/N
Total
, where R represents the
frequency of events having a magnitude greater than or equal M, N
M>M
0
is the
number of events greater than magnitude M
0
, and N
Total
is the total number of
events. To find N(M
0
) (i.e., total number of events with magnitude larger than
M
0
) one can multiply R(M) and N
Total
resulting in the following equation:
N
M>M
0
(t) =


Z
V
G(ΔP
p
(x,y,z)
)ρ
(x,y,z)
dv


R(M
0
) (3.3)
3.3 Results
WeapplyourmodeltotheParadoxValley1996–2005saltwaterdisposalproject
and the KTB 2004–2005 induced seismicity experiment.
3.3.1 Paradox Valley 1996-2005
The Paradox Valley Unit (PVU), is a salinity control project located in South-
westColorado, whichdisposesofbrineinadeepinjectionwelltoreducetheamount
of highly saline groundwater entering the Dolores River. The PVU project started
large scale injections in 1996. It injects 14,000 to 16,000 ft deep into Precambrian
and Paleozoic rocks (Yeck et al., 2015). Until now, approximately 8 million cubic
38
meters of fluid has been injected through this project (Yeck et al., 2015). A seismic
network installed in 1985 has been dedicated to monitoring seismicity in the area.
Between 1985 and July 1996, only 12 tectonic earthquakes were detected within
35 km of the well while this number increases to thousands as the injection starts
(Ellsworth, 2013; Ake et al., 2005). A M
W
3.8 in 2000, a M
W
3.6 in 2004, and
a M
W
4 earthquake in 2013 are among the largest events that have happened
in the vicinity of Paradox Valley area. A series of injection tests for reservoir
characterizationwereperformedbetween1991-1995andafterthisperiod, thebrine
injection started in 1996 (Ake et al., 2005). In 1999, earthquake rates increased to
its maximum level and several M≥ 3 earthquakes persuaded scientists to introduce
a new procedure in 2000 that involved periodic 20 day shutdowns and a 33%
reduction in the injection rates; Employing these measures after 2000 reduced
earthquake rates drastically (Figure 3, left). We use the reservoir parameters
reported in the literature for Paradox Valley, presented in the Appendix (see Table
A1 and Figure A1), to model Paradox Valley injection as a radial, unbounded
reservoir (Ake et al., 2005).
The Cumulative number of M≥ 0 earthquakes in Paradox Valley versus time is
plotted in Figure 3 (right). Initially, the seismicity rate is high but decreases with
the start of the 33% rate reduction and 20 days shut off practices start around
2000. Our "Unbounded Geomechanical Probabilistic Induced Seismicity Model"
(Unbounded GPIS Model), precisely captures this behavior (Figure 3, right). In
the case of model prediction for a "Bounded Geomechanical Probabilistic Induced
Seismicity Model" (Bounded GPIS Model), we observe that the induced number
of earthquakes is much higher and does not fit the observed trend in Paradox
Valley. Unbounded reservoir behavior is consistent with previous geological models
(Bremkamp and Harr, 1988), and studies for PVU (King et al., 2014; Block et al.,
39
2015) which shows confinement from northern and southern parts, while reservoir
isolation is not observed on the sides. The modeled pressure behavior shown in
Figure 3 (left), is consistent with inferred downhole pressure data from Ake et al.
(2005).
Figure 3.3: (left) Injection rate history from late-1996 to mid-2004 in Paradox Val-
ley, Colorado, along unbounded model pressure behavior and calculated pressure
data from Ake et al. (2005). (right) Cumulative number of observed M≥ 0 seismic
events from late-1996 to mid-2004, model predictions for bounded, and unbounded
reservoirs in Paradox Valley, Colorado.
3.3.2 KTB 2004-2005
The German continental deep drilling program (KTB) is located in southeast-
ern Germany near several tectonic units (Emmermann and Lauterjung, 1997). A
pilot and main well, reaching as deep as 4 km and 9.1 km deep have been drilled on
40
the site. Three main fluid injection and production experiments were performed
in the field between 1994 and 2005. In this study, we focus on the last series of
induced seismicity experiments between 2002 and 2005 (Shapiro et al., 2006a).
Initially, between 2002 and 2003, a total volume of∼ 22300 m
3
of brine was
produced from the open-hole section of the KTB pilot well between 3850-4000 m
depth. During the production phase of the project, seismicity was absent. Start-
ing in June 2004, over a period of 11 months,∼ 84600 m
3
of water was injected
into the pilot well. The injection rate was almost constant with an average rate
of 200 l/min and 10 MPa of wellhead pressure (Figure 4, left). We simulated
expected wellhead pressure due to injection and as shown in Figure 4 (left) there
is a good match between the simulated results and wellhead pressure data which
shows that the corresponding parameters for pressure simulation are representa-
tive. We define net injected volume, as the injection volume subtracted by 22300
m
3
produced brine from the reservoir in 2003. We observe that∼ 110 days from
the onset of injection, net injected volume reaches zero. During the first 110 days
of injection, the previously produced saltwater from the reservoir is restored, and
the net positive injection starts afterward. During the first 110 days of the reser-
voir replenishment period, no seismicity is recorded whereas after the positive net
injection starts (after 110 days), seismicity is observed (Figure 4).
The pilot wellbore at the KTB site is at∼ 1.6 km distance to the projected
positionoftheFranconianlineamentwhichisconsistentwiththecalculatedbound-
ary distance from the wellbore based on the well-test results of Stober and Bucher
(2005). It has been shown in the literature that Franconian lineament is imper-
meable (Stober and Bucher, 2005) and may act as a sealing reservoir boundary.
A permeability estimate of 2× 10
−15
m
2
and a reservoir radius of∼ 1.6 km, both
based on pressure transient analysis were calculated (Stober and Bucher, 2005).
41
We use the reservoir parameters reported in the literature for KTB, presented in
the supporting information (see Table S2 and Figure S1) to model KTB injection
as a radial, bounded, cylindrical reservoir (Shapiro et al., 2006a).
We investigate the non-linearity observed in KTB’s induced seismicity response
(Figure 4 right) and try to capture the physical mechanism behind it through our
inducedseismicitymodel. Alineartrendinthecumulativenumberofeventsversus
time is observed initially which changes after∼ 210 days of injection highlight-
ing an acceleration in the seismicity rate. Using our framework, we develop the
KTB induced seismicity model based on a bounded reservoir model that honors
previous pressure transient testing results (Stober and Bucher, 2005) and geology
of the KTB site. Bounded GPIS Model in Figure 4 (right), captures the induced
seismicity behavior observed in the field. The early-time linear growth in the num-
ber of events corresponding to an infinite-acting pressure behavior is followed by a
non-linear growth in the number of events when the overpressure front senses the
reservoir boundaries and the bounded reservoir behavior begins. For comparison,
we present the results from our unbounded induced seismicity simulation denoted
as Unbounded GPIS Model, which clearly shows that the unbounded model cannot
capture the induced seismicity behavior at KTB.
To further document the importance of reservoir boundary conditions, we com-
pare our results to predicted induced seismicity rates from one of the previous
works. We examine the performance of Segall and Lu (2015) induced seismic-
ity model in capturing KTB 2004-2005 seismicity. The Segall and Lu (2015)
model computes seismicity rate as a result of injection in a homogeneous, poroelas-
tic medium, considering full poroelastic coupling and time-dependent earthquake
nucleation. The Rudnicki (1986) poroelastic solution for a point source injection
42
along seismicity rate model relating changes in coulomb stress to changes in seis-
micity rate, are the basis of this model. As is shown in Figure 4 (right) for KTB,
the originalSegall and Lu (2015) model, exceptduring thevery earlytime, behaves
linearly and cannot capture the induced seismicity behavior. The original Segall
and Lu (2015) pore-pressure model is based on a 3D unbounded reservoir model,
We changed their original 3D unbounded, homogeneous, and isotropic pressure
model to a 2D radial, homogeneous, isotropic model based on Carslaw and Jaeger
(1959)andLee et al.(2003), whichisreferredtoas"Bounded Segall and Lu (2015)"
in the following. As is shown in Figure 4 (right), the Bounded Segall and Lu (2015)
model demonstrates a significant improvement in the prediction of the number of
events at KTB.
3.3.3 Induced Seismicity in Oklahoma
Monthly injection volume in the state of Oklahoma, almost doubled from about
80 million barrels/month in 1997 to about 160 million barrels/month in 2013
(Walsh and Zoback, 2015) while the seismicity increased abruptly in 2009 (Fig-
ure 1). This highly non-linear relationship observed between injection rate and
cumulative number of M≥ 2.5 earthquakes in Oklahoma may not be explained by
linear induced seismicity models. Our induced seismicity model introduced in this
work is capable of capturing non-linear induced seismicity behaviors such as KTB
and Paradox Valley, Colorado. The general non-linear trend in observed cumula-
tive induced seismicity in Oklahoma might be justified based on a combination of
injection rate and boundary conditions. Due to the complexity and extended spa-
tial scale of wastewater disposal in Oklahoma, we will reserve the discussions about
Oklahoma for a separate future publication. Intensified seismogenic consequences
due to the effect of bounded injection reservoirs has been hypothesized before for
43
the 2011 Prague earthquake sequence (Keranen et al., 2013) and may explain sig-
nificant time delays between injection activity and seismicity onset, followed by a
rapid acceleration in seismicity rates.
Figure 3.4: (left) Cumulative injected volume at German KTB from 2004 to 2005
for 335 days, original, and simulated wellhead pressure data. (right) Cumulative
numberofseismiceventsduring2004-2005injectionexperimentatKTB,Germany,
along developed bounded and unbounded model predictions, original Segall and
Lu (2015) model, and bounded Segall and Lu (2015) model.
3.4 Conclusions
We developed a model for the probabilistic assessment of induced seismicity
which incorporates geomechanical analysis of state of stress, pore pressure changes
and reservoir boundary conditions. Deviations from linear relationship between
injected volume and cumulative number of seismic events, are expected as a result
of changing injection rates and reservoir boundary conditions. Injection rate and
flow boundary effects may be more pronounced during long-term injections such
as wastewater disposal whereas an approximate linear behavior may be expected
44
during short-term, constant-rate injections such as hydraulic fracturing. We val-
idate our model against seismicity data from the German KTB and the Paradox
Valley saltwater disposal projects. We observe an acceleration in seismicity rate
at KTB despite an almost constant injection rate which can be explained by the
effects of reservoir boundaries on pressure and expected seismicity. Through our
model, we capture the deceleration in induced seismicity in Paradox Valley and
associate it with 33% reduction in injection rate as well as occasional shut-ins of 20
day intervals. The differences in flow boundary conditions can drastically change
the seismic response to injection operations and for a bounded reservoir, induced
seismicity rate increases almost exponentially. Through simulation of different
outcomes corresponding to different injection policies, our proposed model may be
used as a tool for wastewater disposal hazard assessment.
3.5 Appendix
The mathematical solution for pressure increase in an Infinite Acting Reser-
voir, the case in which the well is assumed to be placed in a porous medium of
infinite radial extent, is presented in Equation A3.1 (Lee et al., 2003). Due to
the extended geology of permeable layers in Paradox Valley, specially expanding
in North-West and South-East directions (Yeck et al., 2015), along the extended
pattern of observed induced seismicity in these directions (Bremkamp and Harr,
1988), we describe Paradox Valley injections by an infinite acting reservoir model.
ΔP
p
(r,t) =
qBη
4πkh
Ei

−
r
2
φηc
t
4kt
!
(A3.1)
45
f(ΔP
p−req
) =
1
q
2πσ
2
ΔP
p−req
e
−(ΔP
p−req
−μ
ΔP
p−req
)
2
2σ
2
ΔP
p−req
(A3.2)
F (ΔP
p
) =
1
2
(1 +erf(
ΔP
p
−μ
ΔP
p−req
√
2σ
ΔP
p−req
)) (A3.3)
G(ΔP
p
) =
1
2
(erf(
ΔP
p
−μ
ΔP
p−req
√
2σ
ΔP
p−req
−erf(
0−μ
ΔP
p−req
√
2σ
ΔP
p−req
)) (A3.4)
To get an insight into the number of events greater than a certain magnitude,
we assume a Gutenberg-Richter type frequency-magnitude relationship for induced
earthquakes with an exponent, b∼ 1. We define R(M
0
) as N
M(M>M
0
)/N
Total
,
where R represents the frequency of events having a magnitude greater than or
equal M, N
M>M
0
is the number of events greater than magnitude M
0
, and N
Total
is the total number of events. In the case of a Gutenberg-Richter frequency-
magnitude relationship for induced earthquakes:
R(M
0
) = 10
−bM
0
(A3.5)
where b commonly known as b-value is a constant. N
M>M
0
represents the
number of earthquakes with magnitudes greater than or equal to M
0
. N
M>M
0
formulation is general 3D representations; any specific geometry, analytical or
numerical solution can be used for pore-pressure.
N
M>M
0
(t) =


Z
V
1
2
(erf(
ΔP
p
−μ
ΔP
p−req
√
2σ
ΔP
p−req
−erf(
0−μ
ΔP
p−req
√
2σ
ΔP
p−req
))ρ
(x,y,z)
dv


10
−bM
0
(A3.6)
46
Table 3.1: Nominal Paradox Valley Parameters
Parameter and Variable Value and Unit
Permeability (k) 4.5× 10
−15
m
2
Fluid viscosity (η) 4× 10
−4
Pas
Total compressibility (C
t
) 7.25× 10
−10
Pa
−1
Direct velocity strengthening magnitude (a) 0.003
Nominal friction (μ) 0.7
Effective normal stress (¯ σ) 16.67 MPa
Background stressing rate ( ˙ τ
0
) 0.001 MPa/yr
Height (h) 280 m
Porosity (φ) 0.05
Hypothetical external radius (r
e
) 10000 m
Wellbore radius (r
w
) 0.152 m
Mean required pressure increase for shear failure (μ
ΔP
p−req
) 3.2 MPa
Standard deviation in required pressure for failure (σ
ΔP
p−req
) 2 MPa
Initial reservoir pressure (P
initial
) 43.64 MPa
where ΔP
p
is the increase in fluid pressure, t is time, r is the radius from
the injection point, k is formation permeability, φ is formation porosity, η is fluid
viscosity,c
t
is total compressibility of the porous medium,h is reservoir height, and
B is formation volume factor for injected wastewater assumed to be one. Table
A3.1 shows the input parameters for pore-pressure simulation used in Equation
A3.1 for Paradox Valley, Colorado.
The approximate mathematical solution for pressure increase in Bounded
Cylindrical Reservoir - the case in which well is assumed to be placed in the cen-
ter of a cylindrical area with no flow across the exterior boundary, is presented in
EquationA3.7(Lee et al.,2003). Duetotheboundedgeologyofpermeablegeologic
structure at KTB, bounded by Franconian Lineament known to be impermeable
(Stober and Bucher, 2005), we describe KTB injections by a bounded cylindrical
reservoir model.
47
Table 3.2: Nominal KTB Parameters
Parameter and Variable Value and Unit
Permeability (k) 5× 10
−15
m
2
Fluid viscosity (η) 10
−3
Pas
Injection rate (q) 3.02× 10
−3
m
3
/s
Total compressibility (C
t
) 7.25× 10
−10
Pa
−1
Direct velocity strengthening magnitude (a) 0.003
Nominal friction (μ) 0.65
Effective normal stress (¯ σ) 16.67 MPa
Background stressing rate ( ˙ τ
0
) 0.001 MPa/yr
Height (h) 150 m
Porosity (φ) 0.04
External reservoir radius (r
e
) 1650 m
Wellbore radius (r
w
) 0.152 m
Mean required pressure increase for shear failure (μ
ΔP
p−req
) 1.44 MPa
Standard deviation in required pressure for failure (σ
ΔP
p−req
) 0.48 MPa
Minimum horizontal stress (Sh
min
) 80 MPa
Maximum horizontal stress (SH
Max
) 174.35 MPa
Standard deviation in minimum horizontal stress (σ
Sh
min
) 0.3 MPa
Standard deviation in maximum horizontal stress (σ
SH
Max
) 0.5 MPa
ΔP
p
(r,t) =
qBη
2πkh
"
2
r
2
eD
− 1

r
2
D
4
+t
D
!
−
r
2
eD
ln(r
D
)
r
2
eD
− 1
−
3r
4
eD
− 4r
4
eD
ln(r
D
)− 2r
2
eD
− 1
4(r
2
eD
− 1)
2
#
(A3.7)
where r
eD
=r/rw, r
D
=re/rw, t
D
=kt/(φηc
t
r
2
w
), rw is well radius and re is
external radius. All the other parameters are defined in Equation A3.1.
Table A3.2 shows the input parameters for pore-pressure simulation used in
Equation A3.7 for KTB, Germany.
Figure A3.1 shows the normal distribution for KTB’s minimum total stress,
S
3
, KTB’s maximum total stress, S
1
, and required pore-pressure for shear failure,
ΔP
p−req
for KTB and Paradox Valley. S
1
and S
3
distributions are based on the
reported range for stress at KTB. ΔP
p−req
distribution for KTB is calculated based
48
on Equation 1 in the main text. ΔP
p−req
distribution for Paradox Valley is based
on the reported ranges of pressure increase in Paradox Valley (Ake et al., 2005).
Figure3.5: KTBandParadoxValleyparameters(A)NormaldistributionofKTB’s
minimum total stress, S
3
. (B) KTB’s maximum total stress, S
1
. Required pore-
pressure for shear failure, ΔP
p−req
, for KTB (C) and Paradox Valley (D).
49
In the case of KTB, the injection depth is 13123 ft and the reported pressure is
wellhead pressure which is 1450 psi (Shapiro et al., 2006a). The downhole pressure
including a hydrostatic gradient would be 7355 psi which yields a gradient of 0.56
psi/ft. Brudy et al. (1997) report a gradient of 0.88 psi/ft for minimum stress
at KTB which would be the fracture gradient. Since injection gradient is lower
than fracture gradient for KTB, we would not expect the reservoir to experience
fracturing. Therefore, we use a homogeneous radial pressure model to capture the
induced seismicity response.
In Paradox Valley, the injection depth is 14107 ft and downhole pressure of
11603 psi is reported (Ake et al., 2005). The resulting injection gradient is 0.822
psi/ftwhichisgreaterthan0.69psi/ftreportedfracturegradientinParadoxValley
(Ake et al., 2005). Therefore, we conclude that fracturing might have been possible
in Paradox Valley wastewater disposal.
King et al. (2016) study the far-field reservoir pressurization in Paradox Valley,
Colorado. They identify 4 main injection and shut-in cycles during the injections
from 1993 to 2016, namely phase I to IV. Pressure transient testing was used by
utilizingF.A.S.T.WellTestsoftware, WellTest toanalyzethepressuredatainorder
to find a simple idealized model that captures the pressure in the wellbore and the
field. After comparison of several different models, based on the pressure and
pressure derivative analysis, King et al. (2016) conclude that a radial symmetric
model can best represent the pressure behavior.
Now we present the pressure data for phase I and II of the injection. Figure
A3.2 shows the recorded and modeled values of pressure change divided by flow
rate versus time during phase I and II cycles after July 10th, 1997 and July 26th,
1999. Derivative of the same parameters is also plotted which shows the signature
of a radial flow. Although we are not modeling every fault and fracture in the
50
domain with the use of a radial flow solution, however the data shows that use of
a radial flow model captures the overall behavior of the system to a great degree.
Figure 3.6: Recorded and modeled values of pressure change divided by flow rate
versus time during phase I (top) and II (bottom) cycles after July 10th, 1997 and
July 26th, 1999. Derivative of the same parameters is also plotted which shows
the signature of a radial flow. Adapted from King et al. (2016).
Figure A3.3, from King et al. (2016), shows the recorded and modeled values
of pressure change divided by flow rate versus time during phase III and IV cycles
after January 8th, 2001 and January 14th, 2007, respectively. Derivative of the
same parameters is also plotted which shows the signature of a radial flow. The
51
orange line shows the emergence of radial flow which manifest later in the injec-
tion process. This observation is consistent with our radial pressure model of the
injection in Paradox Valley. Prats et al. (1961) show that for large times and far
wellbore pressure behavior for a fractured well, the pressure behavior is almost
identical to a radial reservoir without any fractures. Due to the extended nature
of wastewater disposal and affected areas, we use a pseudo-radial pressure solution
for assessing the pore-pressure in the reservoir. The pseudo-radial pressure solu-
tion and radial pressure solutions are identical in the mathematical form and are
used in this dissertation.
52
Figure 3.7: Recorded and modeled values of pressure change divided by flow rate
versus time during phase III (top) and IV (bottom) cycles after January 8th, 2001
and January 14th, 2007. Derivative of the same parameters is also plotted which
shows the signature of a radial flow. The orange line shows the emergence of radial
flow (Adapted from King and Block, 2016).
Figure 3.8: Schematic diagram of radial reservoir model and expected seismicity
in the case of bounded and unbounded reservoirs
53
Chapter 4
Probabilistic Bounds on
Maximum Magnitude Induced
Seismicity (MMIS)
4.1 Introduction
The possible magnitude of induced seismicity in subsurface injection applica-
tions is one of the main questions that scientists would like to answer. There
have been efforts to come up with models to predict the maximum magnitude
induced seismicity as a function of the injected volume. However, current models
in the literature lack inclusion of important parameters such as flow rate, reservoir
permeability, flow parameters, as well as reservoir boundary conditions. Here we
develop a model to predict the number of events above any certain magnitude. The
model considers a randomly distributed network of faults that undergoes frictional
failure due to overpressure-induced effective stress changes. The model includes
the effect of flow properties and flow boundary conditions on over-pressure, as
well as the effect of initial stress regime on failure. Theoretical predictions are
compared to observations of induced seismicity in 18 injection induced seismicity
cases. We observe a linear relationship between the log of number of events and
log of injected volume except for the cases in bounded reservoirs and non-constant
flow rates where non-linear behavior is observed. Our model captures an increase
54
in the number and risk of higher magnitude events with an increase in injection
rate, injected volume, and presence of reservoir boundaries. Through application
of our model, injection well operators may be able to calibrate the model based on
site-specific seismicity observations and set operational policies such as reducing
the injection flow rate in order to reduce the risk of larger seismic events.
Historically it is known that subsurface injection might create earthquakes
(Evans, 1966). The Rocky Mountain Arsenal and Rangely, Colorado earthquakes
are all examples of a fluid induced seismicity (Healy et al., 1968; Herrmann et al.,
1981). Starting in 2009, the Central and Eastern U.S. started experiencing an
increase in seismicity rate from background rate. This attracted the attention of
scientists to answer the key questions about the magnitude and number of these
earthquakes with relation to injection activities (Healy et al., 1968; Frohlich, 2012;
Horton,2012;Ellsworth,2013;Keranen et al.,2013;Kim,2013;Zhang et al.,2013).
First, let us review some of the models in the literature for number and mag-
nitude of earthquakes. Shapiro et al. (2010a) using seismogenic index (SI), an
parameter independent of injection parameters, proposed a linear relationship
between the logarithm of the number of induced earthquakes and the logarithm
of the injected fluid volume. McGarr (2014), based on the analysis of several
case histories of induced seismicity concluded that the maximum magnitude and
the maximum seismic moment appeared to have an upper bound proportional to
the total volume of injected fluid. McGarr (2014) proposed that the maximum
seismic moment is limited to the volume of injected liquid times the modulus of
rigidity. Elst et al. (2016) argued that the size of the largest earthquake may not
only be limited to the size of the reservoir affected by injection but in the case
of faults oriented favorably with respect to the tectonic stress field rupture may
55
extend beyond the affected area and the limit of the maximum magnitude earth-
quake might be controlled by the tectonics and not the affected area of the fault.
Elst et al. (2016) showed that injection process may only control the number of
earthquakes triggered but the magnitude of the events are controlled by tectonics
and follow sampling statistics of the Gutenberg-Richter distribution for tectonic
earthquakes. Elst et al. (2016) introduced a probabilistic view to the deterministic
cap proposed by McGarr (2014) with the possibility of sampling higher magni-
tude events as the injection volume increases. Elst et al. (2016) calculated the
total number of events as a function of injected volume by using Shapiro et al.
(2010a) seismogenic index which assumes a linear relationship between logarithm
of injected volume and the logarithm of total number of seismic events.
In Chapter 3, We showed that the cumulative number of earthquakes versus
injected volume for Oklahoma 1997–2014 (Walsh and Zoback, 2015; Weingarten
et al., 2015), KTB 2004–2005 injection experiment (Shapiro et al., 2006a), and
Paradox Valley 1996–2004 saltwater disposal (Ake et al., 2005; King et al., 2014)
show significant deviation from linearity. By developing a physics-based and prob-
abilistic model for cumulative number of seismic events, we showed that these
non-linearity might be related to the injection rate changes or reservoir bound-
aries. We concluded that, reservoir flow parameters, boundary conditions and
flow rate are among the key parameters that control induced seismicity. Current
induced seismicity magnitude models in the literature do not consider these key
parameters.
In this Chapter we focus on including key parameters such as injection rate,
initial state of stress, injection volume and reservoir characteristics to waste-water
disposal induced earthquake magnitude model which is currently missing in the
literature. In the Chapter three, we introduced a probabilistic induced seismicity
56
framework that integrates the physical processes of pressure diffusion in the reser-
voir and frictional failure along faults with the initial state of stress and reservoir
boundary conditions. We conclude that a physics-based, probabilistic framework
of induced seismicity magnitude can successfully capture the site-specific relation-
ships among injection parameters, reservoir parameters, and expected seismicity
magnitudes. Wecomparetheresultsofdifferentinjectionratescenariosoninduced
seismicity. We observe that given the same volume of injected fluid, lower rate
injections show lower number of induced seismicity as well as lower MMIS. We
perform a sensitivity analysis on the effect of permeability, reservoir size and flow
rate on MMIS. We observe that changes in any of these parameters results in sub-
stantial changes in MMIS output of the reservoir. Therefore, we conclude that
comparing very different reservoirs with different reservoir, geomechanical, and
operational characteristics in order to find a universal maximum magnitude seis-
micity model that fits all the reservoirs despite their very different characteristics
may not be practical.
4.2 Method
We derive a mathematical model to predict the possible magnitude of seismic
events induced by injection-related overpressure. There are five key elements of the
model, (1) A random spatial distribution of faults in the reservoir, (2) the distri-
bution of pressure change required for failure on these faults, (3) the overpressure
caused by fluid injection, (4) the resulting probability of inducing an earthquake
at each point in the domain, which is used to compute the total number of pos-
sible induced events, and (5) the most probable maximum magnitude, M
max
and
associated confidence interval. Below we describe these elements of the model.
57
4.2.1 Fault Distribution and Required Overpressure for
Failure
We show that the required pressure increase for shear failure will have a normal
distribution given initial normal stress distribution for stresses shown in Figure 4.1
(A).Wefindtheprobabilityofactivatingpre-stressedfaultsduetooverpressure, as
the integral of probability density function (PDF) of ΔP
p−req
, from zero to ΔP
p
,
which is defined as the difference in the cumulative density functions (CDF) at
those points. We refer to this integral as G(ΔP
p
), which presents the probability
of activating a fault due to a pressure increase of ΔP
p
as shown in Figure 4.1
(C), the area beneath the probability curve. The probability of inducing a seismic
eventissmallerfartherawayfromthewellboreduetosmallerinducedoverpressure,
ΔP
p
, corresponding to a smaller area beneath the PDF (Figure 4.1. C, left), and
larger closer to the wellbore due to the higher ΔP
p
(Figure 4.1. C. right). This
probability is integrated over the volume of the domain and multiplied by the total
number of faults present in the reservoir in order to calculate the total number of
possible seismic event.
Thecorrespondingnumberofinducedseismiceventsversustimefortwocasesof
identical constant injection rates with different bounded and unbounded reservoir
boundary conditions are illustrated in Figure 4.1 (B). At the early injection times,
the presence of outer reservoir boundaries has no effect on the resolved pressure
changes. As the injection continues, the reservoir pressure starts to be influenced
by the boundaries. The flow is governed by a pseudo-steady flow behavior. Figure
4.1 (B) shows that for identical injection rates, a bounded reservoir will eventually
experience an acceleration in seismicity while both reservoirs will share the same
initial behavior. The observed higher number of induced earthquakes in Figure 4.1
(B) is due to higher reservoir pressure increase in a bounded reservoir compared
58
to an unbounded reservoir after the common initial infinite-acting period. Sealing
faults, pinch-outs and low permeability layers such as shale are among the possible
reservoir boundaries that strongly inhibit fluid flow (Manzocchi et al., 2008). The
absenceoftheseelementscreatesunboundedreservoirswherethepressurebehavior
is described as infinite-acting (Figure 4.1, D).
Figure 4.1: (A) Mohr-Coulomb failure envelope along the initial and ultimate
Mohr circles in blue and red respectively. The normal distribution for S
1
and S
3
are shown. The distribution for required pore-pressure for shear failure is shown
in blue. (B) The number of events versus time in two cases of bounded versus
unbounded reservoir. (C) The probability of inducing an earthquake is equal the
area beneath its ΔP
p−req
distribution shown for far-wellbore and near-wellbore.
(D) The possible reservoir boundaries.
The main question is how the number of events (N) is related to the injection
activities. Shapiro et al. (2010a) suggested that a linear relationship between loga-
rithm of injected volume and logarithm of injected volume may be present. One of
59
theshortcomingofShapiro et al.(2010a)approachistheabsenceofkeyparameters
such as injection flow rate, initial state of stress and reservoir parameters. These
parameters have been shown to have important effects on induced seismicity. We
introduce N
M>M
0
as the number of events greater than magnitude M
0
, which is
integrated over study volume:
N
M>M
0
(t) =


Z
V
1
2
(erf(
ΔP
p
−μ
ΔP
p−req
√
2σ
ΔP
p−req
)−erf(
0−μ
ΔP
p−req
√
2σ
ΔP
p−req
))ρdv


10
−bM
0
(4.1)
where μ
ΔP
p−req
, σ
ΔP
p−req
and σ
2
ΔP
p−req
are the mean, standard deviation and
variance of required pore-pressure increase for shear failure distribution, ρ is the
fault density,b (commonly known as b-value) is a constant anderf(x) is the error
function.
Elst et al. (2016) studied the maximum magnitude possible induced seismic-
ity as a function of injected volume. They show that with the assumption of
Gutenberg-Richter distribution for induced seismicity magnitude distribution, the
most probable maximum magnitude to be observed in a sample size of N is:
ˆ
M
max
=M
0
+
1
b
log
10
N (4.2)
where b is the slope of the power-law in Gutenberg-Richter distribution, M
0
is
a reference magnitude, andN is the number of earthquakes observed greater than
magnitude M
0
. Elst et al. (2016) derived that the 90% confidence level for the
maximum magnitude falls between -0.48 and +1.3 magnitude units of the
ˆ
M
max
.
We use their,
ˆ
M
max
relationship in Equation 2 and by inserting N
M>M
0
from last
chapter, we find the following equation for
ˆ
M
max
as:
60
ˆ
M
max
=
1
b
log
10


Z
V
1
2
(erf(
ΔP
p
−μ
ΔP
p−req
√
2σ
ΔP
p−req
)−erf(
0−μ
ΔP
p−req
√
2σ
ΔP
p−req
))ρdv


(4.3)
As it can be observed from Equation 3,
ˆ
M
max
is a function of pressure and
required increase in pore-pressure for shear failure. In other words,
ˆ
M
max
is a
function of all the reservoir flow parameters as well as the initial stress state of
the reservoir. Next, in the sensitivity analysis section we study the total num-
ber of events greater than magnitude M
0
as well as the most probable maximum
magnitude induced seismicity.
4.3 Results
In Chapter three, we showed that injection flow rate, reservoir permeability
and reservoir size are among the key parameters in injection induced seismicity.
Here we perform sensitivity analysis for these three main parameters: (1) reservoir
size (2) permeability (3) flow rate and (4) periodic injections. We observe that
non-linear behavior in the log-log plot for the number of events versus the injected
volume is due to the boundary effects. This non-linear behavior is also seen in
ˆ
M
max
. The reservoir and injection parameters used in this Chapter for the base-
case are presented in Table 1.
4.3.1 Permeability
The plot of the logarithm of the number of events greater than M
0
= 2 versus
the logarithm of the injected volume is shown in Figure 4.2 (left). It is observed
that a linear relationship in the log-log plot of the number of events greater than
61
M
0
= 2 versus the injected volume is present. Looking at the most probable
maximum magnitude seismicity, we observe a linear relationship between log of
injected volume and the
ˆ
M
max
. We observe that as the reservoir permeability
increases the total number of events as well their magnitude, increases. This is
due to the fact that in this case study, the faults are critically stressed and are on
the verge of failure. Therefore, they only will need minimal pressure increase for
shear failure. As the permeability increases, the area affected by the pore-pressure
expands, which results in a greater number of events and a greater
ˆ
M
max
.
Figure 4.2: Sensitivity analysis for the effect of permeability on N
total
and
ˆ
M
max
(left) linear log-log plot for the number of events greater than M
0
= 2 versus the
injected volume (right)
ˆ
M
max
plotted versus the logarithm of the injected volume
showing a linear behavioir.
4.3.2 Reservoir Size
Figure 4.3 shows the log-log plot of the number of events greater than M
0
= 2
(left) and
ˆ
M
max
(right), both versus the logarithm of the injected volume, demon-
strating a non-linear behavior. It is observed in Figure 4.3, (left) that the linear
relationship between the logarithm of the injected volume and the logarithm of the
number of events is not always present as suggested by Shapiro et al. (2010a). We
observe that the presence of reservoir boundaries is the main reason for deviation
62
from the linear behavior. As it can be observed in Figure 4.3, as the reservoir size
decreases, the deviation from linear behavior increases, where the most non-linear
log-log and semi-log behavior is related to the smallest reservoir. As the reservoir
size increases and the effects of boundary conditions diminish we observe more of
a linear behavior which is due to the non-existence of boundaries.
Figure 4.3: Sensitivity analysis for the effect of reservoir size on N
total
and
ˆ
M
max
(left) Non-linear log-log plot for number of events greater than M
0
= 2 versus
injected volume (right)
ˆ
M
max
plotted versus log of injected volume showing a non-
linear behavior.
4.3.3 Flow Rate
Figure 4.4 shows the log-log plot of the number of events greater than M
0
= 2
(left) and
ˆ
M
max
(right) both versus the logarithm of the injected volume showing
a linear behavior. A linear relationship in the log-log plot of the number of events
greater than M
0
= 2 versus the injected volume is observed. Looking at the most
probable maximum magnitude seismicity, we observe a linear relationship between
the logarithm of injected volume and the
ˆ
M
max
. We observe that as the injection
rate increases, the total number of events, as well as their magnitude, increases.
63
Figure 4.4: Sensitivity analysis for the effect of injection rate on N
total
and
ˆ
M
max
(left) linear log-log plot for number of events greater than M
0
= 2 versus injected
volume (right)
ˆ
M
max
plotted versus log of injected volume showing a linear behav-
ior.
4.3.4 Net Injected Volume
One of the important factors on induced seismicity behavior is the effect of
reservoir injection history. Production from a reservoir reduces the initial pore-
pressure in the reservoir which counteracts the effects of injection in a reservoir.
Walsh and Zoback (2015) observe that enhanced oil recovery (EOR) in Oklahoma
has the lowest chance of inducing earthquakes. This is due to the effects of simul-
taneous injection and production which does not allow the pore-pressure increase
to go beyond a certain limit and to expand beyond certain points. There will
be an steady state for pressure which does not promote any seismic events. This
behavior is justified by Equation 2 where the total number of events is a function
of pressure increase.
4.3.5 Periodic Injections
In this section, we analyze the effects of different injection scenarios on induced
seismicity response. We study three different injection rate scenarios for the inter-
val of 120 days. These three scenarios compromise of 1) a constant flow rate of 159
64
m
3
/day 2) a periodic injection of 159m
3
/day with the intervals of 10 days of injec-
tion and 10 days of no injection, and finally 3) a reduced rate to half which equals
79.5m
3
/day. The injected volume for Case 1 and 3 are the same through the 120
days of injection in this simulation. However, as it can be observed in Figure 6
(left), the magnitude of
ˆ
M
max
and N
total
, number of events greater than M
0
= 2
in the Case of lower flow rate injection is smaller compared to the same volume
injection with higher periodic injection. Both cases injected the same amount of
volume into the reservoir. However, lower rate injection reduced the seismicity risk
more than the periodic injection. For the sake of comparison, a constant injection
of 159m
3
/day is simulated and as it is expected due to the double injected volume
at the end of 120 days of injection, a higher
ˆ
M
max
and N
total
, number of events
greater than M
0
= 2 is expected. Based on our observations, we conclude that
in order to reduce the total number of induced earthquakes and their magnitudes,
lower injection rates would be more beneficial rather than periodic injection rate
reductions.
Figure 4.5: Sensitivity analysis for the effect of periodic injection on N
total
,
ˆ
M
max
andP
p
(left) plot of pressure versus time for three cases of injections with different
flow rates of 159 m
3
/day (2) periodic injection of 159 m
3
/day with the intervals
of 10 days of injection and 10 days of no injection and finally (3) reduced rate to
half which equals 79.5 m
3
/day. (right) semi-log plot for number of events greater
than M
0
= 2 versus time on the right axis and
ˆ
M
max
plotted versus time on the
left axis.
65
4.4 Possibility of Universal Maximum Magni-
tude Induced Seismicity Model
We discussed the importance of reservoir permeability, injection rate, reservoir
boundaries and injection schedule on induced seismicity outcome. We showed that
these parameters have substantial effects on MMIS. McGarr (2014), and Elst et al.
(2016)plottedthemaximummagnitudeinducedseismicityagainstthelogarithmof
the injected volume in order to find a universal relationship between the maximum
possible magnitude seismicity and the injection volumes. The cases that have
been plotted are very different in terms of the permeability structure, reservoir
boundaries, injection rate, and injection schedule. As a result, we conclude that
comparing these cases without considering the effects of all the mentioned key
parameters may result in inaccurate conclusions about the maximum magnitude
possible induced seismicity. We recommend to study each case on its own due to
their very different characteristics. For the individual wastewater disposal cases
it would be a more accurate approach to consider the effects of injection rate,
boundary conditions, permeability structure, initial state of stress, net injected
volume, and injection schedule, as these factors have important effects on the
induced seismicity behavior.
4.5 Conclusions
We introduced a new model for the most probable maximum magnitude seis-
micity due to wastewater injection. Our model is a physics-probabilistic-based
66
model which accounts for initial state of stress, reservoir flow and injection charac-
teristics. The inclusion of these important parameters in previous induced seismic-
ity models in the literature is absent. By performing a sensitivity analysis for the
effect of reservoir permeability, reservoir boundary conditions, injection rate, and
injection schedule we observe that reservoir boundaries cause a non-linear behavior
on the log-log plot of the number of induced seismicity events as well as semi-log
plot of the most probable maximum magnitude induced seismicity. This results in
increasing the number and magnitude of the events. We observe that increasing
theflowrateandreservoirpermeabilitywillpromotelargerareasaffectedbywaste-
water disposal. In the case of critically stressed faults close to failure an increase
in the number and magnitude of possible induced seismicity is observed. Given
the same volume of waste-water injection, we observe that a reduced injection rate
results in fewer numbers of induced seismicity and smaller events compared to
the same volume injected periodically with intervals of shut-in but with a double
injection rate. Using our induced seismicity model, injection well operators can
calibrate the model based on induced seismicity data to plan different injection
and mitigation scenarios for induced seismicity in order to reduce the possibility
of inducing a large earthquake.
Appendix
We introduce N
M>M
0
as the number of events greater than magnitude M
0
,
which is integrated over study volume:
67
Table 4.1: KTB Parameters
Parameter and Variable Value and Unit
Permeability (k) 5× 10
−15
m
2
Fluid viscosity (η) 10
−3
Pas
Injection rate (q) 3.02× 10
−3
m
3
/s
Total compressibility (C
t
) 7.25× 10
−10
Pa
−1
Nominal friction (μ) 0.65
Effective normal stress (¯ σ) 16.67 MPa
Height (h) 150 m
Porosity (φ) 0.04
External reservoir radius (r
e
) 1650 m
Wellbore radius (r
w
) 0.152 m
Mean required pressure increase for shear failure (μ
ΔP
p−req
) 1.44 MPa
Standard deviation in required pressure for failure (σ
ΔP
p−req
) 0.48 MPa
Minimum horizontal stress (Sh
min
) 80 MPa
Maximum horizontal stress (SH
Max
) 174.35 MPa
Standard deviation in minimum horizontal stress (σ
Sh
min
) 0.3 MPa
Standard deviation in maximum horizontal stress (σ
SH
Max
) 0.5 MPa
N
M>M
0
(t) =


Z
V
1
2
(erf(
ΔP
p
−μ
ΔP
p−req
√
2σ
ΔP
p−req
)−erf(
0−μ
ΔP
p−req
√
2σ
ΔP
p−req
))ρdv


10
−bM
0
(A4.1)
where μ
ΔP
p−req
, σ
ΔP
p−req
and σ
2
ΔP
p−req
are the mean, standard deviation and
variance of required pore-pressure increase for shear failure distribution, ρ is fault
density, b commonly known as b-value is a constant, erf(x) is error function
formulated as:
erf(x) =
1
√
π
Z
x
−x
e
−t
2
dt (A4.2)
Elst et al. (2016) showed that with the assumption of Gutenberg-Richter distri-
bution for induced seismicity magnitude distribution, the most probable maximum
magnitude to be observed in a sample size of N is:
68
M
max
=M
0
+
1
b
log
10
N
M>M
0
(A4.3)
where b, M
0
and N
M>M
0
have the same definition as before.
M
max
=M
0
+
1
b
log
10


Z
V
1
2
(erf(
ΔP
p
−μ
ΔP
p−req
√
2σ
ΔP
p−req
)−erf(
0−μ
ΔP
p−req
√
2σ
ΔP
p−req
))ρdv


10
−bM
0
(A4.4)
69
Chapter 5
A Probabilistic, Geomechanics
based Interpretation of
Stimulated Rock Volume in
Hydraulic Fracturing
5.1 Introduction
In order to asses and monitor the performance of hydraulic fracturing, micro-
seismicmonitoringisextensivelyused. Theclassicviewofasingleplanarhydraulic
fracture is challenged by the production and microseismic observations from low
permeability shale reservoirs, where the existence of a pre-existing natural network
of fracture is the key to production. Stimulated rock volume (SRV) which refers
to a portion of microseismic cloud associated with fracturing, is used to examine
the performance of hydraulic fracturing. In this study, we identify the key param-
eters that promote a larger SRV and a successful hydraulic fracturing. We use
linear elastic fracture mechanics (LEFM) to illustrate near-field and far-field stress
behavior at the crack tip and surrounding the fracture, respectively. Using Mohr-
Coulomb shear failure criteria, we define proximity to failure (PF) as a parameter
to capture the probability of shear failure on a network of natural fractures. We
use PF along the stress perturbation around the hydraulic fracture from LEFM to
70
simulate the possible SRV of a Frac-job in a typical shale reservoir in four different
settings. 1. High deviatoric stress and fracture pressure, 2. High deviatoric stress
and low fracture pressure, 3. Low deviatoric stress and fracture pressure with high
initial reservoir pressure, and 4. Low deviatoric stress and fracture pressure with
low initial reservoir pressure. We observe that these four different cases which
vary only in deviatoric stress, fracture pressure, and initial reservoir pore-pressure
can drastically yield very different SRV signatures despite their very similar main
hydraulic fracture. We identify that the initial reservoir state of stress, fracture
pressure, reservoir pressure, and pre-existing fractures are among the key param-
eters that determine SRV and the ultimate performance of a hydraulic fracturing
treatment. Knowledge about these parameters is an essential step in the design
and monitoring of fracturing stages, spacing and pumping schedule to reach better
oil and gas recovery from a reservoir.
In recent decades, as the technology of horizontal wells advanced, there is a
boost in the application of hydraulic fracturing in the oil and gas industry. Thou-
sands of hydraulic fracture treatments are conducted each year in a wide range of
geological formations varying from low permeability gas fields (shale formations),
weakly consolidated offshore sediments such as those in the Gulf of Mexico, soft
coal beds for methane extraction and naturally fractured reservoirs (Adachi et al.,
2007). The goal of microseismic monitoring is to find the location of a created
hydraulic fracture and activated network of fractures. Hydraulic fracturing has
become a common practice in the oil and gas industry to the point that, of the
producing wells drilled in North America since the 1950s, about 70% of gas wells
and 50% of oil wells have been hydraulically fractured (Valko and Economides,
1995). Figures 5.1 and 5.2 show the complex network of hydraulic fracture and
71
existing natural fractures. Stimulated rock volume (SRV) referring to the micro-
seismic cloud associated with fracturing is used to examine the performance of
hydraulic fracturing.
Figure 5.1: Hydraulic fracturing in a horizontal wellbore in a shale layer.
Although SRV is used extensively in a hydraulic fracturing treatment design
and evaluation, a study to investigate the effects of different key parameters on
SRV behavior is absent in the literature. In this study, we use LEFM near-field
and far-field solutions to evaluate the extent of SRV accompanied by a single stage
of a planar hydraulic fracture. First, a brief review of LEFM, including near-tip
solution is presented [27, 28, 29, 30]. Then a typical shear failure zone at the tip
of fracture is presented [26]. Due to dependence of shear failure on minimum and
maximum reservoir principal stresses, a single parameter describing shear failure
72
Figure 5.2: Top: Interaction of hydraulic fracture and natural network of fractures
(from Huang et al. 2014). Bottom: SRV around hydraulic fracture stages.
and proximity to failure is absent in the literature. To describe proximity to fail-
ure, we define a single parameter named proximity to failure (PF) as the shortest
distance between the failure envelope and material’s state of stress. We present
a full stress distribution solution for a mode I fracture where a constant pressure
inside the fracture propagates the crack, is presented [31]. Based on the Anderson
faulting theory (Anderson, 1951), a typical reservoir stress state with typical Bar-
nett shale mechanical properties (Table 5.1) is used for fracture propagation and
SRV simulations.
We compare the size and shape of SRV for four different settings: 1. High
deviatoric stress and fracture pressure, 2. High deviatoric stress and low fracture
73
Figure 5.3: (A) Mohr-Coulomb failure envelope along initial and ultimate Mohr
circles in blue and red respectively. Normal distribution for S
1
and S
3
are shown.
Distribution for required proximity to failure is shown in blue. (B) The probability
of inducing an earthquake is equal the area beneath its ΔPF
req
distribution shown
for far-wellbore and near-wellbore.
pressure, 3. Low deviatoric stress and fracture pressure with high initial reser-
voir pressure, and 4. Low deviatoric stress and fracture pressure with low initial
reservoir pressure. We observe very different SRV activated as a result of same
hydraulic fracture but with different reservoir geomechanical conditions. Finally,
we conclude that reservoir state of stress, pore-pressure, injection pressure, and
presence of existing faults and fractures are among the key parameters pertaining
to SRV and observed microseismicity. We show that SRV is highly dependent on
deviatoric stress, fracture pressure, reservoir pressure, and presence of pre-existing
74
faults or fractures. Therefore, any effort on design and diagnosis of hydraulic
fracturing should consider these key elements.
5.2 Linear Elastic Fracture Mechanics
Most models developed for hydraulic fracturing are based on linear elastic frac-
ture mechanics (LEFM), tensile fracture (Mode I fracture), lubrication theory for
fluid flow through fracture, and Carter’s leak-off model for flow through porous
medium (Adachi et al. 2007; Daneshy, 1973; Perkins and Kern, 1961; Geertsma
and de Klerk, 1969). These hydraulic fracturing models are tailored for hard rock
formations with brittle linear elastic behavior and low permeability reservoirs.
LEFM models are based on stress singularity and strain energy dissipation in a
very small zone at the fracture tip.
Linear elastic fracture mechanics (LEFM) is used in fracture toughness test-
ing of rocks (Ouchterlony, 1982; Ingraffea et al., 1984; Labuz et al., 1987), as
well as the design and evaluation of hydraulic fractures (Thiercelin, 1989). Sev-
eral studies have been devoted to the effects of nonlinear processes at the crack
tip including plasticity in hydraulic fracturing applications. These studies have
shown that LEFM based models predict a longer fracture with less width than the
estimated values using pressure build-up test and production data (Ouchterlony,
1989). While, LEFM is applicable to perfectly brittle material, which have a linear
elastic behavior until they fail, rocks usually show a quasi-brittle behavior which
means that they exhibit pre-peak non-linearity and strain softening (Desroches
et al., 1994; Papanastasiou, 1997).
75
5.2.1 Near-Tip Hydraulic Fracture Stress Solution
At the crack tip in geomaterial testing, a zone of micro-cracking is observed
(Labuz et al., 1987). In metals, this zone corresponds to plastic yielding of metals
and in rocks it corresponds to tensile or shear failure of the rock. It has been shown
that depending on the size of the process zone, instead of a LEFM type fracture
behaviorwemighthaveabluntfracture. ThiscannotbesolelydescribedbyLEFM,
and non-linear fracture mechanics (NLFM) might also be needed (Papanastasiou,
1997).
Closed-form solutions for stress distribution around the crack-tip is based on
LEFM (Westergaard, 1939; Sneddon, 1946; Irwin, 1957). Solutions for the stress
distribution around the crack tip (near tip solutions) for a mode I fracture are as
follows,
S
xx
=
−K
I
√
2Πr
[cos
Θ
2
(1−sin
Θ
2
sin
3Θ
2
)] (5.1)
S
yy
=
−K
I
√
2Πr
[cos
Θ
2
(1 +sin
Θ
2
sin
3Θ
2
)] (5.2)
S
xy
=
−K
I
√
2Πr
[cos
Θ
2
sin
Θ
2
cos
3Θ
2
)] (5.3)
where S
xx
is the total compressive stress in x direction, S
yy
is the total com-
pressive stress in y direction, S
xy
is the shear stress on xy plane, and K
I
is the
mode I fracture intensity factor. The definition of the coordinate system, Θ
1
, r
1
,
and stress sign conventions are presented in Figure 5.4.
Based on the near-tip solution that was presented, it can be observed that
LEFM describes a 1/r singularity at the crack tip which is proportional to the
76
Figure 5.4: Fracture geometry and geometric parameters used in the near-tip solu-
tion as well as stress sign convention.
mode I fracture intensity factor. Therefore, near tip solutions are only valid in
the vicinity of the crack tip where 1/r singularity dominates the stress field. This
region is defined as the singularity-dominated zone (Irwin, 1957). However, in
reality there cannot be a singularity in stress. Instead, a zone of failure, called
process zone, develops at the crack tip. As the size of the process zone at the crack
tip increases the elastic stress analysis becomes increasingly inaccurate. Simple
corrections to the LEFM solution based on a moderate size process zone have
been proposed. Both process and singularity-dominated zones are schematically
shown in Figure 5.5. Figure 5.6 shows maximum, S
1
and minimum, S
3
principal
stress concentration at the crack tip for a fracture half-length of 3 ft, S
1
of 2000
psi and S
3
of 1000 psi with a fracture pressure of 1100 psi. Cooler colors at the
tip of the fracture show a tensile regime. Surrounding the face of the fracture, the
compressive stress regime in hot colors are dominant.
5.2.2 Full-Field Hydraulic Fracture Stress Solution
Near-field solution is only valid in the singularity dominated zone with the
assumption of a small process zone and does not include the far-field stresses. The
77
Figure 5.5: Process zone shown in red and singularity-dominated zone shown in
green, at the crack tip.
Figure 5.6: Maximum, S
1
and minimum, S
3
principal stress concentration at the
crack-tip shown on left and right respectively for a fracture half-length of 3 ft.,
S
1
of 2000 psi and S
3
of 1000 psi with fracture pressure of 1100 psi. Cooler colors
at the tip of the fracture show a tensile regime while surrounding the face of the
fracture, compressive stress regime in hot colors are dominant.
full-field stress solution is valid everywhere in the domain and is a function of stress
boundary conditions and fracture geometry. A full stress distribution solution for
a mode I fracture where a constant pressure inside the fracture propagates the
crack, is presented (Pollard, 1987):
78
S
xx
=S
r
xx
− (P
Frac
−S
r
yy
)[rR
−1
cos(Θ− Γ)− 1−a
2
rR
−3
sinΘsin3Γ] (5.4)
S
yy
=S
r
yy
− (P
Frac
−S
r
yy
)[rR
−1
cos(Θ− Γ)− 1 +a
2
rR
−3
sinΘsin3Γ] (5.5)
S
xy
=S
r
xy
− (P
Frac
−S
r
yy
)[a
2
rR
−3
sinΘcos3Γ] (5.6)
WhereS
xx
,S
yy
andS
xy
have the same definition as before,S
r
xx
andS
r
yy
are the
totalremotestressesinxandy directionrespectively,S
r
xy
istheremoteshearstress
whichisequaltozerointhiscase,P
Frac
isthefracturingpressureinsidethefracture
which is reported to be about 200 psi greater than the minimum vertical stress
acting on the face of the fracture for Barnett shale (Vermylen, 2011). R =
q
(r
1
r
2
)
and r
1
, r
2
, r, Θ
1
, Θ
2
, Θ as well as fracture geometry are shown in Figure 5.7.
Figure5.7: Fracturegeometryandgeometricparametersusedinthefull-fieldstress
solution for hydraulic fracturing.
Usingafull-fieldstresssolution, atypicalshearfailurezoneatthetipoffracture
foracaseofS
1
of8040psi,S
3
of1740psi, fracturepressureof2240psi, andfracture
half-length of 100 ft. is shown in Figure 5.8. As it can be observed in Figure 5.5,
79
the size of shear failure zone for a constant failure envelope in a 100 ft. fracture
is roughly 2 ft. wide and 8 ft. long. The full stress field is used in the upcoming
sections to find the size and shape of SRV in typical shale fracturing applications.
In the next section failure criteria for shear failure is presented which is then used
to investigate the SRV.
Figure 5.8: A typical shear failure zone at the tip of fracture for a case of S
1
of
8040 psi, S
3
of 1740 psi, fracture pressure of 2240 psi and fracture half-length of
100 ft. The size of shear failure zone for a constant failure envelope in a 100 ft
fracture is roughly 2 ft wide and 8 ft long.
5.3 Shear Failure Criteria
Shear failure is induced by effective stress changes on faults as a result of
changes in reservoir pressure due to leak-off and total principal stress change due
to stress concentration around the fracture. We present the Mohr-Coulomb failure
criteria for critically stressed faults at the verge of failure (after (Zoback, 2010)):
S
1
−P
p
= (S
3
−P
p
)n (5.7)
80
wheren = (
q
1 +μ
2
f
+μ
f
)
2
,μ
f
isfaultcoefficientoffriction,P
p
ispore-pressure,
S
1
and S
3
are maximum and minimum total principal stresses, respectively. The
state of stress at each point in the reservoir can be described by S
1
and S
3
. To
find proximity to shear failure at each point in the material, we find the minimum
verticaldistancefromfailureenvelope foreachpointintheS
1
-S
3
domainandname
it as change in proximity to failure (ΔPF):
ΔPF =
(S
1
−S
3
n−P
p
(1−n))
√
1 +n
2
(5.8)
We set the sign convention for ΔPF. Negative ΔPF relates to the state of
material where it needs to be driven to failure through higher pore-pressure, higher
deviatoric stress, ora combinationof both whilepositive ΔPF relates tothe points
in the material where shear failure already has happened and those points are
beyond failure.
The state of stress on the population of faults in the reservoir can be described
byS
1−Initial
andS
3−Initial
which are distributions of maximum and minimum prin-
cipalstressesonthefaultsinthereservoir. Tofindtherequiredchangeinproximity
to failure at each fault in order to shear it, we find the minimum vertical distance
from failure envelope for stress state in theS
1−Initial
-S
3−Initial
domain and name it
as required change in proximity to failure. ΔPF
req
. ΔPF
req
is the distribution of
required change in proximity to failure on the network of natural fractures, faults,
joints and fissures around the hydraulic fracture:
ΔPF
req
=
(S
1−Initial
−S
3−Initial
n−P
p−Initial
(1−n))
√
1 +n
2
(5.9)
81
where P
p−Initial
is initial reservoir pore-pressure, S
1−Initial
and S
3−Initial
are
maximum and minimum total principal initial reservoir principal stress distribu-
tions acting on discontinuities, respectively.
ΔPF
req
is a linear combination ofS
1−Initial
,S
3−Initial
, andP
p−Initial
. Moos and
Zoback (1990), and Zoback (2010) among others show that stress distribution in
fault networks can be described by a normal distribution. The principal stresses,
S
3−Initial
andS
1−Initial
, that drive failure on individual faults within the network of
reservoir-bounded faults are assumed to follow a normal distribution. In this case,
based on probability theory, the linear combination of several normal distributions
remainsnormallydistributed. Therefore, givenanormaldistributionformaximum
and minimum stresses, ΔPF
req
, the required proximity to failure, can also be
described by a normal distribution.
Figure 5.7 (A) shows the distributions ofS
1−Initial
,S
3−Initial
, and ΔPF
req
along
with the Mohr circles before and after hydraulic fracturing with the shear failure
envelope. As it is shown in Figure 5.7, by the changes in S
1
, S
3
and change in
radius of Mohr circle, both pore-pressure increase and change in stress around the
fracture are the mechanisms that drive shear failure.
5.4 Probability of an Event
Based on probability theory, one can find the probability of activating pre-
stressed faults due to changes in pore-pressure and effective stress, as the inte-
gral of probability density function (PDF) for ΔPF
req
from zero to ΔPF. This
is defined as the difference in the cumulative density functions (CDF) at those
points. In this study, we are only concerned with positive values of ΔPF
req
, which
82
Figure 5.9: (A) Mohr-Coulomb failure envelope along initial and ultimate Mohr
circles in blue and red respectively. Normal distribution forS
1−Initial
andS
3−Initial
are shown. Distribution for required pore-pressure for shear failure is shown in
blue. (B) The probability of inducing an earthquake is equal the area beneath
its ΔPF
req
distribution shown for far-wellbore and near-wellbore. (D) Possible
reservoir boundaries.
correspond to an increase in proximity to failure due to leak-off and stress concen-
tration in hydraulic fracturing. Therefore, we only focus on the positive part of
the ΔPF
req
distribution from zero to ΔPF. We refer to this integral asG(ΔPF ),
which presents the probability of activating a fault due to an increase in proximity
in failure of ΔPF, shown in Figure 2 (C), as the area beneath the probability
curve. Specifically, the integral can be expressed by:
G(ΔPF ) =
Z
ΔPF
0
f(ΔPF
req
)dΔPF =F (ΔPF )−F (0) (5.10)
83
where f and F are the PDF, and the CDF, respectively. The probability of
inducing a seismic event is smaller the farther away from the hydraulic fracture due
to a smaller induced pressure and stress concentration, ΔPF. This corresponds
to a smaller area beneath the PDF (Figure 5.7, B, left), and larger, closer to the
hydraulicfracture due tothe higher ΔPF (Figure 5.7, B,right). It should benoted
that the required proximity to failure change for failure can have any arbitrary
probability density function and Equation (5.10) holds true for any PDF.
5.5 Sensitivity Analysis for SRV Shape and Size
Microseismic events are created due to shear slippage/failure around the
hydraulic fractures (Pearson 1982; Rutlede and Philips 2009; Warpinski 2004).
A combination of pore-pressure increase due to fracturing fluid leak-off plus stress
concentration at the tip of hydraulic fracture and surrounding areas, results in
shearfailureandmicroseismicactivity. ThereforeLeak-offoffracfluidandmechan-
ical stress concentration at the crack tip are the main 2 mechanisms leading to
microseismicity. Due to a very low matrix permeability in shale reservoirs, here
we mainly focus on the mechanical stress perturbation due to hydraulic fracture
propagation.
We simulate the growth of fracture from injection point to 80 ft of fracture half-
length and observe the size and shape of SRV from the starting point to the 80 ft.
of half-length of fracture. We compare the size and shape of SRV for four different
settings of: 1. High deviatoric stress and fracture pressure, 2. High deviatoric
stress and low fracture pressure, 3. Low deviatoric stress and fracture pressure
with high initial reservoir pressure, and 4. Low deviatoric stress and fracture
pressure with low initial reservoir pressure. Table 5.1 shows the input parameters
84
Table 5.1: Nominal Hydraulic Fracture Simulation Parameters
Parameter and Variable Value and Unit
Nominal friction (μ) 0.6
Mean required increase in proximity to failure (μ
ΔPFreq
) 1000 psi
Standard deviation in required pressure increase for failure (σ
ΔPFreq
) 1000 psi
Case(1) Maximum principal stress (S
1
) 6000 psi
Case(1) Minimum principal stress (S
3
) 4000 psi
Case(1) Fracture pressure (P
Frac
) 4500 psi
Case(1) Reservoir pore-pressure (P
Initial
) 3000 psi
Case(2) Maximum principal stress (S
1
) 6000 psi
Case(2) Minimum principal stress (S
3
) 4000 psi
Case(2) Fracture pressure (P
Frac
) 4100 psi
Case(2) Reservoir pore-pressure (P
Initial
) 3000 psi
Case(3) Maximum principal stress (S
1
) 6000 psi
Case(3) Minimum principal stress (S
3
) 5000 psi
Case(3) Fracture pressure (P
Frac
) 5100 psi
Case(3) Reservoir pore-pressure (P
Initial
) 4000 psi
Case(4) Maximum principal stress (S
1
) 6000 psi
Case(4) Minimum principal stress (S
3
) 5000 psi
Case(4) Fracture pressure (P
Frac
) 5100 psi
Case(4) Reservoir pore-pressure (P
Initial
) 3000 psi
for four hydraulic fracture simulation case studies used in this paper.Next, we
present the results of SRV simulation for four sensitivity analysis cases.
5.5.1 Case (1). High deviatoric stress and fracture pres-
sure
First, we present the Case (1) with high deviatoric stress and fracture pressure.
The simulation results for Case (1) is shown in Figure 5.10. As it is presented
in Table 5.1, the deviatoric stress and fracture pressure is relatively higher that
the other three cases in this study. Higher deviatoric stress results in more prone
to shear failure areas for the reservoir. Therefore, a smaller amount of pressure
perturbation will result in shear failure, while higher injection pressure results in
85
higher stress concentration in the reservoir which results in more extended SRV
for this case. Figure 5.10 shows the ΔPF parameter at the vicinity of a hydraulic
fracture. The areas where the shear potential parameter is positive are prone to
shear failure while the areas where the shear potential parameter is negative are
much less prone to failure or more stable. The hydraulic fracture tip due to high
stress concentration at the crack tips, show a higher proximity to failure shown
with warmer colors on the colormap of proximity to failure. Higher proximity to
failure and larger area of positive shear potential function results in larger shear
failure zone at the tip of a fracture and subsequently larger SRV and microseismic
cloud associated with the fracture.
Figure 5.10: Case (1). High deviatoric stress and fracture pressure. The change in
proximity to failure is plotted. Warmer colors refer to higher proximity to failure
resulting in SRV while cooler colored areas are less prone to shear failure and SRV.
86
5.5.2 Case (2). High deviatoric stress and low fracture
pressure
Now let us look at Case (2) with High deviatoric stress and low fracture pres-
sure. The simulation result for Case (2) is shown in Figure 5.11. As it is presented
in Table 5.1, the only difference between Case (1) and (2) is the lower fracture
pressure in Case (2). High deviatoric stress results in more prone to shear failure
areas for the reservoir hence a smaller amount of pressure perturbation will result
in shear failure, while lower injection pressure results in lower stress concentration
in the reservoir which results in less extended SRV for this case. Figure 5.11 shows
the ΔPF parameter at the vicinity of a hydraulic fracture. The hydraulic fracture
tip due to high stress concentration at the crack tips, show a higher proximity to
failure shown with warmer colors in the colormap of proximity to failure. Lower
proximity to failure and smaller areas of positive proximity to failure results in
a smaller shear failure zone compared to Case (1) at the tip of the fracture and
subsequently smaller SRV and microseismic cloud associated with the fracture.
5.5.3 Case (3). Low deviatoric stress and fracture pressure
with high initial reservoir pressure
Case (3) is devoted to Low deviatoric stress and fracture pressure with high
initial reservoir pressure. The simulation result for Case (3) is shown in Figure
5.12. As it is presented in Table 5.1, we have a Lower deviatoric stress and fracture
pressure for this case, however, we have an high initial reservoir pressure which
results in a reservoir that is at the verge of shear failure. Higher Reservoir pressure
results in more prone to shear failure environment for the reservoir and therefore
a smaller amount of stress perturbation will result in shear failure. However, the
87
Figure 5.11: Case (2). High deviatoric stress and low fracture pressure. The
change in proximity to failure in plotted. Warmer colors refer to higher proximity
to failure resulting in SRV while cooler colored areas are less prone to shear failure
and SRV.
amountofinjectionpressureislowerinthiscasecomparedtoCase(1)andtherefore
as Figure 5.12. shows the proximity to failure observed in this case is lower and
the extent of SRV and microseismic cloud associated with fracture is smaller.
5.5.4 Case (4). Low deviatoric stress and fracture pressure
with low initial reservoir pressure
Case (4) depicts SRV behavior for a Low deviatoric stress and fracture pressure
with low initial reservoir pressure. The simulation result for Case (4) is shown in
Figure 5.13. As it is presented in Table 5.1, a lower deviatoric stress and fracture
pressure with lower initial reservoir pressure is the characteristic of this case which
makes it harder for the reservoir faults and fractures to reach shear failure. Less
prone to shear failure environment results in smaller SRV and microseismicity and
88
Figure 5.12: Case (3). Low deviatoric stress and fracture pressure with high initial
reservoir pressure. The change in proximity to failure in plotted. Warmer colors
refer to higher proximity to failure resulting in SRV while cooler colored areas are
less prone to shear failure and SRV.
as it can be seen, dominant presence of cooler colors and negative proximity to
failure is observed in Figure 5.13.
5.6 Conclusion
Through simulation of different cases we observed that reservoir state of stress,
pore-pressure, injection pressure and presence of existing faults and fractures are
among the key parameters pertaining to SRV and observed microseismicity.
We observe very different SRVs activated as a result of same hydraulic fracture
but with different reservoir geomechanical conditions. Reservoir deviatoric stress,
fracture pressure and initial reservoir pore-pressure are among the key parameters
pertaining to SRV and observed microseismicity. Higher deviatoric stress shows to
have the largest effect in promoting shear failure in the reservoir and a larger SRV.
Higher reservoir pore-pressure and injection pressure along with higher deviatoric
89
Figure 5.13: Case (4). Low deviatoric stress and fracture pressure with a low initial
reservoir pressure. The change in proximity to failure in plotted. Warmer colors
refer to higher proximity to failure resulting in SRV while cooler colored areas are
less prone to shear failure and SRV.
stress create a suitable environment for shear failure and extension of SRV and
productive SRV.
TheidealcombinationformaximizingSRVandproductiveSRVishighfracture
pressure, high initial reservoir pressure, high stress difference between minimum
and maximum reservoir stress. Presence of natural fractures, faults, joints and
geologic discontinuities are among the most important contributors to SRV which
eventually enables the wellbore to reach father away and produce more. One of
the key completion strategies should focus on creating larger SRVs with smaller
fracture spacing to maximize production from the pay-zone.
90
Chapter 6
Conclusions and Future Research
6.1 Conclusions
The overall motivation of this work was to develop an induced seismicity model
that incorporates physics of induced seismicity along current statistical knowl-
edge of earthquakes. Due to the large possible spatial extent of effective stress
change associated with wastewater disposal, the uncertainty involved in the hydro-
geological and fault characteristics of the study area requires the introduction of a
probabilistic method which considers a range of possible input parameters. To our
knowledge, our probabilistic-physics-based induced seismicity model presented in
Chapter 3, is the first and only induced seismicity model that includes probability
as well as physics of induced seismicity. In the next paragraphs, we briefly present
the highlights of the conducted research as they were discussed in Chapter 1.
Mohr-Coulomb shear failure criteria in different applications, namely wastewa-
ter disposal, hydraulic fracturing, production, and geothermal projects are studied,
and the differences in shear failure are highlighted.
We presented a deterministic approach to induced seismicity, in which based
on the knowledge of the state of stress, fault friction factor and trajectory, the
amount of pressure increase required for shear failure is calculated. The determin-
istic approach is suitable for fault stability analysis where a known fault is subject
to stress changes. We study induced seismicity in two field cases in California,
and Illinois basin. Simulating the pressure increase due to injection, we examined
91
the effect of permeability structure, operational parameters, and reservoir geom-
etry. We observed that injection layer width and permeability-contrast relative
to the crystalline basement influence the induced seismicity potential. Employing
the deterministic approach, we studied the difference between induced seismic-
ity response in waste water disposal versus secondary, and tertiary oil recovery.
Our results suggested that net injected fluid volume (i.e., that difference between
injected and produced volumes) is one of the key factors controlling induced seis-
micity potential. We observed that secondary and tertiary oil recovery, as well
as hydraulic fracturing, have less induced seismicity potential than waste water
disposal as a result of less net injected fluid volumes.
Accurate description of the relationship between injection parameters and
expected seismicity rates is a key component to seismic hazard assessment. We
developed a model for the probabilistic assessment of induced seismicity which
incorporates geomechanical analysis of the state of stress, pressure changes and
reservoir boundary conditions. We observe that a linear relationship between
injected volume and cumulative number of seismic events may not always hold
due to injection rate and flow boundary effects. Injection rate and flow bound-
ary effects may be more pronounced in long-term injections such as wastewater
disposal whereas a linear behavior may be expected in short-term, constant-rate
injections such as hydraulic fracturing. We validated our model against seismicity
data from German KTB and Paradox Valley saltwater disposal project, Colorado.
We observed that acceleration in seismicity rate at KTB despite an almost con-
stant injection rate can be explained due to the effects of reservoir boundaries on
pressure and expected seismicity. Through our model, we captured the deceler-
ation in induced seismicity in Paradox Valley, Colorado and associate it mainly
with 33% reduction in injection rate as well as 20 days occasional intervals of
92
shut-in. The differences in flow boundary conditions can drastically change the
seismic response to injection operations and for a bounded reservoir, induced seis-
micity rate increases exponentially. The proposed model may provide a pathway to
assess the expected seismic response to fluid injection as a function of operational
parameters and reservoir boundary conditions.
We introduced a new model for the most probable maximum magnitude seis-
micity due to wastewater injection. Our model is a physics-probabilistic-based
model which accounts for initial state of stress, reservoir flow and injection char-
acteristics. The inclusion of these important parameters in previous induced seis-
micity models in the literature is absent. By performing sensitivity analysis for the
effect of reservoir permeability, reservoir boundary conditions, injection rate and
injection schedule we observe that reservoir boundaries cause a non-linear behav-
ior on the log-log plot of number of induced seismicity events as well as semi-log
plot of the most probable maximum magnitude induced seismicity which results
in increasing the number and magnitude of the events.
We observed that increasing the flow rate and reservoir permeability will pro-
mote larger areas affected by waste-water disposal and in the case of critically
stressed faults close to failure results in an increase in the number and magnitude
of possible induced seismicity. Given the same volume of waste-water injection, we
observe that a reduced injection rate will results in less number of induced seis-
micity and smaller events compared to the same volume injected periodically with
intervals of shut-in but with a double injection rate. By using our induced seis-
micity model, injection well operators get the opportunity to calibrate the model
based on induced seismicity data and to easily try out and plan different injection
scenarios and mitigation scenarios for induced seismicity in order to reduce the
possibility of inducing a large earthquake.
93
Through simulation of different hydraulic fracture propagation we observed
that reservoir state of stress, pore-pressure, injection pressure, and presence of
existing faults and fractures are among the key parameters pertaining to SRV and
observed microseismicity. We observed very different SRVs activated as a result
of same hydraulic fracture but with different reservoir geomechanical conditions.
Reservoir deviatoric stress, fracture pressure and initial reservoir pore-pressure are
amongthekeyparameterspertainingtoSRVandobservedmicroseismicity. Higher
deviatoric stress shows to have highest effect in promoting shear failure in the
reservoir and a larger SRV. Higher reservoir pore-pressure and injection pressure
along with higher deviatoric stress create a suitable environment for shear failure
and extension of SRV and productive SRV. The ideal combination for maximizing
SRV and productive SRV is high fracture pressure, high initial reservoir pressure,
high stress difference between minimum and maximum reservoir stress. Presence
of natural fractures, faults, joints and geologic discontinuities are among the most
important contributors to SRV which eventually enables the wellbore to reach
father away and produce more. One of the key completion strategies should focus
on creating larger SRVs with smaller fracture spacing to maximize production from
the pay-zone.
6.2 Future Research
Based on the results obtained using the simulator that we developed in this
dissertation, we captured the effects of pore-pressure increase and stress concen-
tration on induced seismicity. We used a non-coupled solution for pore-pressure.
Next generation of induced seismicity models can include the poroelasticity and
coupled solution in order to capture this physics, too. An study on the effect
94
of similar waste-water disposal in a gas reservoir and liquid reservoir would be
interesting in the light of poroelasticity.
Geothermal reservoirs experience induced seismicity and microseismicity.
Developing a coupled thermal simulator which can capture the effect of tem-
perature changes on pore-pressure and stress change would give realistic results
in terms of reservoir and risk management for geothermal applications (Jha and
Juanes, 2014).
Datagatheringwouldbeanotherimportantfrontininducedseismicityresearch
as reduces uncertainty in the model results. Systematic reservoir characteriza-
tion techniques, including core analysis, porosimetry, permeability measurement,
knowledge about the state of stress, location of faults and the geology of the target
formation would be a great addition to any wastewater disposal practice.
95
Reference List
Adachi, J., E. Siebrits, A. Peirce, and J. Desroches (2007), Computer simulation
of hydraulic fractures, International Journal of Rock Mechanics and Mining Sci-
ences, 44(5), 739–757.
Ake, J., K. Mahrer, D. O’Connell, and L. Block (2005), Deep-injection and closely
monitored induced seismicity at Paradox Valley, Colorado, Bulletin of the Seis-
mological Society of America, 95(2), 664–683.
Albright, J. N., C. F. Pearson, et al. (1982), Acoustic emissions as a tool for
hydraulic fracture location: Experience at the Fenton Hill hot dry rock site,
Society of Petroleum Engineers Journal, 22(04), 523–530.
Anderson, E. M. (1951), The dynamics of faulting and dyke formation with appli-
cations to Britain, Hafner Pub. Co.
Block, L. V., C. K. Wood, W. L. Yeck, and V. M. King (2014), The 24 january
2013 ml 4.4 earthquake near paradox, colorado, and its relation to deep well
injection, Seismological Research Letters, 85(3), 609–624.
Block, L. V., C. K. Wood, W. L. Yeck, and V. M. King (2015), Induced seis-
micity constraints on subsurface geological structure, paradox valley, colorado,
Geophysical Journal International, 200(2), 1170–1193.
96
Bremkamp, W., and C. Harr (1988), Area of least resistance to fluid movement and
pressure rise, Paradox Valley Unit, Salt Brine Injection Project, Bedrock, Col-
orado, a report prepared for the US Bureau of Reclamation, Denver, Colorado,
39.
Brown, W., and C. Frohlich (2013), Investigating the cause of the 17 may 2012 m
4.8 earthquake near timpson, east texas (abstr.), Seismol. Res. Lett, 84, 374.
Brudy, M., M. Zoback, K. Fuchs, F. Rummel, and J. Baumgärtner (1997), Estima-
tion of the complete stress tensor to 8 km depth in the ktb scientific drill holes:
Implications for crustal strength, Journal of Geophysical Research: Solid Earth,
102(B8), 18,453–18,475.
Carslaw, H. S., and J. C. Jaeger (1959), Conduction of heat in solids, Oxford:
Clarendon Press, 1959, 2nd ed.
Commission, U. N. R., et al. (2012), Technical report on central and eastern
united states seismic source characterization for nuclear facilities: Washington,
dcandpaloalto, california, Electric Power Research Institute, US Department of
Energy, and US Nuclear Regulatory Commission joint research report, variously
paginated, available at http://pbadupws. nrc. gov/docs/ML1204/ML12048A804.
pdf.
Cornell, C. A. (1968), Engineering seismic risk analysis, Bulletin of the Seismolog-
ical Society of America, 58(5), 1583–1606.
Culham, W., et al. (1974), Pressure buildup equations for spherical flow regime
problems, Society of Petroleum Engineers Journal, 14(06), 545–555.
97
Dempsey, D., and J. Suckale (2016), Collective properties of injection-induced
earthquake sequences: 1. model description and directivity bias, Journal of Geo-
physical Research: Solid Earth.
Desroches, J., E. Detournay, B. Lenoach, P. Papanastasiou, J. Pearson,
M. Thiercelin, and A. Cheng (1994), The crack tip region in hydraulic frac-
turing, in Proceedings of the Royal Society of London A: Mathematical, Physical
and Engineering Sciences, vol. 447, pp. 39–48, The Royal Society.
Ellsworth, W. L. (2013), Injection-induced earthquakes, Science, 341(6142),
1225,942.
Elst, N. J., M. T. Page, D. A. Weiser, T. H. Goebel, and S. M. Hosseini (2016),
Induced earthquake magnitudes are as large as (statistically) expected, Journal
of Geophysical Research: Solid Earth, 121(6), 4575–4590.
Emmermann, R., and J. Lauterjung (1997), The German continental deep drilling
program KTB: overview and major results, Journal of Geophysical Research:
Solid Earth, 102(B8), 18,179–18,201.
Evans, D.M.(1966), Thedenierareaearthquakes&idtheRockyMountainArsenal
disposal well, Mountain Geologist, 3.
Frohlich, C. (2012), Two-year survey comparing earthquake activity and injection-
well locations in the Barnett Shale, Texas, Proceedings of the National Academy
of Sciences, 109(35), 13,934–13,938.
Frohlich, C., W. Ellsworth, W. A. Brown, M. Brunt, J. Luetgert, T. MacDonald,
and S. Walter (2014), The 17 may 2012 m4. 8 earthquake near Timpson, East
Texas: An event possibly triggered by fluid injection, Journal of Geophysical
Research: Solid Earth, 119(1), 581–593.
98
Goebel, T., S. Hosseini, F. Cappa, E. Hauksson, J. Ampuero, F. Aminzadeh, and
J. Saleeby (2016), Wastewater disposal and earthquake swarm activity at the
Southern end of the Central Valley, California, Geophysical Research Letters,
43(3), 1092–1099.
Gutenberg, B., and C. F. Richter (1944), Frequency of earthquakes in California,
Bulletin of the Seismological Society of America, 34(4), 185–188.
Hall, H. N., et al. (1953), Compressibility of reservoir rocks, Journal of Petroleum
Technology, 5(01), 17–19.
Healy, J., W. Rubey, D. Griggs, and C. Raleigh (1968), The Denver earthquakes,
Science, 161(3848), 1301–1310.
Herrmann, R. B., S.-K. Park, and C.-Y. Wang (1981), The Denver earthquakes of
1967-1968, Bulletin of the Seismological Society of America, 71(3), 731–745.
Hornbach, M. J., H. R. DeShon, W. L. Ellsworth, B. W. Stump, C. Hayward,
C. Frohlich, H. R. Oldham, J. E. Olson, M. B. Magnani, C. Brokaw, et al.
(2015), Causal factors for seismicity near Azle, Texas, Nature communications,
6.
Horton, S. (2012), Disposal of hydrofracking waste fluid by injection into subsur-
face aquifers triggers earthquake swarm in central Arkansas with potential for
damaging earthquake, Seismological Research Letters, 83(2), 250–260.
Hosseini, S. M. (2012), Hydraulic fracture mechanism in unconsolidated forma-
tions, M.Sc. Thesis, The University of Texas at Austin.
99
Hosseini, S. M., F. Aminzadeh, et al. (2016), Effects of the earth characteristics
on induced seismicity potential, in SPE Western Regional Meeting, Society of
Petroleum Engineers.
Hosseini, S. M., et al. (2013), On the linear elastic fracture mechanics application
in Barnett shale hydraulic fracturing, in 47th US rock mechanics/geomechanics
symposium, American Rock Mechanics Association.
Hsieh, P. A., and J. D. Bredehoeft (1981), A reservoir analysis of the Denver
earthquakes: A case of induced seismicity, Journal of Geophysical Research:
Solid Earth, 86(B2), 903–920.
Hubbert, M. K., and W. W. Rubey (1959), Role of fluid pressure in mechanics of
overthrust faulting i. mechanics of fluid-filled porous solids and its application
to overthrust faulting, Geological Society of America Bulletin, 70(2), 115–166.
Ingraffea, A. R., K. L. Gunsallus, J. F. Beech, and P. P. Nelson (1984), A short-
rod based system for fracture toughness testing of rock, in Chevron-notched
specimens: testing and stress analysis, ASTM International.
Irwin, G. R. (1957), Analysis of stresses and strains near the end of a crack travers-
ing a plate, Journal of Applied Mechanics, 24(3), 361–364..
Jha, B., and R. Juanes (2014), Coupled multiphase flow and poromechanics: A
computational model of pore pressure effects on fault slip and earthquake trig-
gering, Water Resources Research, 50(5), 3776–3808.
Keranen, K. M., H. M. Savage, G. A. Abers, and E. S. Cochran (2013), Potentially
inducedearthquakesinOklahoma,USA:Linksbetweenwastewaterinjectionand
the 2011 mw 5.7 earthquake sequence, Geology, 41(6), 699–702.
100
Kim, W.-Y. (2013), Induced seismicity associated with fluid injection into a deep
wellinYoungstown, Ohio, Journal of Geophysical Research: Solid Earth, 118(7),
3506–3518.
King, V. M., L. V. Block, W. L. Yeck, C. K. Wood, and S. A. Derouin (2014),
Geological structure of the paradox valley region, colorado, and relationship to
seismicity induced by deep well injection, Journal of Geophysical Research: Solid
Earth, 119(6), 4955–4978.
King, V. M., L. V. Block, and C. K. Wood (2016), Pressure/flow modeling and
induced seismicity resulting from two decades of high-pressure deep-well brine
injection, Paradox Valley, Colorado, Geophysics, 81(5), B119–B134.
Krogstad, S., K.-A.Lie, O.Møyner, H.M.Nilsen, X.Raynaud, B.Skaflestad, etal.
(2015), Mrst-ad–an open-source framework for rapid prototyping and evaluation
of reservoir simulation problems, in SPE reservoir simulation symposium, Soci-
ety of Petroleum Engineers.
Kümpel, H.-J., J. Erzinger, and S. A. Shapiro (2006), Two massive hydraulic tests
completed in deep KTB pilot hole, Scientific Drilling, 3, 40–42.
Labuz, J., S. Shah, and C. Dowding (1987), The fracture process zone in granite:
evidence and effect, in International Journal of Rock Mechanics and Mining
Sciences & Geomechanics Abstracts, vol. 24, pp. 235–246, Elsevier.
Langenbruch, C., and M. D. Zoback (2016), How will induced seismicity in Okla-
homa respond to decreased saltwater injection rates?, Science advances, 2(11),
e1601,542.
Lee, J., J. B. Rollins, and J. P. Spivey (2003), Pressure transient testing, vol. 9,
Henry L. Doherty Memorial Fund of Aime Society of Petroleum.
101
Manzocchi, T., A. Heath, B. Palananthakumar, C. Childs, and J. Walsh (2008),
Faults in conventional flow simulation models: a consideration of representa-
tional assumptions and geological uncertainties, Petroleum Geoscience, 14(1),
91–110.
Matthews, C. S., and D. G. Russell (1967), Pressure buildup and flow tests in wells,
vol. 1, Society of Petroleum Engineers of AIME Dallas, TX.
Matthews, J. D., J. N. Carter, K. D. Stephen, R. W. Zimmerman, A. Skorstad,
T. Manzocchi, and J. A. Howell (2008), Assessing the effect of geological uncer-
tainty on recovery estimates in shallow-marine reservoirs: the application of
reservoir engineering to the SAIGUP project, Petroleum Geoscience, 14(1), 35–
44.
Mayerhofer, M. J., E. Lolon, N. R. Warpinski, C. L. Cipolla, D. W. Walser, C. M.
Rightmire, et al. (2010), What is stimulated reservoir volume?, SPE Production
& Operations, 25(01), 89–98.
McGarr, A. (2014), Maximum magnitude earthquakes induced by fluid injection,
Journal of Geophysical Research: solid earth, 119(2), 1008–1019.
McGuire, R. K. (2004), Seismic hazard and risk analysis, Earthquake Engineering
Research Institute.
McNamara, D. E., H. M. Benz, R. B. Herrmann, E. A. Bergman, P. Earle, A. Hol-
land, R. Baldwin, and A. Gassner (2015), Earthquake hypocenters and focal
mechanisms in Central Oklahoma reveal a complex system of reactivated sub-
surface strike-slip faulting, Geophysical Research Letters, 42(8), 2742–2749.
Moos, D., and M. D. Zoback (1990), Utilization of observations of well bore failure
to constrain the orientation and magnitude of crustal stresses: application to
102
continental, deep sea drilling project, and ocean drilling program boreholes,
Journal of Geophysical Research: Solid Earth, 95(B6), 9305–9325.
Olson, J. E., J. Holder, M. Hosseini, et al. (2011), Soft rock fracturing geometry
and failure mode in lab experiments, in SPE Hydraulic Fracturing Technology
Conference, Society of Petroleum Engineers.
Ouchterlony, F. (1982), Review of fracture toughness testing of rock, SM archives,
7, 131–211.
Ouchterlony, F. (1989), On the background to the formulae and accuracy of
rock fracture toughness measurements using ISRM standard core specimens, in
International Journal of Rock Mechanics and Mining Sciences & Geomechanics
Abstracts, vol. 26, pp. 13–23, Elsevier.
Papanastasiou, P. (1997), The influence of plasticity in hydraulic fracturing, Inter-
national Journal of Fracture, 84(1), 61–79.
Person, M., A. Banerjee, J. Rupp, C. Medina, P. Lichtner, C. Gable, R. Pawar,
M.Celia, J.McIntosh, andV.Bense(2010), Assessmentofbasin-scalehydrologic
impacts of CO2 sequestration, Illinois basin, International Journal of Green-
house Gas Control, 4(5), 840–854.
Petersen, M. D., M. P. Moschetti, P. M. Powers, C. S. Mueller, K. M. Haller, A. D.
Frankel, Y. Zeng, S. Rezaeian, S. C. Harmsen, O. S. Boyd, et al. (2015), The
2014 united states national seismic hazard model, Earthquake Spectra, 31(S1),
S1–S30.
Petersen, M. D., C. S. Mueller, M. P. Moschetti, S. M. Hoover, A. L. Llenos, W. L.
Ellsworth, A. J. Michael, J. L. Rubinstein, A. F. McGarr, and K. S. Rukstales
103
(2016), 2016 one-year seismic hazard forecast for the central and eastern united
states from induced and natural earthquakes, Tech. rep., US Geological Survey.
Pollard, D. D. (1987), The theoretical displacements and stresses near fractures
in rock: with applications to faults, joints, veins, dikes, and solution surfaces.,
Fracture mechanics of rock, pp. 277–349.
Prats, M., et al. (1961), Effect of vertical fractures on reservoir behavior-
incompressible fluid case, Society of Petroleum Engineers Journal, 1(02), 105–
118.
Raleigh, C., J. Healy, and J. Bredehoeft (1976a), An experiment in earthquake
control at Rangely, Colorado, work (Fig. Ib), 108(52), 30.
Report, N. (2013), Induced seismicity potential in energy technologies, National
Academies Press.
Rubinstein, J. L., and A. B. Mahani (2015), Myths and facts on wastewater injec-
tion, hydraulic fracturing, enhanced oil recovery, and induced seismicity, Seis-
mological Research Letters, 86(4), 1060–1067.
Rudnicki, J. W. (1986), Fluid mass sources and point forces in linear elastic diffu-
sive solids, Mechanics of Materials, 5(4), 383–393.
Scheirer, A. H. (2007), Petroleum systems and geologic assessment of oil and gas in
the San Joaquin Basin Province, California, Tech. rep., U.S. Geological Survey.
Schmidt, R., and C. Huddle (1977), Effect of confining pressure on fracture tough-
ness of Indiana limestone, in International Journal of Rock Mechanics and Min-
ing Sciences & Geomechanics Abstracts, vol. 14, pp. 289–293, Elsevier.
104
Schmidt, R., and T. Lutz (1979), K ic and j ic of westerly granite?effects of thick-
ness and in-plane dimensions, in Fracture mechanics applied to brittle materials,
ASTM International.
Schultz,R.,G.Atkinson,D.Eaton,Y.Gu,andH.Kao(2018),Hydraulicfracturing
volumeisassociatedwithinducedearthquakeproductivityintheDuvernayplay,
Science, 359(6373), 304–308.
Segall, P., and S. Lu (2015), Injection-induced seismicity: Poroelastic and earth-
quake nucleation effects, Journal of Geophysical Research: Solid Earth, 120(7),
5082–5103.
Shapiro, S., J. Kummerow, C. Dinske, G. Asch, E. Rothert, J. Erzinger, H.-J.
Kümpel, and R. Kind (2006a), Fluid induced seismicity guided by a continental
fault: Injection experiment of 2004/2005 at the German deep drilling site (ktb),
Geophysical Research Letters, 33(1).
Shapiro, S., C. Dinske, and J. Kummerow (2007), Probability of a given-magnitude
earthquake induced by a fluid injection, Geophysical research letters, 34(22).
Shapiro, S. A., E. Huenges, and G. Borm (1997), Estimating the crust permeabil-
ity from fluid-injection-induced seismic emission at the KTB site, Geophysical
Journal International, 131(2), F15–F18.
Shapiro, S. A., C. Dinske, C. Langenbruch, and F. Wenzel (2010a), Seismogenic
index and magnitude probability of earthquakes induced during reservoir fluid
stimulations, The Leading Edge, 29(3), 304–309.
Shapiro, S. A., C. Dinske, C. Langenbruch, and F. Wenzel (2010b), Seismogenic
index and magnitude probability of earthquakes induced during reservoir fluid
stimulations, The Leading Edge, 29(3), 304–309.
105
Sneddon, I. (1946), The distribution of stress in the neighbourhood of a crack in
an elastic solid, in Proceedings of the Royal Society of London A: Mathematical,
Physical and Engineering Sciences, vol. 187, pp. 229–260, The Royal Society.
Sone, H., M. Zoback, et al. (2010), Strength, creep and frictional properties of gas
shalereservoirrocks, in44th US Rock Mechanics Symposium and 5th US-Canada
Rock Mechanics Symposium, American Rock Mechanics Association.
Stober, I., and K. Bucher (2005), The upper continental crust, an aquifer and
its fluid: hydaulic and chemical data from 4 km depth in fractured crystalline
basement rocks at the ktb test site, Geofluids, 5(1), 8–19.
Thiercelin, M. (1989), Fracture toughness and hydraulic fracturing, in Inter-
national Journal of Rock Mechanics and Mining Sciences & Geomechanics
Abstracts, vol. 26, pp. 177–183, Elsevier.
Townend, J., and M. D. Zoback (2000), How faulting keeps the crust strong, Geol-
ogy, 28(5), 399–402.
Valko, P., and M. J. Economides (1995), Hydraulic fracture mechanics, vol. 318,
Wiley Chichester.
Vermylen, J. P. (2011), Geomechanical studies of the Barnett shale, Texas, USA,
Stanford University.
Walsh, F. R., and M. D. Zoback (2015), Oklahoma’s recent earthquakes and salt-
water disposal, Science advances, 1(5), e1500,195.
Walsh, F. R., and M. D. Zoback (2016), Probabilistic assessment of potential
fault slip related to injection-induced earthquakes: Application to north-central
Oklahoma, USA, Geology, 44(12), 991–994.
106
Warpinski, N., S. Wolhart, C. Wright, et al. (2001), Analysis and prediction of
microseismicity induced by hydraulic fracturing, in SPE Annual Technical Con-
ference and Exhibition, Society of Petroleum Engineers.
Weingarten, M., S. Ge, J. W. Godt, B. A. Bekins, and J. L. Rubinstein (2015),
High-rate injection is associated with the increase in U.S. Mid-continent seis-
micity, Science, 348(6241), 1336–1340.
WellTest, F. (), Version 2012, Calgary, Alberta: Fekete Associates Inc.
Wesson, R. L., and C. Nicholson (1987), Earthquake hazard associated with deep
well injection; a report to the U.S. Environmental Protection Agency, Tech. rep.,
US Geological Survey,.
Westergaard, H. M. (1939), Bearing pressures and cracks, Journal of applied
mechanics, 6(2), 49–53.
Williams, M. (1957), The bending stress distribution at the base of a stationary
crack, J Appl Mech, 24, 109–14.
Yeck, W. L., L. V. Block, C. K. Wood, and V. M. King (2015), Maximum mag-
nitude estimations of induced earthquakes at Paradox Valley, Colorado, from
cumulative injection volume and geometry of seismicity clusters, Geophysical
Journal International, 200(1), 322–336.
Zhang, Y., M. Person, J. Rupp, K. Ellett, M. A. Celia, C. W. Gable, B. Bowen,
J.Evans, K.Bandilla, P.Mozley, etal.(2013), Hydrogeologiccontrolsoninduced
seismicity in crystalline basement rocks due to fluid injection into basal reser-
voirs, Groundwater, 51(4), 525–538.
Zoback, M. D. (2010), Reservoir geomechanics, Cambridge University Press.
107
Zoback, M. D., and H.-P. Harjes (1997), Injection-induced earthquakes and crustal
stress at 9 km depth at the KTB deep drilling site, Germany, Journal of Geo-
physical Research: Solid Earth, 102(B8), 18,477–18,491.
108

Abstract (if available)

Abstract The substantial increase in earthquake rates in the Central and Eastern United States starting in 2009, has raised the question of a possible causal relationship between wastewater disposal and earthquakes. This dissertation focuses on the study of induced seismicity through probabilistic and deterministic modeling of involved physical mechanisms and validating theoretically predicted seismicity rate changes with observations. An accurate description of such rate changes is a key component in seismic hazard assessment at the proximity of geothermal reservoirs, carbon sequestration, hydraulic fracturing and wastewater disposal sites. This dissertation meets the need for a representative model and validated simulator for injection-induced seismicity. ❧ We categorize current induced seismicity models into two broad classes of deterministic and probabilistic models. Deterministic models simulate pressure for fault stability of known faults while probabilistic models focus on the probability of inducing earthquakes on a random distribution of faults. Both classes of existing models do not consider uncertainty in faults stresses and possible reservoir boundary conditions due to different geological settings. In this dissertation, in addition to deterministic induced seismicity modeling, we develop a probabilistic-physics-based induced seismicity model that captures the effects of uncertainty in faults stresses, reservoir properties and boundary conditions. ❧ We developed a probabilistic-physics-based model which considers a randomly distributed network of faults that undergo frictional failure due to pressure changes. We include uncertainty in faults stresses through parameterization of required pressure increase for failure, based on probability distributions of tectonic stresses using Mohr-Coulomb failure criterion. Reservoir boundary effects are explored by including pressure solutions for early-time (i.e., infinite-acting) and late-time (i.e., pseudo-steady-state) periods. We sample the probability density function for required pressure-increase for failure through simulated pressure in the domain and calculate the probability of inducing an earthquake at each point. To obtain the total number of possible earthquakes, we integrate the calculated probability over the study domain and multiply it by the fault density. Using frequency-magnitude distribution laws for induced earthquakes (e.g., Gutenberg-Richter law), we calculate the cumulative number of earthquakes above a specific magnitude. Using deterministic numerical simulation of pressure, we study several induced seismicity cases in the U.S. and Europe. ❧ We show that the differences in flow boundary conditions can drastically change the seismic response to injection operations. We observe an exponential seismicity rate increase in the case of a bounded reservoir versus a linear increase for an unbounded reservoir. Our deterministic modeling results suggest that net injected fluid volume and boundary conditions are the primary key factor, controlling induced seismicity. We conclude that secondary and tertiary oil recovery may experience less induced seismicity than wastewater disposal, likely as a result of less net injected fluid volumes. ❧ We observe that a linear relationship between injected volume and cumulative number of seismic events may not always hold due to injection rate and flow boundary effects. Injection rate and flow boundary effects may be more pronounced in long-term injections such as wastewater disposal whereas a linear behavior may be expected in short-term, constant-rate injections such as hydraulic fracturing. We validate our model against seismicity data from German KTB and Paradox Valley saltwater disposal project, Colorado. The differences in flow boundary conditions can drastically change the seismic response to injection operations and for a bounded reservoir, induced seismicity rate increases exponentially. ❧ Maximum magnitude induced seismicity (MMIS) is one of the critical parameters in earthquake hazard estimation. Using the developed induced seismicity model in this dissertation, we study the relationship between operational injection parameters and MMIS. Current MMIS models in the literature do not include the effects of injection rates, permeability, and reservoir boundaries. We developed a probabilistic-physics-based induced seismicity model for MMIS, which includes the operational and reservoir parameters. Our model is a physics-probabilistic-based model which accounts for initial state of stress, reservoir flow and injection characteristics. We conclude that reservoir permeability, injection rate, and the boundary conditions can drastically change expected MMIS, therefore some of the conclusions and models in the literature based on very different cases of wastewater disposal with very different characteristics may not hold. Instead, we suggest studying each case on its own due to differences in all the key parameters surrounding MMIS. ❧ We study induced seismicity and microseismicity due to hydraulic fracturing. Besides the pressure increase due to diffusive processes, stress shadow around the fracture is capable of changing the effective stress in favor of shear failure in preexisting faults. We developed a simulator to propagate a hydraulic fracture through a porous medium and capture shear failure among a random distribution of faults in the domain to study statistics of induced seismicity associated with hydraulic fracturing. We observed very different stimulated rock volume (SRV) activated as a result of same hydraulic fracture but with different reservoir geomechanical conditions. Reservoir deviatoric stress, fracture pressure, existing faults and initial reservoir pore-pressure are among the key parameters pertaining to SRV and observed microseismicity. The ideal combination for maximizing SRV and productive SRV is high fracture pressure, high initial reservoir pressure, high stress difference between minimum and maximum reservoir stress. One of the key completion strategies should focus on creating larger SRVs with smaller fracture spacing to maximize production from the pay zone.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Reactivation of multiple faults in oilfields with injection and production

PDF

Continuum modeling of reservoir permeability enhancement and rock degradation during pressurized injection

$An extended finite element method based modeling of hydraulic fracturing$

PDF

An extended finite element method based modeling of hydraulic fracturing

PDF

Integrated reservoir characterization for unconventional reservoirs using seismic, microseismic and well log data

PDF

A coupled multiphase flow-geomechanics framework for modeling material, kinematic and contact nonlinearities in faulted reservoirs

$Integrated workflow for characterizing and modeling fracture network in unconventional reservoirs using microseismic data$

PDF

Integrated workflow for characterizing and modeling fracture network in unconventional reservoirs using microseismic data

$A coupled poromechanics-damage mechanics approach to model fracturing in multiphase flow$

PDF

A coupled poromechanics-damage mechanics approach to model fracturing in multiphase flow

PDF

Dynamics of CO₂ leakage through faults

PDF

Integration of multi-physics data into dynamic reservoir model with ensemble Kalman filter

PDF

Modeling and simulation of complex recovery processes

PDF

Geothermal resource assessment and reservoir modeling with an application to Kavaklidere geothermal field, Turkey

PDF

Stress and deformation analysis on fluid-exposed reservoir rocks

$Synergistic coupling between geomechanics, flow, and transport in fractured porous media: applications in hydraulic fracturing and fluid mixing$

PDF

Synergistic coupling between geomechanics, flow, and transport in fractured porous media: applications in hydraulic fracturing and fluid mixing

PDF

A study of diffusive mass transfer in tight dual-porosity systems (unconventional)

$Optimization of CO2 storage efficiency under geomechanical risks using coupled flow-geomechanics-fracturing model$

PDF

Optimization of CO2 storage efficiency under geomechanical risks using coupled flow-geomechanics-fracturing model

Asset Metadata

Creator Hosseini, Seyed Mehran (author)

Core Title Assessing induced seismicity rate increase based on deterministic and probabilistic modeling

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Petroleum Engineering

Publication Date 02/12/2018

Defense Date 01/10/2018

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag earthquake,fault stability,induced seismicity,injection,OAI-PMH Harvest,reservoir,shear failure,wastewater disposal

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Aminzadeh, Fred (committee chair), Ershaghi, Iraj (committee member), Sammis, Charles (committee member)

Creator Email mehranh@utexas.edu,seyedmeh@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c40-473245

Unique identifier UC11266697

Identifier etd-HosseiniSe-6027.pdf (filename),usctheses-c40-473245 (legacy record id)

Legacy Identifier etd-HosseiniSe-6027.pdf

Dmrecord 473245

Document Type Dissertation

Rights Hosseini, Seyed Mehran

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA