## 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, *M*_{max}, 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 *M*_{max} and fluid injection or production parameters is a topic of active research (*8*–*13*), previously proposed models to estimate *M*_{max} 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, *M*_{max}, 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, *D*_{c}. 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} + μ_{d} ⋅ *p*, where Δτ_{0} = τ_{0} − μ_{d}σ_{0} is the uniform background stress drop.

Fluid injection into the reservoir increases its pore pressure, which reduces the fault’s frictional strength μ_{s}(σ_{0}′ − *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, *G*_{c} = 1/2σ′(μ_{s} − μ_{d})*D*_{c}. 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*) = *G*_{c} and apply an additional condition, *d*(*G* − *G*_{c})/*dR* < 0, to identify stable equilibrium position. Finally, we use a circular crack model with uniform stress drop to compute the moment magnitude *M*_{w} 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 *c*_{t}, 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, *r*_{e}, is the only reservoir parameter affecting the maximum rupture size: The area of the largest arrested rupture increases with *r*_{e} (Fig. 2B).

The effect of fault parameters on arrested rupture area is illustrated in Fig. 2 (B and C). Fault parameters (μ_{s}, μ_{d}, *D*_{c}, τ_{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 *D*_{c} (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, , 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 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)(1)where Δ*p* is the average increase in pore pressure inside the reservoir, *A*_{p} 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 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 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* as(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* = *A*_{p}⋅*h*, where *h* is the reservoir thickness, and Eqs. 1 and 2 yield(3)

We note that, for example, a greater stress drop Δτ_{0} leads to smaller values of , 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 is unbiased and accurate within an uncertainty of ±0.5. More generally, the term *F* = μ_{d} ⋅ Δ*p* ⋅ *A*_{p} 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 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 is interpreted as a spontaneously arrested rupture. Because moment can be accommodated by aseismic slip or by multiple earthquakes on the same fault, is a conservative upper magnitude bound for arrested ruptures. Conversely, an earthquake with magnitude above our predicted 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 lead to a more stable situation.

Figure 4 reveals a striking consistency between the relation 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 , with γ = 1.5 × 10^{10}. 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 *M*_{max}, only the maximum size of arrested ruptures . Ruptures can grow larger if they are runaway. It is nevertheless informative to compare it with previous estimates of *M*_{max} of injection-induced earthquakes.

McGarr (*9*) proposed a linear relation between *M*_{0 max} and Δ*V*, which differs significantly from our physics-based prediction . 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 with the estimate 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 , 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 . The range of γ values in Fig. 4 yields a reasonable range of values of Σ. Because estimates maximum magnitude and our 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 *M*_{max} 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} ⋅ Δ*p* ⋅ *A*_{p}. 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 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 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 , but reservoir thickness *h* in the parameter γ is replaced by the ratio of volume of the reservoir and the area of the intersection, *V*/*A*_{p}. For a special case of a vertical fault intersecting a cylindrical reservoir at distance *D* from the reservoir center, .

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 is(4)where *r*_{D} = *r*/*r*_{w} is the dimensionless radius (*r* being the radius and *r*_{w} being the wellbore radius)(5)is the dimensionless time (*t* being the time, *k* and φ being the permeability and porosity of the reservoir, respectively, and *c*_{t} being the compressibility of the fluid), and(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, *p*_{i} being the initial pore pressure, and *p* being the current pore pressure). The initial condition is(7)the inner boundary condition (that is, constant injection rate) is given as(8)whereas the outer boundary condition (that is, no-flow boundary) is given as(9)where *r*_{eD} = *r*_{e}/*r*_{w} is dimensionless radius of the reservoir. The general solution for *r*_{D} ≤ *r*_{eD} and *t*_{D} ≥ 0 is(10)where *J*_{k} and *Y*_{k} are the *k*th-order Bessel functions of the first and second kind, respectively(11)and α_{n}(*r*_{eD}) is the *n*th root of the equation(12)

The last term in the *p*_{D} 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) cracks(13)

(4) The details of weakening inside the process zone are ignored, and the rupture criterion is based on the fracture toughness *K*_{c}, related to the slip-weakening fracture energy by(14)

Although we considered here a slip-weakening model with constant *G*_{c}, a scale-dependent *G*_{c}(*R*) can be incorporated too. The conditions for rupture arrest are(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, η ⋅ *K*_{0}(*R*) = *K*_{c}, 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 (∂(*K*_{0} − *K*_{c})/∂*R* > 0).

The solution of the problem (*15*) was found by numerically integrating *K*_{0}(*R*) (Eq. 13) with adaptive sampling of *R* and *r* to achieve high accuracy of the rupture arrest radius, *R*_{arr}, for various stress drop distributions. The rupture arrest area was then obtained as . 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 *K*_{0}(*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 *K*_{0} does not reach *K*_{c}. If the perturbation is large enough, it initiates a runaway rupture because, once *K*_{0} becomes greater than *K*_{c}, 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, *K*_{0} = *K*_{c} on the decreasing leg of *K*_{0}.

We obtained further physical insight by developing a simplified analytical estimate. For this purpose, we assumed that *K*_{c} is constant and the rupture arrest area is much larger than the fluid-pressurized area *A*_{p}. Following the study of Galis *et al*. (*17*), we approximated stress drop by a point load superimposed on the background uniform stress drop(16)where(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 (*30*). After inserting Eq. 16 into Eq. 13, we obtained(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, *R*_{max − arr}. The first term of Eq. 18 decreases with *R*, whereas the second term increases; hence, *K*_{0}(*R*) has a single minimum. Runaway ruptures are possible only if this minimum exceeds *K*_{c} because, then, the Griffith’s equilibrium criterion cannot be satisfied for any *R*. The rupture radius at the minimum of *K*_{0} gives(19)

The corresponding rupture arrest area is .

### 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 *r* ≤ *R*_{i} and Δτ(*r*) = Δτ_{0} for *r* > *R*_{i}. Galis *et al*. (*17*) performed dynamic rupture simulations on a 30 × 15–km fault with μ_{s} = 0.6778, μ_{d} = 0.525, *d*_{c} = 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 *R*_{i}, 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 or . Galis *et al*. (*17*) introduced an adjustable factor η multiplying *K*_{0} 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 *K*_{0}(*R*) plot in Fig. 5B, η provides a correction for the *y* axis, whereas ϕ provides a correction for the *x* axis.

Next, we analyzed results for stress drop with a Gaussian perturbation, Δτ(*r*) = Δτ_{0} + μ_{d}σ_{0}′(1 − σ^{min}) ⋅ exp(−*r*^{2}/(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, *A*_{max}. 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 *M*_{w} is defined by (*32*)(20)where *M*_{0} 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*)(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 moment(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 km^{2}, the approximation to *M*_{w} 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, . 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}, μ_{d}/μ_{s}, 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 *D*_{c}, 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, , by numerically solving the rupture arrest conditions accounting completely for the finite size of the reservoir. Correlation plots between model parameters and are shown in fig. S6.

Figure S7 reveals good agreement between the point-load approximation () and the finite-reservoir calculation (). 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 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 from . 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 is consistent with 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 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 *M*_{w} from ruptured area *A*_{arr}.

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

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

fig. S7. Comparison of (derived for a point-load approximation of a reservoir) with (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.

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.

- Copyright © 2017 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).