U.S. Department of the Interior

Bureau of Reclamation May 2016

Technical Memorandum TM-86-68330-2013-12

The ML 4.4 Earthquake of January 

24, 2013, Near Paradox, Colorado, 

and Implications for Near-term 

Injection Operations

Colorado Basin Salinity Control Project,

Paradox Valley Unit, Colorado

Upper Colorado Region
Mission Statements

The mission of the Department of the Interior is to protect and 

provide access to our Nation’s natural and cultural heritage and 

honor our trust responsibilities to Indian Tribes and our 

commitments to island communities.

The mission of the Bureau of Reclamation is to manage, develop, 

and protect water and related resources in an environmentally and 

economically sound manner in the interest of the American public.
BUREAU OF RECLAMATION 

Technical Service Center, Denver, Colorado 

Seismotectonics and Geophysics Group, 86-68330 

Technical Memorandum TM-86-68330-2013-12 

The ML 4.4 Earthquake of January 24, 

2013, Near Paradox, Colorado, and 

Implications for Near-term Injection 

Operations 

Colorado Basin Salinity Control Project, 

Paradox Valley Unit, Colorado 

Upper Colorado Region 

Prepared by: 

c-1267/6 

Christopher K. Wood, Geophysicist Date 

Lisa V. Block, Geophysicist Date 

Vdnelss )ng, eophysicist Date 

William L. Yeck Geophysicist Date 

Peer Review Certification 

This report has been reviewed and is believed to be in accordance with the service 

agreement and standards of the profession. 

Peer reviewed by: 

5126-/2cii( 

Pertgcheng Liu, Geophysicist Date 

Technical Memorandum TM-68330-2011-39

iii

Contents

Page

Executive Summary..........................................................................................................1

1 Introduction.................................................................................................................3

2 Description of the Earthquake ..................................................................................7

2.1 Earthquake Magnitude ....................................................................................7

2.1.1 Magnitude Scales......................................................................................9

2.2 Location ...........................................................................................................9

2.3 Focal Mechanism ..........................................................................................11

2.4 Nearby Seismicity..........................................................................................12

2.5 Strong Ground Motions .................................................................................18

2.6 Ground Motion Effects and Standards ..........................................................26

2.7 Triggering Mechanism ...................................................................................30

3 Seismic Hazard Factors ...........................................................................................33

3.1 Geographical Distribution of Seismicity.........................................................33

3.1.1 Relation of Seismicity to Inhabited Areas in Paradox Valley ..................33

3.1.2 Volumetric Evolution of Seismicity..........................................................38

3.2 Seismicity Rates and Magnitude Distribution................................................43

3.2.1 Seismicity Rates......................................................................................44

3.2.2 Gutenberg-Richter Magnitude Distribution .............................................45

3.2.3 Seismogenic Index..................................................................................63

3.3 Maximum Magnitude of Induced Earthquakes..............................................68

3.3.1 Geometrical scaling of maximum magnitude..........................................68

3.3.2 Fluid volume scaling of maximum magnitude.........................................73

3.3.3 Effect of Pore Pressure on Earthquake Magnitude ................................76

3.3.4 Maximum magnitude discussion.............................................................78

3.4 Short-term Earthquake Recurrence Estimates..............................................79

3.5 Discussion .....................................................................................................83

4 Injection Pressures and Flow Rates.......................................................................85

4.1 Injection History .............................................................................................85

4.1.1 Phase I (July 22, 1996 – July 25, 1999) .................................................85

4.1.2 Phase II (July 26, 1999 – June 22, 2000)...............................................87

4.1.3 Phase III (June 23, 2000 – January 6, 2002)..........................................87

4.1.4 Phase IV (January 7, 2002 – present) ....................................................87

4.2 Previous Work ...............................................................................................88

4.3 Fitting of Pressure/Flow Data to Idealized Models........................................89

4.4 Results of Pressure/Flow Model Fitting.........................................................95

4.4.1 Changes in Model Parameters over Time ..............................................95

4.4.2 Model Pressure-History Plots and Type Curves...................................100

4.5 Evaluation of Pressures for Alternative Injection Scenarios........................110

4.6 Pressure-Flow Modeling Discussion ...........................................................116

5 Discussion...............................................................................................................117

6 Conclusions ............................................................................................................123

6.1 Analysis of the January 24, 2013 ML 4.4 Earthquake .................................123
iv

6.2 Potential Future Seismicity..........................................................................123

6.2.1 Short-Term Seismicity...........................................................................123

6.2.2 Seismic Hazard Factors........................................................................124

6.3 Injection Operations.....................................................................................127

References.....................................................................................................................129

Appendix A Earthquake Catalog ............................................................................ A-1

A.1 Earthquake Location Methods..................................................................... A-2

A.2 Description of Data Columns....................................................................... A-3

Tables

Page

Table 2-1. Comparison of local and moment magnitudes for the 3 largest PVU-induced 

earthquakes. .......................................................................................................................8

Table 2-2. Observed ground motions and distances for the Paradox Community Center 

(PVCC), the injection well (PVPP), and the surface treatment facility (PVEF).................20

Table 2-3. Input parameters used for estimating ground motions from the NGA equations. 

See Kaklamanos and others (2011) for a description of the parameters. ........................21

Table 2-4. Summary of current ground vibration limits from blasting set by the Office of 

Surface Mining, Reclamation and Enforcement, U. S. Department of the Interior; taken 

from Title 30, Subchapter K, of the Code of Federal Regulations, dated July 1, 2012. ...28

Table 3-1. Estimated probability of occurrence of earthquakes over the next 3 years for 

the near-well area. ............................................................................................................82

Table 3-2. Estimated probability of occurrence of earthquakes over the next 3 years for 

the NW cluster area. .........................................................................................................83

Figures

Page

Figure 1-1. Location of the Jan. 24, 2013 ML 4.4 earthquake relative to the Paradox 

Valley Unit (PVU) injection well (yellow star), PVU extraction well field (PVEF), and the 

town of Paradox, Colorado (PVCC). ...................................................................................3

Figure 1-2. PVU-induced earthquakes having well-constrained relative hypocenters, 

plotted as a function of distance from the PVU injection well and depth (relative to the 

ground surface elevation at the wellhead) ..........................................................................4

Figure 1-3. Correlation between injection flow rates and shallow seismicity (< 8.5 km 

depth) recorded by the Paradox Valley Seismic Network (PVSN) .....................................5

Figure 1-4. Location of Paradox Valley seismograph stations. ..........................................6
Technical Memorandum TM-68330-2011-39

v

Figure 2-1. Example of automatic spectral fitting procedure used to estimate seismic 

moment, showing selected displacement amplitude spectra for P, S, and S-wave coda at 

station PV01 (42 km distance) and PV07 (33 km distance) ...............................................8

Figure 2-2. Location of PVU-induced seismicity in the near-well region (within 5 km of the 

injection well) and northwest cluster (6 to 8 km from the well). ........................................10

Figure 2-3. First motions and computed focal mechanism of the January 2013 event, and 

map showing the distribution of previously induced earthquakes in the vicinity of the 

January 2013 event and the locations of the January 2013 earthquake and aftershocks 

observed through July 2013..............................................................................................12

Figure 2-4. Map showing the distribution of foreshocks and aftershocks of the Jan. 24, 

2013 event ........................................................................................................................15

Figure 2-5. Occurrence of earthquakes within 1.0 km of the January 24, 2013 event 

epicenter, plotted as a function of date and earthquake magnitude.................................16

Figure 2-6. Epicenters of earthquakes with well-constrained relative locations in the near 

vicinity of the January 2013 ML 4.4 earthquake. (a) Earthquakes having similar strike-slip 

focal mechanisms are indicated by the filled red circles, and those having similar oblique 

normal focal mechanisms are indicated by the filled green circles. (b) Earthquakes colorcoded by year of occurrence and with symbol size scaled by event magnitude. .............17

Figure 2-7. Acceleration ground motion time histories recorded by strong-motion 

instruments located at the community center in the town of Paradox, Colorado (PVCC). 

Traces are offset for clarity................................................................................................18

Figure 2-8. Acceleration ground motion time histories recorded by strong-motion 

instruments located near the PVU injection well and pumping plant (PVPP)...................19

Figure 2-9: Acceleration ground motions recorded by strong-motion instruments located 

at the PVU surface treatment facility and extraction-well field (PVEF).............................19

Figure 2-10. Acceleration response spectra for the Paradox Community Center (PVCC), 

and comparison with NGA results.....................................................................................22

Figure 2-11. Velocity response spectra for Paradox Community Center (PVCC), and 

comparison with NGA results............................................................................................22

Figure 2-12. Acceleration response spectra for the injection well (PVPP), and 

comparison with NGA results............................................................................................23

Figure 2-13. Velocity response spectra for the injection well (PVPP), and comparison 

with NGA results. ..............................................................................................................23

Figure 2-14. Acceleration response spectra for the surface-treatment facility (PVEF), and 

comparison with NGA results............................................................................................24

Figure 2-15. Velocity response spectra for surface-treatment facility (PVEF), and 

comparison with NGA results............................................................................................24

Figure 2-16. Comparison of observed PHA values versus distance, with NGA results. ..25

Figure 2-17. Comparison of observed 1.0-second spectral acceleration values versus 

distance, with NGA results................................................................................................25

Figure 2-18. Comparison of observed peak velocity values versus distance, with NGA 

results................................................................................................................................26

Figure 2-19. Comparison of measured particle velocities at station PVCC (Town of 

Paradox, community center) and the blasting standard of 0.75 in/sec peak particle 

velocity (orange dashed lines) ..........................................................................................29
vi

Figure 2-20. Comparison of measured particle velocities at station PVPP (near the 

injection well-head) and the blasting standard of 0.75 in/sec peak particle velocity (orange 

dashed lines).....................................................................................................................29

Figure 2-21. Comparison of measured particle velocities at station PVEF (Surface 

Treatment Facility and brine extraction) and the blasting standard of 0.75 in/sec peak 

particle velocity (orange dashed lines)..............................................................................30

Figure 2-22: Injection downhole pressure data averaged over daily, 6-month, 18-month, 

and 30-month time periods (top), and occurrence of induced seismicity as a function of 

time and magnitude within 5 km of the injection well (middle) and at distances of 5 to 10 

km from the well (bottom)..................................................................................................32

Figure 3-1. Maps showing the geographical distribution of shallow seismicity recorded in 

the Paradox Valley area from 1991 through 2000: (a) injection tests, 1991-1995 (b) 

continuous injection, 1996-2000 .......................................................................................35

Figure 3-2. Maps showing the geographical distribution of shallow seismicity recorded in 

the Paradox Valley area from 2001 through 2012: (a) continuous injection, 2001-2008 (b) 

continuous injection, 2009-2012 .......................................................................................36

Figure 3-3. Occurrence of PVU-induced seismicity as a function of date, distance from 

the town of Paradox, and earthquake magnitude .............................................................37

Figure 3-4. Example of minimum ellipsoid fit to hypocenters (grey dots) from the near-well 

cluster................................................................................................................................40

Figure 3-5. Map showing the surface projections of the minimum ellipsoids fitting the 

near-well and NW clusters ................................................................................................41

Figure 3-6. Evolution of seismicity ellipsoids for alternative selections of hypocenters 

representing the near-well and the NW clusters...............................................................42

Figure 3-7. Seismicity rates as a function of time for earthquakes occurring at distances 

of less than 5 km, 5 to 10 km, and greater than 10 km from the injection well.................45

Figure 3-8. B-values estimated by two maximum-likelihood methods for the earthquakes

..........................................................................................................................................49

Figure 3-9. B-values estimated by three methods for the earthquakes induced in the 

near-well region (within 5 km of the wellhead) during continuous injection operations 

(1996 – 2012)....................................................................................................................50

Figure 3-10. B-values estimated by three methods for the earthquakes induced in the 

northwest (NW) cluster during continuous injection operations (1996 – 2012) ................51

Figure 3-11. (a) The effect on near-well b-value estimates when changing the maximum 

magnitude (Mmax = 3.5, 4.5, 5.5, 7.0) with constant Mmin (0.8), following Bender’s 

(1983) maximum-likelihood calculations. (b) The effect on near-well b-value estimates 

when changing the minimum magnitude (Mmin = 0.8, 0.9, 1.0, 1.1) with constant Mmax 

(4.5). (c) The median b-value calculated for the near-well region with Mmin = 0.8, 0.9, 

1.0, 1.1 and Mmax = 6. (d) The effect of bin window length for fixed Mmax and Mmin ...55

Figure 3-12. Cumulative frequency-magnitude plots for two year windows of near-well 

events with maximum-likelihood solution following Bender (1983) (Mmin=0.8, Mmax = 

4.3) ....................................................................................................................................58

Figure 3-13. Cumulative frequency-magnitude plots for two year windows of north-west 

cluster events with maximum-likelihood solution ollowing Bender (1983) (Mmin=0.8, 

Mmax = 4.3) ......................................................................................................................61
Technical Memorandum TM-68330-2011-39

vii

Figure 3-14. Cumulative seismogenic Index calculated for the near-well region, for 

a range of minimum magnitudes ranging between 0.5 and 1.5...............................64

Figure 3-15. Near-well instantaneous (upper) and cumulative (middle) seismogenic index 

for a range of values between 0.7 and 1.1..............................................................66

Figure 3-16. Northwest cluster instantaneous (upper) and cumulative (middle) 

seismogenic index for a range of values between 0.8 and 1.1...............................67

Figure 3-17. Estimated maximum magnitude for the near-well area, computed using the 

approach of Shapiro and others (2011) for a circular crack and various values of static 

stress drop ........................................................................................................................70

Figure 3-18. Estimated maximum magnitude for the NW-cluster area, computed using the 

approach of Shapiro and others (2011) for a circular crack and various values of static 

stress drop ........................................................................................................................71

Figure 3-19. Comparison of various measures of fault rupture radius versus moment 

magnitude for an assumed circular rupture ......................................................................72

Figure 3-20. Maximum induced earthquake magnitude as a function of (log) cumulative 

injected fluid volume (lower plot).......................................................................................74

Figure 3-21. Residuals of the least-squares fit to the maximum magnitude versus the log 

of cumulative injected volume, plotted against the quantiles of standard normal.............75

Figure 3-22. Maximum magnitude of induced earthquakes from 39 different sites, from 

Nicol and others (2011).....................................................................................................76

Figure 3-23 Hydraulic overpressures inside a simulated fault zone as a function of radial 

distance from the injection well at the time of shut-in (dashed line) and 6 days after shutin, when the largest-magnitude simulated event occurred (solid line)..............................78

Figure 3-24. Magnitude recurrence data for 2009 through 2012 for the near-well area, 

with calculated Gutenberg-Richter relation assuming Mmin = 1.0 and Mmax = 4.0 through 

6.0, in 0.5 magnitude increments......................................................................................80

Figure 3-25. Magnitude recurrence data for 2009 through 2012 for the NW cluster area, 

with calculated Gutenberg-Richter relation assuming Mmin = 1.0 and Mmax = 4.0 through 

6.0, in 0.5 magnitude increments......................................................................................81

Figure 4-1. Daily average injection flow rate (top), daily average surface injection 

pressure (middle), and daily average downhole pressure at 14,100 ft (4.3 km) depth 

(bottom) during PVU injection operations. ........................................................................86

Figure 4-2. Effective permeability (a), dimensionless wellbore storage constant (b), and 

skin due to damage (c) versus time for pressure build-up cylces in Phases II-IV. Dashed 

lines indicate boundaries of the injection phases. ............................................................97

Figure 4-3. Scatter-plot matrix showing pairwise correlations between permeability (k), 

skin due to damage (sd), dimensionless wellbore storage constant (CD), pressure 

increase, and average flow rate for all pressure build-up cycles in Phases II-IV .............98

Figure 4-4. Effective permeability (a) and dimensionless wellbore storage constant (b) 

versus time for all pressure build-up cycles during Phases II-IV, fixing the value of sd to -

4.35. Dashed lines indicate boundaries of the injection phases.......................................99

Figure 4-5. Scatter-plot matrix showing pairwise correlations between permeability (k), 

dimensionless wellbore storage constant (CD), pressure increase, and average flow rate 

for all pressure build-up cycles in Phases II-IV, holding sd fixed at -4.35. ......................100

Σ(t)

Mmin

Mmin

Mmin
viii

Figure 4-6. Recorded (blue dots) and modeled (red line) downhole pressures and percent 

error (dashed line) for an example Phase I cycle, beginning July 10, 1997, allowing all 3 

parameters to vary (top) and fixing sd at -4.35 (bottom). ................................................101

Figure 4-7. Recorded (blue dots) and modeled (red line) values for change in pressure 

(Δp) divided by flow rate (q), and the recorded (red triangles) and modeled (blue lines) 

values for the time derivative of Δp/q for an example Phase I cycle, beginning July 10, 

1997, allowing all 3 parameters to vary (top) and fixing sd at -4.35 (bottom). ................102

Figure 4-8. Recorded (blue dots) and modeled (red line) downhole pressures and percent 

error (dashed line) for an example Phase II cycle, beginning July 26, 1999, allowing all 3 

parameters to vary (top) and fixing sd at -4.35 (bottom). ................................................104

Figure 4-9. Recorded (blue dots) and modeled (red line) values for change in pressure 

(Δp) divided by flow rate (q), and the recorded (red triangles) and modeled (blue lines) 

values for the time derivative of Δp/q for an example Phase II cycle, beginning July 26, 

1999, allowing all 3 parameters to vary (top) and fixing sd at -4.35 (bottom). ................105

Figure 4-10. Recorded (blue dots) and modeled (red line) downhole pressures and 

percent error (dashed line) for an example Phase III cycle, beginning January 8, 2001, 

allowing all 3 parameters to vary (top) and fixing sd at -4.35 (bottom)............................106

Figure 4-11. Recorded (blue dots) and modeled (red line) values for change in pressure 

(Δp) divided by flow rate (q), and the recorded (red triangles) and modeled (blue lines) 

values for the time derivative of Δp/q for an example Phase III cycle, beginning January 

8, 2001, allowing all 3 parameters to vary (top) and fixing sd at -4.35 (bottom). ............107

Figure 4-12. Recorded (blue dots) and modeled (red line) downhole pressures and 

percent error (dashed line) for an example Phase IV cycle, beginning January 14, 2007, 

allowing all 3 parameters to vary (top) and fixing sd at -4.35 (bottom)............................108

Figure 4-13. Recorded (blue dots) and modeled (red line) values for change in pressure 

(Δp) divided by flow rate (q), and the recorded (red triangles) and modeled (blue lines) 

values for the time derivative of Δp/q for an example Phase IV cycle, beginning January 

14, 2007, allowing all 3 parameters to vary (top) and fixing sd at -4.35 (bottom). ..........109

Figure 4-14. Recorded (open squares) and modeled (red lines) downhole pressures for 

the time period from April 2009 to March 2013...............................................................110

Figure 4-15. Maximum and average pressures for a variety of 1-year scenarios with 

varying shut-in schedules................................................................................................111

Figure 4-16. Pressure increases predicted by the simple Horner (1950) model for the 

cases of a 6-month injection cycle with either a 20-day shut and an injection rate of 230 

gpm, or no shut in and an injection rate of 205 gpm.......................................................113

Figure 4-17. Maximum and average pressures for a variety of 1-year scenarios with 

varying flow rates, with a shut-in occurring for the remainder of the year once a fixed 

volume of 108 gallons is reached. ...................................................................................114

Figure 4-18. Maximum pressure vs. injection volume for a variety of injection scenarios 

with varying shut-in schedules and flow rates. In each data series, the data points 

indicate flow rates from left to right of 156 gpm (1.750” plungers), 180 gpm (1.875” 

plungers), 204 gpm (2.000” plungers), and 230 gpm (2.125” plungers).........................115

Figure 4-19. Maximum pressures obtained over one year of injection for 3 different 

permeabilities and 3 flow rates, with 2 evenly spaced 20-day shut-ins per year............116

Figure A-1. Geographical extent of PVSN earthquake catalog. ..................................... A-1
1

Executive Summary

A widely-felt earthquake occurred in the northern Paradox Valley area on January 

24, 2013, at 4:46:39 UTC time (January 23, 2013, at 9:46:39 pm local time). The 

earthquake had an estimated local magnitude of ML 4.4. Using data from a local 

20-station seismic network, we determined the epicenter was 8.2 km (5.1 mi) 

northwest of the Paradox Valley Unit (PVU) injection well and 5.6 km (3.5 mi) 

southwest of the town of Paradox, Colorado (Figure 1-1). The focal depth was

found to be 4.4 km below local ground surface, which is equivalent to a depth of 

4.1 km below the PVU injection wellhead.

Our analysis indicates that this earthquake was induced by long-term fluid 

injection into the PVU injection well. This event was slightly larger than the ML

4.3 earthquake of May 27, 2000, making it the largest PVU-induced event to date. 

However, it was located nearly four times farther away from the injection well 

than any previous PVU-induced earthquake of comparable magnitude. The 

occurrence of a strongly-felt earthquake several kilometers from the injection well 

is cause for greater concern than for an earthquake of similar size, but located 

nearer to the well. That is because the triggering mechanisms for events occurring 

far from the well are poorly understood, the time required to affect seismicity 

changes in this area by varying injection operations is potentially long (e.g., 

several years), and the effectiveness of changes in injection operations to reduce 

the seismic hazard from induced earthquakes at these distances may not be 

apparent for several years.

Many local residents experienced strong ground shaking partly because the

January 24th, 2013 earthquake occurred much closer to northern Paradox Valley

than previous M≥3.5 induced earthquakes. At the town of Paradox, a Reclamation 

strong motion instrument recorded a peak horizontal acceleration of 0.29 g. Even 

accounting for the close distance to the January earthquake, ground motions were 

unusually large. For frequencies above 2.5 Hz, the horizontal accelerations 

recorded in the town of Paradox were nearly two standard deviations higher than 

the mean value estimated by standard ground motion prediction equations for an 

earthquake of this magnitude and distance. 

This technical memorandum (TM) presents an analysis of the January 24th, 2013 

earthquake, its relation to previously-induced seismicity, and an evaluation of the 

ground shaking that it produced (Section 2). Factors that influence the seismic 

risk associated with deep-well injection at PVU are then reviewed (Section 3). 

Results from analysis of PVU injection pressure-flow data, and suggestions for 

operational changes in injection operations are presented (Section 4). Discussion 

and conclusions are presented in Sections 5 and 6. 
2

Results presented in this TM are primarily from analyses in progress before the 

January 24, 3013 earthquake, or that were done in the weeks immediately 

following it. The purpose of these analyses was to fulfill requirements of the 

Emergency Action Plan for operation of the PVU injection well, and to develop 

alternatives for operations to minimize the potential for additional felt 

earthquakes. Hence, this TM mostly contains analyses that could be completed 

quickly after the earthquake. Following completion of the initial draft of this 

report in April, 2013, more detailed analyses were conducted and submitted for 

publication in peer-reviewed scientific journals (Block and others, 2014; King and 

others, 2014; Yeck and others, 2014; Block and others, 2015), and results have 

been incorporated into the current draft of this TM where appropriate. Detailed

analysis of a few remaining topics are in progress, including pressure-flow 

modeling and seismic hazard analysis, and are expected to be documented in 

peer-reviewed journal articles.
3

1 Introduction

An ML 4.4 earthquake occurred in the Paradox Valley area on January 24, 2013, at 

4:46:39 UTC time (January 23, 2013, at 9:46:39 pm local time). The earthquake 

occurred near the northwest end of the valley (Figure 1-1), and was felt strongly 

by local residents.

Figure 1-1. Location of the Jan. 24, 2013 ML 4.4 earthquake relative to the Paradox Valley 

Unit (PVU) injection well (yellow star), PVU extraction well field (PVEF), and the town of 

Paradox, Colorado (PVCC).
4

Available data indicate that this earthquake was induced by fluid injection at the 

Bureau of Reclamation’s (Reclamation) Paradox Valley Unit (PVU) brine 

disposal well. The earthquake epicenter is within the previously-identified 

northwest (NW) cluster of induced seismicity, a persistent zone of induced 

earthquakes 6 to 8 km northwest of the injection well (Figure 1-1). The NW 

cluster has been seismically active since mid-1997, approximately 1 year after the 

start of continuous fluid injection (July, 1996). The January 24th earthquake 

occurred 4.4 km below the local ground surface (4.1 km below the elevation of 

the injection wellhead), consistent with depths of previously-induced events and 

with the depths of injection target formations (Figure 1-2). Natural earthquakes in 

the region tend to occur substantially deeper.

Figure 1-2. PVU-induced earthquakes having well-constrained relative hypocenters, plotted 

as a function of distance from the PVU injection well and depth (relative to the ground 

surface elevation at the wellhead). Each circle represents a single earthquake, with the width 

of the circle scaled by the event magnitude. The January 2013 earthquake is shown in red, 

and the perforated interval in the injection well (4.3 to 4.8 km depth) is indicated by the blue 

circles on the left side of the graph.

Since 1996, Reclamation has been disposing of brine in a deep injection well 

located just southwest of Paradox Valley, Colorado, as part of the PVU salinity 

control project. This period of continuous injection was preceded by a series of 

injection tests conducted between 1991 and 1995. The PVU injection well 

disposes of brine that would otherwise enter the Dolores River, a tributary of the 

Colorado River. This project is authorized by the Colorado River Basin Salinity 

Control Act, passed by the U.S. Congress in 1974, and is an important 

component of the Colorado River Basin Salinity Control Program (CRBSCP), an 

effort to improve water quality in the Colorado River system. Further information 

on the Paradox Valley Unit can be found in Block and others (2012) or at the 

Reclamation project website (see http://www.usbr.gov/projects, under CRBSCP, 

Paradox Valley Unit - Title II).
5

Six years of pre-injection monitoring recorded no seismicity within 18 km of the 

injection well. The data recorded during this period include just 1 local

earthquake, which occurred almost 19 km from the PVU injection well. Within 

days of the start of the initial injection test in 1991, seismic events were detected 

in the immediate vicinity of the injection well. During subsequent injection tests 

and continuous injection operations, induced earthquakes occurred at increasing 

distances from the injection well. The close correlation between injection 

operations and induced seismicity at PVU has been well established, as is 

illustrated by Figure 1-3.

Figure 1-3. Correlation between injection flow rates and shallow seismicity (< 8.5 km depth) 

recorded by the Paradox Valley Seismic Network (PVSN). The network has been 

continuously operational since 1985, 6 years prior to injection, and no seismicity was 

observed within 18 km of the well prior to injection. The large green dot on the right side of 

the figure represents the January 24, 2013 earthquake. It is located within a group of 

induced earthquakes identified as the northwest (NW) cluster.

Microseismic monitoring of the Paradox Valley area began in 1983, with 

installation of the first stations of the Paradox Valley Seismic Network (PVSN). 

Network installation was timed to provide a pre-injection baseline of naturallyoccurring background seismicity, and stations were located to optimize 

monitoring of any earthquakes that might be induced in the vicinity of the well. 

The network has been continuously operated since 1985, and has been expanded 

and updated over the years to respond to changing seismicity patterns. The current 
6

configuration consists of a surface array of 20 digitally-recorded broad-band, 3-

component sensors. A map showing the locations of the stations is shown in

Figure 1-4. The January 24th earthquake was recorded by the broad-band network, 

as well as by 3 digital strong motion instruments operated by Reclamation and 

located at the Paradox Community Center (PVCC), the PVU injection well-head 

and pumping plant (PVPP), and the PVU extraction field and surface treatment 

facility (PVEF).

Figure 1-4. Location of Paradox Valley seismograph stations.
7

2 Description of the Earthquake

2.1 Earthquake Magnitude

We calculated a moment magnitude of MW 4.0 for the January 24th earthquake, 

using only PVSN broad-band data. Seismic moment, and hence moment 

magnitude, were estimated from the long-period displacement spectra at all PVSN 

stations having good-quality data (Figure 2-1). An 2 ω spectral model is assumed 

{Brune, 1970 #2}, including correction for high-frequency attenuation 

{Anderson, 1984 #110}, and model parameters are estimated from the spectra 

using standard methods {Andrews, 1986 #111;Boore, 1986 #112;De Natale, 1987 

#113;Snoke, 1987 #114;Hough, 1988 #115;Anderson, 1991 #116;Shearer, 2006 

#117;Kilb, 2012 #118}. Although the seismic waveforms from the PVSN 

stations closest to the epicenter were clipped, and therefore not usable in this 

analysis, useful data were obtained from 13 of the more distant stations. 

Assuming a simple Brune (1970, 1971) circular rupture model, we also 

determined a static stress drop of 20 bars and rupture radius of 0.6 km. The 

moment magnitude we computed is consistent with magnitudes computed by the 

United States Geological Survey’s (USGS) National Earthquake Information 

Center using data from regional seismic stations and a variety of magnitude 

scales, which include: moment magnitude, MW, 3.9, local magnitude, ML, 4.4, and 

body wave magnitude, mb 3.8 (Harley Benz, personal communication). These 

magnitude scales are described at the end of this section. 

Because earthquakes previously induced by PVU fluid injection with magnitude 

of 3.5 or greater have been reported using the local magnitude scale, ML, for 

consistency we use the ML of January 24th earthquake when comparing the size of 

this event to the historical PVSN seismicity data. The previous largest PVUinduced earthquake was the ML 4.3 earthquake of May 27, 2000. The January 24, 

2013 earthquake, with ML 4.4, is slightly larger than the May 2000 event, making 

it the largest PVU-induced earthquake to date. A comparison of the local and 

moment magnitudes of the 3 largest PVU-induced earthquakes recorded to date is 

provided in Table 2-1.
8

Table 2-1. Comparison of local and moment magnitudes for the 3 largest PVU-induced 

earthquakes. The source for each magnitude estimate is given in parentheses: UUSS = 

University of Utah Seismograph Stations; SLU = Saint Louis University; USGS = U. S. 

Geological Survey; USBR = Bureau of Reclamation.

Earthquake 

Date (UTC)

Local Magnitude 

(ML)

Moment Magnitude 

(MW)

5/27/2000 4.3 (UUSS) 3.8 (SLU)

11/7/2004 4.0 (USGS) 3.6 (SLU)

1/24/2013 4.4 (USGS) 4.0 (USBR)

3.9 (USGS) 

Figure 2-1. Example of automatic spectral fitting procedure used to estimate seismic moment

(see text), showing selected displacement amplitude spectra for P, S, and S-wave coda at 

station PV01 (42 km distance) and PV07 (33 km distance). Dotted lines indicate pre-event 

noise, and dashed lines indicate the modeled fit to the spectra. Arrows on each spectrum 

indicate the fitted corner frequencies. Carat symbols indicate the limits of the frequency 

bands used for fitting the spectra. The spectra are also corrected for a site attenuation 

parameter kappa {Anderson, 1984 #110}, which ranges in value from 0.01 to 0.05. The 

procedure automatically rejects clipped or otherwise bad data, as is demonstrated by the 

second trace from the top.
9

2.1.1 Magnitude Scales 

Multiple magnitude scales are used to quantify the size of earthquakes. Each 

measures different characteristics of the seismic waveforms, which can result in 

different numerical values between the scales and inconsistencies between 

different regions (Chung and Bernreuter, 1981). The duration magnitude scale

(MD), which uses the logarithm of the duration of the seismic signals, is

commonly used for computing magnitudes of local earthquakes recorded by 

microearthquake networks because it is easy to apply and provides consistent 

results (e.g., Lee and Stewart, 1981). Duration magnitudes work well for the 

majority of the induced seismic earthquakes recorded by the Paradox Valley 

Seismic Network (PVSN) because they are small to moderate in size. Although 

the duration magnitude scale generally provides consistent results between 

earthquakes within a given network, it has certain drawbacks: (1) it may not be 

calibrated well between networks, especially for smaller-magnitude events; and, 

(2) the duration magnitude scale may saturate for larger earthquakes (M>3) 

recorded by PVSN due to truncation of long-duration records. Hence, other scales 

must be used to compare the largest PVU-induced events. The local magnitude 

scale (Richter, 1935), based on computing the logarithm of the peak amplitude of 

the seismic signal (corrected for a standard distance and instrumentation type), 

has been reported for some larger PVU events by the USGS and the University of 

Utah Seismograph Stations (UUSS). Other magnitude scales occasionally 

reported by USGS include the surface-wave (MS) and body-wave (mb) 

magnitudes. As seismic instrumentation has improved in recent years, use of the 

moment magnitude scale (M or MW) (Hanks and Kanamori, 1979) has become 

more common.

2.2 Location

We computed a hypocenter for the January 2013 earthquake using P- and S-wave 

arrival times from PVSN and a local 3-dimensional (3-D) velocity model we had 

developed previously. The earthquake epicenter is at latitude 38.32087 deg. N, 

longitude 108.98408 deg. W, which is 8.2 km (5.1 mi) northwest of the PVU

injection well and 5.6 km (3.5 mi) southwest of the town of Paradox, Colorado 

(Figure 1-1). The focal depth is 4.4 km below local ground surface, which is 

equivalent to a depth of 4.1 km below the PVU injection wellhead.

The January 24, 2013 event is the only induced earthquake with ML 3.5 or greater 

(ML 3.5+) to occur at a radial distance greater than about 2 km from the injection 

well (Figure 2-2). The previous four PVU-induced ML 3.5+ events occurred in a 

narrow band between 1.6 and 2.2 km from the injection well (Figure 2-2). At a 

distance of 8.2 km from the injection well, the January earthquake epicenter is

nearly four times farther from the well than that of any other PVU-induced 

earthquake of comparable magnitude. In addition, the January 24th earthquake is 

only the second earthquake with duration magnitude of MD 3.0 or greater to occur 

more than ~2 km from the injection well. The only other MD 3.0+ earthquake at a 
10

relatively large radial distance is an MD 3.3 event that occurred within the NW 

cluster in June, 2002, about 6.6 km from the well (Figure 2-2).

Figure 2-2. Location of PVU-induced seismicity in the near-well region (within 5 km of the 

injection well) and northwest cluster (6 to 8 km from the well).
11

2.3 Focal Mechanism

First-motion data from PVSN stations indicate that the earthquake was produced 

by strike-slip faulting on either of two conjugate steeply-dipping fault planes: (1) 

a fault striking N13°W, and dipping 79° ENE, or (2) a fault striking N78°E, and 

dipping 84° SSE. The directions of compression (P-axis) and extension (T-axis) 

are N58°W and N32°E, respectively. The strike-slip mechanism is consistent with 

the previous finding of Ake and others (2005), who found that strike-slip 

mechanisms accounted for 89% of 1345 well-determined mechanisms and

determined a mean P-axis of N64°W-N67°W.

First-motion data alone do not determine which conjugate plane is the true fault 

plane. However, precise relative locations of earthquakes in the vicinity of the 

January 24th event show a 1.5-km-long lineation oriented approximately N78o

E, 

consistent with the ENE-striking conjugate plane. Several aftershocks of the

January 24th earthquake also occurred along this lineation (Figure 2-3). Hence, the 

earthquake locations suggest that the second conjugate fault orientation, N 78oE, 

is the rupture plane of the January 24th earthquake.

The January 24th earthquake did not occur on any of the deep faults that were 

mapped during previous PVU geophysical investigations. Most of the faults 

mapped during these early studies trend northwest (~N55°W), approximately 

parallel to Paradox Valley (Block and others, 2012). Because this orientation is 

close to the estimated direction of maximum horizontal stress (Ake and others, 

2005), the major basement faults are not optimally oriented for shear slip. The 

majority of the induced earthquakes observed at PVU occur on unmapped, 

previously aseismic faults. Because these faults may not have much net vertical 

offset, they may be difficult to resolve on deep seismic reflection data.
12

Figure 2-3. First motions and computed focal mechanism of the January 2013 event, and 

map showing the distribution of previously induced earthquakes in the vicinity of the 

January 2013 event (a-quality epicenters: black dots; b-quality epicenters: gray dots) and 

the locations of the January 2013 earthquake and aftershocks observed through July 2013 

(a-quality epicenters: magenta circles; b-quality epicenters: green circles). Except for the 

main event, the a-quality epicenters were determined using a relative location routine and 

precise time differences from waveform cross-correlations. The main event was tied into the 

relative location using time differences from high-quality manual arrival time picks. The bquality epicenters were computed using manual arrival time picks and a local 3-D velocity 

model.

2.4 Nearby Seismicity 

In the 6 months prior to the January 2013 event, PVSN recorded about 30 induced 

seismic events within 1 km of the mainshock epicenter (Figure 2-4). These 

foreshocks had duration magnitudes less than MD 2.0, with most having MD ≤ 1.0. 

We did not observe any increase in maximum event magnitude or any substantial 

change in the rate of MD ≥ 0.5 events in the weeks or months prior to the January 

2013 earthquake. In contrast, a large increase in the rate of detected events with 

magnitude MD 0.0 to 0.5 occurred during the third quarter of 2012, with 5 times 

more MD 0.0-0.5 events recorded during this 3-month period than in any previous 
13

quarter (Figure 2-5). However, the baseline period to determine seismicity rates of 

MD 0.0 to 0.5 events was only about a year and a half at the time of the January 

2013 event, because events this small were not reliably detected in this area until 

after the installation of two additional nearby seismic stations in July, 2011. Rates 

of MD 0.0 to 0.5 events returned to normal levels during the fourth quarter of 

2012. Hence, we have limited data to indicate whether the observed increased rate 

of MD 0.0 to 0.5 events a few months prior to the January 2013 earthquake is 

anomalous.

We found that earthquakes in the vicinity of the January 2013 event can be 

grouped into two distinct faulting types, based on an analysis of 35 events with 

well-constrained relative hypocenters: (1) strike-slip faulting consistent with that 

of the January 2013 earthquake (average mechanism from 20 events is N76°E, dip 

76°SSE; standard deviations of strike and dip are about 5°), and (2) an oblique 

normal focal mechanism (average mechanism from 15 events is (a) strike N42°W, 

dip 71°WSW, or (b) strike N77°E, dip 35° NNW, with standard deviations of 

strike and dip less than 7°). Epicenters of these earthquakes are plotted by focal 

mechanism type in Figure 2-6a. All analyzed events with strike-slip mechanisms 

are distributed along the inferred east-northeast trending fault segment interpreted 

to have ruptured during the January 2013 earthquake. Events with the obliquenormal focal mechanism form a tightly-spaced cluster about 250 to 450 meters 

(m) north of the strike-slip fault plane. Most of the foreshocks detected in the 6

months prior to the January 2013 earthquake occurred near this cluster, and 

relatively few foreshocks occurred near the interpreted mainshock fault plane 

(Figure 2-6a, squares and octagons). However, because of their small magnitude, 

most of these foreshocks do not have a sufficient number of clear first-motions to 

compute robust focal mechanisms.

The area surrounding the inferred rupture of the January 2013 earthquake has 

been seismically active since at least June 1997, about 1 year after the start of 

continuous PVU fluid injection, when PVSN recorded a cluster of earthquakes 

near the eastern end of the fault segment (Figure 2-6b). In October 1998, an 

earthquake was first detected near the center of the fault segment, approximately 

500 to 700 m west of the initial cluster of events. By August 1999, seismicity was 

occurring near the western end of the fault segment, to within 100 m of the 

January 2013 event epicenter (Figure 2-6b). The seismicity pattern delineating the 

western end of the fault segment has not grown in size since 2000, suggesting that 

the segment is truncated in this direction, either by a major northwest-trending 

normal fault or by pinch-out of the primary injection target formation due to 

erosion (Block and others, 2012); (King and others). The seismicity pattern 

delineating the eastern limit of the fault segment has continued to expand slowly, 

and several aftershocks of the January 2013 event have extended seismicity on the 

fault segment to the east by about 100 m. 

The entire fault segment, as inferred from pre-2013 seismicity, appears to have 

ruptured during the January 2013 earthquake. Aftershocks of the January 2013 
14

earthquake occur near both ends of the pre-2013 seismicity zone (Figure 2-4 and 

Figure 2-6), and the inferred fault segment length determined from the seismicity 

distribution is consistent with the rupture dimension expected for the moment 

magnitude of the January 2013 earthquake, as will be discussed in Section 3.3.1. 

Although small earthquakes have occurred since 1999-2000 along most of the 

fault segment that ruptured in January 2013, there was a 13-year delay before the 

occurrence of an earthquake which finally ruptured the entire segment. Only two 

earlier earthquakes of MD ≥ 2.5 occurred on portions of this fault segment: an MD

2.6 event occurred near the center of the fault segment in October 2003, and an 

MD 2.7 event occurred near the east end of the fault segment in July 2008 (Figure 

2-6b). Assuming a similar stress drop as for the January 2013 earthquake (~20 

bars), and a circular crack model, these earthquakes each ruptured less than a 300-

m portion of the fault segment, or less than about 20% of the total fault segment 

length. One possible explanation for the 13-year delay prior to rupture of the 

entire fault segment is that a pore pressure threshold must be exceeded over most 

of the fault surface in order for rupture to propagate. We have previously 

observed that a pore pressure threshold may need to be reached to induce largermagnitude earthquakes in the near-well region, based on correlation of near-well 

larger-magnitude earthquakes and higher injection pressures (Block and Wood, 

2009; 2010). Correlation between larger-magnitude induced earthquakes and an 

increase in pore pressure to a threshold level over most of the fault surface has 

also been suggested by numerical studies (Garagash and Germanovich, 2012). 
15

Figure 2-4. Map showing the distribution of foreshocks and aftershocks of the Jan. 24, 2013 

event. The a-quality hypocenters were determined by a high-precision relative event 

location technique; the b-quality events did not have sufficient data to be included in the 

relative location and therefore have less precise locations. The July, 2012 through July, 2013 

events are color-coded by month of occurrence and include both a- and b-quality 

hypocenters.
16

Figure 2-5. Occurrence of earthquakes within 1.0 km of the January 24, 2013 event 

epicenter, plotted as a function of date and earthquake magnitude. The area of each circle 

in scaled by the number of events recorded in a given quarter-year and magnitude range.
17

Figure 2-6. Epicenters of earthquakes with well-constrained relative locations in the near 

vicinity of the January 2013 ML 4.4 earthquake. (a) Earthquakes having similar strike-slip 

focal mechanisms are indicated by the filled red circles, and those having similar oblique 

normal focal mechanisms are indicated by the filled green circles. A typical focal sphere plot 

for each type of mechanism is shown. Gray circles indicate events that lack well-defined 

focal mechanisms. Earthquakes occurring during the 6-month period before the January 

2013 event are identified by the open squares and octagons, as indicated in the legend. (b) 

Earthquakes color-coded by year of occurrence and with symbol size scaled by event 

magnitude.
18

2.5 Strong Ground Motions

Strong-motion records were obtained from the Paradox Community Center 

(station code PVCC), the PVU injection wellhead and pumping plant (station 

code PVPP), and from the surface treatment facility and extraction well field 

(station code PVEF). The acceleration time histories are shown in Figure 2-7

(PVCC), Figure 2-8 (PVPP), and Figure 2-9 (PVEF), and a summary of the 

observed peak ground motions is provided in Table 2-2. Acceleration and velocity 

response spectra were also computed for each component of motion at each 

station, and these are shown in Figure 2-10 through Figure 2-15.

Figure 2-7. Acceleration ground motion time histories recorded by strong-motion 

instruments located at the community center in the town of Paradox, Colorado (PVCC). 

Traces are offset for clarity. Recorded peak horizontal acceleration (PHA) at PVCC was 

0.28 g.

PVCC

-600

-400

-200

0

200

400

600

13.0 14.0 15.0 16.0 17.0 18.0 19.0 20.0 21.0 22.0 23.0

Time (sec.)

Acc. (cm/s/s)

E

N

Z
19

Figure 2-8. Acceleration ground motion time histories recorded by strong-motion 

instruments located near the PVU injection well and pumping plant (PVPP). Traces are 

offset for clarity. Recorded peak horizontal acceleration (PHA) at PVPP was 0.15 g.

Figure 2-9: Acceleration ground motions recorded by strong-motion instruments located at

the PVU surface treatment facility and extraction-well field (PVEF). Traces are offset for 

clarity. Recorded peak horizontal acceleration (PHA) at PVEF was 0.08 g.

PVPP

-600

-400

-200

0

200

400

600

13.0 14.0 15.0 16.0 17.0 18.0 19.0 20.0 21.0 22.0 23.0

Time (sec.)

Acc. (cm/s/s)

E

N

Z

PVEF

-600

-400

-200

0

200

400

600

14.0 15.0 16.0 17.0 18.0 19.0 20.0 21.0 22.0 23.0 24.0

Time (sec.)

Acc. (cm/s/s)

E

N

Z
20

Table 2-2. Observed ground motions and distances for the Paradox Community Center 

(PVCC), the injection well (PVPP), and the surface treatment facility (PVEF).

Site Distance from 

Epicenter (km)

PHA 

(g)

1.0-sec. 

SA (g)

PHV 

(cm/s)

PVCC 5.6 0.293 0.011 6.1

PVPP 8.2 0.153 0.006 3.5

PVEF 11.3 0.077 0.006 2.0

Strong-motion instruments are designed to make on-scale recordings of strong 

ground shaking, and were installed in the Paradox Valley area to evaluate the 

likely engineering effects of shaking on structures. In contrast, the 20-station 

PVSN is designed to detect the weak motions from small earthquakes, and the 

high-gain instrumentation tends to saturate for large earthquakes. The strong 

motion instruments are 24-bit Kinemetrics model K2 data loggers, with either 

internal Episensor or external FBA-23 accelerometers, each with a nominal fullscale range of 2 g (PVCC and PVPP) or 1g (PVEF). (g is the value of 

gravitational acceleration on the surface of the Earth, equal to ~980 cm/s2 at sea 

level.) The instruments record ground acceleration in 3 mutually perpendicular 

directions (East-West, North-South, and vertical). The first strong-motion 

instrument was installed near the injection wellhead (station code PVPP) in 

January, 1998. In November, 1999, a second strong-motion instrument was 

installed near the surface treatment facility (station code PVEF). A third strongmotion instrument was installed in the town of Paradox, Colorado (station code 

PVCC) following an ML 4.0 induced earthquake on November 7, 2004.

Residents of the nearby community of Paradox, and in surrounding areas,

reported feeling strong shaking from the January, 2013 earthquake. These reports 

are consistent with recordings obtained from the strong motion instruments. These 

data show peak horizontal accelerations (PHA) of 287 cm/s2 (0.29 g) in the town 

of Paradox, 150 cm/s2 (0.15 g) near the wellhead, and 75 cm/s2 (0.08 g) at the 

surface treatment facility.

To get a sense of whether ground motions recorded at Paradox Valley during the 

January, 2013 earthquake differed substantially from what would be expected for 

tectonic earthquakes at similar magnitudes and distances, we have compared the 

recorded horizontal ground motions to ground motions estimated using the Next 

Generation Attenuation (NGA) empirical ground-motion prediction equations

(Power and others, 2008). We use four of the NGA equations (Abrahamson and 

Silva, 2008; Boore and Atkinson, 2008; Campbell and Bozorgnia, 2008; Chiou 

and Youngs, 2008), and compute a mean value by assigning equal weight to each 

NGA equation used. To use the NGA equations, the following site properties 

must be assumed: (1) VS30, defined as the time-averaged shear-wave velocity in 

the upper 30 m of soil; and, (2) Z1.0 and Z2.5, defined as the depths to where shearwave velocities reach 1.0 and 2.5 km/s, respectively. We have assumed site VS30

values corresponding to moderate to soft soil sites, and computed Z1.0 and Z2.5

values using the relations in Abrahamson and Silva (2008) and Campbell and
21

Bozorgnia (1987), following the suggestions of Kaklamanos and others (2011). 

Assumed properties required for the NGA relations are shown in Table 2-3. The 

NGA values are for a randomly oriented horizontal component of motion, but we 

have not corrected the data to obtain the random component.

Table 2-3. Input parameters used for estimating ground motions from the NGA equations. 

See Kaklamanos and others (2011) for a description of the parameters.

Site RJB

(km)

ZTOR

(km)

RRUPT

(km)

Dip 

(deg)

Rake

(deg)

Width

(km)

VS30

(m/s)

Z1.0

(m)

Z2.5

(m)

PVCC 4.9

4.0

6.3

80 5 0.6

300 330 1700

PVPP 7.1 8.1 530 170 1100

PVEF 10.0 10.8 400 220 1300

In general, the observed short-period (< 0.4-second) ground motions from the 

January 24th earthquake are substantially greater than the NGA relations predict, 

given our assumed site properties. The mean NGA acceleration and velocity 

response spectra are shown for each spectral period in Figure 2-10 through Figure 

2-15, and compared to the observed response spectra. For periods greater than 1.0 

second, the observed results are consistent with the NGA relations. The recorded 

high-frequency (> 2.5 Hz) recorded ground motions are substantially higher than 

would be predicted by the NGA equations. This effect is greatest at PVCC, but is 

observed at all stations. Comparison of the recorded peak horizontal accelerations, 

1.0-second spectral accelerations, and peak horizontal velocities are shown in 

Figure 2-16 through Figure 2-18. The recorded PHA at PVCC is nearly two 

standard deviations above the NGA mean PHA (Figure 2-16), and peak horizontal 

velocity (PHV) is more than 1 standard deviation above the NGA mean PHV

(Figure 2-18). At PVEF, the PHA and PHV values are also above the mean NGA 

values, but are within 1 standard deviation of the mean. In contrast, the 1.0-

second spectral acceleration values at all sites are within 1 standard deviation of 

the NGA mean (Figure 2-17).
22

Figure 2-10. Acceleration response spectra for the Paradox Community Center (PVCC), and 

comparison with NGA results.

Figure 2-11. Velocity response spectra for Paradox Community Center (PVCC), and 

comparison with NGA results.

0.01

0.1

1

10

100

1000

0.01 0.10 1.00 10.00

Acceleration (cm/s/s)

Period (sec.)

PVCC - Acceleration Response Spectra

E

N

Z

Mean NGA

0.010

0.100

1.000

10.000

100.000

0.01 0.10 1.00 10.00

Velocity (cm/s)

Period (sec.)

PVCC - Velocity Response Spectra

E

N

Z

Mean NGA
23

Figure 2-12. Acceleration response spectra for the injection well (PVPP), and comparison 

with NGA results.

Figure 2-13. Velocity response spectra for the injection well (PVPP), and comparison with 

NGA results.

0.01

0.1

1

10

100

1000

0.01 0.10 1.00 10.00

Acceleration (cm/s/s)

Period (sec.)

PVPP - Acceleration Response Spectra

E

N

Z

Mean NGA

0.010

0.100

1.000

10.000

0.01 0.10 1.00 10.00

Velocity (cm/s)

Period (sec.)

PVPP - Velocity Response Spectra

E

N

Z

Mean NGA
24

Figure 2-14. Acceleration response spectra for the surface-treatment facility (PVEF), and 

comparison with NGA results.

Figure 2-15. Velocity response spectra for surface-treatment facility (PVEF), and 

comparison with NGA results.

0.01

0.1

1

10

100

1000

0.01 0.10 1.00 10.00

Acceleration (cm/s/s)

Period (sec.)

PVEF - Acceleration Response Spectra

E

N

Z

Mean NGA

0.01

0.1

1

10

0.01 0.10 1.00 10.00

Velocity (cm/s)

Period (sec.)

PVEF - Velocity Response Spectra

E

N

Z

Mean NGA
25

Figure 2-16. Comparison of observed PHA values versus distance, with NGA results.

Figure 2-17. Comparison of observed 1.0-second spectral acceleration values versus distance, 

with NGA results.

0.1

1

10

100

1000

0.1 1 10 100 1000

Acceleration (cm/s/s)

Joyner-Boore Distance (km)

NGA - Predicted Peak Horizontal Acceleration

mean + 2σ

mean + 1σ

Mean NGA

mean - 1σ

mean - 2σ

Recorded PVPP

PVEF

PVCC

0.1

1

10

100

0.1 1 10 100 1000

Acceleration (cm/s/s)

Joyner-Boore Distance (km)

NGA - Predicted 1.0-sec. Spectral Acceleration

mean + 2σ

mean + 1σ

Mean NGA

mean - 1σ

mean - 2σ

Recorded

PVPP PVEF

PVCC
26

Figure 2-18. Comparison of observed peak velocity values versus distance, with NGA results.

The strong-motion records also can be used to evaluate the level of shaking in the 

town of Paradox by comparing them to vibration control standards used by the 

mining and construction industries to minimize shaking damage to residential 

structures. This comparison is done in the following section.

2.6 Ground Motion Effects and Standards

We know of no regulatory standards for control of shaking from induced 

seismicity applicable to the Paradox Valley area. However, well-established 

standards exist for the control of vibrations during blasting and construction 

activities, which are applicable to residential structures. According to an induced 

seismicity protocol recently published by the U. S. Department of Energy (DOE) 

(Majer and others, 2012), “the seismic energy from these activities would be 

similar to that from induced seismic events (in frequency bandwidth and range) 

and thus be applicable to induced seismicity cases”. This protocol specifically 

addresses seismicity induced by fluid injection during geothermal energy 

development. However, we believe that the recommendations in this protocol are 

directly applicable to seismicity induced by fluid injection at the Paradox Valley 

Unit because the physical mechanisms for triggering the induced earthquakes, and 

therefore their characteristics, are essentially the same. According to the DOE 

protocol, of the various standards used for blasting and construction activities, 

those established for the regulation of ground vibrations from surface mine 

blasting are the most applicable to injection-induced seismicity. 

0.01

0.1

1

10

100

0.1 1 10 100 1000

Velocity (cm/s)

Joyner-Boore Distance (km)

NGA - Predicted Peak Horizontal Velocity

mean + 2σ

mean + 1σ

Mean NGA

mean - 1σ

mean - 2σ

Recorded

PVPP

PVEF

PVCC
27

The Office of Surface Mining (OSM), U.S. Department of the Interior, has 

performed extensive studies of ground vibrations and building damage caused by 

surface blasting. These studies began as early as the 1930s. A report published by 

OSM in 1971 summarizes the results from a 10-year program conducted to study 

the effects of vibrations from surface blasting (Nicholls and others, 1971). One of 

the main goals of this program was to establish the relationship between the 

magnitude of ground vibration and the resulting damage to structures. The report 

recommends that peak particle velocity of ground motion, as measured in any of 3 

mutually perpendicular directions, be used as the measurement of vibration from 

blasting for correlation to damage of residential structures. The report states that, 

based on analysis of data from several earlier studies, particle velocity more 

closely relates to structural damage than measurements of displacement or 

acceleration. The report concludes that a peak particle velocity (PPV) limit of 2.0 

in/sec (5.1 cm/sec) will provide protection from damage to residential structures 

in the majority of cases (Nicholls and others, 1971). The report also states that, “If 

the observed particle velocity exceeds 2.0 in/sec in any of the 3 orthogonal 

components, there is a reasonable probability that damage will occur to residential 

structures.” 

Additional studies indicated that the potential for damage to structures varied 

considerably depending on the frequency range of the ground vibrations. This 

variation is due to the natural frequency response of buildings, which can amplify 

the particle velocities. Results published by the Office of Surface Mining in 1980 

indicate that the natural resonance frequencies for typical one- to two-story 

residential structures range from about 4 to 24 Hz, and the maximum 

amplification factors measured for these structures occurred for ground vibrations 

between 5 and 12 Hz (Siskind and others, 1980). As a result of the new studies, 

OSM reduced the maximum peak particle velocity limit for ground vibrations 

under 40 Hz. For ground vibrations in this frequency range, OSM concluded that 

a PPV limit of 0.5 in/sec (1.3 cm/sec) will provide protection from blast damage 

for older homes with plaster walls. Modern drywall (gypsum board) interiorwalled homes were judged to be safe at somewhat higher vibration levels, up to 

0.75 in/sec (1.9 cm/sec) (Siskind and others, 1980). The report left the 

recommended limit for ground vibrations above 40 Hz at the previous level of 2.0 

in/sec. The frequency content of ground vibrations generated by induced seismic 

events fall within the lower frequency range (< 40 Hz), and therefore the lower 

PPV limits would be applicable. 

The current OSM standards for ground vibrations from surface mine blasting are 

stated in Title 30 of the Code of Federal Regulations (2012) and summarized in

Table 2-4. The regulations state that the maximum ground vibration shall not 

exceed the limits given in Table 2-4 “at the location of any dwelling, public 

building, school, church, or community or institutional building”. The regulations 

also state that the particle velocity should be measured in 3 mutually 

perpendicular directions and that the vibration limits apply to each of the 3 

measurements. Three separate peak particle velocity limits are listed in the 
28

Federal Code, depending on the distance of the structure from the blasting site. 

The limits are stricter for farther distances, because the ground motions become 

more dominated by lower frequencies at greater distances from the source. As 

discussed above, the lower frequencies are more damaging to structures than the 

higher frequencies. The PPV limit listed for distances greater than 5000 ft (0.75 

in/sec; 1.91 cm/sec) is the most applicable to ground vibrations from induced 

earthquakes. The OSM limits were adopted by the State of Colorado for the 

regulation of ground vibrations from blasting associated with coal mining in 1980 

(Alt and others, 1980).

Table 2-4. Summary of current ground vibration limits from blasting set by the Office of 

Surface Mining, Reclamation and Enforcement, U. S. Department of the Interior; taken 

from Title 30, Subchapter K, of the Code of Federal Regulations, dated July 1, 2012.

Distance from Blasting Site 

(ft)

Peak Particle Velocity 

(in/sec) (cm/sec)

0 – 300 1.25 3.18

301 – 5000 1.00 2.54

5001 and higher 0.75 1.91

In order to compare the ground vibrations measured by the PVSN strong motion 

instruments to the vibration limits given by the OSM regulations, the acceleration 

time histories were integrated to produce velocity records. The velocity data were 

then corrected for instrument response. The final velocity records are presented in 

Figure 2-19 (PVCC), Figure 2-20 (PVPP), and Figure 2-21 (PVEF).

The peak particle velocity measured in the town of Paradox for January 24, 2013 

earthquake was 6.2 cm/sec (2.4 in/sec, East component). This ground vibration 

level is more than 3 times the 1.9 cm/sec (0.75 in/sec) limit set in the OSM 

regulations for distances from the source exceeding 5000 ft. A comparison of the 

peak particle velocity records measured in the town of Paradox and these 

standards is shown in Figure 2-19. Based on this comparison, it is possible that 

minor damage may have occurred to structures in and near the town of Paradox 

from ground vibrations produced by the earthquake.

The peak particle velocity measured at the surface treatment facility is 1.9 cm/sec 

(0.78 in/sec, East component), slightly exceeding the OSM safe vibration limit

(Figure 2-21). Based on the OSM standards, damage would not be expected to 

have occurred at this location from the earthquake.
29

Figure 2-19. Comparison of measured particle velocities at station PVCC (Town of Paradox, 

community center) and the blasting standard of 0.75 in/sec peak particle velocity (orange 

dashed lines). Recorded ground motions reached a maximum peak particle velocity of 2.4 

in/sec, more than 3 times the standard.

Figure 2-20. Comparison of measured particle velocities at station PVPP (near the injection 

well-head) and the blasting standard of 0.75 in/sec peak particle velocity (orange dashed 

lines). Recorded ground motions reached a maximum peak particle velocity of 1.4 in/sec, 

exceeding the standard.

PVCC - Measured Particle Velocities

-2.5

-2.0

-1.5

-1.0

-0.5

0.0

0.5

1.0

1.5

2.0

2.5

13.0 14.0 15.0 16.0 17.0 18.0 19.0 20.0 21.0 22.0 23.0

Time (sec.)

Vel. (in/s)

E

N

Z

PVPP - Measured Particle Velocities

-2.5

-2.0

-1.5

-1.0

-0.5

0.0

0.5

1.0

1.5

2.0

2.5

13.0 14.0 15.0 16.0 17.0 18.0 19.0 20.0 21.0 22.0 23.0

Time (sec.)

Vel. (in/s)

E

N

Z
30

Figure 2-21. Comparison of measured particle velocities at station PVEF (Surface 

Treatment Facility and brine extraction) and the blasting standard of 0.75 in/sec peak 

particle velocity (orange dashed lines). Recorded ground motions reached a maximum peak 

particle velocity of 0.78 in/sec, slightly exceeding the standard.

2.7 Triggering Mechanism

During fluid injection, earthquakes may be induced by increases in pore pressure, 

by redistribution of stress (from accommodation of the injected fluid into the rock 

or from the occurrence of previously-induced earthquakes), or by cooling and 

shrinking of the rock matrix. The latter thermodynamic effect is important in 

geothermal areas but is unlikely to be a dominant factor at PVU, except possibly 

in the near vicinity of the injection well. Changes in pore pressure and the 

redistribution of stress are likely the major triggering mechanisms for earthquakes 

induced by fluid injection at PVU.

The triggering mechanism for the January 24, 2013 event in the NW cluster is not 

known. Previous analysis shows a general correlation between the occurrence of 

induced earthquakes in the near-well region (within 4 to 5 km of the injection 

well) and long-term average injection pressures (Block and Wood, 2009), as is 

shown in Figure 2-22. This correlation suggests that pore pressure increase is a 

dominant factor contributing to the rate and magnitude of induced seismicity in 

the near-well region. The same analysis, however, does not show a correlation 

between injection pressure and seismicity for earthquakes occurring at distances 

greater than 5 km, such as those in the NW cluster (Figure 2-22, lower plot). 

PVEF - Measured Particle Velocities

-2.5

-2.0

-1.5

-1.0

-0.5

0.0

0.5

1.0

1.5

2.0

2.5

14.0 15.0 16.0 17.0 18.0 19.0 20.0 21.0 22.0 23.0 24.0

Time (sec.)

Vel. (in/s)

E

N

Z
31

Roeloffs and Denlinger (2009) suggested that seismicity was initially triggered in 

the NW cluster by stress redistribution. This hypothesis stems from the fact that 

the seismicity in the NW cluster began only 1 year after the start of continuous 

injection operations, too soon for significant pore pressure changes to have 

propagated 6 km from the well, based on results from an axisymmetric fluid flow 

model. Mahrer and others (2008) proposed a different possibility - a northwesttrending fault zone allows for relatively rapid propagation of pore pressure from 

the vicinity of the injection well to the NW cluster (a feature the radiallysymmetric fluid flow model of Roeloffs and Denlinger would not have accounted 

for). In this case the earthquakes may be primarily driven by pore pressure 

changes. It is also feasible that the triggering mechanisms for earthquakes 

occurring in the NW cluster have changed over time, with stress redistribution 

being the primary cause of the initial seismicity and pore pressure increase 

becoming an increasingly important factor at some later date as the pressure front 

from injection has expanded over time. 
32

Figure 2-22: Injection downhole pressure data averaged over daily, 6-month, 18-month, and 

30-month time periods (top), and occurrence of induced seismicity as a function of time and 

magnitude within 5 km of the injection well (middle) and at distances of 5 to 10 km from the 

well (bottom). In the seismicity plots, the area of each circle is scaled by the number of 

events in a given quarter-year and magnitude range. The low seismicity rate in the smaller 

magnitude bins from about mid-2005 to mid-2007 in the bottom plot is believed to be due to 

an unusually large number of offline stations.
33

3 Seismic Hazard Factors

Factors controlling seismic hazard include seismicity rates, distribution of 

earthquake magnitudes, and the distances from the primary areas of seismicity to 

residences or other sites of concern for strong ground shaking. For the induced 

seismicity at Paradox Valley, some of these factors appear to change substantially 

over time periods of several years or less. Such variability makes it difficult to 

estimate long-term seismic hazard, since it will not be sufficient either to 

determine parameters at a single point in time, or to characterize hazard using the 

occurrence of a particular earthquake. In order to evaluate seismic hazard arising 

from long-term fluid injection at PVU, we examine how parameters normally 

used to characterize hazard have changed over the life of the project, and assess

how they may change in the future.

In this section we review what is currently known about key factors controlling 

seismic hazard arising from injection operations at PVU, and examine how these

have varied over the history of the project. We also discuss important factors for 

which we have little or no data.

3.1 Geographical Distribution of Seismicity

To first approximation, common measures of earthquake shaking (e.g., peak 

horizontal acceleration or velocity) depend on the magnitude of the earthquake 

and the distance between the causative fault and a site of interest. For a given 

magnitude, the expected ground motions (and associated seismic hazard) 

generally are higher if earthquakes are occurring closer to a site. For induced 

earthquakes, maximum magnitude may be related to the volume of rock 

stimulated by injection, and the stimulated volume may be inferred by the extent 

of the cloud of induced seismicity. An analysis of the spatial evolution of induced 

seismicity in the Paradox area therefore can provide important constraints for

evaluating hazards. Data collected by the Paradox Valley Seismic Network

provides the basis for such an analysis.

3.1.1 Relation of Seismicity to Inhabited Areas in Paradox Valley

The geographical distribution of induced seismicity at PVU has evloved over 

time, with clusters of seismicity expanding from the near-well region, during the 

earliest phases of injection, to locations closer to inhabited areas in Paradox 

Valley. Since the inhabited areas are spread out over a large region, the issue of 

the distances between potential seismic sources to sites of concern is complex. 

Below we provide a brief review of the geographical evolution of PVU-induced 

seismicity over time. We then provide a simplified look at the source-to-site 
34

distance parameter by examining the evolution of closest distances from induced 

earthquakes to the local population center, the town of Paradox in northern 

Paradox Valley.

3.1.1.1 Geographical Evolution of Seismicity

The geographical distribution of induced seismicity over time is presented in the 

series of maps shown in Figure 3-1 and Figure 3-2. For reference, we have 

delineated and named four key seismicity areas: (1) Near-well: By the end of the 

injection tests in 1995, earthquakes were occurring 3 to 4 km from the injection 

well (Figure 3-1a). This area of induced seismicity immediately surrounding the 

injection well is referred to here as the “near-well” region. (2) Northwest (NW) 

cluster: In mid-1997, about a year after the start of continuous injection, 

earthquakes began occurring 6 to 8 km northwest of the injection well (Figure 

3-1b). This cluster of induced seismicity is called the “northwest (NW) cluster”. 

(3) Northern valley: In mid-2000, earthquakes were first detected 12 to 14 km 

from the injection well, along the northern edge of Paradox Valley (Figure 3-1b). 

Several distinct clusters of earthquakes have occurred along the northern edges of 

the valley since 2000 (Figure 3-2a). The earthquakes occurring in all of these 

clusters are referred to as “northern valley events”. (4) Southeast (SE) cluster: 

An earthquake was first detected about 6 km southeast of the injection well in 

2004 (Figure 3-2a), but the seismicity rate in this area markedly increased 

beginning in 2010 (Figure 3-2b). This group of earthquakes is referred to as the 

“southeast (SE) cluster”. In recent years, a few isolated earthquakes have been 

detected in previously aseismic areas, including in the center of Paradox Valley 

(Figure 3-2b).

The geographical expansion of seismicity has varied over time. The rate of 

expansion was highest during the later injection tests (1994-1995) and the early 

phase of continuous injection (1996-2000). However in mid-2000, when the PVU 

injection flow rate was decreased by about 33% in response to a magnitude 4.3 

induced earthquake, the geographical expansion of seismicity slowed greatly, and 

the spatial extent of the near-well and NW-cluster areas largely stabilized. While 

new clusters of seismicity appeared in the northern-valley area in the years 

following the reduction in flow rate, seismicity rates there were very low. This 

relatively stable period lasted for about a decade.

Beginning in 2009, the spatial distribution and rates of seismicity within 10 km of 

the injection well began to change. Five induced earthquakes of magnitude 2.5 or 

greater occurred in 2009. This was the highest annual rate of M 2.5+ earthquakes 

in 5 years, and the second highest rate since the injection flow rate was reduced in 

mid-2000. A distinct group of earthquakes developed in 2010 about 6 km 

southeast of the injection well, the SE cluster. While the first detected SE cluster 

event occurred in 2004, with 2 additional events recorded in 2008 and 2009, 

beginning in 2010 the SE cluster experienced increasing rates of seismicity: 11 

events in 2010, 17 events in 2011, and 25 events in 2012.
35

Figure 3-1. Maps showing the geographical distribution of shallow seismicity recorded in the 

Paradox Valley area from 1991 through 2000: (a) injection tests, 1991-1995 (b) continuous 

injection, 1996-2000. All detected earthquakes locating less than 8.5 km deep (relative to the 

ground surface elevation at the injection wellhead) are included.

a)

b)
36

Figure 3-2. Maps showing the geographical distribution of shallow seismicity recorded in the 

Paradox Valley area from 2001 through 2012: (a) continuous injection, 2001-2008 (b) 

continuous injection, 2009-2012. All detected earthquakes locating less than 8.5 km deep 

(relative to the ground surface elevation at the injection wellhead) are included.

Seismicity rates within the northern-valley area have also changed in recent years. 

During each year from 2000 (when the northern valley seismicity was first 

detected) to 2009, between 2 and 33 earthquakes were recorded. In 2010, the rate 

increased markedly: 557 northern-valley earthquakes were recorded, with the 

a)

b)
37

majority occurring in a single swarm lasting just 16 days. Northern-valley 

seismicity rates remained elevated during 2011, with 113 earthquakes recorded, 

but declined back to pre-2010 rates during 2012, with just 10 events recorded.

Beginning in mid-2010 and continuing through the present, several shallow 

earthquakes have been detected beneath the floor of Paradox Valley. Three 

earthquakes were detected in 2010; 2 in 2011; and 5 in 2012. Although the total 

number of events is small (10 events in 3 years), no earthquakes were detected 

beneath the floor of Paradox Valley in the 25 years of seismic monitoring prior to 

2010. The renewed spatial expansion of seismicity and increased seismicity rates 

in recent years may be related to a trend of increasing injection pressures at the 

well, which in turn is related to the long-term average injection flow rates.

3.1.1.2 Epicentral Distance Evolution of Seismicity

Because most of the local population near the injection well resides in the 

northern part of Paradox Valley, with the town of Paradox being the main 

population center, the distance of PVU-induced earthquakes from Paradox is an 

important seismic hazard factor. To evaluate how the distance of PVU-induced 

seismicity from the town of Paradox has changed over time, we plot the 

seismicity as a function of date and distance from Paradox in Figure 3-3. In 

addition, the maps in Figure 3-1 and Figure 3-2 include circles of equal distance 

from Paradox, for distances of 5, 10 and 15 km.

Figure 3-3. Occurrence of PVU-induced seismicity as a function of date, distance from the 

town of Paradox, and earthquake magnitude. All recorded events locating less than 8.5 km 

deep (relative to the ground surface elevation at the injection wellhead) are included. Each 

circle represents a single earthquake, with its diameter scaled by event magnitude.
38

Epicentral distances between the town of Paradox and the nearest induced 

earthquakes have changed markedly over time, as shown by the plot of distance 

versus time (Figure 3-3), and by epicenter maps (Figure 3-1 and Figure 3-2). The 

injection well is located at a distance of 10 km from Paradox. During the earliest 

injection tests (1991-1994), the seismicity was restricted to within about 1 km of 

the injection well, and therefore the epicentral distance between the induced 

seismicity and the town of Paradox ranged from slightly less than 9 km to 11 km 

(Figure 3-3). During the later injection tests (1994-1995), the induced seismicity 

cloud expanded away from the injection well, with earthquakes occurring 

between 7.5 and 11.5 km from Paradox. During continuous injection operations, 

as the boundary of the seismicity cloud has progressively expanded away from the 

well, earthquakes have occurred at both smaller and larger distances from 

Paradox, with shallow, potentially-induced earthquakes now being detected as 

close as 2 km from Paradox, and as far as 17 km.

The earthquakes occurring closest to the town of Paradox include events in the 

center of Paradox Valley (with epicentral distances as small as 1.9 km), events in 

the NW cluster (at distances of about 3.5 to 7 km), and events in some of the 

northern-valley clusters (at distances of about 4 to 7 km). To date, of the

earthquake clusters occurring near the town of Paradox, only the NW cluster has 

produced M>3 earthquakes. As discussed in the following sections, however, both 

the rates and magnitude distributions of PVU-induced seismicity have varied 

significantly over time. Hence, we cannot exclude the possibility that future M>3

earthquakes will occur in the other earthquake clusters near the town of Paradox, 

including either the central or northern-valley areas.

3.1.2 Volumetric Evolution of Seismicity 

Understanding the temporal and spatial evolution of brine in the subsurface is 

important for understanding how earthquakes are induced, and how earthquake 

hazards may be affected. For example, whether injected brine can be expected to 

be smoothly distributed within the rock mass near the well bore, or concentrated

along discrete fracture networks opened up during injection, may have a 

significant effect on where future earthquakes are likely to occur and what 

magnitudes can be expected.

Brine migration can be inferred using several methods, including monitoring 

wells, leveling lines, borehole tiltmeters (Jahr and others, 2007; Jahr and others, 

2008), gravimetry, seismic tomography, and seismicity. Of these methods, only 

induced seismicity data is readily available at PVU. Monitoring wells (e.g., new 

wells or re-entered abandoned wildcat wells) have not been available at PVU. 

Geochemical analysis of brine obtained from a wildcat well drilled in 2008 about 

10 km northwest of the PVU well, and which penetrated the Leadville formation,

indicated that PVU brine had not migrated to that location at the time of drilling. 

To our knowledge, no geodetic leveling lines are available in the area around the 

injection well, nor are there any borehole tiltmeters.
39

The relationship between induced seismicity and the migration of injected brine is 

poorly understood at PVU. Simple numerical models have been used to estimate 

the time and distance variations in pressure away from the wellbore (e.g., 

Roeloffs and Denlinger, 2009), but such models are difficult to validate in the 

absence of monitoring wells or other direct measurement. Following Shapiro and 

others (1997; 2003) we simply assume that the state of stress in the rock 

surrounding the injection well is close to critical, and that small changes in that

stress state resulting from diffusion of pore pressure is the direct cause of most

induced earthquakes. Specifically, we assume that induced earthquakes occur in 

response to a decrease in the effective normal stress on pre-existing fracture 

surfaces resulting from the diffusion of pore-pressure transients produced by 

injection (e.g., Shapiro and others, 1997; Shapiro and others, 2002; Shapiro and 

others, 2005). Fracture initiation is assumed to be adequately described by a 

Coulomb fracture criterion (Jaeger, 1969), and we assume that the observed 

seismicity is the result of frictional failure due to shearing (tensile failure events 

may also occur, but have not been observed in the seismicity data).

Following Shapiro and others (2011), we use the volume enclosed by the 

observed seismicity cloud as a working approximation to the total volume of rock 

stimulated by injection, although this may be an over-estimation if the seismicity 

is not homogeneously distributed. The stimulated volume is important because its 

size may be related to the maximum magnitudes of induced seismicity (Shapiro 

and others, 2011).
40

Figure 3-4. Example of minimum ellipsoid fit to hypocenters (grey dots) from the near-well 

cluster. 

We model the volume of the seismicity cloud by determining the volume of the 

minimum ellipsoid enclosing all hypocenters within an identified cluster (Figure 

3-4). Only high-quality, relatively-relocated hypocenters are used for this 

analysis. To determine minimum ellipsoids, we first selected well-defined, dense 

clusters of events, and ignored isolated hypocenters. This process is necessarily 

somewhat subjective, so alternative selections of hypocenters were considered, as 

shown in Figure 3-5. The final results do not appear to be strongly dependent on 

the hypocenter selections used.
41

Figure 3-5. Map showing the surface projections of the minimum ellipsoids fitting the nearwell and NW clusters. The near-well data is modeled with two data sets: (1) all events within 

the near-well area (blue); and (2) events selected to be in the highest-density areas (purple). 

The northwest cluster events are treated in three data sets: (1) all events (green); (2) NW 

cluster - north sub-cluster (magenta); and, (3) NW cluster – south sub-cluster (orange).

Selected clusters of events for this analysis are shown in map view on Figure 3-5. 

Two clusters were defined for the near-well seismicity: (1) all events located 

within the near-well region; and, (2) selected events that don’t include small subclusters in the southwest and southeast directions. For the NW cluster, two sub-
42

clusters were defined. In addition, all events within these two sub-clusters were 

considered as a single cluster. Once the clusters were defined, we computed 

minimum ellipsoids containing the cumulative seismicity, by year.

We estimate the evolution over time of the stimulated volumes of rock containing 

the near-well and NW clusters of induced seismicity. These volumes are 

determined by finding the minimum ellipsoids containing the cumulative 

hypocenters, by year, within each cluster. A plot of the volumes of the minimum 

ellipsoids is shown in Figure 3-6. This plot also shows the cumulative volume of 

brine injected.

Figure 3-6. Evolution of seismicity ellipsoids for alternative selections of hypocenters

representing the near-well and the NW clusters (colors correspond to the ellipsoids shown in 

map view on Figure 3-5). Cumulative volume of injected brine is also shown.

The results shown in Figure 3-6 indicate that the near-well and NW-cluster 

seismicity clouds initially grew rapidly in volume (as measured by the computed 

ellipsoids), followed by much longer periods of very slow growth. Expansion of 

the near-well seismicity volume occurred most rapidly during the injection tests 

(1991 – 1995), slowed during the early phase of continuous injection (1996-

2000), and then nearly halted after the injection flow rate was decreased in 2000. 

Beginning in about 2011, the near-well seismicity volume began to again grow 

slowly. The volume trend for the NW cluster (and its sub-clusters) is similar to 

that for the near-well region, although delayed in time due to the distance of the 

earthquakes from the injection well.

Seismicity Ellipsoid Volume vs. Injected Volume

1.E+04

1.E+05

1.E+06

1.E+07

1.E+08

1.E+09

1.E+10

1.E+11

1.E+12

199219931994199519961997199819992000

2001

2002

2003

2004

2005

2006

200720082009201020112012

2013

Date

Volume (m^3)

Near-well (ext.) Vol

Near-well Volume

NW-cluster Volume

NW-cluster N+S Vol

NW-cluster S Volume

NW-cluster N Volume

Cumulative Injected Volume
43

During the periods of time following the initial rapid expansion of volume, 

earthquakes continued to occur within each of the seismicity clouds. These events 

have served essentially to “fill-in” the enclosing volume without significantly 

increasing its overall size. After the initial expansion, earthquakes were either 

occurring in the same locations as previous events, on different parts of 

previously-active faults, or on newly-active faults elsewhere within the same

volume. Hence, although the enclosing volume of stimulated rock appears to have 

stabilized, the seismicity distribution occurring within this volume has become

more complex over time rather than tending towards uniformity.

The volume of the NW cluster seismicity cloud is much less than the near-well 

seismicity cloud. Considering the NW cluster as a single ellipsoid, its volume is 

about 7% of the volume of the near-well seismicity volume. Considering the NW 

cluster as two sub-clusters, the sum of the volumes of the sub-clusters is about 3% 

of the near-well volume. One hypothesis for the cause of the NW cluster

seismicity is that it represents the terminus of a high-permeability zone (or 

perhaps a series of interconnected zones) that connects with the well (Ake and 

others, 2005).

3.2 Seismicity Rates and Magnitude Distribution

Seismic hazard is determined by the rates at which earthquakes capable of 

producing damage occur in proximity to sites of interest. These rates can be 

expressed in terms of an overall rate of occurrence of earthquakes equal to or 

greater than some minimum magnitude Mmin , and the magnitude probability 

distribution function, which defines the relative likelihood of earthquakes 

occurring at different magnitudes ( ≥ Mmin ). This separation is useful in seismic 

hazard analysis of naturally-occurring earthquakes because larger-magnitude 

earthquakes tend to occur much less frequently then smaller-magnitude quakes, 

making their rates difficult or impossible to determine from direct observations of 

historical seismicity. In contrast, smaller-magnitude earthquakes tend to occur 

much more frequently, with the relative proportion of small to large-magnitude 

earthquakes typically following an exponential or “Gutenberg-Richter” 

distribution. The parameters defining the magnitude distribution can be 

determined from the smaller-magnitude earthquakes, thus providing a means for 

indirectly estimating the rates of the infrequent, larger-magnitude events.

The magnitude distribution of induced seismicity may not necessarily follow that 

of naturally-occurring earthquakes, since the overall seismicity rate depends on 

the injection history in a complicated way. For induced seismicity, both the 

occurrence rate and the magnitude probability distribution may have a strong time 

dependence. In the following sections, we examine the observed seismicity rates 

and magnitude distributions of induced seismicity at PVU.
44

3.2.1 Seismicity Rates

The cumulative seismicity rate (over all magnitudes) is a critical parameter in 

seismic hazard investigations. Given a specified probability distribution of 

earthquake magnitudes, the higher the total rate of seismicity, the more frequently 

that large, potentially damaging earthquakes are expected to occur.

The rate of PVU-induced seismicity for 0.5 MD ≥ has varied significantly over 

time, as shown in Figure 3-7. This plot shows the number of shallow earthquakes 

recorded per quarter during PVU injection operations, all of which we assume to 

be induced. Only earthquakes with magnitude 0.5 or greater are included in the 

plot. Seismicity rates are shown for earthquakes occurring at 3 distance ranges 

from the injection well: up to 5 km, 5 to 10 km, and greater than 10 km. The 

number of 0.5 MD ≥ earthquakes per quarter has varied from 0 to more than 250. 

In general, the highest seismicity rates occurred early during injection operations 

(prior to 2000), with lower rates occurring in the last 12 years. However, a very 

high seismicity rate (over 220 events per quarter) occurred relatively recently, 

when a swarm of seismic activity occurred in the northern valley area during the 

last quarter of 2010.

The temporal patterns in seismicity rate also differ from region to region. The 

near-well region (≤ 5 km distance) shows a large variation in rate. Seismicity rates 

in this region also show some correlation with injection parameters (Block and 

Wood, 2009; Block and Wood, 2010). Seismicity rates at distances of 5 to 10 km 

from the well show much less temporal variation in rate than the near-well region. 

Seismicity rates in these areas also show much less correlation with injection 

operations. Shallow, likely-induced earthquakes occurring at greater than 10 km 

distance from the well have historically had very low seismicity rates. However, 

rates spiked in late 2010 and early 2011 when several swarms of higher seismic 

activity occurred. 
45

Figure 3-7. Seismicity rates as a function of time for earthquakes occurring at distances of 

less than 5 km, 5 to 10 km, and greater than 10 km from the injection well. Shallow (<= 8.5 

km depth), likely-induced earthquakes with magnitude of M 0.5 or greater are included.

3.2.2 Gutenberg-Richter Magnitude Distribution

In order to estimate the frequency of occurrence of large-magnitude earthquakes,

given an estimate of the cumulative seismicity rate, the earthquake probability

distribution over magnitude must be determined. This distribution quantifies the 

relative numbers of earthquakes occurring as a function of magnitude. Earthquake 

frequency-magnitude distributions for naturally-occurring earthquakes are often

characterized using the Gutenberg-Richter (GR) relation (Gutenberg and Richter, 

1954), which can be expressed as (e.g., Bender, 1983):

[ ]

( ) [ ] /2 /2

10 / 2

10 10 10 / 2 i

a bM

Cumulative i

bM bM a bM

Incremental i

for M M M

for M M M

λ

λ

−

∆ −∆ −

= ∈ ±∆

= − ⋅ ∈ ±∆ (3.1)

where λCumulative represents the mean rate at which earthquakes of magnitude M or 

greater occur, λIncremental represents the mean rate at which earthquakes within 

magnitude bin Mi of width ΔM occur, and a and b are constants. The GR relation 

has been applied to induced seismicity in many areas (e.g., Shapiro and others, 

2007; Shapiro and Dinske, 2009; Shapiro and others, 2010; Bachmann and others, 

2011; Shapiro and others, 2011; Convertito and others, 2012; Barth and others, 

2013; Mena and others, 2013). Unlike tectonic earthquake applications, these 

applications of the GR relation to induced seismicity have mostly been done 

without declustering of the earthquake catalog.

For tectonic earthquakes, the simple GR relation typically is modified to account 

for a minimum magnitude of significance (or detection) , and a maximum Mmin
46

magnitude of occurrence Mmax (e.g., Cornell and Vanmarcke, 1969; Anderson, 

1979; Molnar, 1979; Berrill and Davis, 1980; Anderson and Luco, 1983):

( )

( )

( ) [ ]

min max

max

min

min max

max

/2 /2

min max

10 10 10

10 10 10

0

10 10 10 / 2

0 ,

i

a bM bM

a

a bM bM

Cumulative

bM bM a bM

i

Incremental

for M M

for M M M

for M M

for M M M

for M M M M

λ

λ

λ

− −

− −

∆ −∆ −

 ≡ ⋅ − ≤ 

 = ⋅ −  <≤



 > 





 − ⋅ ∈ ±∆ = 

 < >

(3.2)

where λa represents the cumulative event rate for M M≥ min .

The corresponding complimentary cumulative probability for magnitude,

( ) C P M , is therefore a truncated exponential distribution. This and its associated 

probability density function ( ) mp M can be written in terms of the GR parameters 

as follows (Bender, 1983):

( )

( )

( )

( )

min

max min

min

max min

min

min max

max

min max

1

1 ( )

1

0

( ) ( ) 1

0

M M

C

M M

M M

C

M M m

for M M

e P M for M M M

e

for M M

e dP M for M M M p M e dM

otherwise

β

−β −

−β −

−β −

−β −

 ≤



 − =  <≤

 −

 > 



 <≤ =− =  − 



(3.3)

whereβ = b ln10 . With this definition, the mean earthquake rates in Equation 

(3.2) can be rewritten in terms of the GR magnitude probability densities:

( ) /2 /2

( )

( )

C

Cumulative a

m m

Incremental a m

P M

e e

p M

β β

λ λ

λ λ

β

∆ −∆

= 

− =

(3.4)

B-values for earthquakes observed within the near-well and NW-cluster areas 

were estimated using three methods. The first method simply computes b-values 

from the logarithm of the ratios of the cumulative number of earthquakes at two 

magnitudes M M 1 2 < , where M1 and M2 are chosen from inspection of the data 

so that the GR is approximately valid for M M≥ 1 . Similar approaches have been 

used elsewhere (e.g., Frohlich, 1989; Frohlich and Davis, 1993). Neglecting the 

effect of , this method yields: Mmax
47

b N MN M ≈ log10   Cumulative ( 1 ) Cumulative ( 2 )   (3.5)

where the total number of events of magnitude M or greater occurring in time T

is given by N M Nm M m M T Cumulative ( ) = ≥= ≥⋅ ( ) λa ( ) . Because of the limited 

range of magnitudes available in the PVSN data, we chose 2 1 M M= +1.0 for a 

range of M1 values greater than an assumed of 0.5. The time windows used 

to compute N M( ) are the same for each magnitude, and various time windows 

were tried.

The second method uses the maximum-likelihood estimate for the simple 

Gutterberg-Richter relation, which neglects the effect of (e.g., Aki, 1965; 

Utsu, 1966):

b eMM = − log10 ( min ) (3.6)

where M is the average of the observed magnitudes 

1

N

i

i

M MN =

= ∑ for 

M M i ≥ min .

Lastly, we employ a maximum-likelihood method which corrects for both 

magnitude binning and (Bender, 1983):

1

( 1)

1 1

n n

i

n

i

q nq i k

qq N =

− − = − − ∑ (3.7)

where q m = −∆ exp( ) β , n is the number of magnitude bin intervals, and i k is the 

number of earthquakes in a magnitude bin. The minimum and maximum 

magnitudes are implicit in this formulation, and the b-values are solved for

numerically.

The error in estimated b-value can be estimated using two methods that do not 

have significant biases from binned magnitudes (Marzocchi and Sandri, 2003). 

Following the work of Shi and Bolt (1982):

( ) ( )

( )

2

2

1

log 10 1

N

i

b

i

M M

b

N N

σ =

− = − ∑ (3.8)

Following the work of Tinti and Mulargia (1987):

Mmin

Mmax

Mmax
48

min

1 ; 1

log(10) b

p m p

m Np M M

σ

− ∆   = = +   ∆   − (3.9)

When using the simple ratio method, we computed b-values using 5 sets of 

minimum magnitude pairs ( , M2 ): (0.7, 1.7); (0.8, 1.8); (0.9, 1.9); (1.0, 2.0); 

and (1.1, 2.1). All 5 of these sets were used for analyzing data from the near-well 

region. The smallest magnitude pair (0.7, 1.7) was omitted when analyzing the 

NW-cluster data because of probable slightly higher event detection threshold for 

areas farther from the well. For both maximum-likelihood methods, values 

of 0.7, 0.8, 0.9, 1.0, and 1.1 were used for the near-well region, and values of 0.8, 

0.9, 1.0, and 1.1 were used for the NW cluster. For the purpose of comparing 

results from the corrected maximum likelihood method to those from the other 

two methods, Mmax was initially set to a value of 10 to have minimal effect.

(Results using different values of Mmax are presented later in this section.)

B-values were computed using each of the three methods using running time 

windows over the near-well and NW-cluster data, incorporating data from the 

beginning of continuous injection (July, 1996) through the end of 2012. A 

window length of 1.5 years was used for the near-well region. Because the 

seismicity rate is much lower in the NW cluster than in the near-well region, a 

much longer time window, 3 years, was used to obtain robust b-value estimates 

for that region. B-value estimates are not presented for time periods when the 

seismicity rates are too low to obtain robust estimates.

The two maximum-likelihood methods were also used to compute b-values for 

the seismicity data recorded during the injection tests (1991 to 1995). The ratio 

method was not used for the injection test data because for the earliest injection 

tests the number of events induced was too low to produce reliable results using 

this technique. Seismic data from the first year of injection testing (July, 1991 to 

May, 1992) were combined into one data set to obtain enough data points for 

reliable estimation. The data from each of the last three injection tests (October, 

1993 to April, 1995) were analyzed separately. A minimum magnitude value of 

0.5 was used for these computations, to incorporate more data points.

M1

Mmin
49

Figure 3-8. B-values estimated by two maximum-likelihood methods for the earthquakes 

induced during the injection tests (1991 – 1995) (upper plot). The lower plot shows the 

induced earthquakes as a function of date and magnitude; the area of each circle is scaled by 

the number of earthquakes occurring in a given quarter-year and magnitude bin.
50

Figure 3-9. B-values estimated by three methods for the earthquakes induced in the nearwell region (within 5 km of the wellhead) during continuous injection operations (1996 –

2012). The lower plot shows the induced earthquakes as a function of date and magnitude; 

the area of each circle is scaled by the number of earthquakes occurring in a given quarteryear and magnitude bin.
51

Figure 3-10. B-values estimated by three methods for the earthquakes induced in the 

northwest (NW) cluster during continuous injection operations (1996 – 2012). The lower plot 

shows the induced earthquakes as a function of date and magnitude; the area of each circle 

is scaled by the number of earthquakes occurring in a given quarter-year and magnitude 

bin. 
52

The b-values computed for the injection test earthquake data are presented in 

Figure 3-8, while the results for the period of continuous injection are shown in 

Figure 3-9 (near-well region) and Figure 3-10 (NW cluster). The uppermost plot 

in each figure compares the results from the three estimation methods used. The 

three methods produce fairly consistent results, and show that the b-value changes 

systematically over time.

For the near-well region, b-values generally decrease from about 1.3 to 1.5 at the 

beginning of the injection tests (Figure 3-8) to about 0.6 to 0.7 at the end of 2000, 

after continuous injection had been underway for several years (Figure 3-9). 

From 2001 through 2010, the near-well b-value estimates are fairly constant, 

fluctuating mainly between 0.6 and 0.8. Beginning in early 2011, however, the bvalue estimates show an abrupt increase to about 0.95 to 1.0 by the end of 2011. 

Since 2011, seismicity rates in the near-well region have been too low to compute 

robust b-value estimates. The abrupt increase in b-value in 2011, if valid, would

represent a marked change in the long-term trend in the near-well region. It may 

be too soon to definitively reach such a conclusion, however. Additional data 

gathered over the coming months and years is needed to confirm this change in bvalue trend.

B-value estimates for the NW cluster (Figure 3-10) also show temporal variations. 

As observed for the near-well region, b-values in the NW cluster also initially 

decreased, from about 0.7 to 0.8 when earthquakes first began occurring in the 

NW cluster, to about 0.55 to 0.65 by the end of 1999. B-values then gradually 

increased, to a maximum of about 0.7 to 1.0 by the end of 2004. Since then, bvalues have generally decreased, to about 0.5 to 0.7 in 2011. Unlike the near-well 

region, no abrupt change in b-value is observed for recent years, although because 

of the longer time window used for the analysis (3 years), any shorter-term 

temporal variations would be smoothed. 

Similar temporal variations in b-value have been reported elsewhere. For 

example, Lei and others (2008) reported significant temporal variations in 

computed b-values for induced seismicity recorded over a period of several years. 

Their data show decreasing b-value over the first 5 to 6 years, from values of 

about 1.3 – 1.5 to 0.7 - 0.9, followed by a relatively stable period of about 6 to 7 

years, similar to the trend seen in the PVU near-well region. Their hypothesis, 

supported by results of epidemic-type aftershock sequence (ETAS) modeling, is 

that the variation in b-value is related to variation in the relative contributions of 

fluid-driven (pore pressure-induced) seismicity and stress-triggered aftershock 

seismicity. According to this model, as the stress level is reduced by induced 

earthquakes, the ability to trigger aftershocks declines and therefore fluid-driven 

seismicity becomes more dominant as injection progresses. Variation in the 

relative contributions of pore pressure-induced and stress-induced earthquakes 

may also be an important source of the temporal changes in b-values observed at 

PVU, but to date no analyses have been completed to support the hypothesis for 

this site. 
53

The value of the maximum possible magnitude, , can have a substantial

effect on the computed maximum-likelihood b-value (Figure 3-11a). The trends in 

b-value, however, remain largely the same. Lower b-values are calculated when 

assumed maximum magnitudes are small (i.e. < 4), however similar b-values are 

computed for max M ≥ 4.5. A larger difference in the computed b-value occurs 

when the minimum magnitude is changed (Figure 3-11b). This is either the result 

of an incomplete earthquake catalog, or the fact that Paradox Valley earthquakes 

may not follow Gutenberg-Richter statistics at lower magnitudes. Further

evaluation of the PVSN event detection threshold as a function of time is needed

to determine which factor is the primary cause. Due to the difference in b-values 

from distinct and inputs, a median b-value was also computed using a 

large value of (6.0), and = 0.8, 0.9, 1.0, and 1.1 (Figure 3-11c). Lastly,

the effect of bin window length was investigated (Figure 3-11d). The general 

trend of decreasing near-well b-value with time holds for all bin windows, and b 

values remain similar regardless of window length.

Mmax

Mmin Mmax

Mmax Mmin
54

a)

b)
55

Figure 3-11. (a) The effect on near-well b-value estimates when changing the maximum 

magnitude (Mmax = 3.5, 4.5, 5.5, 7.0) with constant Mmin (0.8), following Bender’s (1983) 

maximum-likelihood calculations. The gray lines show the full range of b-values calculated 

c)

d)
56

for a range of Mmin and Mmax values. (b) The effect on near-well b-value estimates when 

changing the minimum magnitude (Mmin = 0.8, 0.9, 1.0, 1.1) with constant Mmax (4.5). (c) 

The median b-value calculated for the near-well region with Mmin = 0.8, 0.9, 1.0, 1.1 and 

Mmax = 6. This median value was used in seismogenic index calculations. (d) The effect of 

bin window length for fixed Mmax and Mmin. In all plots, b-values with errors > 0.1 are 

removed.

In order to further investigate how closely the Paradox Valley catalog follows 

Gutenberg-Richter statistics, we generated cumulative magnitude-frequency plots 

for each time interval for which a- and b-values were calculated. The cumulative 

data for each time interval were then compared with the corresponding maximumlikelihood magnitude-frequency curves (Figure 3-12 and Figure 3-13).

While the maximum-likelihood curves fit the data over intermediate magnitude 

ranges, there is a distinct roll-off in the frequency of low and high magnitude 

earthquakes. The roll off of low magnitude earthquakes varies over time, likely

due either to changing event detection capabilities or a discrepancy in the 

occurrence of low-magnitude earthquakes from Gutenberg-Richter statistics. The 

roll off of high-magnitude earthquakes is likely due to either an actual 

effect, or potential saturation of the duration-magnitude scale above about MD 3.0. 

In either case, the effect of the high-magnitude roll off has a small effect in 

maximum-likelihood b-value calculations since greater weight is given to the 

lower-magnitude event counts.

Regardless of the root causes of the magnitude roll-offs, it is clear that the shape 

of the frequency-magnitude cumulative distribution has evolved over the lifetime 

of the PVU. While the slope of these plots (b-value) has clearly evolved, so has 

the maximum and minimum magnitude roll-offs. It is also clear that the 

magnitude-frequency distributions are complex, and have varied over time.

Mmax
57
58
59
60

Figure 3-12. Cumulative frequency-magnitude plots for two year windows of near-well 

events with maximum-likelihood solution (dashed black line) following Bender (1983) 

(Mmin=0.8, Mmax = 4.3). The b-value for each fit is shown in the upper right-hand corner of 

the plot.
61
62
63

Figure 3-13. Cumulative frequency-magnitude plots for two year windows of north-west 

cluster events with maximum-likelihood solution (dashed black line) following Bender (1983) 

(Mmin=0.8, Mmax = 4.3). The b-value for each fit is shown in the upper right-hand corner of 

the plot.

3.2.3 Seismogenic Index

The seismogenic index (Shapiro and others, 2010; Shapiro and others, 2011)

provides a decomposition of the cumulative event rate of earthquakes with

M M≥ min , denoted NM t a min ( , ), into two parts: (1) a component that depends 

primarily on the tectonic properties of the injected region; and, (2) a component 

that depends primarily on the volume injected over some time interval. This 

decomposition can potentially provide useful estimates of future induced 

seismicity rates, based on an assumed injection scenario, at least for those 

injection projects where there is good correlation between the number of induced 

events during a time interval versus the volume injected during the same interval.

The seismogenic index Σ(t)represents the first part of this decomposition, and is 

designed to depend primarily on properties of the injection site, including: the 

concentration of pre-existing cracks, the poroelastic storage coefficient, and the 

maximum critical pore-pressure needed to initiate faulting (Shapiro and others, 

2010). The seismogenic index is important to determine for PVU because it 

provides an indication about whether properties of the injection site have changed

over time, and whether the rate of induced earthquakes is well-correlated with the 

rate of injection. Both of these factors affect whether the rate of future induced 

earthquakes can be reasonably inferred from injection rates. Determination of the 
64

seismogenic index also allows for comparison of PVU induced seismicity with 

that of other injection sites, after correcting for differing injection volumes.

Following Shapiro and others (2010; 2011), a time-dependent, cumulative

seismogenic index can be defined as:

( ) ( ) ( )

( ) ( )

10

10 10

log ,

log , log

a min

a min

t N M t Qt

N M t Qt

Σ =    

= − (3.10)

where t is the elapsed time since the start of injection, NM t a min ( , )is the 

cumulative number of events of magnitude or larger at time , as defined in 

Equation (3.2), andQ t( ) is the cumulative injected volume (in units of m3).

Figure 3-14. Cumulative seismogenic Index calculated for the near-well region, for a 

range of minimum magnitudes ranging between 0.5 and 1.5. Two definitions of 

cumulative seismogenic index are used: (1) uncorrected for (“simple”), as defined by 

Equation (3.10) – dashed lines; and, (2) corrected for – solid lines, as defined by 

Equation (3.11) using the time-dependent b-values shown in Figure 3-11c. Because the 

seismogenic index corrected for depends on the estimated b-value, the uncertainty in 

the estimated b-value maps into the computed seismogenic index, leading to considerable 

fluctuation in the estimate. Time periods where b-values have high errors (e.g. the large bvalue maximum in 2008) have been excluded.

Equation (3.10) decomposes the logarithm of the total event rate into two 

components by simply subtracting off the cumulative injected volume from the 

total number of events. Using the simple definition from Equation (3.10), the 

cumulative seismogenic index was estimated for the near-well region using 

Mmin t

Σ(t)

Mmin

Mmin

Mmin

Mmin
65

possible values of ranging between 0.5 to 1.5. These estimates are shown as 

the smooth dashed lines in Figure 3-14, and indicate that the cumulative 

seismogenic index at PVU increased after the start of continuous injection, 

reached a peak around 2000, and has gradually decreased since that time.

The overall level of the simple seismogenic index defined by Equation (3.10)

depends on the value selected for , as illustrated in Figure 3-14. This 

dependence can be removed by subtracting a factor − ⋅ bM t M ( min min , ) from 

Equation (3.10) to obtain a second definition of seismogenic index that is 

corrected for (Shapiro and others, 2010):

Σ = (t N M t Qt bM t M ) log , 10   a min ( ) ( ) + ⋅ ( min min , )   (3.11)

where bM t ( min , ) is the time-dependent b-value (calculated for a specified ). 

While this second definition attempts to correct for the dependence inherent 

in the first definition, there is a practical difficulty in computing it: since b is not 

known exactly, an estimated b-value must be used, which may have significant 

uncertainty that maps into the computed seismogenic index value.

The cumulative seismogenic index defined by Equation (3.11) was calculated for 

the near-well region using the same values of (0.5 to 1.5) as before, which 

are shown as the solid lines in Figure 3-14. The resulting curves, while grouped 

closer together, show considerably more fluctuation over time, which is due to 

uncertainty in the estimated b-value. The overall trends of the cumulative 

seismogenic index are similar, independent of which definition of seismogenic 

index is used. Figure 3-14 shows that cumulativeΣ has decreased gradually since 

about 2000. There is a suggestion that the corrected seismogenic index defined by

Equation (3.11) has slightly increased since 2011, but this trend results only from 

a recent change in the estimated b-value, and will need to be confirmed by 

analyzing data collected over the next several years. The uncorrected curves do 

not show such an increase.

The cumulative seismogenic index defined by Equation (3.11) becomes 

insensitive to changes over time because it depends on accumulated quantities –

both the injected volume and the number of induced events. At PVU, however, 

the record of induced seismicity response to injection volume appears to show 

significant variations over time (Block and Wood, 2010), which are not reflected 

in the cumulative seismogenic index. To examine possible time variations in 

seismogenic index, we therefore define an instantaneous (and b-value corrected) 

seismogenic index, Σ ∆ ′(t t ; ) , for a lag-time ∆t , as follows:

Mmin

Mmin

Mmin

Mmin

Mmin

Mmin
66

( ) ( ) ( )

( ) ( ) ( )

( )

( ) ( )

10

10

, , ; log , ;

, ; log , ;

;

a min a min

min min

a min

min min

NM t NM t t

t t bM t t M

Qt Qt t

NM tt

bM t t M

Qt t

  − −∆

Σ ∆= ′   + ∆⋅

− −∆  

  ∆ ∆ =   + ∆⋅

∆ ∆  

(3.12)

where bM t t ( min , ;∆ ) is the average b-value over the time interval [tt t , + ∆ ]. Using 

a time window of 1.5 years, the instantaneous seismogenic index computed for 

the near-well region is shown in Figure 3-15, along with the cumulative 

seismogenic index, and for comparison, the logarithms of the constituent event 

and flow rates, ∆ ∆∆ NM tt t a min ( , ; ) and ∆ ∆∆ Qt t t ( ; ) , respectively. 

Figure 3-15. Near-well instantaneous (upper) and cumulative (middle) seismogenic index for 

a range of values between 0.7 and 1.1. Each plotted point represents a centered 1.5-

year moving average, with points spaced every 90 days. Data points removed for 2007-2008 

time interval having a large b-value uncertainty. Lower plot shows a comparison of the 

logarithm of the average injection rates and event rates, which are well-correlated only until 

about 2002.

Mmin

Near-well

-4.2

-3.8

-3.4

-3.0

-2.6

-2.2

Jan-96 Jan-97 Jan-98 Jan-99 Jan-00 Jan-01 Jan-02 Jan-03 Jan-04 Jan-05 Jan-06 Jan-07 Jan-08 Jan-09 Jan-10 Jan-11 Jan-12 Jan-13

Inst. Seismogenic Index

ML corr 1.1

ML corr 1.0

ML corr 0.9

ML corr 0.8

ML corr 0.7

median

-4.2

-3.8

-3.4

-3.0

-2.6

-2.2

Jan-96 Jan-97 Jan-98 Jan-99 Jan-00 Jan-01 Jan-02 Jan-03 Jan-04 Jan-05 Jan-06 Jan-07 Jan-08 Jan-09 Jan-10 Jan-11 Jan-12 Jan-13

Cum. Seismogenic Index

ML corr 1.1

ML corr 1.0

ML corr 0.9

ML corr 0.8

ML corr 0.7

median

5.3

5.4

5.5

5.6

5.7

5.8

5.9

Jan-96 Jan-97 Jan-98 Jan-99 Jan-00 Jan-01 Jan-02 Jan-03 Jan-04 Jan-05 Jan-06 Jan-07 Jan-08 Jan-09 Jan-10 Jan-11 Jan-12 Jan-13

Log10 Flow Rate (m^3/yr)

0.0

0.5

1.0

1.5

2.0

2.5

3.0

Log10 Event Rate (#/yr)

dQ/dt

dN/dt - Mmin 0.7

dN/dt - Mmin 0.8

dN/dt - Mmin 0.9

dN/dt - Mmin 1.0

dN/dt - Mmin 1.1
67

The near-well instantaneous seismogenic index shown in Figure 3-15 suggests 

that the seismicity response to injection flow rates has varied considerably over 

time at PVU. After about 2000, the number of events produced per unit volume 

injected decreased substantially, until about 2002, when it again increased. After 

2004, there was another decrease, followed by further fluctuations. The observed 

large variability over time in seismicity response to injection rates may be due to 

changing reservoir properties, to a system that is highly sensitive to small changes 

in conditions, or to a dependence of seismicity rates on other injection parameters, 

such as pressure. In any case, the large variability makes it very difficult to 

directly use injection flow rates as the sole indicator of induced seismicity rates.

Figure 3-16. Northwest cluster instantaneous (upper) and cumulative (middle) seismogenic 

index for a range of values between 0.8 and 1.1. Each plotted point represents a

centered 3-year moving average, with points spaced every 90 days. Lower plot shows a 

comparison of the logarithm of the average injection rates and event rates, which are only 

generally correlated. The 1.5-year-average injection rate from Figure 3-15 is also shown for 

comparison.

The instantaneous seismogenic index for the NW cluster is shown in Figure 3-16, 

and is calculated based on the total volume of fluid injected. This plot suggests 

that the seismicity response to injection at the NW cluster was fairly stable until 

Mmin

-4.6

-4.4

-4.2

-4.0

-3.8

-3.6

-3.4

-3.2

Jan-96 Jan-97 Jan-98 Jan-99 Jan-00 Jan-01 Jan-02 Jan-03 Jan-04 Jan-05 Jan-06 Jan-07 Jan-08 Jan-09 Jan-10 Jan-11 Jan-12 Jan-13

Cum. Seismogenic Index

ML corr 1.1

ML corr 1.0

ML corr 0.9

ML corr 0.8

median

5.3

5.4

5.5

5.6

5.7

5.8

5.9

Jan-96 Jan-97 Jan-98 Jan-99 Jan-00 Jan-01 Jan-02 Jan-03 Jan-04 Jan-05 Jan-06 Jan-07 Jan-08 Jan-09 Jan-10 Jan-11 Jan-12 Jan-13

Log10 Flow Rate (m^3/yr)

0

0.3

0.6

0.9

1.2

1.5

1.8

Log10 Event Rate (#/yr)

dQ/dt - 3yr

dQ/dt - 1.5 yr

dN/dT - Mmin 0.8

dN/dT - Mmin 0.9

dN/dT - Mmin 1.0

dN/dT - Mmin 1.1

Northwest Cluster

-4.6

-4.4

-4.2

-4.0

-3.8

-3.6

-3.4

-3.2

Jan-96 Jan-97 Jan-98 Jan-99 Jan-00 Jan-01 Jan-02 Jan-03 Jan-04 Jan-05 Jan-06 Jan-07 Jan-08 Jan-09 Jan-10 Jan-11 Jan-12 Jan-13

Inst. Seismogenic Index

ML corr 1.1

ML corr 1.0

ML corr 0.9

ML corr 0.8

median
68

about 2005, when it decreased and began to fluctuate. It decreased again from 

about 2008-2011, and has been increasing since 2011. The proportion of the total 

injected fluid influencing the NW cluster is not known, nor is it known if that 

proportion has changed over time. Like the near-well region, the seismicity 

response to injection rate has been variable over time, and may not provide a 

reliable indicator of future seismicity rates given an assumed injection rate.

3.3 Maximum Magnitude of Induced Earthquakes

The maximum possible magnitude of induced earthquakes at Paradox Valley has 

a significant effect on the likely seismic hazard resulting from injection 

operations, since strong ground shaking tends to increase with increasing 

earthquake magnitude. Maximum magnitude in turn depends on several key 

factors, such as the total volume and average rate of injected brine, the location 

and characteristics of preexisting fractures, and the far-field ambient stress. We 

examine several methods that have been used for estimating the maximum 

magnitude of induced earthquakes at other injection sites.

3.3.1 Geometrical scaling of maximum magnitude

Shapiro and others (2011) observed that the maximum magnitudes of induced 

seismicity appear to be related to the size of the observed cloud of induced 

seismicity. In their simplified model, the observed seismicity cloud is assumed to

provide a working approximation to the current total volume of rock stimulated 

by injection. Potential fault rupture surfaces for induced earthquakes are assumed 

to be represented by a set of pre-existing circular disks with random orientation, 

and whose centers are randomly distributed throughout the rock volume. 

Earthquakes are assumed to be possible on any of these pre-existing fault surfaces 

once the stimulated volume has grown in size such that an entire fault surface is 

contained within the stimulated volume. This model is consistent with the 

assumption that pore pressure perturbations must extend over most of a potential 

fault surface in order for rupture to propagate. Shapiro and others (2011) analyzed 

induced seismicity from a number of fluid stimulation sites and found that the 

geometrical scale of the fluid-stimulated volume seems to control the magnitude 

of the largest observed earthquake. In particular, they found that the length of the 

minimum principal axis of the stimulated volume appeared to control maximum 

rupture size.

For a fault rupture size X (diameter) of a pre-existing fault, Shapiro and others

(2011) provide a relationship between earthquake magnitude M and static stress 

drop ∆σ of the earthquake as follows:

( ) 2

10 10 10 MX C = + ∆σ − log log log 1.5 6.07 − (3.13)
69

In this relation,C is a constant on the order of one, and all units are metric.

For the case of a circular crack of radius (which appears to be a reasonable 

model for induced earthquakes occurring at PVU) the static stress drop is 

related to seismic moment via (e.g., Brune, 1970, 1971; Kanamori 

and Anderson, 1975). Using the Hanks and Kanamori (1979) definition of 

moment magnitude, 10 0 log 1.5 9.05 M M = + , equation (3.13) may be rewritten

10 10 ( ) 2

10

log log 7 16 9.05 log 1.5

M R

∆σ − −

= + (3.14)

where R X = 2 , and the units of the radius and stress drop are metric. This 

relation indicates that magnitude scales with the logarithm of rupture area, for the 

case of constant stress drop, consistent with empirical magnitude versus rupture 

area relations for natural earthquakes (e.g., Wells and Coppersmith, 1994; Hanks 

and Bakun, 2002; Ellsworth, 2003; Shaw, 2009).

Plots of estimated maximum magnitude versus time, for the near-well and NWcluster regions, are shown in Figure 3-17 and Figure 3-18, respectively. 

Maximum magnitude at a given time is estimated using Equation (3.14),

assuming a value of R t( ) equal to the semi-minor axis Amin of the minimum 

ellipsoid for the seismicity cloud (Shapiro and others, 2011) as estimated in 

Section 3.1.2. As discussed in the next section, to account for uncertainty in stress 

drops in the near-well area we have considered potential static stress drops 

between 1 and 20 MPa (10 to 200 bars).

R

∆σ

0

3

7

16

M

R

∆σ =
70

Figure 3-17. Estimated maximum magnitude for the near-well area, computed using the 

approach of Shapiro and others (2011) for a circular crack and various values of static stress 

drop. The radius Amin of the semi-minor axis of the ellipsoid is shown by the blue curve. 

Observed maximum magnitudes are shown by red dots.

Considering just the near-well area, and assuming a typical static stress drop of 5 

MPa (50 bars), the computed value of maximum magnitude reached MW 5.2 in

2001, and then remained relatively constant since that time (Figure 3-17). In 

contrast, the largest observed magnitude in the near-well area was the ML 4.3 (MW

~3.8) earthquake of May 27, 2000. Considering just the NW cluster (Figure 3-18), 

the computed value of maximum magnitude reached MW 4.4 in 2000, and then 

remained relatively constant since (assuming a static stress drop of about 50 bars).

The largest observed magnitude in the NW cluster was the MW 4.0 earthquake of 

January 24, 2013. For both the near-well and NW cluster events, the estimated 

maximum magnitude is significantly higher than the maximum magnitude

observed to date, with the greatest discrepancy in the near-well area. It is not 

known if this discrepancy is the result of Equation (3.14) not being applicable to 

PVU seismicity, or if it is more a matter of insufficient time for the maximum 

magnitude earthquake to occur.

Near-well (selected) Results

1

2

3

4

5

6

Jan-92

Jan-93

Jan-94Jan-95

Jan-96

Jan-97

Jan-98Jan-99Jan-00

Jan-01

Jan-02Jan-03Jan-04

Jan-05Jan-06

Jan-07

Jan-08

Jan-09

Jan-10Jan-11Jan-12Jan-13

Date

Magnitude

0

1

2

3

4

5

Radius (km)

200 bars

100 bars

50 bars

20 bars

10 bars

Observed Mmax

Amin (km)
71

Figure 3-18. Estimated maximum magnitude for the NW-cluster area, computed using the 

approach of Shapiro and others (2011) for a circular crack and various values of static stress 

drop. The radius Amin of the semi-minor axis of the ellipsoid is shown by the blue curve. 

Observed maximum magnitudes are shown by red dots.

To use Equation (3.14) to estimate the maximum magnitude of an earthquake 

generated on a circular fault completely enclosed in a stimulated volume, an 

assumption must be made about the range of possible stress drops. Using PVSN 

broad-band data (Section 2.1), the estimated static stress drop for the January, 

2013 earthquake was 4.5 MPa (45 bars), and its estimated radius was 0.54 km. 

Analysis of natural seismicity used to determine magnitude versus area scaling 

relations indicates that earthquake stress drops are relatively constant, from the 

smallest to the largest magnitude events (Shaw, 2009). For induced seismicity, 

however, it is not clear whether the assumption of constant stress drop scaling is 

appropriate. For example, Goertz-Allmann and others (2011) found an increase in 

stress drop with radial distance from the injection site, and also found that stress 

drop correlates with the pore pressure perturbations due to the injection.

We have estimated fault rupture radius versus moment magnitude using Equation 

(3.14) for static stress drops of 2, 5, and 10 MPa (20, 50 and 100 bars). We have 

also estimated fault rupture radius (median and ±1 σ) using several empirical

area-versus-magnitude relations: WC94-RA: (Wells and Coppersmith, 1994), 

Els03-RA: (Ellsworth, 2003), Shw09-RA: (Shaw, 2009), as well as the mean of 

these estimates (equal weight). In using these area-versus-magnitude relations, we 

have assumed a circular fault geometry. We also estimated fault radius using one 

NW Cluster (All) Results

1

2

3

4

5

Jan-98Jan-99Jan-00Jan-01Jan-02Jan-03Jan-04Jan-05Jan-06Jan-07Jan-08Jan-09Jan-10Jan-11Jan-12Jan-13

Date

Magnitude

0

1

2

3

4

Radius (km)

200 bars

100 bars

50 bars

20 bars

10 bars

Observed Mmax

Amin (km)
72

half of the rupture-length-at-depth parameter from WC94-RLD: (Wells and 

Coppersmith, 1994). Plots of these estimates are shown in Figure 3-19, which 

suggests a considerable range in allowable rupture radius for a specified 

magnitude. For example, for a M 5.0 earthquake, allowable radius values using

these parameters range between about 1.1 and 2.1 km. This would suggest that the 

minimum principal axis of the stimulated volume would need to exceed twice this 

range (i.e., about 2.2 to 4.2 km) in order for an M 5.0 induced earthquake to 

occur. This condition appears to be met for the near-well cluster of induced 

seismicity, but not for the NW cluster.

Figure 3-19. Comparison of various measures of fault rupture radius versus moment 

magnitude for an assumed circular rupture. Solid and dashed lines represent radius 

estimates using empirical area-versus-magnitude relations – see text for explanation. Dotted 

lines represent radius estimates using Equation (3.14). The red diamond symbol indicates 

the PVSN estimate of M and R for the January, 2013 earthquake. Heavy line patterns 

indicate: (1) likely range of magnitudes for a maximum radius of 0.75 km (based on the ~1.5

km minimum dimension of the NW cluster); and, (2) likely range of radius values for a M 

5.0 earthquake.

For the NW cluster of induced seismicity, the minimum ellipsoid analysis 

presented in Section 3.1.2 (see Figure 3-5 and Figure 3-6) showed that the 

dimension of the minimum axis of the seismicity cloud is about 0.7 km (using the 

north sub-cluster) to about 1.5 km (using the entire NW cluster). To fit a circular 

fault within this volume, the radius of the fault would need to be about 0.8 km or 

less, suggesting a maximum magnitude between about MW 4.1 and 4.7. This range 

Rupture Radius vs. Magnitude

0.0

1.0

2.0

3.0

4.0

3.5 4.0 4.5 5.0 5.5

Moment Magnitude

Radius (km)

Avg. RA (+1 σ)

Avg. RA

Avg. RA (-1 σ)

Shw 09-RA

WC94-RA

ELS03-RA

WC94-RLD

2013 EQ

20 bars

50 bars

100 bars

Range for Xmin=1.5

Range for M=5
73

exceeds the M 3.9-4.1 moment magnitude estimate of the January, 2013 

earthquake, but suggests that a significantly larger earthquake (e.g., M 5.0) in the 

NW cluster would not be expected without a corresponding increase in the size of 

the stimulated volume to approximately 2.4 km or greater. The dimension of the 

NW cluster has been fairly stable since 2004 (Figure 3-5). Application of these 

estimates to Paradox Valley has considerable uncertainty, however, because the

model is very simple and relies on a number of assumptions.

3.3.2 Fluid volume scaling of maximum magnitude

A correlation between the cumulative injected fluid volume and maximummagnitude induced earthquake has been observed at other injection sites (McGarr, 

1976; McGarr and others, 2002; Nicol and others, 2011). Such a correlation is 

also observed at PVU, as shown in Figure 3-20. Each point on the lower plot in 

Figure 3-20 represents the occurrence of an earthquake having a magnitude larger 

than all previously-induced earthquakes. All earthquakes within 10 km of the 

injection well were included in this analysis, but only the last point on the plot 

(representing the January 24, 2013 ML 4.4 event) is from an earthquake occurring 

more than 2.2 km from the injection well.

The PVU maximum observed earthquake magnitudes approximately fit a 

logarithmic distribution as a function of cumulative injected volume. The least 

squares fit made to the observed data points is included in Figure 3-20 (solid line), 

as well as the 95% confidence intervals (dashed lines). Figure 3-21 shows a 

quantile-quantile (QQ) plot of the residuals of the least-squares fit. The residuals 

approximate a straight line when plotted against the quantiles of a standard 

normal distribution, indicating the residuals are normally distributed, and thus that 

the fit to the data is appropriate.

Since even for naturally occurring earthquakes, the larger-magnitude events occur 

less frequently than the smaller-magnitude events, the longer the observation 

time, the more likely it is that a large-magnitude event will be observed. 

However, for a given time window, the magnitude distribution of earthquakes 

occurring within that window is usually assumed to be time-independent for 

naturally-occurring events. In other words, the likelihood of a large-magnitude 

earthquake occurring early during the observation window is the same as the 

likelihood that it will be observed late in the window. This is not the case for the 

PVU-induced seismicity. Rather, the steady increase in earthquake magnitude 

over time (represented by cumulative volume in Figure 3-20) indicates that during 

PVU fluid injection, subsurface conditions are changing such that largermagnitude induced earthquakes become more likely as injection progresses.

The correlation shown in Figure 3-20 indicates a current estimated maximum 

magnitude of approximately ML 4.2, based on the cumulative volume of fluid 

injected to date. This estimate is only slightly less than the observed maximum 

magnitude of ML 4.4. The range of uncertainty in the estimated maximum 

magnitude is about ML 3.7 to 4.8, based on the 95% confidence intervals. 
74

Figure 3-20. Maximum induced earthquake magnitude as a function of (log) cumulative 

injected fluid volume (lower plot). The solid line is a least squares linear fit to the data. The 

dashed lines represent 95% confidence intervals. The upper plot shows the daily average 

injection flow rate as a function of (log) cumulative volume.

Because the volume-magnitude data are plotted in log-linear space in Figure 3-20, 

time periods early in injection operations have a much greater impact on the plot 

than later time periods. For example, data from the 5 years of injection tests 

(1991-1995) are plotted over more than two-thirds of the volume (x) axis, while 

the subsequent 16 years of continuous injection only cover the last one-third. The 

last 12 years of injection, since the flow rate was reduced in mid-2000, are 

represented by only the last two data points on the plot. These last two points 

represent earthquakes of nearly identical magnitude (the ML 4.3 event of May, 

2000 and the ML 4.4 event of January, 2013). Because the time span since the 

injection flow rate was reduced in mid-2000 is represented by such a short 

segment of the volume-magnitude plot, it is impossible to interpret whether the 

long-term reduction in flow rate implemented in 2000 has changed the nature of 

the volume-magnitude correlation. We cannot determine whether the magnitudes 

are likely to increase at the same rate as a function of injected volume as they did 

prior to 2000 or whether the rate of magnitude increase has been reduced since 

2000.

. 
75

Figure 3-21. Residuals of the least-squares fit to the maximum magnitude versus the log of 

cumulative injected volume, plotted against the quantiles of standard normal. The linearity 

of the residuals indicates that they are approximately normally distributed, confirming the 

validity of the fit shown in Figure 3-20.

A similar relation between maximum magnitude and total fluid volume injected 

can be seen in Figure 3-22, which summarizes the case histories presented in

Nicol and others (2011) and includes a wide range of data from injection sites, as 

well as extraction and hot dry rock sites. Based on the cumulative volume injected 

to date at Paradox Valley (7.5×106 m3), the trend line for data incorporating a 

variety of sites would suggest a mean maximum magnitude of MW 4.8 for induced 

earthquakes, with fairly broad 5% and 95% confidence limits of MW 3.4 and 6.2, 

respectively. This value of maximum magnitude is significantly greater than that 

predicted from the Paradox data alone, and the data point for Paradox is also 

considerably below the trend line. 
76

Figure 3-22. Maximum magnitude of induced earthquakes from 39 different sites, from

Nicol and others (2011). Data point for Paradox Valley is indicated (although the symbol was 

plotted for the ML 4.3 earthquake of May 27, 2000 earthquake, at the scale of the figure it 

applies equally to the January, 2013 M 4.0 earthquake (injected volume = 7.5×106 m3

).

There is much less scatter in the Paradox data (Figure 3-20) than in the fit to the 

Nicol and others (2011) data from many sites (Figure 3-22), however there are 

also substantially fewer data points (14 vs. 39). The trend for Paradox data

indicates a significantly smaller maximum magnitude (estimated ML 4.2; observed 

ML 4.4) than does the multi-site data.

3.3.3 Effect of Pore Pressure on Earthquake Magnitude

Pore pressures may also affect the magnitudes of induced seismic events. 

Empirical and theoretical studies suggest that an absolute pore pressure threshold 

may be needed to trigger large-magnitude events, at least in some cases. Studies 

also suggest that pore pressure gradients may play a role in the occurrence of 

large-magnitude induced earthquakes, especially those that occur after injection 

has ceased.

Empirical PVU data show a general correlation between higher long-term average 

injection pressures and the occurrence of larger-magnitude earthquakes in the 

near-well region (within 5 km of the injection well) (Block and Wood, 2009; 

Block and Wood, 2010). This correlation is shown in the upper two plots of 

Figure 2-22. The uppermost plot shows injection pressures averaged over 

different lengths of time, while the middle plot shows the near-well seismicity 

plotted as a function of date (x axis), earthquake magnitude (y axis), and rate (size 

Paradox Valley
77

of bubbles). Since we have no means for measuring in-situ pore pressures 

directly, the long-term averaging of the injection pressures is a simplistic method 

of accounting for the time delay and amplitude modulation of pressure variations 

at the injection well as the pressures propagate away from the well. To correlate 

pressure trends with the seismicity data, the pressure averaging must be done over 

longer time windows for later periods. The seismicity may take progressively 

longer to respond to changes in injection pressures both because the seismicity 

occurs at greater distances from the well for later time periods and because the 

size and complexity of the reservoir increases over time. While acknowledging 

that this method is very simplistic, the general correlation between relatively high 

long-term average injection pressures and increased rates and magnitudes of 

induced earthquakes in the near-well region suggests that pore pressure affects 

both the rate and magnitudes of induced seismic events. The lack of near-well 

induced earthquakes with magnitude greater than 2.5 during periods of time when 

the long-term average injection pressures are low suggests that a pore pressure 

threshold must be exceeded to trigger earthquakes of magnitude 2.5 or greater. 

Recently, Keranen and others (2013) suggested a similar link between the 

occurrence of large earthquakes in Oklahoma they interpret as injection-induced 

events and a rise in injection pressure (with a 5-year time delay). The empirical 

correlation between injection pressures and seismicity for PVU is only observed 

for earthquakes occurring within 5 km of the PVU injection well; the seismicity 

occurring at greater distances does not show the same trends (Figure 2-22, lower 

plot).

Numerical simulations also indicate a correlation between pore pressure and 

earthquake magnitude. Baisch and others (2010) performed finite element 

modeling of pressures and induced seismicity for an injection simulation. In their 

modeling, the largest-magnitude induced earthquakes occurred after the end of 

injection. While pressures near the simulated wellbore were falling during this 

time period, pressures were rising along the outer rim of the stimulated area, 

where the largest-magnitude events occurred. The investigators attribute both the 

time of occurrence (post-injection) and the location of the largest-magnitude 

events to shallow spatial gradients of pressure. The post-shut-in pressure rise 

along the outer rim of the stimulated area has a shallower spatial gradient than the 

pressure rise that occurs during injection (Figure 3-23). Because of the shallow 

pressure gradient, many neighboring patches in their model are at a similar level 

of stress criticality simultaneously. Therefore, only a small amount of stress 

diffusion is sufficient to cause overcritical conditions over a large area, increasing 

the likelihood for large-magnitude events. As the authors note, the results of their 

modeling are consistent with the observation that for many cases of fluid 

injection, the largest-magnitude events occur during the shut-in period after 

injection has ceased. 
78

Figure 3-23 Hydraulic overpressures inside a simulated fault zone as a function of radial 

distance from the injection well at the time of shut-in (dashed line) and 6 days after shut-in, 

when the largest-magnitude simulated event occurred (solid line). The shaded area marks 

the patches that slipped during the largest-magnitude event. Figure taken from Baisch and 

others (2010).

Further insights into the relationship of pore pressure and the occurrence of largemagnitude induced earthquakes come from the work of Garagash and 

Germanovich (2012). These investigators numerically modeled a pre-stressed 

fault with a locally peaked, diffusively spreading pore pressure field. They 

studied the nucleation of dynamic rupture, the rupture run-out distance, and 

rupture arrest. Their modeling indicates that the largest run-out distances of 

dynamic slip – corresponding to the largest earthquakes – occur when the pore 

pressure increases just enough to activate slip on faults that have low shear prestress. These are the faults that are not optimally oriented in the background stress 

field for shear failure. These faults require a larger rise in pore pressure, a larger 

pressurized region, and a longer nucleation time for slip to occur than the more 

favorably-oriented faults. Although counterintuitive, the results of this modeling 

also indicate that as pressure continues to increase to higher values beyond the 

minimum required for slip activation, the likelihood of inducing large-magnitude 

earthquakes actually diminishes. This finding is related to mechanisms that act to 

arrest the nucleation of dynamic rupture under these conditions. In this scenario, 

then, the occurrence of large-magnitude earthquakes could cease at very high pore 

pressures and then later resume if injection operations are altered and pore 

pressures decline. 

3.3.4 Maximum magnitude discussion

We have examined only a few of the models presented in the literature. Other 

models of maximum magnitude include the hypothesis that larger events occur in 
79

places with fully developed, thick faults with cataclasite fault cores. McClure and 

Horne (2012) developed a degree-of-fault-development factor based on categories 

of fault development observed on wellbore cores, surface outcrops, and imaging 

logs. We have not attempted to implement their classification system for Paradox 

Valley.

3.4 Short-term Earthquake Recurrence Estimates

To estimate seismic hazards arising from PVU induced seismicity, the rate of 

occurrence of earthquakes with magnitude large enough to potentially cause 

damage must be determined. Such damage, in turn, depends on earthquake 

magnitude, distance between the earthquake and the site of interest, faulting type 

and orientation, site conditions, the fragility of structures throughout the area, and 

other factors. Insufficient data exists to reliably estimate all (or even most) of 

these factors. Over the life of the project, both the rates and spatial distribution of

induced seismicity have varied considerably, even under conditions of relatively 

constant average injection rates. Currently, induced seismicity rates do not appear 

to quickly respond to changes in average injection rate, and we know of no

models capable of reliably estimating long-term induced seismicity patterns from 

historical well-head pressure or injection rate data.

To make a rough estimate of the likelihood of potentially-damaging earthquakes 

occurring over the next several years, we therefore make the assumption that both 

injection operations and short-term seismicity rates will not vary greatly from 

those of the recent past. We select a subset of the observed induced seismicity 

data, from 2009-2012, to characterize current short-term seismicity rates for both 

the near-well and NW cluster areas. Gutenberg-Richter a- and b-values are then 

calculated from this data using the maximum likelihood method, as described in 

Section 3.2.2. These results are used to compute a short-term cumulative 

seismicity rate λa (m M≥ min ) assuming no change in average injection rate:

( ) ( ) min max

10 min 10 log log 10 10 bM bM λa mM a − − ≥ =+ − (3.15)

If instead injection rates may change, then an alternative approach is to use the

instantaneous seismogenic index, Σ ∆ ′(t t ; ) , defined in Equation (3.12), and 

computed in Section 3.2.3 for previous time intervals. If the seismogenic index is 

assumed to remain relatively constant in the short term, and also assumed to be 

independent of the injection rate, then the short-term seismicity rate under 

conditions of a change in injection rate can be computed as follows:

( ) ( )

( ) ( )

10

10

log , ; ,

; log ( ; ) , ;

a min a min

min min

mM tt NM t t

t t qt t b M t t M

λ ≥ ∆ =∆ ∆

=Σ ∆ + ∆ − ∆ ⋅ ′ ′

(3.16)
80

where qt t Qt t t (; ) ; ∆ =∆ ∆ ∆ ( ′ ′ ) is the average injection rate assumed for the next 

time interval ∆t′ . However, if the long-term average rate of injection does not 

change significantly over the next several years of injection, then these results 

would not differ from those estimated using Equation (3.15).

Once the cumulative seismicity rate has been estimated, using 

either Equation (3.15) or Equation (3.16), the cumulative rate of earthquakes with 

magnitude m M≥ can be obtained using Equation (3.2) as follows:

max

max min

min

( )

( ) min max

max

1 10 ( ) 1 10

0

a

bM M

a bM M

for M M

m M for M M M

for M M

λ

λ λ

− −

− −

 ≤



 − ≥ = ⋅  << − 

 ≥ 

(3.17)

where the symbols are the same as in Equation (3.2).

Figure 3-24. Magnitude recurrence data for 2009 through 2012 for the near-well area, with 

calculated Gutenberg-Richter relation assuming Mmin = 1.0 and Mmax = 4.0 through 6.0, in 0.5 

magnitude increments.

λa (m M≥ min )

Near-well Seismicity, 2009-2012

0.01

0.1

1

10

100

0 1 2 3 4 5

Magnitude

Cumulative Event Rate (#/year)

Data

Mmax 6.0

Mmax 5.5

Mmax 5.0

Mmax 4.5

Mmax 4.0
81

Figure 3-24 shows the magnitude-recurrence data from the near-well area for the 

four-year period 2009-2012, along with Gutenberg-Richter fits determined using 

the maximum-likelihood method and Equation (3.17). Similar data and fits for the 

NW cluster are shown in Figure 3-25. For both areas, a value of Mmin = 1.0 was 

selected, and recurrence estimates were made for Mmax values ranging between 

4.0 and 6.0. Since M ≥ 4.0 earthquakes have already been observed at PVU, a 

choice of Mmax = 4.0 as an upper magnitude limit is inconsistent with the total 

data set, but is shown in these figures because it visually appears to provide a 

better fit to the recent data.

Figure 3-25. Magnitude recurrence data for 2009 through 2012 for the NW cluster area, 

with calculated Gutenberg-Richter relation assuming Mmin = 1.0 and Mmax = 4.0 through 6.0, 

in 0.5 magnitude increments.

The Gutenberg-Richter fits to the recurrence data overestimate the observed 

events rates for magnitudes greater than about M 2.5 for any choice of

min M ≥ 4.0, as illustrated in Figure 3-24 (near-well area) and Figure 3-25 (NW 

cluster). The discrepancy is greatest for the near-well area. This result is 

consistent with the poor fit of Gutenberg-Richter models to PVU induced event 

data seen for most time periods, as discussed in Section 3.2.2. This misfit may be 

due to factors such as potential saturation of the duration magnitude scale used for 

NW Cluster Seismicity, 2009-2012

0.01

0.1

1

10

100

0 1 2 3 4 5

Magnitude

Cumulative Event Rate (#/year)

Data

Mmax 6.0

Mmax 5.5

Mmax 5.0

Mmax 4.5

Mmax 4.0
82

PVU earthquakes, as well as deviation of the PVU data from Gutenberg-Richter 

statistics. Until the detailed nature of the misfit is determined, the GutenbergRichter fits to the data can be treated as a conservative estimate of short-term 

seismicity rates.

If induced seismicity at PVU is assumed to follow a Poisson process, then the 

short-term earthquake rates determined from the Gutenberg-Richter fits to the 

data can be used to estimate a probability of occurrence for earthquakes with a 

magnitude larger than have been observed in the data. For earthquakes following 

a Poisson process with a time-independent seismicity rate λ (m M≥ ) , the 

probability of occurrence of at least one earthquake with magnitude in the 

next ∆t years is given by:

( ) ( ) 1; 1 mM t Pn t e− ≥∆ λ ≥ ∆ =− (3.18)

We select a time period of 3 years for these estimates, since we have little 

confidence that reasonable results can be obtained for longer time periods, given 

the observed variation in seismicity rates over the lifetime of the project.

The seismicity rates determined from Equation (3.17) are used in Equation (3.18)

to estimate the probabilities of one or more earthquakes occurring within the next 

3 years, and having a magnitude greater than 3.0, 3.5, 4.0, 4.5, 5.0, and 5.5. 

Short-term earthquake occurrence probabilities are computed independently for 

both the near-well area (Table 3-1) and the NW cluster (Table 3-2). The estimated 

short-term earthquake probability results indicate a 20 to 40 % probability of 

another M ≥ 4 earthquake occurring within the next 3 years, and a 5 to 10 % 

probability of an M ≥ 5 earthquake in the case that Mmax ≥ 5.0.

Table 3-1. Estimated probability of occurrence of earthquakes over the next 3 years for the 

near-well area.

M

Probability of earthquake m M≥ within the next 3 years

max M = 4.0 max M = 4.5 max M = 5.0 max M = 5.5 max M = 6.0

3.0 0.91 0.92 0.92 0.93 0.93

3.5 0.53 0.63 0.66 0.68 0.68

4.0 0.26 0.34 0.38 0.39

4.5 0.12 0.17 0.19

5.0 0.05 0.08

5.5 0.02

m M≥

∆t
83

Table 3-2. Estimated probability of occurrence of earthquakes over the next 3 years for the 

NW cluster area.

M

Probability of earthquake within the next 3 years

3.0 0.63 0.65 0.66 0.65 0.66

3.5 0.29 0.37 0.40 0.40 0.41

4.0 0.14 0.20 0.21 0.23

4.5 0.07 0.10 0.11

5.0 0.03 0.05

5.5 0.02

Formally, the assumptions of a Poisson process are not met for this data set since 

no attempt was made to decluster the earthquake catalog to remove aftershocks or 

otherwise account for dependent events. The estimated rates and probabilities 

therefore are higher than would be obtained otherwise. In addition, the 

Gutenberg-Richter model used to estimate the rates of larger-magnitude 

earthquakes overestimates the observed rates of seismicity. Finally, the 

probabilities were calculated for the NW cluster and the near-well area, assuming 

that seismicity in the two areas is independent, however since seismicity in both 

areas is ultimately the product of the same injection sequence, the assumption of 

independence may not be warranted. The calculated probabilities therefore are 

order-of-magnitude, and possibly conservative estimates of short-term earthquake 

recurrence.

3.5 Discussion

A thorough analysis of seismic risk posed by induced seismicity at PVU would 

consider the rates at which damaging ground motions are likely to occur in the 

populated surrounding areas, and the likelihood of damage, economic loss, and 

potential life-loss given those loadings. Such an analysis is beyond the scope of 

this report, and many of the required input parameters are not currently known. 

However, some of the necessary inputs required for such an analysis have been 

discussed in this section.

A probabilistic seismic hazard analysis (PSHA) estimates the average rate 

λ (a A > ) of exceeding a specified ground motion A given an average 

earthquake rate ν (EQ) and a conditional probability of exceeding ground 

motion given the earthquake (Cornell, 1968):

( ) ( ) ( )

, ( ) ( ) ,

Pr |

, Pr | , m r

M R

a A EQ a A EQ

m r a A m r dmdr

λ ν >= >

=ν > ∫∫ (3.19)

m M≥

max M = 4.0 max M = 4.5 max M = 5.0 max M = 5.5 max M = 6.0

A
84

In order to make the calculations which are implicit in Equation (3.19) we need to 

know the average seismicity rates, maximum magnitudes, and the magnitude 

distributions for the areas of induced seismicity. As shown earlier in this section, 

these parameters have changed substantially at PVU over the lifetime of the 

project, making estimates for future behavior highly uncertain. An even bigger 

problem is that the precise relation between the seismicity rate and injection 

parameters that can be controlled (e.g., average flow rate, maximum injection 

pressure, shut-in duration, etc.) is poorly known. Additional required parameters 

include the attenuation of ground motions with distance, site-response properties 

(e.g., VS30 and basin depth), and the density and fragility of residences and other 

structures.

Recent studies have attempted to estimate the probability of earthquake 

occurrence (Bachmann and others, 2011; Barth and others, 2013) and its relation 

to injection parameters in an attempt to devise methods for controlling seismic 

hazards from enhanced geothermal systems (EGS). Other studies of EGS 

injection have estimated seismic hazard for induced seismicity using a standard 

PSHA approach (Convertito and others, 2012; Mena and others, 2013). The longterm behavior of EGS injection likely differs substantially from injection-only 

operations such as PVU, since fluids are withdrawn in EGS systems, and since 

there is a much closer correlation between induced seismicity and recent injection.

Long-term seismic hazard estimates will be difficult to make at PVU without 

developing models that are capable of relating induced seismicity to injection 

operations. Such models have not yet been developed for PVU, although it is 

possible that recent numerical modeling approaches could be used to estimate 

both the pressure response and the seismicity response to injection (e.g., Wang, 

2008; Baisch and others, 2010; Lee and Ghassemi, 2010; Lee and Ghassemi, 

2011; Wang and Ghassemi, 2012). Short-term seismic hazard estimates can be 

made by extrapolating recent seismicity parameters, however even these estimates 

are likely to include substantial uncertainty.
85

4 Injection Pressures and Flow Rates

In order to predict the PVU reservoir response to future injection, it is necessary 

to understand the reservoir characteristics and their evolution over time. We have 

found that industry-standard well-test analysis methods based on idealized models 

of the reservoir provide a relatively simple way to determine these characteristics, 

and can reasonably accurately predict the observed pressure response at the PVU 

wellhead; however, both the determination of the parameters and their 

interpretation in terms of the actual behavior of the reservoir are subject to 

significant uncertainty. This section describes the history of PVU injection, 

examines previous modeling work, reviews the theoretical models used, shows 

the results of the modeling, and uses the models obtained to predict well-head 

pressures under various changes in injection operations.

4.1 Injection History

Between 1991 and 1995, PVU conducted a series of 7 injection tests, an acid 

stimulation test, and a reservoir integrity test. The purpose of these tests was to 

qualify for a permit for long-term injection from EPA. Continuous injection of 

brine began in July 1996, after EPA granted the permit. Since continuous 

injection began, PVU has instituted and maintained three major changes in

injection operations. Each change was invoked to mitigate the potential for 

unacceptable seismicity or to improve injection economics. These injection 

phases are described below. Plots of the daily average injection flow rate, surface 

injection pressure, and downhole pressure at a depth of 14,100 ft (4.3 km) 

throughout the history of PVU injection operations are shown in Figure 4-1.

4.1.1 Phase I (July 22, 1996 – July 25, 1999)

During this initial phase of continuous injection, PVU injected at a nominal flow 

rate of 345 gpm (~1306 l/min), at about 4,950 psi (~34.1 MPa) average surface 

pressure. This corresponds to approximately 11,800 psi (~81.4 MPa) downhole 

pressure at 14,100 ft (4.3 km) depth. To maintain this flow rate, 3 constant-rate 

pumps were used with each operating at 115 gpm. The surface pressure on 

occasion approached the wellhead pressure safety limit of 5,000 psi. At these 

times PVU would shut down 1 or 2 injection pumps, reducing the injection rate 

and allowing the pressure to drop a few hundred psi, before returning to 3-pump 

injection. These shutdowns occurred frequently and lasted for minutes, hours, or a 

few days. Maintenance shutdowns lasted for 1 to 2 weeks and, in mid-1997, a 71-

day shutdown was needed when replacing the operations contractor. The 

shutdowns resulted in an overall average injection rate for phase I of ~300 gpm 
86

(1136 l/min). The injectate during phase I was 70% Paradox Valley Brine (PVB) 

and 30% fresh water.

Due to concerns about corrosion, the decision was made in 1997 to switch from 

using an O2 scavenger to a rust inhibitor, as the O2 scavenger was found to be 

insufficient as a corrosion inhibitor. While the rust inhibitor was substantially 

more effective at inhibiting corrosion, it also increased the risk of precipitation of 

elemental sulfur.

Figure 4-1. Daily average injection flow rate (top), daily average surface injection pressure 

(middle), and daily average downhole pressure at 14,100 ft (4.3 km) depth (bottom) during 

PVU injection operations.
87

4.1.2 Phase II (July 26, 1999 – June 22, 2000)

Following two ML 3.5 induced earthquakes in June and July, 1999, PVU 

augmented injection to include a 20-day shutdown (i.e., a “shut-in”) every 6 

months. Prior to these events, it was noted that the rate of seismicity in the nearwellbore region (i.e., within about a 2-km radius from the wellbore) decreased 

during and following unscheduled maintenance shutdowns and during the 

shutdowns following the injection tests of 1991 through 1995. It was thought that 

the biannual shutdowns might reduce the potential for inducing felt seismicity by 

allowing extra time for the injectate to diffuse from the pressurized fractures and 

faults into the formation rock matrix. When injecting during this phase, the 

injection pressure and flow rate were the same as during Phase I.

4.1.3 Phase III (June 23, 2000 – January 6, 2002)

Immediately following a ML 4.3 earthquake on May 27, 2000, PVU shut down for 

28 days. During this shutdown period, PVU evaluated the existing injection 

strategy and its relationship to induced seismicity. PVU decided to reduce the 

injection flow rate in order to reduce the potential for inducing felt seismicity. On 

June 23, 2000, PVU resumed injection using two pumps rather than alternating 

between 2 and 3 pumps. The biannual 20-day shutdowns were maintained. The 

nominal flow rate during Phase III, while injecting using two pumps, was 230 

gpm (~871 l/min). Accounting for the two 20-day shut-ins per year, the average 

injection flow rate was approximately 205 gpm (776 l/min), a decrease of about 

32% compared to Phase I. The 70:30 ratio of brine to fresh water was maintained.

4.1.4 Phase IV (January 7, 2002 – present)

Beginning with continuous injection operations in 1996, the injectate had been 

diluted to 70% PVB and 30% Dolores River fresh water. A geochemical study 

had predicted that if 100% PVB were injected, it would interact with connate 

fluids and the dolomitized Leadville limestone at downhole (initial) temperatures 

and pressures, and that PVB would then precipitate calcium sulfate, which in turn 

would lead to restricted permeability (Kharaka and others, 1997). During October 

2001, with the decreased injection volume discussed above, the injectate 

concentration concerns were reconsidered. Temperature logging in the injection 

interval documented substantial near-wellbore cooling, indicating that a 

significant displacement of connate brine away from the wellbore had occurred by 

this time, reducing the potential for 100% PVB and connate brine to mix, and 

therefore that if precipitation occurred, it would not be near the wellbore 

perforations where clogging would be a concern. Further discussions indicated 

that, if precipitation occurred, its maximum expected rate would be ~8 tons of 

calcium sulfate per day. To put this amount into perspective, injecting at ~230 

gpm, assuming a density of 9.86 lbs/gal (17% more dense than fresh water), 

results in a daily injection of ~1633 tons. The maximum expected precipitate is 

~0.5% of the daily injection mass.
88

After considering this new information, Reclamation decided to begin injecting 

100% PVB, in order to increase the amount of salt disposed with the reduced 

injection rate initialized in Phase III. Injecting 100% PVB began on January 7, 

2002, following the December-January 20-day shutdown, and has been 

maintained since. The same reduced injection rate as in Phase III (230 gpm) and 

biannual 20-day shutdowns have been maintained. The only noticeable effect of 

the change to 100% PVB injectate has been increasing bottom hole pressure 

because of the increased density of 100% PVB (by about 5%) over the 70% PVB 

to 30% fresh water mix.

4.2 Previous Work

One of the key parameters that defines reservoir behavior is the effective 

permeability. Previous estimates of permeability have varied widely. Drill stem 

tests using a variable flow rate gave an original permeability of 7.97 mD (Harr, 

1989). However, Horner analysis performed at the same time indicated a 

permeability between 1.3 and 1.5 mD (Harr, 1989). Additionally, analysis of core 

samples from the Conoco-Scorup well located approximately 4.6 km to the 

northeast yielded permeabilities ranging from 0.03 to 1.3 mD (Harr, 1989), which 

could suggest significant regional variations in permeability.

A series of 7 injection tests were performed between 1991 and 1995. 

Permeabilities were only calculated for the first 3 tests, which took place between 

July 11, 1991 and June 18, 1992, primarily due to equipment failures in later tests. 

The permeabilities for the first injection and falloff periods were calculated as 4-5 

mD and 2 mD, respectively. The second injection and falloff periods yielded 7.8 

mD and 2.2 mD, respectively. The third falloff period yielded a permeability of 

1.6 mD, while the third injection period was determined to be unanalyzable due to 

numerous flow rate changes caused by equipment failure. The lower 

permeabilities calculated during falloff periods, as compared with injection 

periods, were assumed to indicate fracture closure (Envirocorp Services and 

Technology Inc., 1995).

An early model of the injection test data attempted to fit the data from the seventh 

and final injection test, which took place from August 14, 1994 to April 3, 1995, 

using a fracture network model with 5% porosity and a permeability of 3 mD 

(Envirocorp Services and Technology Inc., 1995). While they obtained reasonable 

fits to the data from this injection test, their model predicted future pressures that 

were significantly lower than observed pressures. While the modeled injection 

scenarios differ from the actual injection scenarios, preventing a direct 

comparison, their model predicted that surface pressure would not exceed 4500 

psi after 10 years of continuous injection at a flow rate of 300 gpm, whereas the 

observed pressures exceeded 4500 psi after less than 4 months of injection at an 

average flow rate of 163 gpm. This suggests that the effective permeabilitythickness and/or porosity-thickness may have been overestimated in this early 

modeling.
89

In 2001, Reclamation funded USGS to perform fluid-flow modeling of the 1996-

2000 pressure and flow data, and a draft report was received in 2009 (Roeloffs 

and Denlinger, 2009). The USGS model consisted of a radially-symmetric finite 

element grid with 20 layers based on stratigraphy, each characterized by four 

parameters (shear modulus, Poisson ratio, Skempton’s coefficient, and hydraulic 

diffusivity). The model was calibrated by trial-and-error matching of the predicted 

flows and pressures with those measured at the wellhead. For the layer 

representing the Leadville formation, a porosity of 10% and a permeability of 28 

mD were assumed for a thickness of 146 meters (479 feet), with significantly 

lower values for other layers. In contrast, Bremkamp and Harr (1988) determined

a porosity in the Leadville of ≥10% for only 2 feet of its 416-foot thickness, using 

sonic logging, or 18 feet, using visual binocular examination of recovered core.

Hydraulic diffusivity (calculated from permeability, formation thicknesses, 

porosity, and shear modulus), was found to be the parameter most strongly 

effecting pressure distribution in the model. The USGS model apparently did not 

incorporate skin factor or wellbore storage parameters.

The USGS model provides a reasonably good fit to the flow rate data for 

approximately the first 500 days, but then the modeled flow rates begin to 

decrease, while the measured flow rates remain approximately constant. Roeloffs 

and Denlinger (2009) suggest that formation permeability may have increased, or 

that the flow may have reached a zone of higher permeability away from the well.

The permeability of intact limestone and dolomite generally varies from 0.01 to 

0.1 mD (Bear, 1972). The permeability value for intact rock is known as primary, 

or matrix, permeability, and is generally only applicable to laboratory samples, as 

rocks over larger scales will contain at least some degree of fracturing. 

Limestones may also have increased permeability due to vugular porosity. The 

permeability due to fracturing and secondary porosity is known as secondary 

permeability and can vary by several orders of magnitude, even exceeding 105 mD 

in some highly fractured locations (Cox and others, 2000). Thus, the estimated 

permeabilities of a few mD, while higher than would be expected to be measured 

under laboratory conditions, are reasonable for a moderately fractured reservoir.

4.3 Fitting of Pressure/Flow Data to Idealized Models

It was unclear initially if complicated numerical models were required to fit the 

pressure-flow data at the wellhead, or instead if simple idealized models could 

suffice. In the absence of monitoring wells or other independent means of 

verifying predicted pressures within the injection reservoir at distances away from 

the wellbore, we were concerned that use of excessively complex models might

result in overfitting of the wellhead data, and yield a correspondingly poor 

predictive capability of the models in evaluating potential injection scenarios. We

therefore tried fitting the wellhead data for individual injection and shut-in cycles 

using several simple, idealized models commonly employed in the well-testing 
90

and reservoir-modeling industry. F.A.S.T. WellTest software v.7.6.0, described in 

Fekete (2012), was used for this analysis.

Biot (1941) developed a model of a fluid-filled porous material based on a 

conceptual model of a coherent solid skeleton and freely moving pore fluid

(Detournay and Cheng, 1993). Under this model, for the case of an isotropic 

applied stress tensor σ , the constitutive equations are

11 12

21 22

a ap

a ap

ε σ

ς σ

= +

= +

(4.1)

where ς is the increment of fluid content, p is the fluid pressure, and ε is the 

volumetric strain (Wang, 2000). This model includes full coupling between the 

solid and fluid components. Substituting physically meaningful constants, the 

most general form of the linear constitutive equations for isotropic material 

response is

11 1

2 69 3

3

ij

ij ij kk ij

kk

p G GK H

p

H R

σ

ε δσ δ

σ

ς

  =− − +    

= +

(4.2)

where K and G are the bulk modulus and shear modulus of the drained elastic 

solid, and H , and R are constitutive constants characterizing the coupling 

between the solid fluid stress and strain (Detournay and Cheng, 1993).

At reservoir depths, the rock compressibility is generally small, and thus a 

simplifying assumption may be made in order to uncouple the solid and fluid 

components. This is a common simplification in reservoir modeling applications, 

and is used in our analysis.

For the uncoupled system, and assuming incompressible rock, the governing 

equation for radial flow in cylindrical coordinates is

2

2

1 t pp p c

r rr k t

∂∂ ∂ φ μ + =

∂∂ ∂ (4.3)

where r is the radial distance from the well in centimeters, t is time in seconds, 

p is the reservoir pressure in atmospheres at distance r and time t , φ is the 

formation porosity as a fraction of bulk volume, k is the formation permeability 

in Darcies, μ is the fluid viscosity in centipoises, and t c is the total

compressibility in volumes per volume per atmosphere (Horner, 1951).

For a point source, the change in pressure resulting from an applied flow is
91

2

0 Ei

4 4

t q r c

p p kh kt

μ φμ

π

  − =    

(4.4)

where 0 p is the initial reservoir pressure in atmospheres at t = 0 , q is a constant 

rate of production of the well (starting at t = 0 ) in cubic centimeters of subsurface 

volume per second, h is the layer thickness in centimeters, Ei ; 0 ( x x ) > is the 

exponential integral (Abramowitz and Stegun, 1964), and other variables are the 

same as in Equation (4.3). This solution assumes an infinite reservoir and an 

infinitely small wellbore radius (Horner, 1951).

For sufficiently small x > 0 , Ei( x)may be accurately approximated by 

Ei log ( x x ) ≈ + γ ( ), where γ is Euler’s constant (Abramowitz and Stegun, 1964). 

Thus if the parameters of the model are homogeneous and constant over time, 

then Equation (4.4) may be approximated as

0 2

4 log . 4 t

q kt p p const

kh r c

μ

π φμ

  − =−   +

 

(4.5)

which shows that for continuous injection q < 0 , the pressure change at the 

wellbore radius wr is asymptotic to log(t) , and with a slope proportional to the 

flow rate q and inversely proportional to the permeability-thickness product kh .

Typical well-testing procedures consider a period of constant production starting 

at time t = 0 , followed by a shut-in period starting at time 0t > 0, then the well

pressure during the shut-in period at time 0t +τ , can be obtained by 

superimposing two equations of the form of Equation (4.4). Using the wellbore 

radius wr for r, this leads to the equation

( )

2 2

0

0

Ei Ei

44 4

wt wt

w

q rc rc

p p kh k t k

μ φμ φμ

π ττ

      − =     −         +  

(4.6)

where w p is the pressure at the wellbore in atmospheres and wr is in centimeters. 

Using the previous approximation forEi( x) at small x > 0 , the pressure change 

during the shut-in period described by Equation (4.6) may be approximated as

0

0 log 4 w

q t

p p kh

μ τ

π τ

  + − =−     (4.7)
92

For continuous injection q < 0 , Equation (4.7) shows that the pressure fall-off at 

the wellbore radius following shut-in is asymptotic to 0 log t τ

τ

  +     , where 0t is the 

injection duration, τ is the time since shut-in, and where the slope is proportional 

to the flow rate and inversely proportional to the permeability-thickness 

product kh . Horner (1951) demonstrates that the error introduced by this 

approximation typically drops to 0.25% within seconds of closing the well.

In dimensionless variables, the governing equation equivalent to Equation (4.3) is

2

2

1 D DD

D DD D

p pp

r rr t

∂ ∂∂

+ =

∂ ∂∂ (4.8)

where, in arbitrary units, D w r rr = is the dimensionless radial distance from the 

well, D 2

w t

kt

t

r c φμ = is dimensionless time, and D ( 0 ) kh

p pp

qμ = − is the dimensionless 

pressure change on the formation side at distance Dr and time Dt (Agarwal and 

others, 1970). Using the previously defined units, the equations for Dr and Dt are 

unchanged, while the equation for D p becomes D 0.0689 ( 0 ) kh

p p p

qμ = − . All 

model properties are assumed to remain constant with time, and pressure is 

calculated based on the observed injection rate, which varies over time.

In solving Equation (4.8) for a finite wellbore, two additional variables, the 

wellbore storage constant and skin factor, are introduced in the inner boundary 

conditions, in order to ensure continuity across the wellbore boundary. In 

dimensionless variables, the boundary conditions are

1

1

D

wD D

D

D D r

dp p C

dt r =

  ∂ − =   ∂  

(4.9)

and

1 D

D

wD D

D r

p p ps

r =

    ∂ = −     ∂    

(4.10)

where wD p is the dimensionless pressure drop within the wellbore, CD is the 

dimensionless wellbore storage constant, s is the dimensionless skin factor, and 

other parameters are the same as in Equation (4.8) (Agarwal and others, 1970). 

All fluid and rock properties are assumed to remain constant.

q
93

This analysis used injection flow rates and surface pressures recorded at the 

wellhead since continuous injection began in 1996. The sample rate of that data is 

non-uniform over time. Injection data prior to 2003 were recorded on an older 

SCADA system in a format that is currently inaccessible. For this period, only 

daily average pressures and flow rates are available. Starting in 2003, a newer 

SCADA system was installed, and that data is available at much higher sample 

rates (up to 1 sample every 2 seconds). We found that mixing data at the two 

sample rates created artificial trends in the data. For example, a large apparent 

decrease in the wellbore storage constant occurred when mixing the daily-average 

with the high-frequency data. While misfits were lower using data with the higher 

sample rate, and the absolute values obtained for parameters were likely more 

accurate, we are more interested in trends in the parameters rather than in their 

particular values. We therefore decided to use the daily average data over the 

entire well operating history for this analysis.

Reservoir temperature is a significant factor affecting the modeling, primarily 

because of its effect on the fluid viscosity, as the temperature of the injectate is 

generally assumed to be equal to the undisturbed reservoir temperature. While this 

assumption is probably adequate for short-term injection testing, the volume of 

fluid that has been injected into the PVU injection well has likely led to cooling a 

significant distance from the well, and thus the reservoir temperature varies in 

both time and space. However, as our model cannot account for temperature 

variations, and as we have no data to constrain temperatures away from the well, a 

single temperature was used for the modeling. A temperature of 37.8 ̊C was used, 

corresponding to the temperature measured 3 days into a shut-in that occurred in 

March 1994 (Subsurface Technology, 2001).

Other input parameters include a reservoir thickness of 100 feet (30.5 meters), 

porosity of 3%, and salinity of 2.6 kg/L after January 8, 2002, when injection of 

100% brine began, and 1.82 kg/L prior to this date, when the injectate was a 

mixture of 70% brine and 30% fresh water. The salinity of the fresh water was 

assumed to be negligible relative to the brine, and therefore a salinity 70% that of 

the brine was used for the mixture. The assumed thickness of 100 feet is 

approximately equal to the thickness of the perforated interval of the upper 

Leadville, as logging has indicated that the middle and lower Leadville 

perforations are covered with fill (Subsurface Technology, 2001), and it is 

unknown to what extent fluid is able to migrate vertically outside of the 

perforated interval. Note that this is significantly smaller than the value of 479 

feet used by Roeloffs and Denlinger (Roeloffs and Denlinger, 2009), which is 

approximately equal to the entire thickness of the Leadville formation. 

Viscosity and fluid compressibility are calculated using the mean temperature and 

salinity. The calculated viscosity is 1.030 cP during injection of 70% brine and 

1.348 cP during injection of 100% brine. The calculated fluid compressibility is 

1.933×10-6 psi-1 during the injection of 70% brine and 1.731×10-6 psi-1 during the 

injection of 100% brine.
94

While the input parameters have a significant effect on the absolute values of the 

permeability, the effect on the relative values is minimal, and thus we do not 

expect our interpretations or the results of the forward modeling performed in 

Section 4.5 to be significantly affected by the choice of input parameters. For 

example, if we had used the value of 479 feet chosen by Roeloffs and Denlinger 

(Roeloffs and Denlinger, 2009) rather than 100 feet for the layer thickness, the 

calculated permeabilities would be reduced by a factor of 4.79, but the ratios of 

permeabilities between cycles would not be affected, nor would the pressures 

obtained in forward modeling. Variations in input parameters in space or time, 

such as a change in the thickness of the Leadville away from the well or reservoir 

cooling over time, could potentially have more significant effects, but our model 

cannot account for changes in these parameters, nor do we have data to constrain 

potential changes. 

Downhole pressures were calculated by adding a constant value to the surface 

pressures, accounting for the weight of the brine and the frictional effect due to 

flow, assuming new pipe conditions with no scale buildup inside the tubing

(Mahrer and others, 2004). The friction term is negligible for the injection rates 

and tubing size used at PVU. It is possible that the friction term becomes more 

significant over time due to scale build-up, but as we have no way of determining 

this, we are unable to incorporate the potential changes in friction into our 

analysis. A constant pressure difference between the wellhead and downhole 

pressure of 6822 psi was used for continuous injection Phases I-III, and 7133 psi 

for injection Phase IV, when injectate density increased as a result of the switch to 

100% brine.

After comparisons of several different models, a vertically stratified, radially 

symmetric model with no-flow boundaries in the vertical direction and infiniteacting boundaries in the horizontal directions was selected. Models with finite noflow or constant pressure boundaries in the x- and/or y-directions were 

considered, as well as radial composite models and models incorporating fractures 

and anisotropic permeability, but these models provided either a worse fit or an 

insignificant improvement in the fit to the data, and so the simplest model that 

adequately fit the data was selected. This model does not incorporate individual 

fractures, and so we assume that the near wellbore region is sufficiently fractured 

so that it can be adequately represented as a porous medium for the purpose of 

this simplified pressure-flow modeling.

There are 3 free parameters in our model: the effective permeability, k, skin due to 

damage, sd, and dimensionless wellbore storage constant, CD. While the entire 

flow rate history is incorporated into the calculation of modeled pressures, the 

parameters are fit individually for each build-up period. This approach is used to 

assess whether these parameters have changed over time, and because there was 

no clear reason for assuming that a single set of parameters would be adequate to 

fit all the data.
95

For Phases II-IV, when two 20-day shut-ins were scheduled every year, a cycle is 

defined to be the build-up period between two scheduled shut-ins. Unscheduled 

11-day and 17-day shut-ins occurred in the cycles beginning July 26, 1999 and 

January 21, 2006, respectively, so these flow periods were each split into two 

cycles for the purpose of analysis. All cycles contain unscheduled shut-ins, which 

range in length from a few hours to 10 days, primarily due to equipment issues. In 

order to analyze the 6-month cycles as a single flow period despite the short shutins, it is necessary to input a small but non-zero flow rate for the short shut-in 

periods. 

For Phase I, cycles are defined as the flow periods between two unscheduled 

multi-day shut-ins, which range in length from 4 to 70 days. Cycles range from 47 

to 345 flow days in length; build-up periods of only a few days are excluded from 

analysis.

Defining cycles in this matter leads to a total of 38 cycles: 7 in Phase I, 3 in Phase 

II, 4 in Phase III, and 24 in Phase IV. 

4.4 Results of Pressure/Flow Model Fitting

Model parameters were fit simultaneously for each injection cycle, using a 

Simplex automatic parameter estimation procedure to find a local minimum in the 

misfit between the data and the model. All parameters for each injection cycle 

were fit independently from the other injection cycles. Section 4.4.1 describes 

how the computed model parameters have evolved over the various phases of 

injection. Section 4.4.2 describes the quality of the model fits to the pressure-flow 

data.

4.4.1 Changes in Model Parameters over Time

In Phase I, values for k range from 9.06 to 10,200 mD, and values for sd range 

from -5.82 to 2000, with a strong positive correlation between these two 

parameters. Values for CD range from 6.27×105 to 1.46×106

. The magnitude of 

variations in k and sd is believed to be an artifact of attempts by the model to fit 

data that includes factors that the model cannot account for. In particular, the 

model seems to be unable to fit periods of frequent flow rate changes.

In Phase II, the variation in model parameters is much smaller than found in 

Phase I, with k ranging from 13.64 to 21.93, sd ranging from 6.05 to -3.99, and CD

ranging from 8.80×105 to 1.37×106

. 

In Phases III and IV, for 27 of the 28 cycles in these phases, effective 

permeabilities range from 10.39 to 29.21 mD, with a mean of 18.03 mD and 

median of 16.74 mD, and skin values range from -5.98 to -1.98, with a mean of -
96

4.35 and median of -4.47. In these phases, neither of these parameters shows any 

clear trends with time, as shown in Figure 4-2.

The remaining cycle was initially a significant outlier, with a permeability of 140

and skin value of 26.31. There are no obvious differences in the build-up curves 

between this cycle and the other cycles, so as in early phases, it is likely that these 

differences are not real, but are due to the method of model fitting. For this 

reason, the permeability and wellbore storage constants was re-fit after fixing the 

value of sd at the mean of the preceding and following cycles. This led to a 

permeability of 20.72, with a minimal increase in the misfit. 

Dimensionless wellbore storage constants for Phases II-IV vary from 9.01×105 to 

1.74×106

, which correspond to dimensional wellbore storage constants of 4.17 to 

8.03 ft3/psi. This is significantly larger than the number expected from 

considering only the actual wellbore volume, suggesting that the fractured area 

around the well may essentially be acting as an enlarged wellbore or an extension 

of the wellbore (Fekete, 2012). The wellbore storage constants are somewhat 

higher prior to 2003, with a decrease of about 21% between the mean values from

September 1999 to January 2002 (the beginning of Phase IV, when the change 

from injection of 70% brine to injection of 100% brine occurred) and the mean 

values beginning with June 2002 (Figure 4-2b). The magnitude of the decrease 

between the mean of the cycles in these two time periods is approximately 

2.94×105

, which is 2.6 times the standard deviation of the wellbore storage 

constants beginning in June 2002, indicating that the decrease is statistically 

significant. It is possible that the decrease in the wellbore storage constant 

indicates a decrease in the perforated interval of the wellbore due to narrowing of 

the wellbore and/or precipitation of elemental sulfur, both of which have been 

observed in past wireline logs (Subsurface Technology, 2001).

There is significant tradeoff between the parameters, leading to concern that the 

resulting parameters may not be robust, as evidenced by the two cycles with 

outlying values of both k and sd. The values of k and sd have a strong positive 

correlation (Figure 4-3). In order to minimize the tradeoff, the values for k and CD

are refit after fixing sd at a value of -4.35. This value for sd corresponds to the 

mean of the calculated values for the cycles in Phase III and Phase IV.
97

Figure 4-2. Effective permeability (a), dimensionless wellbore storage constant (b), and skin 

due to damage (c) versus time for pressure build-up cycles in Phases II-IV. Dashed lines 

indicate boundaries of the injection phases.
98

Figure 4-3. Scatter-plot matrix showing pairwise correlations between permeability (k), skin 

due to damage (sd), dimensionless wellbore storage constant (CD), pressure increase, and 

average flow rate for all pressure build-up cycles in Phases II-IV. Note the strong correlation 

between permeability and skin.

Fixing sd leads to much less variation in the values obtained for k, particularly 

those in Phase I and Phase II. The values for all cycles vary from 14.02 to 27.12, 

with a mean of 18.76 and median of 18.36. With the exception of an outlying 

cycle in 2002, the values are highest during Phases I and II, when the injection 

rate was higher, then lowest during Phase III, and have remained relatively 

constant since Phase IV began in 2002 (Figure 4-4).

The decrease in k between Phase II and Phase III may be an effect of fracture 

closure due to the decrease in flow rate and resultant decrease in reservoir 

pressure. Over time, the pressures have gradually increased to the pressures seen 

in Phase II, even surpassing these pressures in recent years. This increase in 

pressure could cause fractures to re-open and permeabilities to return to the values 

seen in Phase II. However, the permeabilities have remained lower than they were 

in Phase I, which could be an effect of precipitation in the near wellbore area, 

which could limit fluid flow and decrease the average effective permeability.

k (mD)

-6 -5 -4 -3 -2 2500 3000 3500 4000

10 15 20 25

-6 -5 -4 -3 -2

sd

Cd

10^6 1.4*10^6

2500 3000 3500 4000

Pressure Increase

(psi)

10 15 20 25 10^6 1.4*10^6 180 220 260 300

180 220 260 300

Average Flow Rate

(gpm)
99

Figure 4-4. Effective permeability (a) and dimensionless wellbore storage constant (b) versus 

time for all pressure build-up cycles during Phases II-IV, fixing the value of sd to -4.35. 

Dashed lines indicate boundaries of the injection phases.

The dimensionless wellbore storage constant varies from 4.70×105 to 1.89×106

, 

corresponding to a dimensional wellbore storage constant varying from 2.18 to 

8.73 ft3/psi. For Phases II- IV, the wellbore storage constant shows a similar trend 

to the values fit without fixing sd, with the highest values occurring prior to 2002, 

a 25% decrease from the mean values prior to June 2002 and the mean values 

beginning in June 2002, and relatively constant values since 2003. The decrease 

in the mean wellbore storage constant between January, 2002 and June, 2002 is 

3.1 times the standard deviation of the data beginning in June, 2002, again 

indicating that the decrease appears to be statistically significant.

With the skin value fixed, the permeability, wellbore storage constant, pressure 

increase during the cycle, and average flow rate all appear to be independent of 

one another (Figure 4-5). 
100

Figure 4-5. Scatter-plot matrix showing pairwise correlations between permeability (k), 

dimensionless wellbore storage constant (CD), pressure increase, and average flow rate for all

pressure build-up cycles in Phases II-IV, holding sd fixed at -4.35. 

The degree of tradeoff between the permeability and skin value demonstrates that 

we are unable to uniquely resolve the values for both parameters with the current 

dataset. It is possible, though not certain, that injecting at a fixed flow rate for a 

longer time period could allow the observation of boundary effects, which could 

help to resolve these parameters.

4.4.2 Model Pressure-History Plots and Type Curves

For each injection cycle, the free parameters of the model were fit simultaneously, 

using automatic parameter estimation to find a local minimum in the misfit

between the “observed” and modeled pressure at the bottom of the well. Bottomhole pressures (“observed”) were calculated from the measured surface pressures 

by accounting for the weight of the brine in the wellbore, and for frictional effects 

of the flow. Interruptions in injection during an injection cycle were handled 

using superposition to account for the variable rate history. The results of this 

modeling are a series of pressure-history plots and type curves, one for each 

injection cycle, which illustrate the modeled flow regime and how well it fits the 

data.

k (mD)

10^6 1.4*10^6 1.8*10^6 180 220 260 300

16 18 20 22 24

10^6 1.6*10^6

Cd

Pressure Increase

(psi)

2500 3000 3500 4000

16 18 20 22 24

180 220 260 300

2500 3000 3500 4000

Average Flow Rate

(gpm)
101

Pressure-history modeling results from a typical injection cycle in Phase I are 

shown in Figure 4-6. This figure includes results for the case when all 3 model 

parameters are allowed to vary (upper plot), and for the case when sd is held fixed 

(lower plot). These plots suggest that an idealized radial-flow model can provide a 

reasonable fit to this data, with errors of less than about 5 %. The model using 3 

free parameters fits the data significantly better than the model using two free 

parameters. The corresponding type curves (pressure change and derivative) are 

shown in Figure 4-7 and suggest that a radial-flow model fits the data reasonable 

well.

Figure 4-6. Recorded (blue dots) and modeled (red line) downhole pressures and percent 

error (dashed line) for an example Phase I cycle, beginning July 10, 1997, allowing all 3 

parameters to vary (top) and fixing sd at -4.35 (bottom). Recorded downhole pressures are 

calculated by adding a constant value of 6822 psi to measured surface pressures.
102

Figure 4-7. Recorded (blue dots) and modeled (red line) values for change in pressure (Δp) 

divided by flow rate (q), and the recorded (red triangles) and modeled (blue lines) values for 

the time derivative of Δp/q for an example Phase I cycle, beginning July 10, 1997, allowing all 

3 parameters to vary (top) and fixing sd at -4.35 (bottom).

Modeling results for typical injection cycles in the remaining phases are shown as 

follows: Phase II – pressure history (Figure 4-8) and type curve (Figure 4-9); 

Phase III – pressure history (Figure 4-10) and type curve (Figure 4-11); Phase IV 

– pressure history (Figure 4-12) and type curve (Figure 4-13). For each phase, the 

cycle shown is reasonably representative of all cycles in that phase. Unlike the 

Phase I cycle, model results from the later phases do not show a significantly 

greater misfit when sd is fixed rather than allowed to vary.

The fits in the later phases are generally better than in Phase I, particularly 

compared to the time prior to 1998 when there were more variations in flow rate, 
103

with large variations in the estimated parameters. In Phase I, the modeling 

produced large variations in the estimated parameters and errors in the predicted 

pressures at times exceed 5%. In Phase II, the fits are significantly improved with 

errors mostly within 2%. In Phases III and IV, the errors are generally less than 1 

to 2%. While the magnitude of the pressure falloff during the shut-ins is matched 

relatively well, the shapes of the falloff curves are not well-matched in the later 

phases, which suggests an effect that is not accounted for in the model, such as 

closing of fractures during shut-in periods.

Whether or not sd is fixed, the models generally fit the data within 1-2% error for 

the build-up cycles in Phases III and IV, with slightly higher errors in Phase II and 

significantly higher errors in Phase I. Errors tend to be lowest in periods of 

continuous flow, and highest at the beginning of cycles and around mid-cycle 

shut-ins.
104

Figure 4-8. Recorded (blue dots) and modeled (red line) downhole pressures and percent 

error (dashed line) for an example Phase II cycle, beginning July 26, 1999, allowing all 3 

parameters to vary (top) and fixing sd at -4.35 (bottom). Recorded downhole pressures are 

calculated by adding a constant value of 6822 psi to measured surface pressures.
105

Figure 4-9. Recorded (blue dots) and modeled (red line) values for change in pressure (Δp) 

divided by flow rate (q), and the recorded (red triangles) and modeled (blue lines) values for 

the time derivative of Δp/q for an example Phase II cycle, beginning July 26, 1999, allowing 

all 3 parameters to vary (top) and fixing sd at -4.35 (bottom).
106

Figure 4-10. Recorded (blue dots) and modeled (red line) downhole pressures and percent 

error (dashed line) for an example Phase III cycle, beginning January 8, 2001, allowing all 3 

parameters to vary (top) and fixing sd at -4.35 (bottom). Recorded downhole pressures are 

calculated by adding a constant value of 6822 psi to measured surface pressures.
107

Figure 4-11. Recorded (blue dots) and modeled (red line) values for change in pressure (Δp) 

divided by flow rate (q), and the recorded (red triangles) and modeled (blue lines) values for 

the time derivative of Δp/q for an example Phase III cycle, beginning January 8, 2001, 

allowing all 3 parameters to vary (top) and fixing sd at -4.35 (bottom).
108

Figure 4-12. Recorded (blue dots) and modeled (red line) downhole pressures and percent 

error (dashed line) for an example Phase IV cycle, beginning January 14, 2007, allowing all 3 

parameters to vary (top) and fixing sd at -4.35 (bottom). Recorded downhole pressures are 

calculated by adding a constant value of 7133 psi to measured surface pressures.
109

Figure 4-13. Recorded (blue dots) and modeled (red line) values for change in pressure (Δp) 

divided by flow rate (q), and the recorded (red triangles) and modeled (blue lines) values for 

the time derivative of Δp/q for an example Phase IV cycle, beginning January 14, 2007, 

allowing all 3 parameters to vary (top) and fixing sd at -4.35 (bottom).

The earlier portions of the type curves shown above appear to indicate primarily 

the transition from wellbore storage to radial flow, with radial flow (identifiable 

as a flat spot on the derivative curve) emerging around 1000 hours in Phases IIIV, and possibly earlier in Phase I. However, the amount of scatter in the 

derivative for all phases, much of which is related to changes in flow rate, makes 

it difficult to resolve the precise times that these transitions occur. Boundary 

effects, which cause an increase in the derivative following the radial flow period, 

are not observed.

Daily-averaged pressure-flow data were used in these simulations because the 

pre-2003 high-resolution data are presently not in an accessible format. The 

modeling results suggest that daily-averaged data are insufficient to resolve the 

wellbore storage period, and that the high-resolution data will be needed.
110

It has been noted that the rate of increase in the maximum pressure during each 

cycle appears to be greater since 2010 than previously (Figure 4-1), which has led 

to speculation that the pressures are being affected by a reservoir boundary or a 

zone of decreased permeability. Figure 4-14 shows actual and modeled pressures 

for a single set of parameters for the cycles from April 2009 through March 2013. 

The mismatch between the observed and modeled maximum pressures is 

relatively small for all cycles, suggesting that a radially symmetric model with 

infinite-acting boundaries is adequate to model the data at this time. It is likely 

that the increase in maximum pressure in recent years is related to a decrease in 

the number of unscheduled maintenance shut-ins, leading to a decrease in the total 

shut-in time per cycle. 

Figure 4-14. Recorded (open squares) and modeled (red lines) downhole pressures for the 

time period from April 2009 to March 2013. The model uses a single set of input parameters, 

as shown in the bottom left corner. Recorded downhole pressures are calculated by adding a 

constant value of 7133 psi to measured surface pressures.

4.5 Evaluation of Pressures for Alternative Injection 

Scenarios

If we assume that the current reservoir characteristics will remain constant in the 

near future, it is possible to model the maximum pressures for a range of injection 

scenarios. While it is difficult to determine the precise relationship between 

seismic hazard and pressures near or away from the well, a highly simplified 

assumption is that higher maximum pressures at the well correlate with increased 

seismic hazard. Thus, it follows that for a given injection volume per unit time, 

the injection scenario that produces the lowest pressure (maximum and/or 

average) during that time is preferable, subject to logistical constraints. 

To explore the relative differences of various injection scenarios, we use the

radially symmetric model with boundary conditions developed in the preceding 
111

section. The parameters used for the forward modeling are k=18.47, sd=-4.23, and 

CD=1.06×106. These parameters are obtained from the fit to the cycle beginning 

October 19, 2011, allowing k, sd, and CD to all vary simultaneously. The 

parameters from this cycle were chosen because it was the most recent full-length 

cycle with good data for the entire cycle, as the cycle beginning April 16, 2012 

had a period near the beginning of the cycle where data was not recorded 

consistently, degrading the quality of fits, and the most recent cycle was 

substantially shorter than the other cycles, due to the premature shut-in that 

followed the earthquake of January 24, 2013. The entire flow history up to 

January 24, 2013 is used in the modeling, and a shut-in time of 83 days is 

assumed before resuming flow. The starting pressure is adjusted to constrain the 

modeled pressure to match the observed pressure at the end of the 83-day shut-in.

In deciding how to proceed with future injection, there are several factors that 

could be adjusted, including the flow rate and the frequency and duration of shutins. 

The first range of scenarios involved fixing the flow rate and total shut-in duration 

during a 1-year period, and adjusting the duration of individual shut-ins. The flow 

rate was fixed at 230 gpm and the total shut-in time per year was fixed at 40 days, 

in line with the past shut-in schedule. Scenarios varied from one 40-day shut in at 

the end of the year to 40 one-day shut-ins distributed evenly throughout the year. 

Both the maximum and average pressures increase as the shut-in length increases, 

suggesting that a large number of short shut-ins may be more effective than a 

small number of long shut-ins at controlling pressures at the well, as illustrated in 

Figure 4-15. 

Figure 4-15. Maximum and average pressures for a variety of 1-year scenarios with varying

shut-in schedules. Flow rate is fixed at 230 gpm and total shut-in time is fixed at forty days. 

Scenarios range from one 40-day shut-in to 40 one-day shut-ins.

4000

4100

4200

4300

4400

4500

4600

4700

4800

4900

0 5 10 15 20 25 30 35 40

Pressure (psi)

Shut-in Duration (days)

Wellhead Pressure vs. Shut-in Duration for Fixed Total Shut-In Days

Max. pressure

Ave. pressure
112

We consider a number of scenarios to compare the effect of variable shut-in 

duration on pressures. In these scenarios, shut-in durations are varied between 1 

and 40 days, but the total annual number of shut-in days and amount of fluid 

injected is kept fixed. All shut-ins for a particular scenario are evenly-spaced and 

of equal-length. The reason for keeping the same number of total shut-in days for 

all scenarios is to account for the assumption that a certain number of 

maintenance days each year are required for operations. For this idealized case, 

we further assume that maintenance can be carried out in 1-day increments, and 

that 40 days of maintenance per year is a reasonable average total number. The

general results do not depend on the specific increment or the total number of 

shut-in days.

The modeling results shown in Figure 4-15 suggest that shorter, more frequent 

shut-ins result in a lower pressure increase, given a specified yearly total volume 

of injection total and a total number of shut-in days. To understand the basis for 

this result, it is useful to consider the limiting case of an infinite number of shutins of infinitely short duration, and equivalent to injection at a constant rate such 

that the total yearly volume of fluid injected remains the same as the scenarios 

with shut-ins. For a constant injection rate q , injection duration I t , and shut-in 

duration St , injecting at the average injection rate q qt t t = + IIS ( ) produces the 

same total injection volume over the full injection cycle I S t t + . As an example, if 

we consider a typical 6-month injection cycle with a single 20-day shut-in, and a 

constant injection rate of 230 gpm, then the same total volume could be injected 

using an average injection rate of 205 gmp and no shut-ins. Using the simple 

Horner (1951) model described by Equation (4.3), and ignoring wellbore storage 

and skin effects, the pressure increase for these two cases is quite different 

because the model predicts that the change in pressure is proportional to the 

product of injection rate and the logarithm of injection duration. For this model, 

the difference in injection rate for these two cases has a much stronger effect than 

the difference in injection duration, and therefore the scenario with the lowest 

injection rate results in the lowest pressure increase. This effect is illustrated in 

Figure 4-16.
113

Figure 4-16. Pressure increases predicted by the simple Horner (1950) model for the cases of 

a 6-month injection cycle with either a 20-day shut and an injection rate of 230 gpm, or no 

shut in and an injection rate of 205 gpm. In each case the total amount injected over the 6-

month period is the same, but injecting at the lower flow rate produces a lower pressure 

increase even though the duration of injection is longer. This results because the slope of the 

pressure curve is greatest at the start of injection, and continues to decrease with time, and 

because the pressure increase is scaled by flow rate. For this illustration we assumed the 

following model parameters: k = 18 md, h = 100 ft, φ = 0.03, μ = 1.03 cP, and ct = 1.93×10-6

psi-1

. For simplicity, the pressure fall-off curve for the 20-day shut-in scenario is not shown.

This model is for illustration only and is not a prediction of actual PVU surface pressures.

To further investigate the effects of reducing the flow rate on injection pressures, 

several scenarios with varying average flow rates and shut-in lengths are modeled. 

In these scenarios, we assume a constant average yearly injection rate (or 

equivalently, a constant injection volume). The scenarios last for 1 year, with a 

single shut-in occurring after the total injection volume has been reached. The 

injection volume is fixed at 108 gallons, which is near the average volume per 

year injected in recent years. Scenarios range from continuous injection for the 

entire year at a flow rate of 190.26 gpm to injection at a flow rate of 380.52 gpm 

for half the year followed by a shut-in for the remainder of the year. 

The scenarios with a higher flow rate reached significantly higher maximum 

pressures during the year, as illustrated in Figure 4-17. This suggests that for a 

given injection volume, reducing the flow rate is a more effective way of reducing 

3000

3500

4000

4500

5000

0 40 80 120 160 200

Pressure Change (psi)

Time (days)

Wellhead Pressure for 0.5-year Injection Cycle

230 gpm with 20-day shut-in scenario

205 gpm with no shut-in scenario

20-day shut-in

Same cumulative

volume injected
114

the maximum pressures than increasing the shut-in time. The average pressures 

also increase as the flow rate increases, although by a much smaller amount. 

However, this interpretation relies on the assumption that effective permeability is 

independent of flow rate, which may not be valid. As discussed in Section 4.4.1, 

modeled effective permeabilities appear to have been higher during the period 

when flow rate was higher, so it is plausible that if flow rate were to be decreased, 

effective permeability would also decrease. However, the potential magnitude of 

this decrease, if it occurs at all, is unknown.

Figure 4-17. Maximum and average pressures for a variety of 1-year scenarios with varying 

flow rates, with a shut-in occurring for the remainder of the year once a fixed volume of 108

gallons is reached.

Since not all of these scenarios are logistically feasible, scenarios are also created 

using the four available plunger sizes of 2.125”, 2.000”, 1.875”, and 1.750”, 

which correspond to injection rates of 230 gpm, 204 gpm, 180 gpm, and 156 gpm, 

respectively. The values of k and CD obtained by fitting the cycle beginning 

October 19, 2011 with sd fixed at -5.028 were used. Several injection scenarios 

are considered, including shut-ins of 5, 10, 15, and 20 days every 3 or 6 months, 

and no shut-ins.

3500

4000

4500

5000

5500

6000

6500

7000

7500

8000

180 220 260 300 340 380

Pressure (psi)

Flow Rate (gpm)

Wellhead Pressure vs. Flow Rate for Fixed Injection Volume

Max. pressure

Ave. pressure
115

Figure 4-18. Maximum pressure vs. injection volume for a variety of injection scenarios with 

varying shut-in schedules and flow rates. In each data series, the data points indicate flow 

rates from left to right of 156 gpm (1.750” plungers), 180 gpm (1.875” plungers), 204 gpm 

(2.000” plungers), and 230 gpm (2.125” plungers).

Maximum pressures obtained after 1 year were compared for the four plungersize scenarios, and the results are shown in Figure 4-18. For a given injection 

volume, schedules with fewer shut-in days appear to have a smaller maximum 

pressure, confirming the results found in the previous set of scenarios. If we chose 

a maximum pressure not to be exceeded in a year of injection, this graph could be 

used to select the injection scenario that injected the greatest volume while 

staying below that pressure.

To investigate the effects of possible changes in effective permeability, models 

were also created using a fixed injection schedule and varying the effective 

permeability and flow rate, using the four available plunger sizes. The fixed 

schedule was comparable to the past injection schedule, with two 20-day shut-ins 

per year. Three permeabilities are considered: 16.18 mD, 24.19 mD, and 18.47

mD, corresponding to the minimum and maximum permeabilities fit in Phase IV 

with sd fixed at -4.35 and the permeability considered in the previous model, 

respectively. The values for sd and CD from the previous model were maintained.

The maximum pressures obtained over 1 year of injection are shown in 

Figure 4-19, which shows that even with a range of only about 8 mD, the 

effective permeability has a significant effect on the maximum pressure.

3000

3200

3400

3600

3800

4000

4200

4400

4600

4800

5000

6.00E+07 7.00E+07 8.00E+07 9.00E+07 1.00E+08 1.10E+08 1.20E+08 1.30E+08

Maximum Pressure (psi)

Injection Volume (gallons)

Max. Wellhead Pressure After 1 Year Vs. Injection Volume 

No shut-ins

5 days / 6 months

10 days / 6 months

15 days / 6 months

20 days / 6 months

5 days / 3 months

10 days / 3 months

15 days / 3 months

20 days / 3 months
116

Figure 4-19. Maximum pressures obtained over one year of injection for 3 different 

permeabilities and 3 flow rates, with 2 evenly spaced 20-day shut-ins per year.

4.6 Pressure-Flow Modeling Discussion

If we assume both that the current reservoir model is correct and that it will 

remain valid for projecting future system responses, even in the case of significant 

changes to the flow rate and injection schedule, then it follows that injection at a 

lower flow rate, with shut-ins that are significantly shorter than the current 20-day 

shut-ins, would be the most effective scenario evaluated for the reduction of both 

maximum and average wellhead pressures.

However, both these assumptions are subject to significant uncertainty. While a 

model with infinite-acting boundaries is adequate to fit the present data, it is 

possible that boundary effects may become significant in the future, which would 

lead to higher pressures. Also, temperature effects likely play a significant role in 

holding open fractures in the near wellbore region, and it is unknown whether 

these phenomena due to temperature would be sufficient to hold open these 

fractures if a lower injection flow rate were to be used, or whether decreasing the 

injection rate would lead to fracture closure, reducing the effective permeability.

Additionally, even if changes to the flow rate and schedule are successful in 

reducing the pressures at the well, it is unknown how these pressures correlate 

with pressures away from the well, or more importantly, with seismic risk.

2200

2700

3200

3700

4200

4700

5200

5700

15 16 17 18 19 20 21 22 23 24 25

Maximum Pressure (psi)

Permeability (mD)

Max. Wellhead Pressure After 1 Year Vs. Permeability

230 gpm

204 gpm

180 gpm

156 gpm
117

5 Discussion

We interpret the shallow seismicity occurring in the vicinity of the PVU injection 

well as induced by fluid injection and occurring in response to a decrease in the 

effective normal stress on pre-existing fracture surfaces. Fracture initiation is 

assumed to be adequately described by a Coulomb fracture criterion (Jaeger, 

1969), and the observed seismicity is interpreted to be the result of frictional 

failure due to shearing. Focal mechanisms analyzed to date are consistent with 

simple shear failure; tensile-failure events have not been recognized in the 

recorded data (Ake and others, 2005).

During fluid injection, the effective normal stress on pre-existing fractures can

decrease as a result of the following effects: increasing pore pressure, 

redistribution of stress (from accommodation of the injected fluid into the rock or 

from the occurrence of previously induced earthquakes), or cooling and shrinking 

of the rock matrix. The latter thermodynamic effect is important in geothermal 

areas but is unlikely to be a dominant factor at PVU, except possibly in the near 

vicinity of the injection well. Changes in pore pressure and stress redistribution 

are likely the major triggering mechanisms for the earthquakes induced by fluid 

injection at PVU.

Simple correlation of PVU injection pressure and earthquake data suggests that 

pore pressure increase is a dominant factor contributing to the rate and magnitude 

of induced seismicity occurring within 5 km of the injection well, but no clear 

correlation is observed for earthquakes occurring at greater distances, as shown in 

Figure 2-22, and discussed in greater detail in Block and Wood (2009; 2010). 

Since there are no observation wells for measuring in-situ pore pressures directly

at PVU, we use long-term averaging of the injection pressures as a simplistic 

method of accounting for the time delay and amplitude modulation of pressure 

variations at the injection well as the pressures propagate away from the well. To 

correlate pressure trends with the near-well seismicity data, the pressure 

averaging must be done over longer time windows for later periods. The 

seismicity may take progressively longer to respond to changes in injection 

pressures both because more of the seismicity occurs at greater distances from the 

well for later time periods and because the size and complexity of the reservoir 

increases over time. While acknowledging that this method is very simplistic, the 

general correlation between relatively high long-term average injection pressures 

and increased rates and magnitudes of induced earthquakes in the near-well 

region suggests that pore pressure affects both the rate and magnitudes of induced 

seismic events occurring within 5 km of the well (Figure 2-22). The same 

analysis, however, does not show a correlation between injection pressure and the 

rates and magnitudes of earthquakes occurring at distances greater than 5 km, 

such as within the NW cluster (Figure 2-22).
118

Roeloffs and Denlinger (2009) suggested that seismicity was initially triggered in 

the NW cluster by stress redistribution. This hypothesis stems from the fact that 

the seismicity in the NW cluster began only one year after the start of long-term 

injection operations, too soon for significant pore pressure changes to have 

propagated 6 to 8 km from the well, based on results from their axisymmetric 

porous-medium fluid flow model.

Ake and others (2005) suggested an alternative explanation for the early onset of 

the NW-cluster seismicity: one or more northwest-trending faults allow for fluid 

flow, and relatively rapid propagation of pore pressure, from the vicinity of the 

injection well to the NW cluster. The concept of northwest-trending, relatively 

high-permeability fault zones acting as fluid and pore pressure conduits is 

consistent with the mapped geology and the local stress field. A series of deep, 

northwest-trending normal faults was mapped during early PVU geophysical 

investigations (Block and others, 2012; King and others, 2014). The subparallel 

alignment of these faults (~N55°W), the P axes of N64°W-N67°W determined 

from the focal mechanisms of induced events (Ake and others, 2005; Block and 

others, 2015), and the direction of maximum horizontal stress of N84°W 

interpreted from borehole televiewer data (Steve Hickman, personal 

communication, 2003) suggests that the northwest-trending basement fault zones 

may tend to dilate and accommodate fluid flow during PVU injection. In this 

model, earthquakes occurring in the NW cluster, such as the January 2013 event, 

are inferred to be induced primarily by pore pressure changes. 

The lack of correlation between the average injection pressures at the well and the 

rates and magnitudes of earthquakes occurring in the NW cluster, as well as the 

lack of correlation between the near-well and NW-cluster seismicity patterns (as 

shown in Figure 2-22), is difficult to explain using a simple model of relatively 

rapid pore pressure response between the well and the NW cluster. A possible 

explanation for this discrepancy is that, if the basement fault zones are acting as 

conduits for fluid flow and pore pressure propagation to the NW cluster, then 

perhaps these pathways are pressure-sensitive. When the pore pressures in the 

vicinity of the well increase to sufficiently high levels, the conduits for fluid flow 

and pore pressure propagation to the northwest may open, and pore pressures may 

then increase relatively quickly in the NW cluster. When pore pressures in the 

vicinity of the injection well fall to sufficiently low levels, the conduits may close 

and the pressures in the NW cluster may then be relatively isolated from further 

pressure reductions at or near the injection well. Studies of the local geology and 

the spatial and temporal patterns of PVU-induced seismicity indicate that the NW 

cluster may be partially surrounded by impermeable boundaries (King and others, 

2014), and therefore elevated pore pressures in the NW cluster may require an 

unusually long time to decay. This type of pressure-sensitive response may help 

explain the lack of correlation between the seismicity pattern in the NW cluster 

and both the seismicity pattern in the near-well region and the temporal variation 

in average injection pressure. 
119

The occurrence of the January 2013 ML 4.4 earthquake, at a distance from the 

injection well nearly four times greater than the radial distances of previous PVUinduced earthquakes of comparable size, is part of a broader trend of recently 

changing seismicity patterns related to PVU fluid injection. Patterns of PVUinduced seismicity largely stabilized for a decade following a decrease in the 

injection flow rate by one-third in mid-2000. Since 2010, however, the seismicity 

patterns have been changing, with seismicity rates increasing in some areas and 

seismicity occurring in previously aseismic areas. For example, a distinct group of 

earthquakes, now referred to as the SE cluster, developed in 2010 about 6 km 

southeast of the injection well. While the first detected SE cluster event occurred 

in 2004, a total of only 3 SE cluster events were recorded prior to 2010. From 

2010 through 2012, 53 SE cluster events were recorded. Seismicity rates within 

the northern-valley areas have also changed in recent years. The number of 

northern-valley earthquakes recorded each year from 2000 (when the northern 

valley seismicity was first detected) to 2009 ranges from 2 to 33. In 2010, the rate 

increased markedly: 557 northern-valley earthquakes were recorded, with the 

majority occurring in a single swarm lasting just 16 days. Northern-valley 

seismicity rates remained elevated during 2011, with 113 earthquakes recorded, 

but declined back to pre-2010 rates during 2012, with just 10 events recorded. 

Beginning in mid-2010 and continuing through the present, several shallow 

earthquakes have been detected in the center of Paradox Valley. Three 

earthquakes were detected in 2010; two in 2011; and five in 2012. Although the 

total number of events is small (10 in 3 years), no earthquakes were detected 

beneath the floor of Paradox Valley in the 25 years of seismic monitoring prior to 

2010. The renewed spatial expansion of seismicity and increased seismicity rates 

in recent years may be related to the trend of increasing injection pressures at the 

well.

The January 2013 earthquake did not occur on any of the subsurface normal faults 

that were mapped with seismic reflection surveys during early PVU geophysical 

investigations. As discussed above, most of these faults trend close to the 

estimated direction of maximum horizontal stress (Ake and others, 2005), and 

therefore they are not now optimally oriented for shear slip. Instead, most of the 

induced earthquakes observed at PVU, including the January 2013 event, appear 

to occur on favorably oriented pre-existing fractures that have been reactivated by 

injection (Block and others, 2015). Because these faults may not have significant 

vertical offset, they may be difficult to resolve on deep seismic reflection data. 

The distribution of precise relative-relocated hypocenters as well as focalmechanism analyses provides a way to identify these faults (Block and others, 

2014; Block and others, 2015).

A thorough analysis of seismic risk from induced seismicity at PVU would 

consider the rates at which damaging ground motions are likely to occur in the 

populated surrounding areas, and the likelihood of damage, economic loss, and 

potential life-loss given those loadings. Such an analysis is beyond the scope of 
120

this report, and many of the required input parameters are not known with any 

degree of confidence. Necessary inputs to a typical seismic hazard analysis 

include mean seismicity rates, maximum magnitudes, and the magnitude 

distribution for induced seismicity, but these parameters have changed 

substantially at PVU over the lifetime of the project, making estimates for future 

behavior highly uncertain. An even bigger problem is that the precise relation 

between the seismicity rate and injection parameters which can be controlled 

(e.g., average flow rate, maximum injection pressure, shut-in duration, etc.) is 

poorly known. Additional unknown parameters include the attenuation of ground 

motions with distance, site-response properties (e.g., VS30 and basin depth), and 

the density and fragility of residences and other structures.

Long-term seismic hazard estimates will be difficult to make at PVU without 

developing models that are capable of relating induced seismicity to injection 

operations. Such models have not yet been developed for PVU, although it is

possible that recent numerical modeling approaches could be used to estimate 

both the pressure response and the seismicity response to injection. Short-term 

seismic hazard estimates can be made by extrapolating recent seismicity 

parameters, however even these estimates are likely to include substantial 

uncertainty.

Simple pressure-flow modeling suggest some operational steps that could be 

taken to reduce the wellhead pressures observed in recent years, and to potentially 

reduce the rate of occurrence of induced earthquakes. If we assume that this 

model is valid and can provide useful predictions of future pressure response to 

injection, then it follows that injection at a lower flow rate, with shut-ins that are 

significantly shorter than the current 20-day shut-ins, would be the most effective 

scenario evaluated for the reduction of both maximum and average wellhead 

pressures. 

However, the model assumptions are subject to significant uncertainty. While a 

model with infinite-acting boundaries is adequate to fit the present data, it is 

possible that boundary effects may become significant in the future, which would 

lead to higher pressures. The model assumes that reservoir permeability is 

independent of wellhead pressure, but it is possible that reduced injection pressure 

would result in reduced effective permeability, which could substantially reduce 

the amount of pressure reduction obtained from reducing the flow rates. 

Temperature effects may also play a significant role in holding open fractures in 

the near wellbore region, and it is unknown whether these phenomena due to 

temperature would be sufficient to hold open these fractures if a lower injection 

flow rate were to be used, or whether decreasing the injection rate would lead to 

fracture closure, reducing the effective permeability. Additionally, even if 

changes to the flow rate and schedule are successful in reducing the pressures at 

the well, it is unknown how these pressures correlate with pressures away from 

the well, or more importantly, with seismic risk
121

Results presented in this TM are primarily from analyses in progress before the 

January 24, 3013 earthquake, or that were done in the weeks immediately 

following it. The purpose of these analyses was to fulfill requirements of the 

Emergency Action Plan for operation of the PVU injection well, and to develop 

alternatives for operations to minimize the potential for additional felt 

earthquakes. Hence, this TM mostly contains analyses that could be completed 

quickly after the earthquake. Following completion of the initial draft of this 

report in April, 2013, more detailed analyses were conducted and submitted for 

publication in peer-reviewed scientific journals (Block and others, 2014; King and 

others, 2014; Yeck and others, 2014; Block and others, 2015), and results have 

been incorporated into the current draft of this TM where appropriate. Detailed 

analysis of a few remaining topics are in progress, including pressure-flow 

modeling and seismic hazard analysis, and are expected to be documented in 

peer-reviewed journal articles.

.

123

6 Conclusions

6.1 Analysis of the January 24, 2013 ML 4.4 Earthquake

Our analysis indicates that the ML 4.4 earthquake that occurred on January 24, 

2013 (UTC; January 23, 2013 MST) was induced by long-term fluid injection into 

PVU injection well #1. This event is slightly larger than the ML 4.3 earthquake of 

May 27, 2000, making it the largest PVU-induced event to date. It locates 8.2 km 

northwest of the injection well, a radial distance nearly four times greater than 

other PVU-induced earthquakes of comparable magnitude. 

Because this event occurred closer to northern Paradox Valley - where much of 

the local population resides - than previous large induced earthquakes, strong 

ground shaking was experienced by many local residents. Ground motion 

recordings from a strong motion instrument located at the Paradox Community 

Center indicate a peak horizontal acceleration of 0.28 g. For frequencies above 

2.5 Hz, the horizontal accelerations recorded in the town of Paradox were nearly 2 

standard deviations higher than the mean value estimated by standard ground 

motion prediction equations for an earthquake of this size and distance. Possible 

reasons for the unusually large ground vibrations in northern Paradox Valley 

produced by this earthquake include: (1) site effects resulting form the local 

geologic structure, consisting of a thick salt section and soft soil layer overlying 

bedrock; (2) unusual diffraction effects related to the “salt wall”; or, (3) decreased 

attenuation of high-frequency ground motions resulting from the shallow focal 

depth.

The peak particle velocity generated by this earthquake in the town of Paradox 

was 6.2 cm/sec (2.4 in/sec, horizontal component). To put this in perspective, this 

ground vibration level is more than 3 times the 1.9 cm/sec (0.75 in/sec) limit set 

in the Department of Interior, Office of Surface Mining (OSM) regulations for

control of vibrations from mine blasting. As a result, minor damage to structures 

located in northern Paradox Valley may have occurred.

6.2 Potential Future Seismicity

6.2.1 Short-Term Seismicity

The apparent insensitivity of NW-cluster seismicity to past injection operations 

suggests that either the response time of this seismicity to injection operations is 

very long (several years), or that the threshold for triggering earthquakes on the 

pre-existing faults in this area is low (earthquakes over a range of magnitudes 

may occur with very little stress or pore pressure perturbation). Hence, we must 
124

conclude that the risk of additional felt earthquakes in the NW cluster, some 

comparable in magnitude to (or even larger than) the January event, will be 

significant regardless of any changes in injection operations that may occur.

6.2.2 Seismic Hazard Factors

Reliably forecasting the long-term behavior of induced seismicity at PVU is 

extremely difficult because of the large uncertainties involved. Here we 

summarize a few of the factors that we have examined in this report that relate to

the long-term seismic hazard.

6.2.2.1 Geographical Distribution of Seismicity

PVU-induced seismicity has expanded spatially over time. Most of the 

earthquakes induced during the injection tests (1991-1995) occurred within 2.5 

km of the injection well. During continuous injection operations (1996-present), 

clusters of induced earthquakes have developed in areas much farther from the 

injection well, with shallow, potentially-induced earthquakes now occurring as far 

as 16 km from the injection well. The geographical expansion of seismicity 

slowed markedly after the injection flow rate was reduced in mid-2000, and 

seismicity rates simultaneously declined. However, over the last 3 years, new 

areas of seismicity have appeared and the seismicity rates in other areas have 

increased. Some of the recent seismicity is occurring closer to the more-populated 

northern end of Paradox Valley than earlier seismicity.

The volume of the stimulated rock mass, as indicated by the clouds of induced 

seismicity, is orders of magnitude larger than the cumulative volume of fluid 

injected. One possible explanation for this, which is consistent with the 

dominance of clear linear patterns within the induced seismicity, is that the 

majority of the injected fluid is traveling through a network of discrete, preexisting fractures and that very little of the fluid is entering the pore space within 

the larger volume of rock surrounding the well. This hypothesis is consistent with 

the natural low permeability of the target reservoirs. An implication of this 

hypothesis is that the pore pressures within these fractures may be substantially 

increased at large distances from the injection well, while pressures within the 

pore spaces in isolated sections of the rock mass not connected to the fracture 

network may be considerably lower. Pore pressure increase within fractures at 

large distances from the injection well may contribute to the occurrence of 

induced seismicity distant from the well, such as that in the NW cluster.

6.2.2.2 Seismicity Rates

Both the overall rate and magnitude distribution of PVU-induced seismicity has 

varied substantially over time. These characteristics also vary spatially, with, for 

example, the near-well region and NW cluster showing distinctly different trends.

The temporal variation in rate and magnitude distribution makes it difficult to 

predict the future behavior of the induced seismicity. Any quantitative measures 
125

of seismicity rates based on current characteristics or historical trends may not be 

extrapolated with confidence very far into the future. Also, because the seismicity 

appears to have a long response time to changes in injection operations, variations 

in the rates or magnitudes of induced earthquakes related to changes in injection 

operations may not be noticeable until several years after those operational 

changes occur.

6.2.2.3 Maximum Earthquake Magnitude

Studies indicate that the maximum magnitude of induced earthquakes may be 

related to several factors, including the cumulative volume of fluid injected, the 

volume of rock stimulated by injection, and both the magnitude of pore pressure 

and the spatial pore pressure gradient. The models considered here for examining 

the question of the maximum magnitude earthquake likely to be induced for a 

given injection site are simplistic, considering the effect of one of these variables 

at a time. An evaluation of the combined effects of multiple parameters is 

complex and beyond the scope of this report.

The maximum magnitudes of earthquakes induced at PVU over time show a 

general correlation with the cumulative volume of injected fluid. This correlation 

indicates that, as of the end of 2012, enough fluid has been injected to trigger a 

maximum earthquake with magnitude between 3.7 and 4.8. This model suggests 

that the maximum earthquake magnitude will increase as additional fluid is 

injected. However, because this relationship is logarithmic, the trend is largely 

determined by the early history of injection operations, including the fluid volume 

injected during the tests conducted between 1991 and 1995. This type of analysis 

does not provide enough resolution to evaluate whether the rate of increase in 

event magnitude as a function of fluid volume injected changed when the 

injection flow rate was decreased in 2000. If we assume that the same trend is 

valid for all periods of time, and if injection were to continue at the same average 

flow rate (as in 2000-2012) for an additional 5 years, then the relationship implies 

a maximum earthquake magnitude of between 3.8 and 4.9 at the end of 2017. 

After 10 years (end of 2022), the maximum magnitude estimate increases to 3.9 -

5.05.

A model that correlates maximum earthquake magnitude with volume of 

stimulated rock (as estimated by the spatial distribution of induced seismicity) 

indicates generally similar results as that above but allows for separate analyses 

for the near-well region and NW cluster. This model indicates a maximum 

earthquake magnitude of about 5.2 for the near-well area and about 4.4 for the 

NW cluster. According to this model, the majority of the risk for generating 

earthquakes of this size is due to the volume of rock mass that was stimulated by 

fluid injected prior to 2001. Since the spatial extent of the induced seismicity in 

the near-well and NW cluster areas has not expanded substantially since 2001, the

corresponding estimates of volume of stimulated rock mass and maximum 

earthquake magnitude have not increased much since that time. This type of 

model suggests that the maximum earthquake magnitude will not increase 
126

substantially beyond these current estimates unless the near-well or NW-cluster 

seismicity clouds expand in the future (either during or after additional injection).

This model provides estimates of maximum earthquake magnitudes due to fluid 

injection that has already occurred at PVU, but it does not predict when these 

events may occur. Because of the potentially long time period over which pore 

pressure diffusion takes place, the largest induced earthquake could occur years 

after the risk has been elevated, and potentially after injection has ceased. 

A general correlation between high average injection pressures and the rate and 

magnitudes of earthquakes induced within 5 km of the PVU injection well has 

been noted previously (Block and Wood, 2009). When long-term average 

injection pressures fall to sufficiently low levels, the occurrence of M 2.5+ nearwell induced earthquakes ceases while smaller-magnitude events continue to 

occur. This analysis implies that the local pore pressure may need to exceed some 

threshold level in order for the largest-magnitude events to occur. A higher pore 

pressure perturbation for larger-magnitude earthquakes is also indicated by 

numerical modeling of pore pressure diffusion on a pre-stressed fault (Garagash 

and Germanovich, 2012). The empirical observations and numerical model imply 

that the risk of inducing large-magnitude earthquakes can be reduced if pore 

pressures are sufficiently lowered. While no such empirical correlation between 

injection pressures and the occurrence of seismicity at distances greater than 5 km 

from the well (such as in the NW cluster) has been observed, the numerical 

modeling nonetheless suggests that a reduction in pore pressure could potentially 

reduce the seismic risk in this area as well. The response of pore pressures in 

fractures located several km from the well to changes in injection operations is 

unknown however, both in terms of the magnitude of pressure change and the 

time delay for changes to occur.

The spatial gradient of the pore pressures may also play a role in triggering largemagnitude earthquakes. Finite element modeling of pressures and induced 

seismicity for an injection simulation shows a correlation between the largestmagnitude induced earthquakes and shallow spatial pressure gradients (Baisch 

and others, 2010). Results from this modeling are consistent with the observation 

that in many instances the largest injection-induced earthquakes occur after 

injection operations cease. Under this situation, pore pressure continues to rise in 

some areas after shut-in but with a shallower spatial gradient than during 

injection. The results are also consistent with the observation that in many cases 

the largest-magnitude events tend to occur near the edges of the stimulated region 

rather than closer to the injection well. Such a pattern is observed at PVU, with 

the largest induced earthquakes occurring near the edges of the near-well and 

NW-cluster regions. The implication of this modeling is that the magnitudes of 

PVU-induced earthquakes could continue to increase after injection operations are 

stopped at some point in the future, and that these large events are most likely to 

occur around the edges of the stimulated regions (including the NW cluster).
127

6.3 Injection Operations

Simple pressure-flow modeling performed using industry-standard well testing 

software does not indicate any measurable change in reservoir properties in recent 

years. Permeability and near-wellbore storage values appear to be stable. This 

modeling implies that the increase in injection pressures (and expanding 

seismicity) may be a natural consequence of long-term injection into a 

(potentially unbounded) reservoir with limited permeability.

The modeling also indicates that the current pattern of injection may not be the 

most efficient in terms of injecting a given volume of fluid with the least amount 

of pressure rise. In model simulations, the current injection scenario, in which 2 

20-day shut-ins occur each year, causes a greater pressure rise than scenarios 

where the same volume of fluid is injected at the same rate but more frequent, 

shorter shut-ins are used. Simulations also suggest that the same volume could be 

injected with even less pressure rise if no shut-ins are used and instead the flow 

rate is decreased. However, because permeability may change substantially when 

flow rate or pressures decline (and fractures potentially close), these model 

simulations have significant uncertainties.

Pressure modeling indicates that keeping the same shut-in schedule and

decreasing the flow rate from 230 gpm to 204 gpm, a decrease of approximately 

11%, would decrease the maximum pressure seen in a 1-year period by about 570 

psi relative to maintaining the same flow rate, assuming no change in 

permeability. The corresponding decrease in the annual average pressure is 

estimated to be 480 psi. Keeping the same flow rate and number of shut-in days 

but changing the shut-in schedule from 2 20-day shut-ins per year to 40 1-day 

shut-ins per year would decrease the maximum pressure by approximately 250 

psi, and decrease the annual average pressure by about 70 psi, with no change in 

the injected volume. We expect that the combined effect of implementing both of 

these changes would be greater than either of the individual effects. 

129

References

Abrahamson, N., and Silva, W., 2008, Summary of the Abrahamson & Silva NGA GroundMotion Relations: Earthquake Spectra, v. 24, no. 1, p. 67-97.

Abramowitz, M., and Stegun, I.A., 1964, Handbook of Mathematical Functions with Formulas, 

Graphs, and Mathematical Tables: Washington, D.C., United States Department of 

Commerce, National Bureau of Standards.

Agarwal, R.G., Al-Hussainy, R., and Jr., H.J.R., 1970, An Investigation of Wellbore Storage 

Constant and Skin Effect in Unsteady Liquid Flow: I. Analytical Treatment: Society of 

Petroleum Engineers Journal, v. 10, no. 3, p. 279-290.

Ake, J., Mahrer, K., O, C.D., and Block, L., 2005, Deep-injection and closely monitored induced 

seismicity at Paradox Valley, Colorado: Bulletin of the Seismological Society of 

America, v. 95, no. 2, p. 664-683.

Aki, K., 1965, Maximum Likelihood Estimate of b in the Formula logN=a-bM and its Confidence 

Limits: Bulletin of the Earthquake Research Institute, v. 43, p. 237-237.

Alt, B., Henry, S., Kraeger-Rovey, C., Paulin, I.J., Ellison, D.T., Javernick, H.J., and Nelligan, 

M.T., 1980, Regulations of the Colorado Mined Land Reclamation Board for Coal 

Mining: Colorado Division of Minerals and Geology, 469 p.

Anderson, J.G., 1979, Estimating the seismicity from geological structure for seismic-risk studies: 

Bulletin of the Seismological Society of America, v. 69, no. 1, p. 135-158.

Anderson, J.G., and Luco, J.E., 1983, Consequences of slip rate constraints on earthquake 

occurrence relations: Bulletin of the Seismological Society of America, v. 73, no. 2, p. 

471-496.

Bachmann, C.E., Wiemer, S., Woessner, J., and Hainzl, S., 2011, Statistical analysis of the 

induced Basel 2006 earthquake sequence: introducing a probability-based monitoring 

approach for Enhanced Geothermal Systems: Geophysical Journal International, v. 186, 

no. 2, p. 793-807.

Baisch, S., Vörös, R., Rothert, E., Stang, H., Jung, R., and Schellschmidt, R., 2010, A numerical 

model for fluid injection induced seismicity at Soultz-sous-Forêts: International Journal 

of Rock Mechanics and Mining Sciences, v. 47, no. 3, p. 405-413.

Barth, A., Wenzel, F., and Langenbruch, C., 2013, Probability of earthquake occurrence and 

magnitude estimation in the post shut-in phase of geothermal projects: Journal of 

Seismology, v. 17, no. 1, p. 5-11.

Bear, J., 1972, Dynamics of Fluids in Porous Media, Dover Publ.

Bender, B., 1983, Maximum likelihood estimation of b values for magnitude grouped data: 

Bulletin of the Seismological Society of America, v. 73, no. 3, p. 831-851.
130

Berrill, J.B., and Davis, R.O., 1980, Maximum entropy and the magnitude distribution: Bulletin of 

the Seismological Society of America, v. 70, no. 5, p. 1823-1831.

Biot, M.A., 1941, General Theory of Three‐Dimensional Consolidation: Journal of Applied 

Physics, v. 12, no. 2, p. 155-164.

Block, L., and Wood, C., 2009, Overview of PVU-Induced Seismicity from 1996 to 2009 and 

Implications for Future Injection Operations: Bureau of Reclamation, 16 pps. p.

Block, L., Yeck, W., King, V., Derouin, S., and Wood, C., 2012, Review of Geologic 

Investigations and Injection Well Site Selection, Paradox Valley Unit, Colorado: Bureau 

of Reclamation, 62 p.

Block, L.V., and Wood, C.K., 2010, Evolving characteristics of seismicity induced by long-term 

fluid injection at Paradox Valley, Colorado, American Geophysical Union Fall Meeting 

abstract #S13B-2000.

Block, L.V., Wood, C.K., Yeck, W.L., and King, V.M., 2014, The 24 January 2013 ML 4.4 

Earthquake near Paradox, Colorado, and Its Relation to Deep Well Injection: 

Seismological Research Letters, v. 85, no. 3, p. 609-624.

Block, L.V., Wood, C.K., Yeck, W.L., and King, V.M., 2015, Induced seismicity constraints on 

subsurface geological structure, Paradox Valley, Colorado: Geophysical Journal 

International, v. 200, no. 2, p. 1170-1193.

Boore, D.M., and Atkinson, G.M., 2008, Ground-Motion Prediction Equations for the Average 

Horizontal Component of PGA, PGV, and 5%-Damped PSA at Spectral Periods between 

0.01 s and 10.0 s: Earthquake Spectra, v. 24, no. 1, p. 99-138.

Bremkamp, W., and Harr, C.L., 1988, Area of Least Resistance to Fluid Movement and Pressure 

Rise, 49 p. p.

Brune, J.N., 1970, Tectonic Stress and the Spectra of Seismic Shear Waves from Earthquakes: J. 

Geophys. Res., v. 75, no. 26, p. 4997-5009.

Brune, J.N., 1971, Correction [to “Tectonic stress and the spectra, of seismic shear waves from 

earthquakes”]: Journal of Geophysical Research, v. 76, no. 20, p. 5002-5002.

Campbell, K.W., and Bozorgnia, Y., 1987, Campbell-Bozorgnia NGA Ground Motion Relations 

for the Geometric Mean Horizontal Component of Peak and Spectral Ground Motion 

Parameters: Pacific Earthquake Engineering Research Center, University of California 

UCB/PEER 2007/02 [Technical report through June 2006].

Campbell, K.W., and Bozorgnia, Y., 2008, NGA Ground Motion Model for the Geometric Mean 

Horizontal Component of PGA, PGV, PGD and 5% Damped Linear Elastic Response 

Spectra for Periods Ranging from 0.01 to 10 s: Earthquake Spectra, v. 24, no. 1, p. 139-

171.

Chiou, B.S.J., and Youngs, R.R., 2008, An NGA Model for the Average Horizontal Component of 

Peak Ground Motion and Response Spectra: Earthquake Spectra, v. 24, no. 1, p. 173-215.

Chung, D.H., and Bernreuter, D.L., 1981, Regional relationships among earthquake magnitude 

scales: Reviews of Geophysics and Space Physics, v. 19, no. 4, p. 649-663.
131

Convertito, V., Maercklin, N., Sharma, N., and Zollo, A., 2012, From Induced Seismicity to Direct 

Time‐Dependent Seismic Hazard: Bulletin of the Seismological Society of America, v. 

102, no. 6, p. 2563-2573.

Cornell, C.A., 1968, Engineering seismic risk analysis: Bulletin of the Seismological Society of 

America, v. 58, no. 5, p. 1583-1606.

Cornell, C.A., and Vanmarcke, E.H., 1969, The Major Influences on Seismic Risk, in Proceedings 

of the 4th World Conference on Earthquake Engineering, Santiago, Chile, p. 69-83.

Cox, D.O., Stinson, S.H., and Stellavato, J.N., 2000, Well Testing in Ultra-High Permeability 

Formations, SPE Annual Technical Conference and Exhibition: Dallas, Texas, Copyright 

2000, Society of Petroleum Engineers Inc.

Detournay, E., and Cheng, A.H.-D., 1993, Fundamentals of Poroelasticity, in Fairhurst, C., ed., 

Comprehensive Rock Engineering: Principles, Practice and Projects, Vol. II, Analysis 

and Design Method: Pergamon Press, p. 113-171.

Ellsworth, W.L., 2003, Magnitude and Area Data for Strike Slip Earthquakes, Working Group on 

California Earthquake Probabilities, Earthquake probabilities in the San Francisco Bay 

region - 2002-2031, U.S. Geological Survey Open-file Report 03-214: U.S. Geological 

Survey, p. 6.

Envirocorp Services and Technology Inc., 1995, Report of Evaluation of Injection Testing for 

Paradox Valley Injection Test No. 1: Durango, CO, Bureau of Reclamation, p. 67 p.

Fekete, 2012, WellTest User Manual (7.6.0 ed.): Calgary, Fekete Associates Inc.

Frohlich, C., 1989, The Nature of Deep-Focus Earthquakes: Annual Review of Earth and 

Planetary Sciences, v. 17, no. 1, p. 227-254.

Frohlich, C., and Davis, S.D., 1993, Teleseismic b values; Or, much ado about 1.0: Journal of 

Geophysical Research: Solid Earth, v. 98, no. B1, p. 631-644.

Garagash, D.I., and Germanovich, L.N., 2012, Nucleation and arrest of dynamic slip on a 

pressurized fault: Journal of Geophysical Research: Solid Earth, v. 117, no. B10, p. 27 

pps.

Goertz-Allmann, B.P., Goertz, A., and Wiemer, S., 2011, Stress drop variations of induced 

earthquakes at the Basel geothermal site: Geophysical Research Letters, v. 38, no. 9, p. 

L09308.

Gutenberg, B., and Richter, C.F., 1954, Seismicity of the Earth and Associated Phenomena (2nd 

ed.): Princeton, New Jersey, Princeton University Press, 310 p.

Hanks, T.C., and Bakun, W.H., 2002, A bilinear source-scaling model for M-log A observations

of continental earthquakes: Bulletin of the Seismological Society of America, v. 92, no. 

5, p. 1841-1846.

Hanks, T.C., and Kanamori, H., 1979, A moment magnitude scale: Journal of Geophysical 

Research: Solid Earth, v. 84, no. B5, p. 2348-2350.

Harr, C., 1989, Proposed Location for Injection Well No. 2: Grand Junction, CO, Bureau of 

Reclamation, p. 33 p.
132

Horner, D.R., 1951, Pressure Build-up in Wells, 3rd World Petroleum Congress: The Hague, 

Netherlands, World Petroleum Congress, p. 25-43.

Jaeger, J.C., 1969, Elasticity, Fracture and Flow with Engineering and Geological Applications 

(Third ed.): London, Chapman & Hall, 268 p.

Jahr, T., Jentzsch, G., Gebauer, A., and Lau, T., 2008, Deformation, seismicity, and fluids: Results 

of the 2004/ 2005 water injection experiment at the KTB/ Germany: Journal of 

Geophysical Research, v. 113, no. B11, p.; Citation B11410.

Jahr, T., Jentzsch, G., Letz, H., and Gebauer, A., 2007, Tilt observations around the KTB-site, 

Germany; monitoring and modelling of fluid induced deformation of the upper crust of 

the Earth, in Tregoning, P., and Rizos, C., eds., Dynamic planet; monitoring and 

understanding a dynamic planet with geodetic and oceanographic tools.: Springer-Verlag. 

Berlin - Heidelberg, Federal Republic of Germany. 2007.

Kaklamanos, J., Baise, L.G., and Boore, D.M., 2011, Estimating Unknown Input Parameters when 

Implementing the NGA Ground-Motion Prediction Equations in Engineering Practice: 

Earthquake Spectra, v. 27, no. 4, p. 1219-1235.

Kanamori, H., and Anderson, D.L., 1975, Theoretical basis of some empirical relations in 

seismology: Bulletin of the Seismological Society of America, v. 65, no. 5, p. 1073-1095.

Keranen, K.M., Savage, H.M., Abers, G.A., and Cochran, E.S., 2013, Potentially induced 

earthquakes in Oklahoma, USA: Links between wastewater injection and the 2011 Mw 

5.7 earthquake sequence: Geology.

Kharaka, Y.K., Ambats, G., Thordsen, J.J., and Davis, R.A., 1997, Deep well injection of brine 

from Paradox Valley, Colorado; potential major precipitation problems remediated by 

nanofiltration: Water Resources Research, v. 33, no. 5, p. 1013-1020.

King, V., Block, L., Yeck, W., Wood, C., and Derouin, S., submitted, Geological Structure of the 

Paradox Valley Region, Colorado, and Relationship to Seismicity Induced by Deep Well 

Injection: Journal of Geophysical Research: Solid Earth.

King, V.M., Block, L.V., Yeck, W.L., Wood, C.K., and Derouin, S.A., 2014, Geological structure 

of the Paradox Valley Region, Colorado, and relationship to seismicity induced by deep 

well injection: Journal of Geophysical Research: Solid Earth, v. 119, no. 6, p. 4955-4978.

Lee, S.H., and Ghassemi, A., 2010, A Three-Dimensional Thermo-Poro-Mechanical Finite 

Element Analysis of a Wellbore On Damage Evolution, 44th U.S. Rock Mechanics 

Symposium and 5th U.S.-Canada Rock Mechanics Symposium: Salt Lake City, Utah, 

American Rock Mechanics Association.

Lee, S.H., and Ghassemi, A., 2011, Three-Dimensional Thermo-Poro-mechanical Modeling of 

Reservoir Stimulation and Induced Microseismicity in Geothermal Reservoir Stanford 

Geothermal Workshop.

Lee, W.H.K., and Stewart, S.W., 1981, Principles and Applications of Microearthquake networks: 

New York, Academic Press, v. Supplement 2.

Lei, X., Yu, G., Ma, S., Wen, X., and Wang, Q., 2008, Earthquakes induced by water injection at 

approximately 3 km depth within the Rongchang gas field, Chongqing, China: Journal of 

Geophysical Research, v. 113, no. B10, p.; Citation B10310.
133

Mahrer, K., Ake, J., O'Connell, D., and Block, L., 2004, 2003 Status Report - Paradox Valley 

Seismic Network, Paradox Valley Unit, Southwestern Colorado: Bureau of Reclamation 

Technical Memorandum No. D8330-2004-04, 126 p.

Mahrer, K., Block, L., and Wood, C., 2008, 2007 Status Report, Paradox Valley Seismic Network, 

Paradox Valley Project, Colorado: Bureau of Reclamation, 96 p.

Majer, E., Nelson, J., Robertson-Tait, A., Savy, J., and Wong, I., 2012, Protocol for Addressing 

Induced Seismicity Associated with Enhanced Geothermal Systems: U. S. Department of 

Energy Geothermal Technologies Program, 46 p.

Marzocchi, W., and Sandri, L., 2003, A review and new insights on the estimation of the b-value 

and its uncertainty: Annals of Geophysics, v. 46, no. 6, p. 1271-1282.

McClure, M., and Horne, R., 2012, The Effect of Fault Zone Development on Induced Seismicity 

in Stanford Geothermal Workshop.

McGarr, A., 1976, Seismic moments and volume changes: Journal of Geophysical Research, v. 

81, no. 8, p. 1487-1494.

McGarr, A., Simpson, D., and Seeber, L., 2002, Case histories of induced and triggered seismicity, 

in William H.K. Lee, H.K.P.C.J., and Carl, K., eds., International Geophysics: Academic 

Press, p. 647-661.

Mena, B., Wiemer, S., and Bachmann, C., 2013, Building Robust Models to Forecast the Induced 

Seismicity Related to Geothermal Reservoir Enhancement: Bulletin of the Seismological 

Society of America, v. 103, no. 1, p. 383-393.

Molnar, P., 1979, Earthquake recurrence intervals and plate tectonics: Bulletin of the 

Seismological Society of America, v. 69, no. 1, p. 115-133.

Nicholls, H.R., Johnson, C.F., and Duvall, W.I., 1971, Blasting Vibrations and Their Effects on 

Structures: U. S. Department of the Interior, Office of Surface Mining, Reclamation and 

Enforcement, 110 p.

Nicol, A., Carne, R., Gerstenberger, M., and Christophersen, A., 2011, Induced seismicity and its 

implications for CO2 storage risk: Energy Procedia, v. 4, no. 0, p. 3699-3706.

Power, M., Chiou, B., Abrahamson, N., Bozorgnia, Y., Shantz, T., and Roblee, C., 2008, An 

Overview of the NGA Project: Earthquake Spectra, v. 24, no. 1, p. 3-21.

Richter, C.F., 1935, An instrumental earthquake magnitude scale: Bulletin of the Seismological 

Society of America, v. 25, no. 1, p. 1-32.

Roeloffs, E., and Denlinger, R., 2009, An Axisymmetric Coupled Flow and Deformation Model 

for Pore Pressure Caused by Brine Injection in Paradox Valley, Colorado: Implications 

for the Mechanisms of Induced Seismicity, Preliminary Peport to the Bureau of 

Reclamation: U.S. Geological Survey, 31 p.

Shapiro, S.A., and Dinske, C., 2009, Scaling of seismicity induced by nonlinear fluid-rock 

interaction: Journal of Geophysical Research: Solid Earth, v. 114, no. B9, p. B09307.

Shapiro, S.A., Dinske, C., and Kummerow, J., 2007, Probability of a given-magnitude earthquake 

induced by a fluid injection: Geophysical Research Letters, v. 34, no. 22, p.; L22314.
134

Shapiro, S.A., Dinske, C., Langenbruch, C., and Wenzel, F., 2010, Seismogenic index and 

magnitude probability of earthquakes induced during reservoir fluid stimulations: The 

Leading Edge, v. 29, no. 3, p. 304-309.

Shapiro, S.A., Huenges, E., and Borm, G., 1997, Estimating the crust permeability from fluidinjection-induced seismic emission at the KTB site: Geophysical Journal International, v. 

131, no. 2, p. 15-18.

Shapiro, S.A., Krüger, O.S., Dinske, C., and Langenbruch, C., 2011, Magnitudes of induced 

earthquakes and geometric scales of fluid-stimulated rock volumes: Geophysics, v. 76, 

no. 6, p. WC55-WC63.

Shapiro, S.A., Patzig, R., Rothert, E., and Rindschwentner, J., 2003, Triggering of Seismicity by 

Pore-pressure Perturbations: Permeability-related Signatures of the Phenomenon: Pure 

and Applied Geophysics, v. 160, no. 5-6, p. 1051-1066.

Shapiro, S.A., Rentsch, S., and Rothert, E., 2005, Characterization of hydraulic properties of rocks 

using probability of fluid-induced microearthquakes: Geophysics, v. 70, no. 2, p. F27-

F33.

Shapiro, S.A., Rothert, E., Rath, V., and Rindschwentner, J., 2002, Characterization of fluid 

transport properties of reservoirs using induced microseismicity, in Shapiro Serge, A., 

and Gurevich, B., eds., Seismic signatures of fluid transport.: Society of Exploration 

Geophysicists. Tulsa, OK, United States. 2002.

Shaw, B.E., 2009, Constant Stress Drop from Small to Great Earthquakes in Magnitude-Area 

Scaling: Bulletin of the Seismological Society of America, v. 99, no. 2A, p. 871-875.

Shi, Y., and Bolt, B.A., 1982, The standard error of the magnitude-frequency b value: Bulletin of 

the Seismological Society of America, v. 72, no. 5, p. 1677-1687.

Siskind, D.E., Stagg, M.S., Kopp, J.W., and Dowding, C.H., 1980, Structure Response and 

Damage Produced by Ground Vibration From Surface Mine Blasting: U. S. Department 

of the Interior, Office of Surface Mining, Reclamation and Enforcement, 83 p.

Subsurface Technology, 2001, Report of Operations: Houston, no. 60D5207, p. 24 p. .

Tinti, S., and Mulargia, F., 1987, Confidence intervals of b values for grouped magnitudes: 

Bulletin of the Seismological Society of America, v. 77, no. 6, p. 2125-2134.

Utsu, T., 1966, A Statistical Significance Test of the Difference in b-value between Two 

Earthquake Groups: Journal of Physics of the Earth, v. 14, no. 2, p. 37-40.

Wang, H., 2000, Theory of linear poroelasticity with applications to geomechanics and 

hydrogeology: Princeton, NJ, Princeton University Press.

Wang, X., and Ghassemi, A., 2012, A 3 D Thermal-poroelastic Model for Geothermal Reservoir 

Stimulation, Stanford Geothermal Workshop.

Wang, Y.C., 2008, Discrete element simulation of hydraulic fracturing and induced seismicity in 

engineered geothermal systems, in Gurgenci, H., and Budd, A., eds., Proceedings of the 

Sir Mark Oliphant international frontiers of science and technology Australian 

geothermal energy conference.: Geoscience Australia. Canberra, A.C.T., Australia. 2008.
135

Wells, D.L., and Coppersmith, K.J., 1994, New empirical relationships among magnitude, rupture 

length, rupture width, rupture area, and surface displacement: Bulletin of the 

Seismological Society of America, v. 84, no. 4, p. 974-1002.

Yeck, W.L., Block, L.V., Wood, C.K., and King, V.M., 2014, Maximum magnitude estimations of 

induced earthquakes at Paradox Valley, Colorado, from cumulative injection volume and 

geometry of seismicity clusters: Geophysical Journal International, v. 200, no. 1, p. 322-

336.

A-1

Appendix A Earthquake Catalog

The earthquake catalog used for the analyses presented in this report is provided 

in the attached file “TM 86-68330-2013-12 earthquake catalog.csv”. This commadelimited ascii file contains a catalog of earthquakes that occurred within 

approximately 10 km of the perimeter of the Paradox Valley Seismic Network 

(PVSN) from 1985 through July, 2013 (Figure A-1). The earthquakes were 

detected and their locations computed primarily using data from PVSN. Data 

from the nearby Ridgway Seismic Network, also operated by Reclamation, 

supplemented the PVSN data for some of the events that occurred from 1985 to 

2007. Identified explosions are excluded from the catalog. 

Figure A-1. Geographical extent of PVSN earthquake catalog.
A-2

A.1 Earthquake Location Methods

This catalog contains one hypocenter for each earthquake. The hypocenter was 

computed using either a relative event location procedure or an absolute event 

location method. The following description of the event location methods is taken 

from Block et al. (2015):

We calculate two sets of hypocenters from the PVSN data. Initially, we 

compute absolute earthquake locations using manually-determined Pwave and S-wave arrival times. These hypocenters are computed using 

three-dimensional (3-D) P-wave and S-wave velocity models developed 

from hypocenter-velocity-station correction inversions of data from 

naturally-occurring and induced local earthquakes and explosions 

recorded by PVSN. Both the hypocenter-velocity inversion and 

earthquake location software we use were developed in-house and 

represent an extension of the work of Block (1991). We subsequently 

compute precise relative event locations from inversion of arrival time 

differences obtained from cross-correlations of P-wave and S-wave 

arrivals extracted from waveforms recorded at the same station for pairs of 

events. No time differences from manual time picks are included, with the 

exception of data from the five events with magnitude of 3.5 or greater. 

Because most of the waveforms for these events are clipped but many of 

the first breaks are clear, we incorporate them into the event relative 

location using differences of high-quality manually picked arrival times. 

For the relative location inversion, we use a double-difference-type 

algorithm. The 3-D velocity models described above are used and remain 

fixed during the relative location inversion. No event clustering is 

performed, so an event may tie to any other event in the dataset. Data from 

event pairs with small separation distances are weighted most strongly in 

the inversion; data from event pairs with larger separation distances (as a 

fraction of the event-to-station distance) are down-weighted or eliminated 

from the inversion to reduce potential bias from unmodeled velocity 

variations in the source regions. Because the absolute earthquake locations 

are not well-constrained by the time difference data alone, we keep the 

locations of a few widely-spaced (> ~3 km apart) events with wellconstrained absolute locations fixed during the relative location inversion. 

An earthquake must tie, either directly or indirectly (possibly through 

multiple event pairs), to an event with a fixed location to be retained in the 

relative location inversion.

Approximately 87% of all induced earthquakes and 93% of induced 

earthquakes with duration magnitude (MD) ≥ 0 occurring within 10 km of 

the PVU injection well are considered well-constrained in the relative 

location inversion. We assign these hypocenters a quality factor of ‘a’. 

For the remaining events, we retain the absolute locations and assign a 

hypocenter quality factor of ‘b’.
A-3

Both the absolute and relative location uncertainties affect the reliability 

of the seismicity patterns produced. When considering the relative 

locations of widely spaced earthquakes, the absolute location uncertainties 

of the events fixed during the hypocenter relative location are the most 

important factor. The standard errors of the absolute locations of these 

fixed events vary from 50 to 100 meters horizontally and 70 to 150 meters 

vertically. When considering the spatial distribution of closely spaced 

events (< ~ 3 km apart), the relative location uncertainties control the 

spatial resolution of the hypocenter patterns. Preliminary analysis of the 

relative location errors indicates that the error of an individual earthquake 

with respect to all other tied events is generally less than 50 meters 

horizontally and 100 meters vertically. Relative location errors of 

earthquakes in close proximity to each other are likely less, but a detailed 

analysis of the relative location errors as a function of event separation 

distance is still in progress. The absolute errors of the events that do not 

tie into the relative location procedure (the b-quality events) have 

significantly higher errors, with average standard errors of approximately 

500 meters horizontally and 800 meters vertically.

A.2 Description of Data Columns

The earthquake catalog file contains the following columns:

Event_ID: Identification number of event

Year: Earthquake origin time year

Month: Earthquake origin time month

Day: Earthquake origin time day

Hour: Earthquake origin time hour

Minute: Earthquake origin time minute

Second: Earthquake origin time second

Latitude_(deg): Hypocenter latitude in decimal degrees North.

Longitude_(deg): Hypocenter longitude in decimal degrees; negative values are 

West

Elevation_(m): Hypocenter elevation in meters with respect to mean sea level

Magnitude: Earthquake magnitude. 

Magnitude_type: Type of magnitude calculation. Magnitudes < M 3.5 are 

duration magnitudes (MD). Magnitudes ≥ M 3.5 are local magnitudes (ML).

Magnitude_source: Organization that computed the magnitude. Duration 

magnitudes were computed by the Bureau of Reclamation (USBR). Local 

magnitudes were obtained from the U.S. Geological Survey (USGS) or the 

University of Utah Seismograph Stations (UUSS).
A-4

Quality: Indicates whether the hypocenter was computed with a relative location 

method (quality = ‘a’) or an absolute location method (quality=’b’)

RMS_residual_(s): Root-mean-square (rms) time residual in seconds. For bquality events, this is the rms of the absolute time residuals. For a-quality events, 

this is the rms of the time difference residuals. 

Nabstimes: Number of arrival times used to compute the absolute location of bquality events. (Zero for a-quality events.)

Neventpairs: Number of event pairs providing time differences for the relative 

location of a-quality events. (Zero for b-quality events.)

Ntimediffs: Number of time differences used to compute the relative location of

a-quality events. (Zero for b-quality events.)

Nstations: Number of stations providing data for the hypocenter calculations, 

either in terms of absolute arrival times (for b-quality events) or time differences 

(for a-quality events).

Maxgap_(deg): Maximum azimuthal gap in ray coverage to the nearest integer 

degree. For the event relative location (a-quality events), this considers data for all 

event pairs.

Min_dist/depth: Horizontal distance to the closest seismic station providing data 

for the hypocenter calculation divided by the focal depth.

Focal_mechanism_category: Focal mechanism category for the events near the 

January 2013 ML 4.4 earthquake for which focal mechanisms were analyzed, as 

described in section 2.4. Categories are: ‘1’ for strike-slip focal mechanisms, ‘2’ 

for oblique normal focal mechanisms, and ‘0’ for events not analyzed.

References:

Block, L.V., 1991. Joint hypocenter-velocity inversion of local earthquake arrival 

time data in two geothermal regions, Ph.D., Massachusetts Institute of 

Technology, Cambridge, MA.

Block, L.V., C. K. Wood, W. L. Yeck, and V. M. King, 2015, Induced seismicity

constraints on subsurface geological structure, Paradox Valley, Colorado, 

Geophysical Journal International, 200, 1172-1195.