Role of Fluid Injection on Earthquake Size in Dynamic Rupture Simulations on Rough Faults

An outstanding question for induced seismicity is whether the volume of injected fluid and/or the spatial extent of the resulting pore pressure and stress perturbations limit rupture size. We simulate ruptures with and without injection‐induced pore pressure perturbations, using 2‐D dynamic rupture simulations on rough faults. Ruptures are not necessarily limited by pressure perturbations when (1) background shear stress is above a critical value, or (2) pore pressure is high. Both conditions depend on fault roughness. Stress heterogeneity from fault roughness primarily determines where ruptures stop; pore pressure has a secondary effect. Ruptures may be limited by fluid volume or pressure perturbation extent when background stress and fault roughness are low, and the maximum pore pressure perturbation is less than 10% of the background effective normal stress. Future work should combine our methodology with simulation of the loading, injection, and nucleation phases to improve understanding of injection‐induced ruptures.


Introduction
An important question in the study of induced seismicity is whether earthquake magnitudes are limited by the volume of injected fluid or some other injection-related parameter (e.g., Baisch et al., 2010;Maurer & Segall, 2018;McGarr, 2014;McGarr & Barbour, 2017;Shapiro et al., 2011Shapiro et al., , 2013, or follow naturally occurring (Gutenberg-Richter) size variability (van der Elst et al., 2016). For example, McGarr and Barbour (2017) propose an upper bound on seismic moment released by induced earthquakes, M max 0 , defined by where G is shear modulus and ΔV is injected volume. The premise of such an approach is that a pore pressure perturbation diffuses through the medium, perturbing the effective stress in a finite volume of crust sufficient to induce and maintain rupture, while stress conditions outside the perturbed region do not allow rupture. To evaluate this hypothesis, we consider the behavior of individual simulated ruptures perturbed by spatially variable pore pressure increases.
Linear elastic fracture mechanics predicts that under uniform background stress conditions and constant fracture energy, a crack introduced to an elastic solid will grow unstably if its length exceeds a critical ©2020. The Authors.
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

10.1029/2020GL088377
Key Points: • Rupture size is not necessarily limited by the volume of injected fluid in earthquake simulations with imposed pore pressure perturbations • Stress heterogeneity arising from geometric roughness may be the primary cause for rupture termination on rough faults • High initial shear stress or pore pressure can trigger a rupture larger than the volume-based magnitude limit on faults with low roughness Supporting Information: • Supporting Information S1 Correspondence to: J. Maurer, jmaurer@mst.edu value a c . Assuming linear slip-weakening friction on a preexisting fault, a c is proportional to the ratio of peak minus residual strength (τ p −τ r ) and the square of the static stress drop Δτ (Andrews, 1976): where G is shear modulus, ν is Poisson's ratio, D c is the slip-weakening distance, and f is a factor related to the geometry of the problem. In this scenario, there are two possibilities: a crack that does not reach half-length a c will naturally self-arrest, while a crack that does will slip indefinitely. Galis et al. (2017) applied this reasoning to fluid-induced earthquakes to estimate the size of the largest self-arresting ruptures for spatially variable peak strength. They considered a stress perturbation due to a pore pressure perturbation in an otherwise-uniform background stress, approximated as a point load. Since background stress is uniform, when the stress is low the localized strength drop provided by the perturbation drives slip into the (unfavorable) stress environment beyond the pressurized zone. If the background shear stress is high enough, the rupture will continue to grow without limit.
Norbeck and Horne (2018) considered quasi-dynamic simulations of induced earthquakes on flat faults with linear slip weakening friction. Based on their simulations, they proposed that induced earthquakes are gov- (f b is the initial background shear to effective normal stress ratio, f D is dynamic friction, and σ 0 ¼ σ 0 − Δp, where σ 0 is the total normal stress and Δp is the pore pressure perturbation.) Events on faults for which f b /f D <1 were limited to the pressurized zone, while f b /f D >1 resulted in runaway ruptures, irrespective of volume injected.
In these studies, the only source of stress heterogeneity is that of the perturbations in pore pressure. However, preexisting stress heterogeneity on faults occurs due to geometric roughness and past fault slip, among other sources.  and  investigated the role of heterogeneity on the size distribution of induced earthquakes on 1-D flat faults using a fracture mechanics approach. They solved the crack equation of motion numerically (Freund, 1990) for a suite of stochastic (fractal) shear stress profiles. Ruptures arrest naturally due to variations in shear stress, and  showed that the distribution of rupture size was controlled by the interaction between the spatial distribution of the pore pressure perturbation and the statistical characteristics of the fractal stress profiles. In their model, stress heterogeneity was imposed as an initial condition, and the rupture size calculation did not account for the potential effects of fault roughness (which influences both shear and normal tractions) and off-fault plasticity. These effects result in fracture energy that cannot be predicted a priori, and higher background stress required for rupture (Dieterich & Smith, 2009;Fang & Dunham, 2013).
In this study, we address these issues and explore the hypothesis that induced earthquakes are limited in size by the magnitude and/or spatial extent of the pore pressure perturbation, in the context of 1-D rough (fractal) faults embedded in a 2-D elasto-visco-plastic medium and obeying a rate-state friction law with strong dynamic weakening (Dunham et al., 2011a(Dunham et al., , 2011b. In contrast to the slip-weakening models discussed above, rate-state friction does not have a well-defined residual strength. However, for strong rate weakening friction, there exists a critical stress level τ pulse , at which self-sustaining rupture on flat faults is just possible (Dunham et al., 2011b;Zheng & Rice, 1998). When the background shear stress is close to τ pulse (referred to here as "low-stress"), ruptures are pulse-like: slip occurs in a narrow pulse just behind the rupture front, and shear strength recovers behind the rupture tip (e.g., Beeler & Tullis, 1996;Cochard & Madariaga, 1994;Zheng & Rice, 1998).
We simulate earthquakes with and without pore pressure and stress perturbations to determine whether rupture size is limited by the volume of injected fluid and/or the spatial extent of the stress changes. Since faults are geometrically rough, we generate several thousand stochastic realizations in order to characterize results statistically. At low background shear stress, one might expect the extent of the stress and pore pressure perturbations to exert some control on rupture lengths. However, we find that events may be larger than the pressurized region even at low stress if the magnitude of the pore pressure perturbation is sufficiently large. Ruptures are not confined when stress is high, consistent with Norbeck and Horne (2018) and Galis et al. (2017). Our results suggest that dynamic effects and in situ stress conditions interact with pore pressure and poroelastic stress perturbations to influence rupture size and that low stress conditions may not be sufficient to guarantee ruptures smaller than an injection-related threshold.

2-D Dynamic Earthquake Simulations
We use the 2-D plane strain rupture dynamics code FDMAP (Dunham et al., 2011a(Dunham et al., , 2011bKozdon et al., 2011bKozdon et al., , 2013; see section 6). The model employs a rate-and-state friction formulation in the slip law form with strong rate weakening on the fault and Drucker-Prager visco-plasticity in the off-fault material (Dunham et al., 2011a;Noda et al., 2009;Rice, 1983). There is no quasi-static nucleation phase; events are artificially initiated by adding a Gaussian shear stress perturbation at the first time step. Once initiated, the rupture process is entirely self-governed. Faults are 1-D self-similar fractal profiles and are oriented such that they lie along the y = 0 line of the model domain; flat faults are on the line exactly while rough faults follow it on average. Roughness, parameterized by amplitude to wavelength ratio α (supporting information, Figure S1), is band-limited, with minimum and maximum wavelengths of 300 m and 60 km. Values of α on natural faults are thought to vary over an order of magnitude or more, ranging from 0.001 or less on mature faults like the San Andreas, up to perhaps 0.01 (e.g., Brodsky et al., 2016;Candela et al., 2009Candela et al., , 2012Fang & Dunham, 2013;Sagy & Brodsky, 2009). The initial stress (including any initial pore pressure) is spatially uniform in the medium; the pore pressure perturbation is spatially variable as described in section 2.3. Resolved tractions on rough faults varies along the fault (see section 2.2), so prior to simulation, the fault profile is shifted such that the least stable part of the fault is located at the origin, where the initiating stress perturbation is applied.

Stress and Slip on Geometrically Rough Faults
Fault roughness provides additional resistance to slip above that of friction, hence rougher faults require higher stress levels for events to propagate (Dieterich & Smith, 2009;Fang & Dunham, 2013). This effect is termed "roughness drag" by Fang and Dunham (2013) and is proportional to slip (s), roughness level (α), and inversely proportional to the minimum roughness wavelength, λ min . In most of our simulations, λ min = 300 m and τ drag is approximately 10 MPa (s/λ min ) (α/10 −3 ) 2 ; however, τ drag increases as λ min decreases (see supporting information). In comparison with the flat-fault simulations (Figure 1), ruptures on rough faults arrest over a wider range of initial background stress ratios and may even arrest and then re-nucleate due to interacting stresses around fault bends (Bruhat et al., 2016).

Pore Pressure Perturbation Models
FDMAP does not model the nucleation phase of rupture; therefore, we run experiments imposing several different pore pressure perturbations as part of the initial conditions. We simulate pore pressure and poroelastic stress changes based on an injector location centered with respect to the fault but offset by 2 km. Events are initiated at the origin, where both the resolved stress ratio (see section 2.1) and the pore pressure are highest. Figures 1a-1c and supporting information Figure S2 show pressure and poroelastic stress changes along the y = 0 line of the model domain for each pore pressure perturbation model.
1. Pressure Model 0 (PM0) is the reference case with no pore pressure perturbation. 2. Pressure Models 1 and 2 (PM1 and PM2; Figures 1a and 1b, respectively) are two realizations of injection into an infinite 2-D (plane strain) poroelastic medium with uniform poroelastic and hydraulic properties, using line source solutions from Rudnicki (1986). We account for the change in total stress from both poroelasticity and pore-pressure in the medium and on the fault. Pressure decays with distance from the origin r as expð−r 2 =4ctÞ , with diffusivity c and time t. (Parameters for the simulations are given in supporting information Tables S1-S2.) The pore pressure perturbation profiles used in our simulations are for 1,000 days of injection with different rates and diffusivities. Peak pore pressure on the y = 0 plane (max Δp) is 2 MPa for PM1 and 19 MPa for PM2 and drops to 10 kPa at 19 km from the origin for PM1 and 12.5 km for PM2 (Figures 1a-1b). 3. In Pressure Model 3 (PM3; Figure 1c), we introduce a high-permeability (k) zone 20 km wide, oriented perpendicular to the fault in the out-of-plane direction and centered at the origin (initiation region), between two symmetric outer regions with low permeability (supporting information Figures S3-S4). We simulate the same volume of injection as in PM1, the only difference being the presence of the high permeability zone. The resulting pressure distribution drops sharply at the boundaries by ∼4 MPa on the y = 0 line, introducing an additional length scale into the problem. We solve numerically for the pressure distribution (Elsworth & Suckale, 2016) (details in the supporting information) and use the pressure to calculate the effective stress in the medium, and ignore poroelastic stress perturbations.

Flat Faults with Strong Rate-Weakening Friction
As a reference, we ran a suite of simulations on flat faults. We show results for PM0, PM2, and PM3 in Figure 1a; note that PM1 ruptures behave qualitatively similar to PM3 but with a smaller effect, so are omitted for clarity. For these simulations, σ 0 ¼ 62 MPa. The stress perturbation required to initiate events results in an slip peak at the origin (see Figures 1e and 1f). Ruptures may arrest immediately or transition to a pulse-like or crack-like rupture mode, depending on the stress ratio f b (Figure 1d).
For PM0 events (solid circles in Figure 1a), there is a narrow transition near τ pulse from self-arresting ruptures to full-fault ruptures, over a range less than 3% of τ b =σ 0 . At low background shear stress (≲0:32σ 0 ) and no pore pressure perturbation, ruptures arrest, while at higher stress ruptures are self-sustaining, consistent with previous work (Dunham et al., 2011a;Gabriel et al., 2012;Zheng & Rice, 1998).
PM2 ruptures initiate, grow, and become full fault at lower levels and over a broader range of background stress ratios than PM0 simulations, due to the decreased strength from additional pore pressure in the nucleation region. Ruptures become self-sustaining at τ b =σ 0 ≈ 0:30, lower than the reference case, even though the stress beyond ±10 km from the origin (L rup /L = 0.33) is very similar to the unperturbed model. That is, the decrease in the pore pressure perturbation results in an increase in fault strength, such that away from the origin the fault is nearly as strong as the unperturbed case. Ruptures are able to propagate through the strong region (once initiated inside the weaker perturbed zone), at stress levels where they could not initiate. This is due partially to the larger shear stress drop in the nucleation region for the perturbed case and partially to the strong dynamic weakening. Ruptures at lower stress may arrest due to the increase in fault strength encountered outside the perturbed region (Figure 1e), consistent with equation (1). Figure 1a) show evidence of arresting due to the spatially variable pore pressure perturbation. Figure 1f shows an example, where the rupture begins to propagate at a constant rate, then dies out upon reaching the edge of the perturbed zone. This is the clearest example of pressure controlling where the rupture stops. At higher background stresses, ruptures grow beyond the pressurized zone to the edge of the computational domain. Thus, for PM3, the increase in frictional strength at the edge of the pressurized region may influence rupture arrest for a small range of stress ratios ∼0.27−0.30.

PM3 ruptures (solid diamonds in
To summarize, flat fault simulations show that (1) pore pressure perturbations leads to rupture at lower shear stress (or larger ruptures) relative to the reference case, and (2) the spatial extent of pore pressure perturbations may limit ruptures in a narrow range of stress conditions, but (3) at high shear stress (τ b =σ 0 > τ pulse ) ruptures are unbounded, consistent with the results of Galis et al. (2017). The question we consider next is how geometric roughness impacts rupture size.

Results on Rough Faults
Results for rough faults at a background effective normal stress of 62 MPa are shown here; results for 126 MPa are shown in the supporting information. For these simulations, α = 0.004−0.012 and f b ∼0.015−0.45. Note that the values of f b are lower than inferred in previous studies of induced seismicity (0.6-0.8; e.g., Walsh & Zoback, 2016), which is because the minimum roughness wavelength in the simulations is much larger than that expected on natural faults (see supporting information). Fault strength at high slip speed depends on fault roughness (due to τ drag ), thus faults with smaller minimum roughness wavelength require higher stress to rupture (see supporting information for more details). Figure 2 shows two example simulations on the same fault with identical parameters, one with no pressure perturbation (PM0) and one with perturbed pressure model PM3. Slip in Figure 2a, without a perturbation, does not extend outside the nucleation region, and therefore is considered an "arrested" rupture, while the simulation with PM3 in Figure 2b ruptures ∼40% of the fault. In this simulation, stress perturbations due to 10.1029/2020GL088377 Geophysical Research Letters fault geometry dominate the initial stress heterogeneity on the fault (10× larger than the pore pressure perturbation). However, the perturbed rupture propagated outside the nucleation region, suggesting that the length scale over which the pressure perturbation acts is an important factor in determining final rupture size. Comparing the initial and final stresses in Figures 2c and 2d show that the PM3 rupture arrests due to encountering low-stress barriers at restraining bends. Supporting information Figures S5 and S6 show additional simulation examples.
In Figure 3, we show summary results for several hundred simulations, illustrating two background stress ratios and roughness levels. The left column in Figure 3 shows empirical frequency-length distributions, while the right column shows frequency-moment distributions. The gray-scale bars at the top left show the spatial extents of pore-pressure perturbation for the different models. Additional event size distributions are shown in supporting information Figures S7-S10. Figure 3 demonstrates the importance of the length scale of the pressure perturbation. Pressure models PM1 and PM3 have the same total injected volume, but PM3 ruptures propagate farther than PM1. The pore pressure perturbation has less of an impact on rupture size at high roughness.
The right column of Figure 3 shows frequency-moment distributions. Moment per unit length in the out-of-plane direction (D), is defined as the product of the shear modulus G with the length-averaged slip s(ξ), where ξ is the arc length along the fault trace of length L: There is a minimum moment imposed by the initiation process of approximately 2 × 10 13 N m /m, while the upper bound on moment corresponds to a full fault rupture (60 km) times a few meters of slip, giving a "full-fault" moment between ∼10 15 and10 16 N m /m, depending on the amount of slip. The injected volume (see supporting information Tables S3-S5 for relevant parameters) is ΔV = 4 × 10 3 m 3 /m for PM1 and PM3, and 2 × 10 4 m 3 /m for PM2. M max 0 from equation 1 is then 2.8 × 10 14 N m/m for PM1 and PM3, and 1.55 × 10 15 N m /m for PM2.
At high background stress (f b = 0.347), all of the moment distributions exceed the hypothesized bounds. At low background stress ratios (f b = 0.282), the distributions tend to tail off well before reaching the hypothesized bounds. At best, PM3 ruptures at low stress (f b = 0.282) arrest close to the magnitude limit theorized by McGarr and Barbour (2017), which may indicate that the pore pressure perturbation may have a secondary role in stopping ruptures when roughness and stress are low (and compare to Figure 1f for PM3 rupture on a flat fault). Even at low background stress, the strong pore pressure perturbations (max Δp∼30% of the background normal stress) of PM2 are sufficient to induce large ruptures greater than the McGarr and Barbour (2017) limit in our simulations (Figure 3, top panel). Figure S11 show perturbed versus nonperturbed moment for several roughness/stress combinations. As with Figure 3, at higher roughness (Figures 4c and 4d), the maximum size of perturbed events is controlled primarily by roughness and background stress, and secondarily by the injection-induced stress perturbation. In particular, for high stress and high roughness, the largest perturbed event (i.e., out of the whole population of events with the same stress conditions and fault roughness) is less than four times larger than the largest nonperturbed event out of the whole population. The perturbation has a stronger impact on rupture size at low roughness. At low stress and roughness (Figure 4a), strongly perturbed events (PM2) tend to be much larger (by more than an order of magnitude in moment) than nonperturbed events, while moderate pressure changes (PM1) result in little difference between perturbed and nonperturbed ruptures.

Discussion
On flat faults, we find empirically (Figure 1) that the criteria for when ruptures exceed the pressurized zone is related to the ratio of the background shear stress and τ pulse : where the bar in σ 0 emphasizes that this is the effective normal stress. Zheng and Rice (1998) showed that faults for which f b ≈ τ pulse =σ 0 could sustain pulse-like ruptures, while Norbeck and Horne (2018) showed that if this criteria is met only locally inside of a pressurized zone, ruptures would be limited by the spatial extent of the zone. Replacing f D in their slip-weakening simulations with τ pulse =σ 0 as a modified criteria, our results qualitatively agree with this conclusion.
In contrast to flat faults, on rough faults (with the parameter ranges, we have considered: 10 −3 < α < 10 −2 , σ 0 ∼ 100 MPa, Δp∼1−10 MPa), the pore pressure perturbation plays a less important role compared to stress perturbations from geometry. Comparison of rupture magnitudes with those predicted by the McGarr and Barbour (2017) relationship indicates that ruptures are not limited by the volume injected; either ruptures arrest due to local high-strength patches, or ruptures exceed the hypothesized boundary. The exception is at low roughness and low background stress, where pressure decay may result in ruptures arresting in some cases (Figure 3, low stress PM3 ruptures; cf. Figure 1f). These results suggest that the role of the perturbation in limiting rupture size is secondary to that of the in situ stress level and heterogeneity.
The results shown in Figure 3 demonstrate that stress heterogeneity arising from fault roughness exerts primary control on stopping ruptures. However, the spatial distribution of the pore pressure perturbation clearly plays an important role. Comparing PM1 with PM3 ruptures, which have identical injected volume, PM3 ruptures can reach larger size than PM1 ruptures regardless of stress and roughness, and can be larger than PM2 ruptures at high roughness. This may be because the higher available stress drop from the perturbation distributed over a smaller region is not able to overcome the resistance to slip of very poorly oriented fault segments. Thus, the pore pressure perturbations does impact rupture size, but not in the simple manner suggested by equation 1. Instead, the preexisting stress state, including both the mean value and the heterogeneity in stress and interactions with the spatial distribution and magnitude of the pore pressure perturbation to impact rupture size.
The results presented in this study demonstrate that the addition of pore pressure to a given background stress state encourages larger ruptures. However, the results do not address whether the pore pressure distributions considered in this study are realistic in natural settings. For example, perhaps events in Figure 3 exceeding the moment limits of equation 1 would have nucleated a smaller event at a lower pore pressure perturbation. While it is possible to reach high pore pressure consistent with PM2 in localized areas around an injector (Häring et al., 2008), this level of perturbation would not be expected at large depths and/or distances from the injector. Thus, care must be taken in interpreting the results. In the simulations presented in this study, the initial stress for both the perturbed and unperturbed models is identical prior to introducing the pore pressure perturbation, and the initiation stress pulse is also the same. These conditions are not reflective of rupture initiation in nature, where pore pressure increases in time and stress and slip velocity co-evolve on the fault interface. However, no events at low stress exceed the hypothesized limits without additional pore pressure, so the artificial initiation alone is not sufficient to produce large events.
Future research should address the limitations of this study and focus on sequence simulations of induced earthquakes that account for nucleation and aseismic slip processes explicitly, and allow rupture to occur naturally, rather than artificially imposing a particular pressure perturbation and comparing rupture size. Simulations that account for both gradual pressure build-up as well as the dynamic effects that occur during rupture are required to fully resolve how stress and frictional strength change throughout the earthquake cycle and determine whether the results presented here are relevant in more realistic scenarios.

Conclusions
We have conducted an extensive set of simulations to explore how injection-induced pore pressure and poroelastic stress changes impact the size of dynamic ruptures on rough faults. We find that rupture size is not limited by injected volume except when roughness, background stress, and the pressure perturbation are all low. Events can grow beyond the pressurized zone and exceed published magnitude limits if τ b > τ pulse or the pore pressure perturbation is large. Higher pore pressure perturbations tend to result in larger ruptures; however, at low background stress and high roughness events never grow as large as published limits. Only in the limited case of low to no roughness and low background stress (τ b ≤ τ pulse ) do events appear to ever be limited in size by the size of the perturbed region. Instead, the results indicate that rupture size is primarily controlled by the in situ stress level and heterogeneity and only secondarily by pressure. This is likely partly due to the stress ratio on geometrically rough faults varying up to 30-70% from the background level for the parameter ranges considered here, compared to 15% or less for the modeled pressure-induced perturbations. Future research is required to determine whether our results hold for naturally nucleated earthquakes, but at present we suggest that, once nucleated by fluid injection, induced earthquakes are not required to stop at the boundaries of the pressurized region.

Data and Resources
The code for FDMAP is available from https://bitbucket.org/ericmdunham/fdmap. Data from the simulations is available from Maurer (2020), last accessed April 13, 2020.