Abstract
The Paradox Valley Unit (PVU), a salinity control project in southwest Colorado, disposes of brine in a single deep injection well. Since the initiation of injection at the PVU in 1991, earthquakes have been repeatedly induced. PVU closely monitors all seismicity in the Paradox Valley region with a dense surface seismic network. A key factor for understanding the seismic hazard from PVU injection is the maximum magnitude earthquake that can be induced. The estimate of maximum magnitude of induced earthquakes is difficult to constrain as, unlike naturally occurring earthquakes, the maximum magnitude of induced earthquakes changes over time and is affected by injection parameters. We investigate temporal variations in maximum magnitudes of induced earthquakes at the PVU using two methods. First, we consider the relationship between the total cumulative injected volume and the history of observed largest earthquakes at the PVU. Second, we explore the relationship between maximum magnitude and the geometry of individual seismicity clusters. Under the assumptions that: (i) elevated pore pressures must be distributed over an entire fault surface to initiate rupture and (ii) the location of induced events delineates volumes of sufficiently high pore-pressure to induce rupture, we calculate the largest allowable vertical penny-shaped faults, and investigate the potential earthquake magnitudes represented by their rupture. Results from both the injection volume and geometrical methods suggest that the PVU has the potential to induce events up to roughly MW 5 in the region directly surrounding the well; however, the largest observed earthquake to date has been about a magnitude unit smaller than this predicted maximum. In the seismicity cluster surrounding the injection well, the maximum potential earthquake size estimated by these methods and the observed maximum magnitudes have remained steady since the mid-2000s. These observations suggest that either these methods overpredict maximum magnitude for this area or that long time delays are required for sufficient pore-pressure diffusion to occur to cause rupture along an entire fault segment. We note that earthquake clusters can initiate and grow rapidly over the course of 1 or 2 yr, thus making it difficult to predict maximum earthquake magnitudes far into the future. The abrupt onset of seismicity with injection indicates that pore-pressure increases near the well have been sufficient to trigger earthquakes under pre-existing tectonic stresses. However, we do not observe remote triggering from large teleseismic earthquakes, which suggests that the stress perturbations generated from those events are too small to trigger rupture, even with the increased pore pressures.
1 INTRODUCTION
1.1 The Paradox Valley Unit
The United States Department of the Interior, Bureau of Reclamation (Reclamation), operates the Paradox Valley Unit (PVU), in southwest Colorado, to reduce the amount of highly saline groundwater entering the Dolores River (Fig. 1). The PVU operates under the Colorado River Basin Salinity Control Project, a multistate collaborative program to lower salinity in the Colorado River. The Dolores River, a tributary of the Colorado River, is a major contributor to the salt load in the Colorado River, naturally picking up 20 000 ktonnes of salt per year as it crosses Paradox Valley (Ake et al.2005). Salt enters the Dolores River from the inflow of highly saline groundwater (∼8 times more saline than sea water). To reduce the inflow of brine, the PVU intercepts shallow brine groundwater at nine extraction wells and disposes of the brine in a single deep injection well at surface pressures up to 35 MPa. PVU conducted injection tests between 1991 and 1995 and subsequently began long-term brine disposal in 1996. To date, roughly 7.7 million cubic metres (m3) of fluid have been injected.
Overview map of Paradox Valley and the current Paradox Valley Seismic Network (PVSN) (yellow triangles), including the epicentres of induced earthquakes (grey circles). The majority of induced earthquakes occur near the injection well, on the southwestern edge of Paradox Valley.
The PVU injection well is immediately southwest of Paradox Valley, a northwest trending valley formed by collapse of a northwest-trending diapiric salt-cored anticline (Cater 1970; Gutiérrez 2004; Trudgill 2011). The injection well penetrates Triassic- through Cambrian-age sedimentary rock layers and granitic Precambrian basement (Fig. 2). With the exception of the Paradox formation, which consists primarily of highly deformed salt layers, the geological units are generally subhorizontal. The basement and overlying sedimentary layers are offset by a series of northwest-trending high-angle normal faults. These deep faults do not extend to the ground surface. Their locations have been mapped using deep seismic reflection and well log data but are only approximately known (King et al.2014). The Mississippian Leadville Formation, a relatively low porosity (<10 per cent) highly fractured carbonate, is the primary injection target formation. Some of the underlying early- to mid-Palaeozoic limestone and sandstone units and the Precambrian basement are considered supplemental reservoirs. The well casing was perforated in several intervals between the top of the Leadville formation (4.3 km depth) and the bottom of the borehole (4.8 km depth). The overlying Paradox salt formation acts as a confining layer.
In conjunction with the PVU injection well, Reclamation operates the Paradox Valley Seismic Network (PVSN), a surface seismic array currently consisting of 20 three-component broad-band seismometers (30-s to 50 Hz bandwidth, 24-bit sampling at 100 Hz), to monitor earthquakes induced by PVU brine injection (Fig. 1). The network has been continuously operational since 1985, 6 yr prior to the start of injection in 1991. The network consisted initially of 10 short-period vertical-component seismometers (1 Hz natural period) sampled using 12 and 16-bit analogue-to-digital converters at 100 Hz; the array was expanded to 16 short-period stations by 1999, including two three-component stations. Since injection began, PVSN has recorded ∼6000 shallow (mostly <6.5 km depth) earthquakes within ∼16 km of the injection well. In contrast, a review of the historical PVSN data files indicates only one local earthquake for the 6-yr pre-injection monitoring period, and it occurred about 19 km from the well. The shallow seismicity observed since the beginning of injection has occurred at increasing distance from the PVU injection well over time (Block et al.2014). Based on the near complete lack of seismicity detected during 6 yr of pre-injection seismic monitoring, the general correlation of the depths of the earthquakes and the depth of injection, and the spatiotemporal evolution of the seismicity, we interpret that most, and possibly all, of the shallow earthquakes recorded since the start of injection were induced by PVU fluid injection. The spatio-temporal evolution of PVU induced earthquakes and the relationship between earthquake locations and geological structures are detailed in recent publications by King et al. (2014) and Block et al. (in press).
About 100 of the events interpreted to be induced have a duration magnitude of Md 2.5 or greater, the approximate threshold for human detection in the Paradox Valley area. The largest events include a local magnitude ML 4.3 (MW 3.8) in 2000, ML 4.1 (MW 3.6) in 2004 and most recently a ML 4.4 (MW 4.0) in 2013 (Block et al.2014). Reclamation has altered injection operations several times to reduce the rate of occurrence of such events. For long-term operational planning it is important to understand the seismic hazard associated with deep-well injection at PVU, especially the maximum magnitude earthquake that could be induced by continued injection. In this paper, we estimate the maximum magnitude earthquake that PVU can induce by considering the relationship between earthquake magnitude, the total volume of brine injected, and the geometry of distinct seismicity clusters.
1.2 Maximum earthquake magnitude
In areas of naturally occurring seismicity, the magnitude-dependent earthquake rate distribution is often considered to be stationary, implying there is no temporal variation in the mean rates of earthquakes of a given magnitude. Previous studies of induced seismicity have used such statistics to evaluate either the probability of inducing earthquakes exceeding a given magnitude or the maximum magnitude (Shapiro et al.2010; Hallo et al.2014). However, for PVU-induced seismicity, both the overall seismicity rates and the Gutenberg–Richter (GR) b-values vary substantially through time and by location (Block & Wood 2010). This adds difficulty to applying GR statistics to estimate the probability of inducing events of a given magnitude. A key parameter for defining the earthquake rate distribution is the maximum magnitude earthquake. Therefore, it is important to understand the variations in maximum earthquake magnitude in time and space to evaluate the seismic hazard posed by the PVU. Here we consider two methods for estimating the maximum magnitude of seismicity induced by the PVU.
The first method to estimate the maximum magnitude of earthquakes induced by fluid injection utilizes the total cumulative volume of fluid injected. Observational studies indicate that the maximum earthquake magnitude is linearly proportional to the logarithm of cumulative volume of fluid injected (McGarr 1976, 2014; McGarr et al.2002; Nicol et al.2011). The seismic and injection history at the PVU are well documented and therefore this relationship is easily explored. We examine the observed maximum magnitude earthquakes as a function of time during PVU fluid injection and compare the observed maximum magnitude values to the cumulative fluid volume injected.
The second method that we consider to estimate maximum magnitude relies on the well-explored relationship between fault rupture size and earthquake magnitude (Brune 1970; Aki 1972; Mark 1977; Hanks & Kanamori 1979; Wyss 1979; Scholz 1982). In the case of tectonic earthquakes, estimating maximum earthquake magnitude from potential fault rupture size often relies on mapped surface faults and earthquake aftershock zones. Induced earthquakes, however, often occur in previously aseismic areas where slip occurs on unmapped subsurface faults, and therefore potential fault rupture size can only be estimated from hypocentre patterns of induced events. In Paradox Valley, there is a rough correlation between the fault rupture size estimated by aftershock lineations (Figs 3 and 4), and earthquake magnitude (Fig. 5).
Earthquake epicentres in regions surrounding the four largest earthquakes induced by the PVU. Earlier events are shown in yellow, aftershocks (within 6 months of the main shock) are shown in blue, and the remaining events are shown in grey. The rupture plane length (red dashed line) is estimated by aftershock locations. In some cases, the relative location uncertainties of the main shocks cause them to appear offset from the assumed ruptured fault plane. Focal mechanisms for the 2013 January and 1999 June events (Reclamation), and the 2000 May and 2004 November events (Saint Louis University Earthquake Center) are shown and agree well with observed aftershock lineations.
Earthquake locations in regions surrounding the next four largest earthquakes induced by the PVU. See Fig. 3 description. Available focal mechanisms are shown (Reclamation).
The relationship of rupture radius and earthquake magnitude for a variety of models including: empirical rupture area-magnitude relationships from Shaw (2009) (Shw09, orange), Ellsworth (2003) (Els03, light green), Wells & Coppersmith (1994) (WC94-RA, light blue), the average of these relationships ±1σ (black dotted lines), the relationship shown in equation (1) with stress drops of 2, 5 and 10 Mpa (blue, green and red dashed lines, respectively), and estimated fault radius using one half of the rupture-length-at-depth parameter from Wells & Coppersmith (1994) (WC94-RLD, purple). Values from earthquakes shown in Figs 2 and 3 are shown as grey or white circles. Duration magnitudes are from Reclamation and local magnitudes are from the NEIC.
Shapiro et al. (2011) observed that at many locations of injection induced earthquakes, including the PVU, large-magnitude induced earthquakes are underrepresented relative to what is predicted by GR statistics. Shapiro therefore concluded that a large portion of a pre-existing fault surface must be stimulated by injection in order to induce an earthquake that ruptures the full fault. We refer to this case as full fault stimulation (FFS). In this case fault slip is mediated by pore-pressure diffusion, as discussed by Garagash & Germanovich (2012). In the alternative case, an earthquake is triggered from stimulating a small portion of a fault plane, and with dynamic rupture then propagating into regions that have not been stimulated by injection. This case predicts a higher rate of large-magnitude events as compared to FFS, and does not appear to be the dominant type of triggering at the PVU. Prior to the onset of injection, Reclamation monitored seismicity for 6 yr and observed only one earthquake in the vicinity of the PVU (Block et al.2014). The near absence of pre-injection seismicity could be due to either a low rate of crustal deformation, a lack of critically stressed faults, or both.
Shapiro et al. (2011) investigated the relationship between the geometric scale of induced seismicity clusters and maximum earthquake magnitudes. They found that the largest observed magnitude is largely controlled by the minimum principal axis (Lmin) of an ellipsoid fitting the seismicity cloud, under the assumption that faults are represented by pre-existing randomly oriented and located circular cracks. In the case of Paradox Valley, multiple distinct seismicity clusters exist with aseismic gaps between them. This implies that effective stresses have been insufficiently altered within the aseismic regions to induce earthquakes. Other factors that may result in aseismic gaps include local heterogeneities in fracture distribution, geology and local stresses. Therefore distinct clusters should be regarded independently when estimating maximum magnitudes.
Ellipsoids may not accurately capture the 3-D geometry of individual seismicity clusters and can overestimate the maximum crack radius within a seismic cloud. In this study, we use a Density-Based Spatial Clustering of Applications with Noise (DBSCAN) algorithm (Ester et al.1996) to determine earthquake clusters and then apply a 3-D convex hull algorithm (Barber et al.1996) to define the shape of each seismic cloud. We estimate the maximum earthquake magnitude from the radius of the largest vertical circular fracture within each cloud rather than from Lmin. We restrict our analysis to vertical fractures because previous focal mechanism analyses of PVU-induced earthquakes indicate that 89 per cent of events studied occur on near-vertical fault planes (Ake et al.2005), which is consistent with the (>M3) largest PVU induced earthquakes for which we have well-constrained focal mechanisms (Figs 3 and 4). We apply this method to seismicity clusters observed at PVU as a function of time and compare the estimated and observed maximum earthquake magnitudes. This geometrical analysis depends on three major assumptions: (1) a critical pore-pressure must affect an entire fault surface to initiate full fault rupture, (2) the location of induced events delineates the volume where sufficient pore-pressure alteration exists to induce full-fault rupture and (3) the fault rupture planes are adequately represented by vertical penny-shaped cracks. We note that if dynamic rupture were to extend outside of the zone of pore-pressure alteration, then maximum magnitude estimations based on geometric constraints would be invalid.
2 DATA
We calculate two sets of hypocentres from the PVSN data. Initially, we compute absolute earthquake locations using manually determined P-wave and S-wave arrival times. The location algorithm uses local 3-D P-wave and S-wave velocity models that we developed from hypocentre-velocity inversions of earthquake and explosion data recorded by PVSN. We subsequently compute precise relative event locations from arrival time differences obtained from time-domain cross-correlations of windowed P-wave and S-wave arrivals extracted from filtered waveforms recorded at the same station for pairs of events. No time differences from manual time picks are included, with the exception of data for 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, they are incorporated into the event relative location using differences of high-quality manually determined arrival times. No event clustering is performed; an event may tie to any other event in the dataset. However, 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 well-constrained absolute locations fixed during the relative location. 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 per cent of earthquakes occurring within 20 km of the PVU injection well are well constrained in the relative event location procedure. 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 m horizontally and 100 m vertically. We use this set of relative hypocentres for the geometric estimated maximum magnitude analysis described below.
Earthquake magnitude estimates are obtained from a variety of sources. Below magnitude 3.5, we use duration magnitude (Md) calculated by Reclamation from PVSN waveform coda. Prior to this work, we recomputed coda durations and duration magnitudes for all PVSN-recorded local earthquakes using a consistent, semi-automated procedure (Block & Wood 2011). Duration magnitudes for pre-2008 events cannot be used for the largest (approximately ML ≥ 3.5) PVU-induced earthquakes due to clipped waveforms. Also, accurate coda durations are not available for many of the PVSN waveforms from these events because of insufficient record length. Therefore, for magnitude 3.5 and above we use local magnitude (ML) reported by the National Earthquake Information Center (NEIC). These include the 1999 ML 3.6, 2000 ML 4.3, 2004 ML 4.1 and 2013 ML 4.4 events. We chose to use ML for larger events as ML has consistently been computed for all large events in our catalogue. Where possible, we also show MW estimates. This includes a MW 4.0 estimate for the 2013 January earthquake (Block et al.2014) as well as a MW 3.8 for the 2000 May earthquake and MW 3.6 for the 2004 November earthquake (Saint Louis University Earthquake Center). Although some of our methodology directly depends on the observed magnitude of earthquakes, and we compare our results to calculated moment magnitudes, slight uncertainties in reported magnitudes should not substantially affect our observations and interpretations.
3 METHODS
3.1 Injection volume versus earthquake magnitude
In order to investigate the relationship between maximum earthquake magnitude and cumulative injected fluid volume, we find each induced earthquake whose magnitude is larger than all previously induced earthquakes and plot it compared to the cumulative injected fluid volume at that time. PVU has maintained a record of daily fluid volume injected since the beginning of the injection tests in 1991. The cumulative fluid volume for the day on which a given earthquake occurred is calculated from these daily totals, and therefore the cumulative fluid volumes used in this analysis are accurate to within a 1-d period (which, at current injection rates, corresponds to about 1090 m3). We calculate the linear least-squares fit between these earthquake magnitudes and the log of the cumulative injected volumes, as well as 95 per cent confidence intervals and 95 per cent prediction intervals.
This analysis does not account for the spatial clustering of earthquakes, but rather utilizes all induced events in the data set. The volume of injected brine affecting individual clusters of seismicity relative to the total volume of injected brine is unknown, and therefore it is not possible to perform this analysis for individual clusters. However, all but the most recent observed maximum magnitude earthquake occurred within 2.2 km of the injection well. In addition, the observed maximum magnitude earthquakes tend to occur at increasing distance from the injection well over time. Hence, the correlation of observed maximum magnitude earthquakes with injected fluid volume is most relevant to the near-well region of induced seismicity and inherently accounts for the injected fluid acting over a spatially expanding region around the injection well over time.
3.2 Geometric maximum magnitude estimations
3.2.1 DBSCAN cluster analysis
Before investigating the geometry of earthquake hypocentres it is necessary to classify individual clusters. As there are large aseismic regions between distinct clusters, implying either a lack of sufficient pore-pressure alteration to induce earthquakes or a lack of fractures favourably oriented for failure in the given stress field, we consider it unrealistic to simultaneously consider all induced events when estimating maximum earthquake magnitudes from hypocentre geometry. In map-view, seismicity is located in visually distinct clusters (Fig. 1). In three dimensions, it is more difficult to visually determine what constitutes an individual cluster. We therefore employ a DBSCAN algorithm to define individual seismicity clusters and remove outliers. The DBSCAN algorithm (Ester et al.1996) relies on point density for determining individual clusters. We chose this clustering algorithm because it is efficient with large datasets, there is no need to predetermine the number of clusters, the algorithm handles irregularly shaped clusters, and it is effective at excluding outliers from clusters. We employed the DBSCAN algorithm from the machine-learning Python module Scikit-learn (Pedregosa et al.2011).
The DBSCAN algorithm requires two inputs when performing cluster analysis: the distance to search for neighbouring points (ε) and the minimum number of points required in order to create an individual cluster (minPts). In the case of Paradox Valley, we find that values of ε = 0.75 km and minPts = 20 replicate well what we would visually classify as clusters. With higher ε and lower minPts (lower required cluster densities), distinct clusters merge while with lower ε and higher minPts, a large number of subclusters form. An example of a DBSCAN selected cluster as compared to surrounding seismicity is seen in Fig. 6(a).
Visual outline of how maximum circular crack radius is calculated for the Nearwell Cluster. (a) First the DBSCAN algorithm clusters earthquakes. Earthquakes included in the cluster appear in blue while outliers appear in black. (b) A convex-hull is fit around the cluster data points. (c) The maximum circular crack is calculated for a 2-D slice along a plane in the convex hull (in this example the X–Z plane).
3.2.2 Convex-hull and maximum circular crack radius analysis
After determining individual seismicity clusters using the DBSCAN algorithm, we next determine the largest vertical circular fault plane that can be contained within the cluster. In order to perform this analysis, we first compute the 3-D convex-hull of the seismicity cluster. A convex-hull is the minimum convex shaped envelope which contains the full set of points. We use the computer program Qhull (Barber et al.1996) to compute convex hulls. The convex hull of an example seismicity cluster is illustrated in Fig. 6(b).
For each convex hull, we find the largest 2-D vertical circular crack enclosed within the hull. In order to compute this we first find the intersection between the hull and a vertical plane, the result of which is a 2-D polygon (Fig. 6c). We use a Monte Carlo approach to find the maximum radius circle enclosed by the 2-D polygon. We first select random points within the polygon. For each point, we calculate the distance to the closest edge of the bonding polygon. The centre point is the point where this distance is largest. An example of the largest circular crack bound in a 2-D slice through the polygon is shown in Fig. 6(c) . We repeat this process for a subset of all possible vertical planes intersecting the cluster, incrementing the strike of the plane by ∼2° and translating the slicing plane through the cluster in 100 m increments in order to find the maximum radius of a circular crack.
In order to track the temporal variation of the clusters, we perform this evaluation for all events up to a given time. The DBSCAN algorithm is only performed on the full data set to ensure that the set of events considered in a given cluster is consistent in time, while the convex hull analysis is performed in year increments. In order to perform this analysis, >3 hypocentres must exist to create a 3-D convex hull.
3.2.3 Magnitude–rupture size relationship
This relation indicates that magnitude scales with the logarithm of rupture area, for the case of constant stress drop. This is consistent with empirical magnitude versus rupture area relations for natural earthquakes (Wells & Coppersmith 1994; Hanks & Bakun 2002; Ellsworth 2003; Shaw 2009). With eq. (1) we estimate earthquake magnitude for a given rupture size and static stress drop. Using PVSN data, the estimated static stress drop for the 2013 January ML 4.4 earthquake was 2 MPa, and its estimated radius was 0.6 km (Block et al.2014). 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 et al. (2011) found an increase in stress drop with radial distance from the injection site and that stress drop correlates with the pore-pressure perturbations due to injection. We therefore calculate earthquake magnitudes for static stress drops of 2, 5 and 10 MPa (Fig. 5).
4 RESULTS
4.1 Injected volume versus maximum magnitude
We find good correlation (coefficient of determination, R2 = 0.92) between maximum earthquake magnitude and the cumulative volume of injected fluid over time (Fig. 7). Each point on Fig. 7(b) represents the occurrence of an earthquake having a magnitude larger than all previously induced earthquakes. All of these earthquakes occurred within 2.2 km of the injection well except for the most recent event, which occurred 8.2 km northwest of the well. 95 per cent confidence intervals give a maximum earthquake magnitude range at the current cumulative injected volume of M 3.8–4.7, while the 95 per cent prediction interval ranges from M 3.3 to 5.2. Quantile–Quantile (QQ) plots of the residuals of the least-squares fit 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.
(a) Flow rate as a function of the cumulative injected volume. The grey line shows the flow rate during the injection tests (1991–1995), while the black line shows the flow rate during long-term injection (1996–2013). (b) Observed maximum magnitude PVU-induced earthquakes as a function of the cumulative injected volume. Least-squares fit shown in black, with 95 per cent confidence interval (black dashed) and 95 per cent prediction interval (grey dashed).
Because maximum earthquake magnitudes correlate with the log of cumulative injected volume, earthquakes during the early injection history of the PVU have a disproportionately large influence on the least-squares fit. Nine of the observed maximum magnitude events, covering more than two-thirds of the correlation plot, occurred during the 5 yr of injection testing (1991–1995), although the fluid injected during this time period accounts for less than 10 per cent of the total fluid volume injected to date. In contrast, the 17 yr of long-term injection (1996–2013), during which more than 90 per cent of the fluid was injected, are represented by only five observed maximum magnitude earthquakes. The rate of magnitude increase was much larger early in the injection history than it is currently. Maximum earthquake magnitudes increased from 0.7 to 2.8 during the 5 yr of injection testing, an increase of 2.1 magnitude units. During the first 4 yr of long-term injection (1996–2000), observed maximum magnitudes increased an additional 1.5 magnitude units. During the last 13 yr (2000–2013), the maximum earthquake magnitude increased by only 0.1 unit. Hence, the correlation shown in Fig. 7 suggests that the growth rate of maximum earthquake magnitude is currently small and is expected to decline further over time. For example, given the best fit line, to raise the estimated maximum magnitude from 4.4 to 5.0, an additional 1.79 × 107 m3 of injected fluid would be required, which is equivalent to a time period of ∼45 yr at the current rate of injection.
The steady increase in earthquake magnitude over time indicates that during PVU fluid injection, subsurface conditions are changing such that larger magnitude induced earthquakes become more likely as injection progresses. This is different from naturally occurring earthquakes, where larger-magnitude events occur less frequently than the smaller-magnitude events, but the magnitude of the largest earthquake is generally assumed to be time invariant. At PVU we observe a time-dependent magnitude–frequency relationship.
4.2 Geometric maximum magnitude results
4.2.1 Earthquake clustering
The DBSCAN algorithm found five distinct earthquake clusters, which agrees well with what we would visually interpret (Fig. 8). The largest cluster occurs in the area directly surrounding the well on the southern side of the valley (Nearwell Cluster). All of the observed maximum magnitude earthquakes shown in Fig. 7 occur within this cluster with the exception of the January 2013 ML 4.4 event, which occurred in a large cluster roughly 7.5 km northwest of the injection well (Northwest Cluster). The Northwest Cluster is the second largest cluster. Across the valley from the injection well, two clusters are identified. The largest cluster (Across Valley #1) occurs nearly 12 km north of the injection well. It is accompanied by a smaller cluster about 5.5 km directly to the west (Across Valley #2 Cluster). The smallest identified cluster occurs directly SE of the injection well (Southeast Cluster), and forms a tight lineation.
Identified earthquake clusters from the DBSCAN algorithm. Excluded events (black dots) are events with poor relative location quality.
Shallow, likely induced earthquakes not belonging to clusters occur sparsely throughout the region (Fig. 8). In addition, a few microclusters were not identified by the DBSCAN algorithm. The isolated earthquakes and microclusters, under the assumptions in this analysis, could not produce large magnitude earthquakes. With the exception of a Md 2.9 earthquake near the Across Valley #2 cluster, earthquakes with Md larger than 2 have not been observed to date in these areas.
The estimated maximum magnitude is controlled by the size of the earthquake cluster and therefore, the parameters we use for the DBSCAN algorithm. We investigated the effect of more inclusive clustering. By reducing the density constraints of the DBSCAN algorithm (with ε = 1.2 km and minPts = 25), we incorporated the largest microcluster, located immediately southwest of the Nearwell Cluster, as well as many of the surrounding events (outliers) into the Nearwell Cluster. Although the change in cluster geometry in plain view was large, our estimated maximum magnitudes for the Nearwell Cluster increased by less than 0.5.
In this analysis we do not report results from the clusters north of Paradox Valley (Across Valley #1 and Across Valley #2). Many of the earthquakes within these clusters, most notably those that occurred prior to the installation of three-component seismic stations in the northern valley area in 2007, have relatively large hypocentre uncertainties, making maximum magnitude estimations using our methodology less reliable. Also, it is unclear what physical mechanisms induce these distant events. Thus far, our analyses indicate that the maximum magnitude from these clusters is small relative to that of the near well region.
4.2.2 Maximum magnitudes
From our geometric constraints, we find that the largest vertical crack possible is in the Nearwell Cluster with a radius of 1.49 km and strike of 99° (Table 1, Fig. 9). Using equation 1, we calculate magnitude from the rupture-radius of this crack for stress drops of 2, 5 and 10 MPa as MW 4.8, 5.0 and 5.2, respectively. The Northwest Cluster supports the second largest crack size, with a radius of 0.69 km and calculated magnitudes for 2, 5 and 10 MPa stress drops of MW 4.1, 4.3 and 4.5, respectively. The Southeast Cluster facilitates a small crack radius (0.06 km), with maximum magnitude estimates ranging from 1.9 to 2.3 for stress drops between 2 and 10 MPa, which is generally near or below the threshold of human detection.
Estimated maximum earthquake magnitude given a 5 MPa stress drop as a function of strike for the Nearwell, Northwest, and Southeast Clusters. Maximum values are shown as white stars. The observed primary and secondary fault plane orientations from Ake et al. (2005) are marked as vertical dashed lines.
Radius and strike of largest vertical crack in each cluster, and computed earthquake magnitudes for given stress drops.
| Cluster | Maximum vertical crack radius (km) | Strike of crack (°) | Magnitude for Δ σ = 2 MPa | Magnitude for Δσ = 5 MPa | Magnitude for Δσ = 10 MPa | Observed maximum magnitude |
|---|---|---|---|---|---|---|
| Nearwell | 1.49 | |$\phantom{0}$|99 | 4.8 | 5.0 | 5.2 | 4.3 ML (3.8 MW) |
| Northwest | 0.69 | 143 | 4.1 | 4.3 | 4.5 | 4.4 ML (4.0 MW) |
| Southeast | 0.06 | |$\phantom{0}$|85 | 1.9 | 2.2 | 2.3 | 1.9 Md |
| Cluster | Maximum vertical crack radius (km) | Strike of crack (°) | Magnitude for Δ σ = 2 MPa | Magnitude for Δσ = 5 MPa | Magnitude for Δσ = 10 MPa | Observed maximum magnitude |
|---|---|---|---|---|---|---|
| Nearwell | 1.49 | |$\phantom{0}$|99 | 4.8 | 5.0 | 5.2 | 4.3 ML (3.8 MW) |
| Northwest | 0.69 | 143 | 4.1 | 4.3 | 4.5 | 4.4 ML (4.0 MW) |
| Southeast | 0.06 | |$\phantom{0}$|85 | 1.9 | 2.2 | 2.3 | 1.9 Md |
Radius and strike of largest vertical crack in each cluster, and computed earthquake magnitudes for given stress drops.
| Cluster | Maximum vertical crack radius (km) | Strike of crack (°) | Magnitude for Δ σ = 2 MPa | Magnitude for Δσ = 5 MPa | Magnitude for Δσ = 10 MPa | Observed maximum magnitude |
|---|---|---|---|---|---|---|
| Nearwell | 1.49 | |$\phantom{0}$|99 | 4.8 | 5.0 | 5.2 | 4.3 ML (3.8 MW) |
| Northwest | 0.69 | 143 | 4.1 | 4.3 | 4.5 | 4.4 ML (4.0 MW) |
| Southeast | 0.06 | |$\phantom{0}$|85 | 1.9 | 2.2 | 2.3 | 1.9 Md |
| Cluster | Maximum vertical crack radius (km) | Strike of crack (°) | Magnitude for Δ σ = 2 MPa | Magnitude for Δσ = 5 MPa | Magnitude for Δσ = 10 MPa | Observed maximum magnitude |
|---|---|---|---|---|---|---|
| Nearwell | 1.49 | |$\phantom{0}$|99 | 4.8 | 5.0 | 5.2 | 4.3 ML (3.8 MW) |
| Northwest | 0.69 | 143 | 4.1 | 4.3 | 4.5 | 4.4 ML (4.0 MW) |
| Southeast | 0.06 | |$\phantom{0}$|85 | 1.9 | 2.2 | 2.3 | 1.9 Md |
The largest observed earthquake in the Nearwell Cluster is a ML 4.3, or MW 3.8, event that occurred 2000 May 27 (Fig. 10). The largest predicted earthquake based on the size of the seismicity cluster for a 5 MPa stress drop is MW 5.0. There is little ellipticity in the convex hull of the fault plane for this earthquake, and therefore the estimated maximum magnitude only varies slightly with strike. These estimated maximum magnitudes are significantly larger than the maximum earthquake magnitude observed to date (MW 3.8). The vertical crack size is partly constrained by the vertical extent of the seismicity cluster. The largest vertical crack strikes at 99°. This agrees well with the preferred fault strike observed in focal mechanisms of 86° (Ake et al.2005).
Maximum crack radius and maximum earthquake magnitude as a function of time for the (a) Nearwell, (b) Northwest, and (c) Southeast Clusters. Observed maximum magnitudes through time are shown as crosses and stars.
When the estimated maximum magnitude for the Nearwell Cluster is computed as a function of time, the temporal variations of predicted and observed maximum magnitudes match well (Fig. 10). The size of the Nearwell Cluster increased quickly over the injection test period (1991–1995), and the observed and estimated maximum magnitudes also increase substantially during the period. The seismicity cluster's growth stalled from 1995 to 1997, matching a time when no maximum magnitude earthquakes were observed. From 1997 to about 2000 the seismicity cluster again increased in size, and maximum magnitude earthquakes were again observed. Since 2000, the seismicity cluster has grown very slowly in size, and therefore the estimated maximum magnitude has not increased substantially. During this time, no maximum magnitude earthquakes have been observed in the Nearwell Cluster.
In the case of the Northwest Cluster, predicted maximum magnitudes (MW 4.3, for 5 MPa stress drop; Table 1 and Fig. 9) match well with the 2013 January 24 ML 4.4 (MW 4.0) observed maximum magnitude earthquake. The crack orientation of the largest predicted event (143°) is misaligned from the strike of the 2013 January event fault plane (78°; Block et al.2014) and from the Ake et al. (2005) most commonly observed fault plane orientation for the Nearwell Cluster (86°). Along the strike of the 2013 earthquake, the predicted maximum magnitude is MW 4.1 given a 5 MPa stress drop and MW 3.8 given a 2 MPa stress drop, the stress drop estimated for this event (Block et al.2014). The time history of the maximum predicted magnitude is similar to that observed in the Nearwell Cluster (Fig. 10), and matches the observed maximum magnitude earthquakes. After seismicity in the cluster initiated in 1997, the largest potential vertical crack radius grew to over 80 per cent of the current maximum vertical crack radius by 1998. From 1998 onward, the seismicity cluster has slowly increased in size.
In the Southeast Cluster, the computed magnitudes are small, generally below the threshold of human detection. The largest earthquake observed in this cluster is an Md 1.9, which agrees well with the computed maximum magnitude of MW 2.2 for a stress drop of 5 MPa. The time history is poorly constrained on this cluster due to the lack of events early in the cluster's history, making it difficult to compute the convex hull. The estimated maximum magnitude occurs on a fault plane oriented at 85°. This agrees extremely well with preliminary focal mechanisms for larger events in the cluster as well as with the primary mechanism observed by Ake et al. (2005) for the Nearwell Cluster, likely due to the tight lineation of the cluster.
5 DISCUSSION
There is general agreement between maximum magnitude estimations based on cumulative injected volume and those based on the geometry of the seismic cloud. The 95 per cent prediction interval from the cumulative volume relationship and the maximum vertical crack size of any orientation in the Nearwell Cluster (for a 5 MPa stress drop) place upper bounds on the maximum magnitude of M5.2 and M5.0, respectively. Maximum magnitudes predicted by the geometric method are smaller for the other clusters studied. Under the assumptions of this study, it is therefore reasonable to assume that the upper bound magnitude that PVU could likely induce, based on injection that has occurred to date, is in the low magnitude 5 range.
In the Nearwell Cluster, the geometric model overestimates the maximum magnitude observed to date. One reason for this may be that a circular crack fault model may not accurately characterize the seismically active faults. The vertical extent of the faults may be constrained by the local geology. Fault height may be constrained by the overlying salt formation and underlying Precambrian units. Fault rupture length may also be limited by the size of the blocks formed by the northwest-trending normal faults, especially since one of the most commonly observed fault plane orientations from focal mechanism analyses (86°; Ake et al.2005) is oblique to the northwesterly strike of these mapped basement faults.
Slow pore-pressure diffusion may prevent the maximum magnitude earthquake from occurring until many years after the size of the cluster has grown sufficiently to incorporate the fault rupture plane. A cluster may be seismically active for many years before pore-pressure increases enough along an extensive fault surface to cause rupture of the entire segment. We have recently reported a 13-yr delay in rupture of a 1.4-km-long fault surface in the Northwest Cluster (Block et al.2014). Although smaller events had occurred along all portions of this fault segment since 2000, the entire fault segment did not rupture in a single event until 2013 January, producing the latest observed maximum magnitude event. This maximum magnitude event (MW 4.0) correlates very well with the maximum magnitude estimate from the geometric method, which has remained fairly constant at a value of about MW 4 (for 5 MPa stress drop) since 2000 (Fig. 10b). This observation for the Northwest Cluster suggests that, although the geometric model currently overestimates observed maximum earthquake magnitude in the Nearwell Cluster, perhaps the predicted maximum magnitude event will occur after sufficient time passes.
Temporal variations in both the observed maximum magnitude earthquakes and those predicted from the fluid volume and geometric models show that the greatest increase in maximum earthquake magnitude occurred early during PVU injection operations, especially during the injection test period (1991–1995). The rate of increase of maximum earthquake magnitude subsequently declined and has been very small since 2000. This trend is likely due to the fact that the calculated magnitude scales with the log (R2) or log(volume), and the rate of increase of R or total volume decreases with an increase in the volume of injectate assuming simple geometries. In addition, changes in injection operations that have been made to reduce the occurrence of larger-magnitude events could partly explain the reduced rate of increase in observed maximum magnitude. Beginning in mid-1999, after the occurrence of two ML 3.5–ML 3.6 earthquakes, PVU implemented two 20-d shut-down periods each year to provide time for near-well pore pressures to decline. In mid-2000, after the occurrence of a ML 4.3 earthquake, the injection flow rate was decreased by one-third. Most recently, after the occurrence of a ML 4.4 earthquake in 2013 January, the effective injection flow rate was decreased by 8 per cent and weekly shut-downs were implemented in place of the bi-annual shutdowns (based on results from pressure-flow modelling).
We observe that the rate of growth of a given cluster can be extremely high during its initiation (Fig. 10). This has the implication that if a new cluster forms in a previously aseismic location, it has the potential to produce felt earthquakes in a relatively short period of time (1–2 yr), according to the geometric maximum magnitude model. Observations at PVU substantiate this prediction, as earthquakes greater than Md 2.5 (the approximate level for human detection in the Paradox Valley area) occurred in the Nearwell Cluster during the injection test period and earthquakes of this magnitude occurred in the Northwest Cluster only about 1 yr after seismicity in the cluster initiated (Fig. 10). The alternative scenario has also been observed, however. Both seismicity rates and maximum earthquake magnitudes were relatively low for more than 5 yr in the Southeast Cluster prior to the observed rates and magnitudes increasing. Differences in the growth rates of individual seismicity clusters may be related to the earthquake triggering mechanics (relative effects of pore-pressure increase versus static stress change, for example).
If the aseismic region between two clusters becomes seismically active, such that the clusters merge, the maximum magnitude earthquake computed with the geometric method can quickly increase. At PVU, however, the earthquake clusters closest to the injection well (Northwest, Nearwell and Southeast Clusters) are misaligned as compared to the apparent strike of the faults that generate larger earthquakes (Figs 3 and 4; Ake et al.2005). Therefore, if one or more of these three clusters merge, the size of fractures optimally oriented for failure would not be expected to change substantially and the observed maximum magnitude earthquake may not greatly increase.
An important limitation of both the injection volume and geometric methods for estimating maximum magnitude is that they are cumulative, and therefore they cannot account for any change in subsurface conditions that would cause the maximum magnitude potential to decrease. During long-term PVU fluid injection, we have observed two distinct time periods during which maximum earthquake magnitude temporarily declined in the near-well area (within 5 km of the injection well). No earthquakes larger than Md 2.5 (Md 2.5+) were observed in the near-well area for more than 2.5 yr, beginning a few months after a 33 per cent decrease in injection flow rate in mid-2000. In contrast, during the 2.5 yr prior to this operational change, more than 35 Md 2.5+ earthquakes were observed, including events up to ML 4.3. Similarly, no Md 2.0+ near-well earthquakes were recorded during a 2.5-yr period beginning a few months after an extended (2+ month) shut-down of injection operations in late 2005 to early 2006. We have correlated these periods of decreased maximum earthquake magnitude with decreased long-term average injection pressures and infer that the lack of larger-magnitude events during these time periods is due to a decrease in pore pressures (Block & Wood 2010; Block et al.2009, 2014). Since the most recent operational changes following the 2013 January ML 4.4 earthquake, maximum injection pressures have been approximately 3 MPa lower than prior to the event. The corresponding decrease in pore pressures likely means that the potential maximum magnitude earthquake for the near-well area of induced seismicity is currently less than indicated by the cumulative models presented here.
Triggering of injection related earthquakes from the transient stresses generated by the seismic waves of remote earthquakes has been linked to the presence of critically loaded faults (van der Elst et al.2013). In Paradox Valley, remote triggering is not apparent (Fig. 11), suggesting that the small stress perturbations from remote earthquakes are insufficient to initiate slip. The lack of pre-injection seismicity combined with the apparent lack of remote triggering suggests that faults are not close to failure in the absence of an increase in pore-pressure over a sufficiently large volume. Dynamic rupture therefore appears to be unlikely outside regions of pore-pressure alteration. If it were the case that dynamic rupture propagates into regions where the pore-pressure has not been altered, our geometric analysis would fail to capture the true maximum magnitude. This case appears unlikely, but it is impossible to completely rule out. The magnitude of such an earthquake would not be limited by injected fluid volume or size of induced seismicity clusters, but instead by the entire pre-existing fault size and the local tectonic stresses.
Histogram of Paradox Valley earthquakes for 10 d prior and after the six largest teleseismic events during PVU operations, following van der Elst et al. (2013). Bottom plot shows sum of all histograms. There is no discernible trend of remote triggering of induced earthquakes at PVU.
Zhang et al. (2013) observed that in many recent cases of induced seismicity, earthquakes occur within crystalline basement. The authors suggest that the stress drops from earthquakes in crystalline basement is often larger than in shallower rock layers. This would suggest that earthquake magnitudes might be larger in crystalline basement relative to those in overlying formations. In the case of Paradox Valley, an early flow profile found that 22 per cent of the injected brine entered the lower perforations, one of which perforates the Precambrian basement (Envirocorp Services and Technology Inc. 1995). It is likely that this percentage has decreased over time due to the fact that the well has slowly filled with sediment, burying the deeper perforations. As of 2001, fill covered all but the top perforation into the upper Leadville. It is unknown what amount of brine is currently entering the Precambrian basement. Within 500 m of the well, assuming the geology matches that at the well; only 4 out of 705 well-constrained earthquakes appear to be in the Precambrian basement. This may suggest that brine is largely not entering the basement at the well. It remains possible that pore-pressure has been significantly altered in the basement away from the well. Due to the complexity of the subsurface structure and the lack of detailed subsurface imaging, the exact amount of induced earthquakes that occur in the Precambrian basement is unknown.
6 CONCLUSIONS
Results from two types of maximum earthquake magnitude analyses indicate that the largest earthquake that could likely be induced by PVU fluid injection, based on injection that has occurred to date, would have a magnitude of ∼M 5 and occur in the near-well region of induced seismicity (within 5 km of the injection well). These analyses are based on the cumulative volume of fluid injected to date and the current size of induced seismicity clusters. These analyses may overestimate the current seismic hazard because the effects are cumulative and do not take into account recent changes that have been made in injection operations to reduce pore pressures, or potential constraints on maximum fault plane rupture size related to the geological structure.
Both observed and predicted maximum earthquake magnitudes indicate that the greatest increase in maximum earthquake magnitude at PVU occurred during the injection test period (1991–1995). The growth rate of maximum earthquake magnitude subsequently declined and has been extremely low since 2000. The maximum earthquake magnitude projected from the cumulative injected fluid volume versus earthquake magnitude relation predicts that the maximum magnitude potential would only increase by 0.6 magnitude units from another ∼45 yr of fluid injection, at current (2014) flow rates.
Maximum magnitudes calculated from the size of seismicity clusters in general overpredict earthquake magnitudes as compared to those observed. The cause of this overprediction is unclear but could be attributed to the fact that an insufficient time has passed for maximum magnitude events to occur, the seismicity cluster volume is heterogeneous and not all locations have sufficient pore-pressure alterations to rupture, or that the magnitude-rupture area relations used in this study are inappropriate (possibly due to geological constraints on fault plane rupture geometries).
An important observation when considering seismic hazard is the fact that seismicity clusters can quickly grow in size. Therefore, there is the potential that large clusters could quickly develop in currently aseismic regions and produce felt events. These analyses are not predictive of such effects, which likely can only be observed through continuous seismic monitoring.
When comparing the two methodologies utilized in this paper, we believe that both are valid and produce plausible estimations of maximum magnitudes at Paradox Valley. Each methodology may be better suited in specific situations. A benefit of maximum magnitude estimations based on the geometry of seismic clusters is the fact that they account for spatial complexities in earthquakes locations, including aseismic gaps between seismic clusters. The method is hindered by the fact that it requires detailed, well-constrained earthquake hypocentres and assumptions about fault geometries (e.g. penny-shaped cracks) and stress drops. Maximum magnitude estimations based on injection parameters benefit from their simplicity and may be more readably applied to areas of induced seismicity that lack dense seismic monitoring and therefore precise earthquake hypocentres. This methodology could be especially pertinent in new areas of induced seismicity. Injection volumes are commonly available from waste-water injection wells and observed maximum magnitudes are more easily available then precise hypocentres. As well, the method relies on maximum magnitude observations and therefore does not require assumptions about fault geometry or stress drops. A disadvantage of the method is the fact that it cannot account for observed spatial complexities. Both methods are hindered by the fact that they are cumulative and thus predicted maximum magnitudes can never decrease. In reality, the combination of these methods is the most informative as both geological structure and injection volume play a role in controlling maximum magnitudes.
This paper summarizes studies done in support of the PVU. We thank Andy Nicholas, Tyler Artichoker and Ed Warner for their support of PVU seismic monitoring and analysis. We thank Mark Meremonte, Glenda Besana-Ostman, David Copeland, Michael Gilliam and Scott Cartwright for their maintenance of PVSN stations and data acquisition systems. We thank Art McGarr and an anonymous reviewer of this paper for their helpful and constructive comments. The reports prepared under contract to Reclamation cited in this paper, as well as annual PVSN seismicity reports, are available at the PVU project website: http://www.usbr.gov/uc/wcao/progact/paradox.
