Research ArticleGEOPHYSICS

Induced seismicity provides insight into why earthquake ruptures stop

See allHide authors and affiliations

Science Advances  20 Dec 2017:
Vol. 3, no. 12, eaap7528
DOI: 10.1126/sciadv.aap7528

Abstract

Injection-induced earthquakes pose a serious seismic hazard but also offer an opportunity to gain insight into earthquake physics. Currently used models relating the maximum magnitude of injection-induced earthquakes to injection parameters do not incorporate rupture physics. We develop theoretical estimates, validated by simulations, of the size of ruptures induced by localized pore-pressure perturbations and propagating on prestressed faults. Our model accounts for ruptures growing beyond the perturbed area and distinguishes self-arrested from runaway ruptures. We develop a theoretical scaling relation between the largest magnitude of self-arrested earthquakes and the injected volume and find it consistent with observed maximum magnitudes of injection-induced earthquakes over a broad range of injected volumes, suggesting that, although runaway ruptures are possible, most injection-induced events so far have been self-arrested ruptures.

INTRODUCTION

The hazard posed by induced and triggered seismicity is a growing concern, notably in the geoenergy industry, both fossil and renewable (1, 2). If viewed as large-scale experiments, anthropogenic earthquake sequences may also provide opportunities to advance our understanding of earthquake processes. The magnitude of the largest expected earthquake, Mmax, is a crucial parameter for seismic hazard analysis. For natural earthquakes, it is estimated from a combination of paleoseismological data, recorded or historical seismicity, fault geometry, and empirical scaling relations (3, 4). The use of dynamic rupture modeling for the assessment of extreme properties of earthquakes is still emergent (5, 6)—notably due to large uncertainties on initial conditions, fault geometry, and rheology—and has been used only in a few studies of injection-induced earthquakes (7). Although the relation between Mmax and fluid injection or production parameters is a topic of active research (813), previously proposed models to estimate Mmax lack a relation between earthquake size and the triggering stress perturbation, consistent with physics-based earthquake models.

Although it has been demonstrated that earthquakes can be initiated by stress perturbations caused by fluid injection and production (1), the conditions controlling the arrest (and hence final size) of induced earthquakes are still poorly understood. A working assumption of several existing models is that induced ruptures remain confined to the volume perturbed by fluid pressure (8, 9) or poroelastic stresses (14). However, in principle, an earthquake nucleated in an area of enhanced stress can continue propagating into the surrounding areas if sufficient stored elastic energy is available there. These earthquakes are often qualified as “triggered” rather than “induced” (15). The concept of a “critically stressed crust” (16), supported by geophysical and geodetic observations, implies that even tectonically stable regions (that is, intraplate regions with low strain rate) are capable of storing significant elastic energy. Statistical models of induced seismicity (12) allow for these triggered earthquakes but do not specify their physics.

Numerical simulations and fracture mechanics provide fundamental insights into the nucleation of dynamic ruptures by localized stresses and their subsequent propagation and arrest beyond the nucleation area. Galis et al. (17) found that the capability of a rupture to propagate along a homogeneous fault depends on the amplitude and area of the stress perturbation. Small or weak perturbations result in arrested ruptures that stop spontaneously at a finite distance from the nucleation area. Large or strong stress perturbations result in runaway ruptures that are only stopped by the boundaries of the fault. These two rupture behaviors are also observed in laboratory experiments of frictional sliding of polymethyl methacrylate blocks driven by a lateral point load (18). At low imposed force, arrested ruptures propagate along the block interface but stop spontaneously before reaching the end of the block. At large enough shear force, a runaway rupture leads to slip of the whole block. Kammer et al. (19) applied fracture mechanics to predict the length of the arrested ruptures as a function of loading force and the transition to runaway ruptures. Garagash and Germanovich (20) identified the two modes of ruptures in their analysis of nucleation and arrest of ruptures due to a line-source injection, resulting in a one-dimensional (1D) pore-pressure perturbation. These theoretical and experimental results suggest that also earthquake ruptures may propagate in the two regimes.

Determining whether earthquakes propagate as arrested or as runaway ruptures would provide physics-based constrains on earthquake size. Although the nucleation conditions of natural earthquakes are poorly known, they may be better constrained for earthquakes initiated by pore-pressure perturbations during fluid injection or withdrawal. Therefore, to gain deeper physical insight into the propagation and arrest of earthquake ruptures, we focus on studying injection-induced earthquakes. Building up on our previous works (17, 19, 21, 22), here, we develop a quantitative model for the effect of pre-existing elastic energy on the size of induced earthquakes based on numerical simulation and fracture mechanics theory of the nucleation and arrest of dynamic rupture. Our model of rupture arrest is qualitatively similar to that of Garagash and Germanovich (20), despite the fact that they analyzed a 1D case, whereas we study a 2D pore-pressure perturbation.

Through our investigations, we show that the observed relation between maximum magnitudes of induced earthquakes, Mmax, and net volume of injected fluid, ΔV, is consistent with predictions of our model. Our work is therefore a first-order step toward the integration of earthquake physics into a theoretical and computational framework for modeling injection-induced and triggered seismicity.

RESULTS

Rupture arrest in 3D dynamic simulations of earthquakes induced by fluid injection

We consider earthquake ruptures nucleated by a stress perturbation localized on a small portion of a fault. To introduce the model through a basic example, we consider a horizontal planar fault intersecting a cylindrical reservoir (Fig. 1). Fluid is injected at a constant rate at the center of the reservoir. Before injection, the fault is loaded by large-scale tectonic stresses, resulting in uniform background shear stress τ0 and effective normal stress σ0′ (positive if compressive). Upon injection, the pore-pressure perturbation p varies in space and time, as given by an axisymmetric analytical solution of the diffusion equation for single-phase (that is, water) fluid flow (see Materials and Methods). The effective normal stress is σ′ = σ0′ − p. Fault slip is governed by linear slip-weakening friction with static and dynamic friction coefficients μs and μd, respectively, and characteristic slip-weakening distance, Dc. These frictional properties are uniform over the fault. The background stress is characterized by a nondimensional strength excess ratio S = (μsσ0′ − τ0)/(τ0 − μdσ0′). The “potential stress drop,” that is, the shear stress in excess of dynamic strength available to drive slip, is Δτ = τ0 − μdσ′ = Δτ0 + μdp, where Δτ0 = τ0 − μdσ0 is the uniform background stress drop.

Fig. 1 Arrested and runaway earthquake ruptures.

The red line and the symbols delineate the rupture arrest area as a function of pore pressure inside the perturbation. The rupture extent is also illustrated by final slip maps. The sketch in the inset depicts the relative position of a fault and a reservoir.

Fluid injection into the reservoir increases its pore pressure, which reduces the fault’s frictional strength μs0′ − p). A rupture nucleates if and where the fault strength drops below the shear stress. The ensuing rupture propagation is modeled with a 3D dynamic rupture simulation method (17). The typical behavior resulting from these simulations is illustrated in Fig. 1. In this example, the pore-pressure perturbation on the fault has a fixed area but increasing amplitude. To estimate an upper bound on rupture size at a given pressure, we assume that the fault is locked until that pressure is reached and then slip is released in a single event. If the pressure is low, ruptures are self-arrested. Their arrest area is significantly larger than the fluid-pressurized area. As the pressure increases, so does the rupture arrest area. When a critical pressure is reached, there is a sharp transition to runaway rupture: The rupture area jumps to a size limited only by the assumed fault boundaries. More generally, we find that whether the rupture becomes arrested or runaway depends on both the area and amplitude of the pore-pressure perturbation.

Fracture-mechanics estimates of the size of arrested ruptures

Although it is possible to conduct a parametric study of the arrest of induced earthquakes based on 3D dynamic rupture simulations, here, we develop a more efficient approach based on fracture mechanics. We validate and calibrate the approach through 3D dynamic simulations and then apply it to assess the influence of reservoir and fault parameters on the size of arrested ruptures.

To estimate the area of arrested ruptures, we apply Griffith’s crack equilibrium criterion. The criterion is essentially an energy balance equating the elastostatic energy release rate due to crack growth, G, and the fracture energy dissipated by unit of crack growth, Gc = 1/2σ′(μs − μd)Dc. To keep the derivation tractable, we approximate the rupture as circular and ignore the differences between shear crack modes (II and III). We apply known integral relations to numerically compute G(R) as a function of rupture radius R for any axisymmetric distribution of stress (21). We determine the equilibrium radii R that satisfy G(R) = Gc and apply an additional condition, d(GGc)/dR < 0, to identify stable equilibrium position. Finally, we use a circular crack model with uniform stress drop to compute the moment magnitude Mw from the rupture area (the derivation of the approach and its validation are presented in Materials and Methods).

Considering now an injection at a constant flow rate into a sealed reservoir, we investigate the dependency of rupture arrest area on reservoir and fault parameters. In Fig. 2A, we compare the temporal evolution of rupture arrest area for a fixed set of fault parameters (with S = 2) and for a representative range of reservoir parameters summarized in table S1. We find that the detailed distribution of pore pressure inside the reservoir is largely irrelevant: Rupture arrest area is controlled by the integral over the spatial extent of pore pressure on the fault. Consequently, reservoir parameters other than its radius do not affect the sizes of arrested ruptures, only their admissible timing. The injection rate q, fluid compressibility ct, porosity φ, and reservoir thickness h have similar effects. At times long enough so that reservoir pressure is saturated, these parameters affect the value and temporal evolution of the average pressure but not the shape of the pressure distribution (see Fig. 3). For example, reducing injection rate q delays the onset of arrested ruptures, whereas increased injection rate leads to earlier onset, yet in both cases, the minimal and maximal rupture arrest areas are identical. On the other hand, permeability k and viscosity μv alter the pressure distribution but not its integral; hence, they do not affect rupture timing. Reservoir radius, re, is the only reservoir parameter affecting the maximum rupture size: The area of the largest arrested rupture increases with re (Fig. 2B).

Fig. 2 Influence of reservoir- and fault-related parameters on the evolution of rupture arrest area.

(A) Varying reservoir-related parameters, (B) varying strength parameter S and reservoir radius re, and (C) varying fault-related parameters. The reference solution (bold blue line) is the same in all three panels. In each panel, the same color depicts the solutions for the same varying parameter. For (A) and (C), thin and bold lines indicate solutions for lower and greater value of the corresponding parameter compared to that of the reference solution. For stresses in (C), the thin and thick lines indicate lower and higher values of S, respectively. Transition to runaway ruptures is marked by the dashed lines. The horizontal gray lines indicate the area of the intersection of the reservoir and the fault; it also corresponds to the minimal arrested area. The horizontal dotted gray line in (A) indicates constant maximum arrested area (in other panels, maximum arrested area varies). Note that in (C), the lines for τ0 and σ are almost identical.

Fig. 3 Comparison of pore-pressure distribution inside sealed cylindrical reservoirs induced at time t = 100 days by steady-flux injection at its center, for varying compressibility ct, viscosity μv, porosity φ, and permeability k.

Higher and lower values of parameters are depicted by thicker and thinner lines of the same color, respectively. The numbers indicate the values of the spatial integral of pressure in megapascals per square kilometer.

The effect of fault parameters on arrested rupture area is illustrated in Fig. 2 (B and C). Fault parameters (μs, μd, Dc, τ0, and σ) influence both the timing of the onset of arrested ruptures and the subsequent temporal evolution of the rupture arrest area. For example, faults with lower background stress τ0 (higher S; Fig. 2B) or larger slip-weakening distance Dc (Fig. 2C) can produce larger arrested ruptures. This is expected because both smaller energy available to drive ruptures and larger fracture energy (larger fault resistance to rupture growth) tend to stabilize a fault.

Magnitude of the largest arrested rupture

The magnitude of the largest arrested rupture, Embedded Image, is an important indicator of the stability of a perturbed fault. A fault capable of producing runaway ruptures is relatively unstable. Noting that the largest arrested ruptures are significantly larger than the fluid-pressurized area, we obtain further insight into what controls Embedded Image by approximating the additional stress drop (induced by the pore-pressure perturbation) as a point load superimposed on the background stress drop. This approximation leads to the following estimate of the maximum arrested moment (see Eqs. 19 and 21)Embedded Image(1)where Δp is the average increase in pore pressure inside the reservoir, Ap is the area of the stress perturbation (here, the area of intersection between the reservoir and the fault), μd is the dynamic friction coefficient, and Δτ0 is the background stress drop. The corresponding moment magnitude Embedded Image is computed using Eq. 20.

The result in Eq. 1 reflects how rupture arrest is controlled by a competition between two sources of elastic energy: injection-induced fluid pressure and tectonic prestress. The contributions of these two sources to the elastic energy available for rupture growth (via the stress intensity factor in Eq. 18) are both positive. However, the energy contributed by injection-induced fluid pressure decays with increasing rupture size, whereas the energy contributed by tectonic prestress increases, thereby creating a trade-off between these two strain-energy sources. At the maximum arrest size, both contributions are comparable. Hence, although arrested ruptures can grow well beyond the reservoir boundaries, the reservoir pressure not only triggers them but also contributes significantly to their energy balance.

To draw the connection between Embedded Image and injection parameters, following McGarr (9), we approximate the average increase in pore pressure after injecting a volume ΔV of fluid into a saturated reservoir of volume V asEmbedded Image(2)where κ is the bulk modulus of the reservoir rock. This approximation assumes an incompressible fluid and neglects the effects of porosity on the ability of a rock formation to accommodate fluid pressure diffusion or distribution. For the scenario of Fig. 1, the volume affected by the stress perturbation is V = Aph, where h is the reservoir thickness, and Eqs. 1 and 2 yieldEmbedded Image(3)

We note that, for example, a greater stress drop Δτ0 leads to smaller values of Embedded Image, reflecting the fact that transition to runaway ruptures will occur earlier for a fault with a greater stress drop than for a fault with a lower stress drop. Comparison to computations accounting for the exact stress drop distribution in ~4250 randomly generated reservoir-fault configurations (see Materials and Methods) shows that the point-load approximation of Embedded Image is unbiased and accurate within an uncertainty of ±0.5. More generally, the term F = μd ⋅ ΔpAp in Eq. 1 can be interpreted as the spatial integral of a localized perturbation of potential stress drop resulting from the intersection of a fault with a fluid-pressurized volume or a poroelastic stress distribution or from pressure diffusion along a permeable fault zone channel.

In Fig. 4, we compare our Embedded Image estimates with the largest magnitudes observed in various injection-induced earthquake sequences as a function of the cumulative net injected fluid volumes in nearby wells. We consider three values of γ corresponding to plausible values of key model parameters (h from 10 to 1000 m, Δτ0 from 0.1 to 10 MPa, κ = 50 GPa, and μd = 0.1). An induced earthquake with magnitude below our predicted Embedded Image is interpreted as a spontaneously arrested rupture. Because moment can be accommodated by aseismic slip or by multiple earthquakes on the same fault, Embedded Image is a conservative upper magnitude bound for arrested ruptures. Conversely, an earthquake with magnitude above our predicted Embedded Image is interpreted as a runaway rupture. In reality, arrested ruptures are likely to stop even before they reach their maximum size because even small or weak barriers can significantly affect their propagation. A runaway rupture needs a stronger and larger barrier to be stopped, typically the end of a fault or a geometrical/geological fault segmentation. In this context, larger values of γ and Embedded Image lead to a more stable situation.

Fig. 4 Comparison of our estimate of magnitude of maximum arrested rupture, Embedded Image, for three values of γ with magnitudes of injection-induced earthquakes over a broad range of injected volumes.

We find that our estimate is equivalent to that of magnitude of maximum possible earthquake by van der Elst et al. (12) for b = 1; therefore, we also indicate the corresponding values of seismogenic index Ap. For completeness, we also include the estimate of maximum possible induced earthquake by McGarr (9). Note that for the events induced during multistage hydrofracking in the Western Canada sedimentary basin (yellow circles), we consider the total volumes for all previous stages and all proximal well completions reported as “maximum volume” by Atkinson et al. (11).

Figure 4 reveals a striking consistency between the relation Embedded Image predicted by our model and data recorded across a very broad range of injected volumes and scales, from laboratory experiments (centimeter scale) to induced seismicity at field scales (hectometers to kilometers), including intermediate-scale natural laboratories (meters to dekameters). Moreover, almost all magnitudes fall below the predicted Embedded Image, with γ = 1.5 × 1010. These two observations suggest that these induced earthquakes may have been arrested ruptures.

DISCUSSION

Our model does not provide an estimate of maximum possible magnitude Mmax, only the maximum size of arrested ruptures Embedded Image. Ruptures can grow larger if they are runaway. It is nevertheless informative to compare it with previous estimates of Mmax of injection-induced earthquakes.

McGarr (9) proposed a linear relation between M0 max and ΔV, which differs significantly from our physics-based prediction Embedded Image. Although for low injected volumes, McGarr’s model largely overpredicts observed maximum magnitudes, several cases from data sets by Buijze et al. (23) and Atkinson et al. (11) exceed his model for large injected volumes. Our model is more consistent with the data as a whole; however, the difference between our model and maximum observed magnitudes increases with injected fluid volume, suggesting that fluid injection can induce even larger events than what has been recorded so far. This difference concerns magnitudes >4 and is therefore significant for seismic hazard analysis.

Comparing Embedded Image with the estimate Embedded Image by van der Elst et al. (12), based on the premise that induced seismicity is Poissonian and follows the Gutenberg-Richter magnitude-frequency distribution, like regular tectonic earthquakes, reveals interesting similarities. Their model predicts Embedded Image, where b is the Gutenberg-Richter exponent. If b = 1, both models predict exponent of 3/2. Equating the prefactors, the seismogenic index Σ of their model can be related to our parameter γ by Embedded Image. The range of γ values in Fig. 4 yields a reasonable range of values of Σ. Because Embedded Image estimates maximum magnitude and our Embedded Image estimates the largest arrested rupture, the similarity between the two estimates is unexpected. For b ≠ 1 (b > 1 is common for induced seismicity), the two models no longer yield consistent results. Using our terminology, van der Elst et al. (12) showed that runaway ruptures that stopped because of natural heterogeneity of tectonic stresses are a viable explanation of the Mmax versus injected volume data. Our model allows both possibilities (the largest observed ruptures could be either runaway or self-arrested), and we show that the latter is viable too. However, our attempts to discriminate the two models by their predictions of the relative timing of the largest earthquake are inconclusive because of the scarcity of data (see Materials and Methods).

Our analysis shows that rupture arrest and the transition to runaway rupture are dominantly controlled by friction parameters and stress state (that is, conditions on the fault before injection), whereas the temporal evolution of pore pressure only influences timing (that is, triggering). Our model also reveals that for the same pore pressure, a larger reservoir can produce larger events than a smaller reservoir. This implies that there is no universal safe pressure limit; what matters is the product F = μd ⋅ ΔpAp. This also explains, conceptually, how low-pressure wastewater disposal into an extended reservoir may induce larger events than hydrofracking operations involving higher pore pressure but over a smaller footprint.

Our model also reveals significant variations in the duration of the relatively stable, arrested rupture regime. This may have important implications for traffic-light systems (TLSs) aiming at operational control of induced earthquake hazard. A TLS that implicitly assumes that magnitudes of induced earthquakes do not change abruptly with time or a TLS with magnitude thresholds higher than Embedded Image could fail to capture an abrupt transition to runaway ruptures. In the framework of our model, a TLS may be particularly challenged in reservoir-fault systems characterized by a short arrested rupture regime.

Although our scaling of Embedded Image with injected volume has been derived for the point-force approximation of the injection scenario captured in Fig. 1, we believe that it is more robust and applicable to a wider range of scenarios. Strictly speaking, the point-force approximation is valid if the pressurized region is small compared to the size of the maximum arrested rupture. However, the verification of the point-force approximation using the finite-reservoir approach revealed that the point-force approximation is robust, as demonstrated by no systematic correlation of residuals with other fault-reservoir parameters (fig. S6), particularly with strength parameter S (from 1.5 to 8), radius of reservoir (from 0.1 to 5 km), or magnitudes (from 4 to 8). The only exceptions are compressibility and porosity, whose effects, as discussed, are neglected in the approximation of the average pore pressure (Eq. 2). Moreover, our approach can be applied to more general injection scenarios than depicted in Fig 1. For example, if an arbitrarily oriented fault intersects a general-shaped reservoir, the resulting scaling remains Embedded Image, but reservoir thickness h in the parameter γ is replaced by the ratio of volume of the reservoir and the area of the intersection, V/Ap. For a special case of a vertical fault intersecting a cylindrical reservoir at distance D from the reservoir center, Embedded Image.

Beyond injection-induced seismicity, our theory for rupture arrest is potentially applicable to faults loaded by localized stresses of other origins, which can be represented by an effective load F, including situations that are relevant for natural seismicity. For instance, F could arise from stochastic stresses (21) or from stress concentrations due to localized fault loading (19). In particular, creep in the deep extension of a fault concentrates stress at the base of the seismogenic zone. In earthquake cycle simulations, this loading nucleates sequences of arrested ruptures at a depth before a runaway rupture. Megathrust earthquakes that remained confined at depth, like the 2015 Gorkha, Nepal earthquake (24, 25), may be examples of these arrested ruptures, a necessary prelude to a much larger event. The belt of background seismicity in Nepal is consistent with concentrated tectonic loading around the bottom of the megathrust seismic coupling zone. The appropriateness of one of our key model assumptions, the Griffith’s fracture criterion, for natural earthquake arrest is supported by microseismicity observations. On strike-slip faults, the criterion predicts that rupture arrest length is larger horizontally (mode II) than vertically (mode III) by a factor 1/(1 − ν), where ν is the Poisson’s ratio. This is consistent with the aspect ratio of rupture areas inferred from aftershocks of microearthquakes in California (26).

Although our approach allows scaling of fracture energy with rupture size, here, we assumed constant fracture energy on each rupture. Considering an empirical model (27) and a thermal pressurization model (28) of fracture energy scaling with slip, which both fit the relevant seismological observations adequately, we find qualitatively different results for rupture arrest size. Thus, further study of earthquake arrest can help advance our understanding of fracture energy scaling.

MATERIALS AND METHODS

Pore pressure in a cylindrical reservoir

We used the general solution of the diffusion equation for single-phase fluid flow for a cylindrical reservoir, with constant injection rate at its center and no-flow boundary condition [Appendix B in the study of Lee et al. (29)]. The dimensionless diffusion equation exploiting axisymmetry isEmbedded Image(4)where rD = r/rw is the dimensionless radius (r being the radius and rw being the wellbore radius)Embedded Image(5)is the dimensionless time (t being the time, k and φ being the permeability and porosity of the reservoir, respectively, and ct being the compressibility of the fluid), andEmbedded Image(6)is the dimensionless pressure (h being the reservoir height, q being the injection rate (negative sign indicates withdrawal), B being the formation volume factor, μv being the fluid viscosity, pi being the initial pore pressure, and p being the current pore pressure). The initial condition isEmbedded Image(7)the inner boundary condition (that is, constant injection rate) is given asEmbedded Image(8)whereas the outer boundary condition (that is, no-flow boundary) is given asEmbedded Image(9)where reD = re/rw is dimensionless radius of the reservoir. The general solution for rDreD and tD ≥ 0 isEmbedded Image(10)where Jk and Yk are the kth-order Bessel functions of the first and second kind, respectivelyEmbedded Image(11)and αn(reD) is the nth root of the equationEmbedded Image(12)

The last term in the pD solution (Eq. 10) is a transient component that vanishes at long times.

Estimation of rupture arrest area

We derived estimates of the rupture arrest area given a spatial distribution of potential stress drop, Δτ, based on the Griffith crack equilibrium criterion and small-scale yielding fracture mechanics. Following the approach in Appendix B in the study of Ripperger et al. (21) and Appendix A in the study of Galis et al. (17), we adopted the following simplifying assumptions:

(1) The rupture is approximately circular, with radius R.

(2) The stress drop distribution is axisymmetric, Δτ(r).

(3) The static stress intensity factor averaged along the crack rim is approximated by the expression for tensile (mode I) cracksEmbedded Image(13)

(4) The details of weakening inside the process zone are ignored, and the rupture criterion is based on the fracture toughness Kc, related to the slip-weakening fracture energy byEmbedded Image(14)

Although we considered here a slip-weakening model with constant Gc, a scale-dependent Gc(R) can be incorporated too. The conditions for rupture arrest areEmbedded Image(15)

The first condition is the Griffith’s crack equilibrium criterion, and the second condition is necessary to find a stable equilibrium at which further perturbation does not lead to further self-accelerating dynamic rupture. Galis et al. (17) considered an adjustable factor η in the Griffith’s criterion, that is, η ⋅ K0(R) = Kc, as a proxy to account for the differences between modes I, II, and III, the departures from circularity, the effect of dynamic overshoot, etc. They found that η = 1 yields results consistent with numerical simulations; therefore, we adopted η = 1. In numerical simulations, the stresses in the hypocentral area were chosen so that the initial condition was near an unstable equilibrium state (∂(K0Kc)/∂R > 0).

The solution of the problem (15) was found by numerically integrating K0(R) (Eq. 13) with adaptive sampling of R and r to achieve high accuracy of the rupture arrest radius, Rarr, for various stress drop distributions. The rupture arrest area was then obtained as Embedded Image. If a solution to Eq. 15 exists, it represents the final radius R of an arrested rupture. Otherwise, rupture is runaway.

We considered stress drop distributions with identical peak amplitude such that τ0 > τs at r = 0 but with varying width (Fig. 5). Such a peak could result from reduced normal stress due to locally increased pore pressure. The stress drop perturbation affects the shape of K0(R) following three cases depending on its width. If the perturbation is too narrow, rupture does not propagate outside the overstressed area because the perturbed K0 does not reach Kc. If the perturbation is large enough, it initiates a runaway rupture because, once K0 becomes greater than Kc, it will remain greater. In the intermediate case, rupture starts propagating outside the overstressed region, but the rupture eventually stops when the stable equilibrium is reached, that is, K0 = Kc on the decreasing leg of K0.

Fig. 5 Fracture mechanics of rupture arrest under nonuniform stress.

Sketch of (A) stress drop distributions Δτ(r) and (B) the corresponding static stress intensity factors K0(R) illustrating different crack behavior (although our approach supports Kc being a function of R, for illustration purposes, Kc is plotted as a constant).

We obtained further physical insight by developing a simplified analytical estimate. For this purpose, we assumed that Kc is constant and the rupture arrest area is much larger than the fluid-pressurized area Ap. Following the study of Galis et al. (17), we approximated stress drop by a point load superimposed on the background uniform stress dropEmbedded Image(16)whereEmbedded Image(17)

Here, Δp is the average increase in pore pressure inside the reservoir. We adopted a definition of Dirac’s distribution in polar coordinates consistent with Embedded Image (30). After inserting Eq. 16 into Eq. 13, we obtainedEmbedded Image(18)

Solving the Griffith’s equilibrium condition (Eq. 15) leads to a quartic equation with too complicated solutions to provide the desired insight. Yet, a compact expression can be derived for the maximum rupture arrest radius, Rmax − arr. The first term of Eq. 18 decreases with R, whereas the second term increases; hence, K0(R) has a single minimum. Runaway ruptures are possible only if this minimum exceeds Kc because, then, the Griffith’s equilibrium criterion cannot be satisfied for any R. The rupture radius at the minimum of K0 givesEmbedded Image(19)

The corresponding rupture arrest area is Embedded Image.

Verification of theoretical estimates of rupture arrest area

We verified our theoretical estimates of rupture arrest area by comparison with 3D dynamic rupture simulations. First, we analyzed results for a step-like distribution of stress drop with Δτ(r) = Δτi for rRi and Δτ(r) = Δτ0 for r > Ri. Galis et al. (17) performed dynamic rupture simulations on a 30 × 15–km fault with μs = 0.6778, μd = 0.525, dc = 0.4 m, and σ = 120 MPa and various values of strength parameter S (obtained by varying τ0) to investigate conditions for runaway ruptures. The fault is large enough so that ruptures that break the entire fault are considered runaway ruptures. Here, we considered the same configurations. Figure 6A shows very good agreement between estimated and simulated rupture arrest areas for all considered values of S and Ri, including the transition to runaway ruptures. The comparison further revealed that our approach worked equally well for circular, elliptical, and square perturbations, although it was based on an assumption of radially symmetric stress drop. To achieve agreement between theory and numerical results, we introduced a correction factor ϕ = 2 such that Embedded Image or Embedded Image. Galis et al. (17) introduced an adjustable factor η multiplying K0 to account for deviations from circular mode I rupture. The factor ϕ introduced here has a different role; it accounts for the effects affecting crack radius R. In terms of the K0(R) plot in Fig. 5B, η provides a correction for the y axis, whereas ϕ provides a correction for the x axis.

Fig. 6 Verification of our theoretical estimates with numerical results.

Comparison of estimated rupture arrest area Aarr and maximum rupture arrest area Amax with rupture arrest area from numerical simulations Embedded Image as functions of overstressed area for (A) a steplike distribution of stress drop and varying strength parameter S (B), a Gaussian distribution of stress drop, and varying minimal normal stress σmin. For a step-like stress drop distribution, we performed simulations with circular, square, and elliptical perturbations marked by circles, squares, and triangles. The transitions from solid to dashed lines indicate transition from arrested to runaway ruptures. The undulated double line indicates interrupted y axis; r.r. indicates runaway ruptures. Lfric = μ ⋅ Dc/(τs0 − τd0) is a characteristic length scale introduced by the slip-weakening process (39).

Next, we analyzed results for stress drop with a Gaussian perturbation, Δτ(r) = Δτ0 + μdσ0′(1 − σmin) ⋅ exp(−r2/(2 χ2)), as would result from an instantaneous pressure point source at r = 0 (31). We adopted the same friction parameters as those of Galis et al. (17), yielding Δτ0 = 6.1 MPa, and considered σmin = 0.5, 0.65, 0.75, and 0.8. Figure 6B shows the overall consistency between estimated and simulated rupture arrest areas as a function of the characteristic area of the stress perturbation, πχ2. Part of the small differences decrease with mesh refinement (fig. S1) and are attributed to the staircase representation of the smooth stress drop distribution in a discrete grid. For the larger σmin values considered, the critical perturbation area at the transition to runaway ruptures was underestimated. This was expected because our theoretical estimates lacked a rupture nucleation criterion to account for the finite size of the initiation area (where the initial stress is higher than the static strength).

Figure 6 also shows very good agreement between theory and numerical results for the maximum area of spontaneously arrested ruptures, Amax. The agreement improved with increasing S for the steplike stress drop and with decreasing σmin for the Gaussian distribution because in this case, the assumptions for the point-load approximation were better satisfied (fig. S2). Although both examples represent spatially variable stress drop and nucleation of ruptures by increased fluid pressure, the results for step-like distributions are more relevant for a sealed reservoir.

Conversion of rupture arrest size to magnitude

Moment magnitude Mw is defined by (32)Embedded Image(20)where M0 is the seismic moment. We derived the seismic moment from the initial stress information using the equation for a circular crack with uniform stress drop (33)Embedded Image(21)

Strictly speaking, the formula should be based on the static stress drop, which differs from the dynamic stress drop Δτ0 due to dynamic overshoot. The average dynamic and static stress drops of arrested ruptures are almost identical (fig. S3), and in highly heterogeneous ruptures, they differ by only 10% on average (34). Because the pressurized region is small compared to the rupture area, we approximated the average stress drop as Δτ0, neglecting the additional stress drop μd ⋅ Δp inside the reservoir. Figure S4 compares magnitudes estimated by this approach with those obtained from dynamic rupture simulations by directly applying the definition of seismic momentEmbedded Image(22)where μ is the shear modulus, and D is the final slip. Rupture simulations were initiated artificially by prescribing an area in which the initial stress slightly exceeded the static frictional strength. We excluded this initiation area from the magnitude computation because it contained artificially large slip. Overall, for rupture areas from 1 to 500 km2, the approximation to Mw using Eq. 21 agrees with the numerical simulation results, with deviations smaller than ~0.1 and insignificant underestimation. It is also naturally consistent with empirical scaling relations based on the self-similar rupture model, in particular, with the relation for stable continental regions by Leonard (35), which assumes a stress drop value (5.84 MPa) that is very similar to the background stress drop in our simulations (6.11 MPa).

Verification of point-load approximation of reservoir

Here, we verified the adequacy of the maximum arrested rupture magnitude estimated under the point-load approximation, Embedded Image. We considered ~4250 randomly generated reservoir-fault configurations. Values of reservoir parameters (that is, viscosity, permeability, compressibility, porosity, well radius, flow rate, height, and radius of the reservoir) were sampled as independent, uniformly distributed variables within ranges reported for common reservoirs (36, 37). The following fault-related parameters were sampled from uniform distributions: μs, μds, depth, and S. Normal background stress was derived from depth assuming a typical density of shallow crust and water (2500 and 1000 kg m−3, respectively). To improve the sampling of lower magnitudes events, we increased the number of models with lower S, reservoir radius, and Dc, the only three parameters that exhibited signs of correlation with magnitude. Histograms of model parameters are shown in fig. S5. For each parameter setting, we computed the exact value of the maximum arrested earthquake magnitude, Embedded Image, by numerically solving the rupture arrest conditions accounting completely for the finite size of the reservoir. Correlation plots between model parameters and Embedded Image are shown in fig. S6.

Figure S7 reveals good agreement between the point-load approximation (Embedded Image) and the finite-reservoir calculation (Embedded Image). Scatter was expected because of the approximation of the average reservoir pore pressure Δp (Eq. 2). In particular, the approximation neglects the effects of compressibility and effects of porosity on the ability of a formation to accommodate fluid. Both assumptions lead to systematic deviation of Δp from the true average pore pressure, which was also manifested by the weak correlation of compressibility and porosity with residuals (fig. S6). The scatter of Embedded Image is mostly confined within ±0.5, with only few points reaching down to −1. This indicates that the influence of compressibility and porosity is not strong enough to cause a significant deviation of Embedded Image from Embedded Image. The approximation of Δp also neglects the effects of viscosity and permeability, but as shown in Fig. 3, these parameters affect the pressure distribution without affecting the average pressure. The nonsymmetric distribution of the orthogonal residuals has mean and median close to zero (− 0.13 and − 0.1, respectively). We concluded that Embedded Image is consistent with Embedded Image with uncertainty ±0.5.

Predictions of the relative timing of the largest earthquake

In the model by van der Elst et al. (12), the probability of occurrence of the largest earthquake is uniformly distributed within a sequence. In our model, the potential size of arrested ruptures grows with ongoing injection; hence, the largest triggered event has a higher chance to occur later in the sequence. However, our model does not make a specific prediction about how strong this tendency is; its effect could be a slight increase of probability with time. Moreover, the induced seismicity sequences considered by van der Elst et al. (12) contain also aftershocks of the largest event, which tends to bias low the rank of occurrence of the largest event.

We applied the two-sample Anderson-Darling (AD) test (38) to assess whether the occurrence rank of the largest event in the recorded induced sequences is drawn from each of the following three distributions: uniform, linearly increasing, and Gaussian (the latter being a proxy for a distribution biased by aftershocks). For the linear and Gaussian distributions, we determined optimal parameters by grid search to minimize the AD test statistic. Comparison of the optimal AD test values showed that all three considered distributions can explain the data with similar level of statistical significance and, therefore, none of them can be rejected (fig. S8). Consequently, our model qualitatively explains this data similarly well as the model proposed by van der Elst et al. (12). However, the current data set is small (only 17 sequences), and a larger data set may facilitate the discrimination between models in the future.

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/3/12/eaap7528/DC1

fig. S1. Comparison of dimensionless rupture arrest area calculated from numerical simulations with grid spacing h = 50 m (circles) and h = 100 m (squares) with our theoretical estimates (bold lines) for varying σmin (indicated by color).

fig. S2. Comparison of stress drop distributions as functions of dimensionless crack radius at the time of Embedded Image for situations from Fig. 6.

fig. S3. Scaling of mean static and dynamic stress drops in results of numerical simulations.

fig. S4. Comparison of various approaches to estimate Mw from ruptured area Aarr.

fig. S5. Distributions of reservoir-fault parameters for all ~4250 configurations used for verification of Embedded Image (point-load approximation) against Embedded Image (finite-reservoir approach).

fig. S6. Distributions of reservoir-fault parameters for all ~4250 configurations used for verification of Embedded Image (point-load approximation) against Embedded Image (finite-reservoir approach).

fig. S7. Comparison of Embedded Image (derived for a point-load approximation of a reservoir) with Embedded Image (derived for a finite reservoir) and corresponding orthogonal residuals.

fig. S8. Evaluation of the probability of occurrence rank of the largest event within a sequence.

table S1. Reservoir and fault parameters used to prepare Fig. 2.

References (40, 41)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

REFERENCES AND NOTES

Acknowledgments: Funding: Research presented in this paper is supported by King Abdullah University of Science and Technology (KAUST) in Thuwal, Saudi Arabia (grants BAS/1339-01-01 and URF/1/2160-01-01) and by the NSF (CAREER award EAR-1151926). Some of the 3D dynamic rupture simulations for verification of our model have been carried out using the KAUST Supercomputing Laboratory. We thank the Agence Nationale de la Recherche through the HYDROSEIS project (Role of fluids and fault HYDROmechanics on SEISmic rupture) under contract ANR-13-JS06-0004-01 for supporting the in situ experiments providing the data (Duboeuf et al.) used in Fig. 4. We also thank S. Goodfellow and L. De Barros for providing their laboratory and in situ data used in Fig. 4. J.P.A. and F.C. thank the Observatoire de la Côte d’Azur for supporting this research. Author contributions: M.G. and J.P.A. developed the main ideas, interpreted the results, and produced the manuscript. P.M.M. contributed to the discussions and interpretations of the results and commented on the manuscript at all stages. F.C. provided the data for Fig. 4, contributed to the discussions of the results, and commented on the manuscript at all stages. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.
View Abstract

Navigate This Article