Understanding rate effects in injection-induced earthquakes

Understanding the physical mechanisms that underpin the link between fluid injection and seismicity is essential in efforts to mitigate the seismic risk associated with subsurface technologies. To that end, here we develop a poroelastic model of earthquake nucleation based on rate-and-state friction in the manner of spring–sliders, and analyze conditions for the emergence of stick-slip frictional instability—the mechanism for earthquakes—by carrying out a linear stability analysis and nonlinear simulations. We find that the likelihood of triggering earthquakes depends largely on the rate of increase in pore pressure rather than its magnitude. Consequently, fluid injection at constant rate acts in the direction of triggering seismic rupture at early times followed by aseismic creep at late times. Our model implies that, for the same cumulative volume of injected fluid, an abrupt high-rate injection protocol is likely to increase the seismic risk whereas a gradual step-up protocol is likely to decrease it.

Derivation of the poroelastic spring-slider equations Frictional evolution Frictional evolution is modeled by the rate-and-state constitutive laws, which are capable of reproducing a wide range of observed seismic and aseismic fault behaviors ranging from preseismic slip and earthquake nucleation to coseismic rupture and earthquake after slip [1,2,3,4]. These laws propose that frictional shear stress τ can be described as where Σ is the effective normal stress (the difference between total normal stress and pore pressure), µ is the coefficient of friction, V is the slider's velocity or slip rate, and Θ is a state variable with the physical interpretation of the fractional contact area that is associated with time dependent creep [5]. It is also related to the age of asperity contacts.
We adopt Ruina's [3] slip law for the coefficient of friction because it fits experimental data at variable normal stress better than Dieterich's [1] aging law [6], whereâ is an experimentally derived parameter, V * is a normalizing velocity, and µ * is a constant appropriate for steady-state at velocity V * . Laboratory experiments on dry rocks show that a step change in normal stress results in a sudden change in the coefficient of friction followed by a displacement-dependent decay back toward the initial steady-state value [7,5,8,9]. Linker and Dieterich [5] interpret this as due to a normal stress effect on the state variable and propose to model the magnitude of the sudden change asαΣ/Σ.
Although the model is based on step test experiments, it captures, at least qualitatively, the pressurization-weakening effect on the coefficient of friction observed from ramping experiments by Olsson [7]. He performed laboratory tests in which the normal stress was increased at constant rate while the load point speed was held constant. He found that shear stress is a function of the normal stress rate. When the normal stress rate was increased by 10 during steady sliding, the rate of increase of shear stress with normal stress (coefficient of friction) decreased by a factor of two-a significant effect.
Here, we combine the proposed state evolution model with the effective stress principle to getΘ where d c is the characteristic sliding distance required to replace the old contact population with a new one,b is a constitutive parameter, andα is a scaling factor. Theoretical and laboratory studies for a sudden change in normal stress show thatα ranges from 0 to µ [5,10,11], but more studies are needed to determine the value ofα for a gradual change in normal stress. From momentum balance of forces acting on the slider, the equations of motion of the system evolution at variable effective normal stress becomė where T = 2π m/k s is the vibration period of the analogous freely slipping system [12].

Poroelastic coupling
To obtain a physical evolution of effective stress on the frictional surface, we couple it with a poroelastic model of pore pressure and rock deformation. Starting with the principle of mass conservation, we specify the change of mass from fluid diffusion to be ∆m diff , mass accumulation due to rock expansion or fluid compressibility ∂ ∂t (ρV f )∆t, and injection source term to beQ∆t. We assume that both the fluid and rock matrix are compressible [13], and so mass balance leads to where the change in fluid mass due to pressure diffusion can be written using Darcy's law where η is fluid dynamic viscosity, k is permeability, and L is the pressure diffusion length.
The mass accumulation term can be expressed as where H is the current height of the slider, ρ 0 is the initial fluid density, c f is fluid compressiblity, H 0 is the initial height of the slider, and W is the position of the piston.
When fluid is injected into a rock that is free to deform in the direction orthogonal to sliding, the addition of mass induces an increase of volume equivalent to where V f,0 is the initial fluid volume. We then derive an expression for rock deformation W from force balance, while using the convention of compression positive, where k n is the normal spring stiffness and Σ 0 is the initial effective stress. We further approximate the mass accumulation term to Note that we consider that the total stress Σ is analogous to overburden stress in the earth, and is therefore constant in time. By substituting Eqs. (7)-(10) into Eq. (6), we find that pore pressure satisfies a diffusion equation that leads to transient behavior at early times and steady-state behavior at late timeṡ where η is fluid dynamic viscosity (η = νρ), Q is the volumetric injection rate per unit area (Q =Q/ρA), and k eff n = (1/k n + c f H 0 /A) −1 is an effective stiffness somewhat equivalent to the uniaxial bulk modulus or the reciprocal of the uniaxial specific storage per diffusion length in a continuum [13]. Since the slider has a unit base area (A = 1), the evolution of the pore pressure as a result of fluid injection followṡ

Governing equations in dimensionless form
The equations describing the dynamic motion of the spring-poroslider system ( Figure. 1B) with an evolving pore pressure, in dimensional form, arė Choosing the following characteristic quantities: , θ c = µ * , and t c = d c /V * , the equations describing the dynamic motion of the system, in dimensionless form, becomė r = t c k eff n , and q = Q/p c . The parameter κ is the normalized shear stiffness, and a, b, α are normalized frictional parameters. The parameter is the normalized oscillation period or ratio of inertial to state-evolution timescales, which may range from 10 −8 to 10 −6 depending on rupture diameter and shear wave speed. The parameter c is the normalized diffusivity or ratio of the pore-pressure to the state-evolution timescales, which may range from 10 −4 to 10 1 depending on reservoir permeability, bulk modulus (storativity), and well-fault distance. The parameter rq is the normalized injection rate, which may range from 10 −5 to 10 −1 depending on injection rate and reservoir size.

Linear stability analysis
The stability of steady frictional sliding to small perturbations in velocity depends on the state of pore pressure. For a constant pore pressure, linear stability analysis leads to the steady-state stability condition by Ruina (1983) [3], which shows that the slider is stable when the dimensionless spring stiffness of the loading system exceeds a critical value given by and it is unstable otherwise. Pore pressure, however, is not constant in time, and its evolution depends on the injection rate and on the poroelastic and hydraulic parameters of the rupture (Eq. (21)). A common approach to stability analysis with time-varying state variables is the use of the quasi-steady-state approximation [14,15,16,17]. Using this approach, we freeze time in the pore pressure solution and then perform a linear stability analysis of the spring-poroslider system at a fixed pore pressure.
The equations describing the quasi-static motion of the spring-poroslider system evolving at variable pore pressure p, in dimensionless form, arė Equations (26) If the real part of the roots λ i are negative for all i, perturbations from the quasisteady-state are damped and the system is stable. If the real part of the roots λ i are positive for some i, then perturbations grow exponentially and the system is unstable. At (λ i ) = 0, we find that the dimensionless critical stiffness is If we express the instability condition (Eq. (29)) in terms of a critical dimensionless injection rate, above which an earthquake is induced, we obtain .

Nonlinear simulations
To validate our analytical instability criterion (Eq. (29)), we simulate the fully dynamic equations of motion of the spring-poroslider system (Eqs. (18)- (21)) with the following initial conditions Analytical vs numerical estimates of critical stiffness As a whole, Supplementary Figure 2 shows that the analytical critical stiffness (Eq. (29)) is in good qualitative agreement with the numerical simulation results. To validate our instability criterion quantitatively, we compare our instability criterion against estimates obtained empirically from the fully dynamic nonlinear simulations ( Supplementary Figure 3). The analytical estimate (blue) works well when the growth rate of perturbations is large compared to the growth rate of the pore pressure. Initially, when pore pressure grows rapidly, estimates differ slightly, but they become indistinguishable at late times, when pore pressure changes relatively slowly. We suspect that the small difference in estimates at early times is due to, at least partially, the use of the QSSA in our analysis, which we discuss in detail in Supplementary Note 6.

Quasi-steady-state approximation
The quasi-steady-state approximation, in general, is an approach to simplify dynamic systems of ordinary differential equations with an initial fast transient, after which some of the dependent variables can be assumed to be in steady-state with regard to the other slowly evolving dependent variables [16]. In particular, the QSSA is a good approach to use in our analysis because it allows us to study the stability of steady frictional sliding to small perturbations in velocity while pore pressure is evolving. The sliding velocity and the velocity-dependent part of the state variable are in steady state with respect to the pore pressure. Here, we analyze the QSSA in the context of singular perturbation theory following the analysis by Segel and Slemrod (1989) [16], identify the small parameter(s) necessary for the validity of the QSSA, and quantify the error associated with it.

Reduced dimensional equations
As shown in Supplementary Note 1, the dynamics of our poroelastic spring-slider model is governed by a system of four coupled nonlinear ODEs. Under quasi-static loading, velocity is the fastest evolving variable of the system and it responds instantaneously (negligible inertia) to small perturbations. Thus we can focus our analysis on a reduced system of ODEs at steady-state velocity V = V 0 , with initial conditions Timescales As a first step in the analysis, we estimate the fast timescale t Θ of the pre-steady-state period and the slow timescale t P for the evolution of pore pressure. To estimate t Θ we make the approximation P ≈ P 0 in Eq. (35). The solution for the state variable becomes ). Subsequently, we take Since the pore pressure evolution is independent of the evolution of the state variable Θ, we estimate t P by solving Eq. (36) with the initial condition P (0) = P 0 to obtain Similarly, we take

Scaled dimensionless equations
During the pre-steady-state, it is reasonable to scale time by t Θ , where the dimensionless time τ is given by with initial conditions where we have defined b =b/µ 0 , α =α/µ 0 , and q = Q/p 0 . After the pre-steady state, the QSSA is assumed to be valid and t P becomes a reasonable timescale. We introduce a new dimensionless scaled timet byt with which the scaled dimensionless governing equations become

Singular perturbation
Approximate solutions can now be obtained by methods of singular perturbation theory [18], for 0 < t Θ /t P

A solution of Eqs. (44)-(45) is obtained of the form
where where Note that Eqs. The results of this analysis remain valid for general initial conditions Θ(0) = Θ i and P (0) = P i [16], where Θ i ranges from 0 to Θ ss and P i ranges from P 0 to P ss . Linearizing the spring-poroslider system about the true steady-state, where thus yields the Ruina (1983) stability condition [3].

QSSA validity conditions
A necessary aspect of the QSSA is that the duration of the pre-steady-state period is much shorter than the characterstic time for the pore pressure evolution. An essential condition for the QSSA to be valid after the pre-steady state is therefore t Θ t P , Note that the initial condition P (0) = P 0 is reasonable for the QSSA only if there is a negligible relative change |∆P/P 0 | in pore pressure during the pre-steady state. We estimate |∆P/P 0 | by An additional condition for the validity of the QSSA is therefore Recall that parameter c is the normalized diffusivity or ratio of the pore pressure to the state evolution timescales and parameter rq is the normalized injection rate, where r = t Θ k eff n and q = Q/P 0 .

Error estimates
This analysis shows that using the QSSA to study the stability of steady frictional sliding to small perturbations in velocity with an evolving pore pressure is justified when conditions Eqs. (61) and (63) are met. In other words, if in a timescale t Θ sliding reaches steady state with a constant pore pressure, then assuming that sliding is in a quasi-steady state with a changing pore pressure is valid when the pore pressure change occurs on a time scale t P that is long compared to t Θ and ∆P | t Θ is small compared to P 0 .
Therefore, we expect that the accuracy of our instability criterion depends on dimensionless parameters c and rq. Here we evaluate the error in the analytical estimate of the critical stiffness required to trigger the first slip event ( Supplementary Figure 4). We indeed find that the error decreases as c or rq decrease. It becomes small (< 15%) when the normalized diffusivity and normalized injection rate reach small values (c ≤ 5 × 10 −2 , rq ≤ 5 × 10 −3 ). It is also interesting to note that the QSSA validity may be extended to instances where c is of order one provided that rq is significantly smaller than one (green curve).

Application to the Denver earthquakes
To bridge the gap between the analysis of the idealized spring-poroslider model and the real world, we express our instability criterion in dimensional form, and identify values of dimensionless parameters c and rq that correspond to real-world settings.
In dimensional form, the critical stiffness from the linear stability analysis is or, equivalently, where the termb−â represents the original velocity weakening effect and the dimensionless termα(d c /V 0 )Ṗ /(Σ − P ) represents an additional weakening effect from fluid pressurization. Note that this pressurization term is maximum at early times and is approximately equal to rq.
The 1960s Denver earthquakes is a good example of a real-world setting, where it is well-documented that injection of wastewater into the fractured Precambrian granite gneiss underneath the Rocky Mountain Arsenal triggered the earthquakes and where injection rate is directly related to the frequency of earthquakes [19,20]. The reservoir spans a depth interval from 3.7 to 7 km below the surface. Experimental data on granite at this depth shows velocity weakening behavior (b−â in the range 0.002 to 0.005, µ 0 = 0.7 to 0.75) [21].
To identify values of dimensionless parameter c that correspond to this setting, we estimate the state evolution timescale t Θ and the pore pressure evolution timescale t p : We find that t Θ ranges from 10 days to 4 months based on field data of the characteristic slip distance and loading rate (d c = 10 −3 to 10 −2 m, V 0 = 10 −9 m s −1 ) [22,4].
The pore pressure evolution time scale t P can be transferred from the poroslider model to field settings, We find that t P is approximately 2 years based on a reservoir analysis of the Denver earthquakes, with transmissivity T = 10 −5 m 2 s −1 , storativity S = 10 −5 , and characteristic In a similar manner, we identify values of dimensionless parameter rq. We translate this quantity from the poroslider model to field settings: and evaluate values based on the reservoir analysis and injection data, with reservoir pressure P 0 = 30 MPa, reservoir width W = 3 × 10 3 m, and field injection rate Q w = 2 to 9 million gal mo −1 [19,23].
Therefore, reasonable estimates of c and rq for this setting would be in the order of 10 −2 to 10 −1 and 10 −3 to 10 −1 , respectively. Note that both estimates are much smaller than one, and thus meet the QSSA validity conditions (Supplementary Figure 5).
Having determined the validity of the QSSA analysis to this setting, we now assess whether pressurization rate effects were likely significant during fluid injection leading to the Denver earthquakes. In dimensionless form, the critical stiffness κ crit is given by Eq. (S29). Prior to fluid injection, the pore pressure is constant and the critical stiffness, is estimated to be 0.005 (b − a = 0.003 to 0.007, σ − p 0 = 1). Shortly following the start of fluid injection, the pore pressure increases rapidly and the dimensionless critical stiffness takes the form: This results in an increase in critical stiffness at early times of around 300% (α = 1, v 0 = 1, b − a = 0.005, rq = 10 −2 ), thus indicating that the additional weakening effect from fluid pressurization is likely significant in this setting.

Phase diagram of injection-induced seismicity
To study the influence of reservoir properties on injection-induced seismicity, we simulate the occurrence of earthquakes as a function of dimensionless injection rate rq = When the normalized diffusivity is higher than one, the dimensionless pressure diffusion time is less than one. Pore pressure reaches steady-state on a very short time scale, and so the rate of change in pore pressure is negligible. To quantify this, we writė It is helpful to express the condition for instability (Eq. (29)) in terms of rq crit so that an earthquake is triggered if the injection rate is higher than a critical value given by We find that the injection rate required to trigger an earthquake is proportional to diffusivity c, which explains the simulation results in Supplementary Figure 6: In the regime of high diffusivity c > 1, q crit ∼ c. Accordingly, earthquakes are more easily triggered when fluid is injected into a low permeability reservoir rock than a high permeability, for a fixed effective stiffness.
Conversely, when the normalized diffusivity is lower than 0.001, the dimensionless pressure diffusion time is more than 100. Pore pressure stays near the initial transientstate throughout the simulation period and so the magnitude of change in pore pressure is negligible. To quantify this, we writeṗ If we express the condition for instability in terms of rq crit , similar to the high diffusivity case above, we find that the injection rate required to trigger an earthquake is This also explains the simulation results in Supplementary Figure 6: In the regime of low diffusivity c < 0.01, q crit ∼ const. Accordingly, in this regime, earthquake triggering is independent of permeability. These results, as a whole, suggest that reservoirs with high hydraulic diffusivity and low stiffness may be safer sites for fluid injection operations compared to sites with low hydraulic diffusivity and high stiffness.  The error is taken as the ratio of the maximum difference between the analytical and numerical estimates of the critical stiffness required to trigger the first slip event to the numerical estimate. The dimensionless parameter c is varied by varying permeability k, and rq is varied by varying injection rate Q.
Supplementary Figure 5: Analytical vs. numerical estimates of critical stiffness κ crit for a case with normalized diffusivity of c = 1 × 10 −2 and normalized injection rate of rq = 1 × 10 −2 . Solid line represents the analytical estimate, and open circles represent numerical ones. Blue represents the combined rate and magnitude effect (setting α = 1), and green represents the magnitude effect (setting α = 0). The black horizontal line represents κ crit prior to fluid injection. The earthquake likelihood is computed as a normalized initial critical stiffness for a given injection duration (κ crit | T /κ crit | T min ). The injection duration is computed as the ratio of a fixed total injection volume to the injection rate (T = V /q).