Analysis of Fluid-Injection-Induced Seismicity Using a Dynamic Sliding Model Incorporating the Rate- and State-Dependent Friction Law

Seismicity events triggered by fluid injection have become a significant concern for various subsurface activities, including geothermal energy development, hydraulic fracturing stimulation, and carbon dioxide sequestration. Presently, with scientific knowledge and monitoring technologies, it is impossible to discriminate between human-induced and natural tectonic earthquakes. However, the injection of fluids at high pressures and volumes can have a certain risk of triggering subsurface seismic events, potentially damaging infrastructure and public safety. Understanding the underlying mechanisms of triggered seismic events is essential for mitigating associated risks and developing effective seismicity management strategies.

Seismicity can be modulated by stress induced by anthropogenic activity in the crust. Foulger et al. (2018) compiled and presented more than 700 possible human-induced earthquake events recorded between 1868 and 2016. Despite the existence of scientific uncertainties, it is presumed that these earthquake events are induced by a wide range of human activities, including the impoundment of water reservoirs; construction of tall buildings; coastal engineering; quarrying; extraction of groundwater, coal, minerals, oil, gas, and geothermal fluids; excavation of tunnels; and waste disposal injection. Recent seismic events related to the oil and gas industry include seismic swarms in Oklahoma, linked to wastewater disposal from hydraulic fracturing operations (Keranen et al. 2014), and Permian Basin earthquakes in Texas, related to the underground disposal of produced water (Skoumal et al. 2020, 2021; Zhai et al. 2021). The two largest triggered earthquakes reported in Oklahoma are the 2011 Mw 5.7 Prague (USGS 2011) and the 2016 Mw 5.8 Pawnee (USGS 2016) events, and two of the five largest earthquakes (Mw 5.2 and 5.4 events) in Texas occurred in the Permian in 2022 (Presley 2023).

In addition to earthquake events linked to wastewater disposal, enhanced geothermal system projects face challenges related to induced seismicity. In 2017, an Mw 5.5 earthquake occurred near an enhanced geothermal system drillsite in Pohang, South Korea, injuring 135 residents and causing more than 75 million USD in direct damage (Lee et al. 2019; Yeo et al. 2020). In addition, gas storage-induced seismicity occurred at the Hutubi Underground Gas Storage site in Xinjiang, China, where the largest earthquake reached Mw 3.0 in 2013, and the epicenter was located at a depth of approximately 4  km (i.e., slightly deeper than the reservoir formation; Zhou et al. 2019). Another case of gas storage-induced seismicity was reported at the Castor Underground Gas Storage project in offshore Spain, which has been considered to have triggered a seismic sequence of more than 1,000 earthquakes, peaking with three Mw 4-scale earthquakes in 2013 (Gaite et al. 2016). These fluid injection-induced seismicity events highlight the importance of understanding the mechanisms underlying the triggered seismicity events and developing reliable models to predict and mitigate induced seismicity.

Possible mechanisms of fluid injection-induced seismicity include the direct fluid pressure effects of injection on the fault and changes in underground stress due to fluid extraction or injection (Ellsworth 2013). Elastic strain energy is released when a fault slips (i.e., an earthquake). A fault slides when the shear stress acting on the fault plane exceeds the contact strength. Increasing the shear stress, reducing the normal stress, and/or elevating the pore pressure can cause the fault to fail, triggering earthquake nucleation. The relative importance of these mechanisms depends on specific injection conditions and geological settings. Recent studies have shown that rate- and state-dependent friction laws play critical roles in controlling the dynamics of fluid-injection-induced seismicity (Alghannam and Juanes 2020). This friction law captures the rate dependence of frictional strength and provides a better representation of the complex dynamics of the fault system.

Various approaches have been proposed to model the dynamics of fluid-injection-induced seismicity. These approaches include empirical seismicity rate models, numerical simulations of coupled-flow geomechanics, and rate and state friction models. Empirical seismicity rate models assume the Gutenberg-Richter law that expresses the relationship between the magnitude and the total number of earthquakes in a particular region in a statistical sense (Shapiro and Dinske 2009). Seismicity rate models may be useful for short-term earthquake forecasting but do not address rupture dynamics. Numerical simulation approaches have been proposed using various discretization methods, such as the FEM (Baisch et al. 2010), boundary element method (Aoki et al. 2019), and discrete element method (Yoon et al. 2013). Numerical models can address the physical mechanisms underlying the relationship between fluid injection and fault slip. However, as most numerical simulation models accommodate a simple friction law, they do not address whether a fault slips seismically or aseismically.

The rate- and state-dependent friction laws are incorporated into numerical analysis to evaluate earthquake cycles (Allison and Dunham 2018; Hirose and Hirahara 2002). In particular, the poroelastic spring-slider model incorporating the rate- and state-dependent friction laws proposed by Alghannam and Juanes (2020) can be used to analyze the conditions and criteria of stick/slip frictional instability observed during the subsurface fluid injection. Their model is a useful prediction tool for determining a safe injection scenario; however, its application is limited to simple pore-pressure changes/distributions under homogeneous stress fields.

In this study, we extended the spring-sliding model with the rate- and state-dependent friction laws proposed by Alghannam and Juanes (2020) to analyze the induced seismicity occurring during various subsurface fluid injection activities. The augmented seismicity analysis model incorporates the rate- and state-dependent friction laws into a 2D finite element model, accounting for the variations in pore pressure and stresses in the formation.

Alghannam and Juanes (2020) proposed a poroelastic spring-slider model to elucidate the mechanism of seismicity triggered by fluid injection. This model represents a fault and reservoir using two springs and a slider that can be infused with fluid, as illustrated in Fig. 1 . When the horizontal spring is pulled by external forces, additional fluid is injected into the slider. This injection process is analogous to the fluid injection into a reservoir.

The governing equations for this model are as follows (Alghannam and Juanes 2020):

The descriptions of the physical quantities in these equations are listed in Table 1 .

The velocity of the slider $V$ represents the velocity of the fault slip. The loading velocity $V0$ represents the strain rate of the surrounding rock. The shear stiffness $ks$ represents the elastic interaction between the fault and the surrounding rock. Eqs. 1 and 2 represent the equations of motion of the system evolution derived from the momentum-balancing forces acting on the slider, whereas Eq. 3 represents the time evolution of the state variable that defines the frictional contact area associated with time-dependent creep. Finally, Eq. 4 is a diffusion equation describing the evolution of the pore pressure.

Eqs. 1 through 4 are nondimensionalized as follows:

The descriptions of the dimensionless physical quantities in these equations are listed in Table 2 .

The conversion relationships between the physical quantities and dimensionless physical quantities are summarized in Table 3 .

In Table 3 , the characteristic time $tc$ and the characteristic shear stress $τc$ were introduced to formulate the conversion relationships. The characteristic quantities are calculated as follows:

This study uses the embedded SIRK method to compute numerical solutions to Eqs. 5 through 8. The embedded SIRK method is a combination of the semi-implicit and embedded Runge-Kutta methods.  Appendix A includes a detailed explanation of these methods. To obtain stable numerical solutions for stiff ordinary differential equations using the explicit Runge-Kutta method, extremely small step sizes must be used, thereby significantly increasing the computation time. To mitigate this challenge, we chose the SIRK method, which possesses a larger absolutely stable region compared to the explicit method. This larger stable region permits the use of larger step sizes, reducing the computation time. Furthermore, the embedded Runge-Kutta method is used to automatically adjust the step sizes.

Eqs. 5 through 8 were solved using the embedded SIRK method. To ensure the validity of the numerical solutions of $u$⁠, $v$⁠, $θ$⁠, and $p$⁠, the calculation results are compared with those obtained through an alternative solution method in MATLAB® (The MathWorks Inc., Massachusetts, USA), specifically the numerical differentiation formulas.

In this validation, we used the initial conditions below, as prescribed by Alghannam and Juanes (2020).

Initially, the slider is sliding at a constant velocity $v|t=0=1$⁠, which is equal to the loading velocity $v0|t=0=1$⁠. Subsequently, the loading velocity suddenly increases to $v0|t>0=1.5$⁠. The other parameters used in this calculation are summarized in Table 4 .

Fig. 2  presents the numerical solutions obtained from the embedded SIRK method, which are depicted in blue, and those obtained from numerical differentiation formulas, which are depicted in orange. A comparison between the blue and orange curves shows that the embedded SIRK method yields accurate calculation results, as illustrated in the figure.

Under this condition, several stick/slip phenomena were observed, with one of them reaching a velocity of $107$⁠. These stick/slip phenomena depend on the parameters and initial conditions. Therefore, we will discuss the dependency of these values on the occurrence of stick/slip phenomena in the next section.

Parametric studies were conducted to fully understand the effects of the model parameters and initial conditions in the governing equations. The results of the parametric studies are shown in Fig. 3 , which illustrates the maximum velocity of the slider (⁠$vmax$⁠). As shown in the figures, the maximum velocity of the slider changes rapidly within the specific boundaries. This observation allows us to separate a rapid sliding region from a slow sliding region using $v=105$ as a cutoff value. In the later discussions, we refer to sliding with $v>105$ as rapid sliding. Mathematically, this rapid sliding is regarded as an unstable condition, as determined through an examination of the eigenvalues of the linearized governing equations.

Fig. 3a  shows that a high injection rate leads to rapid sliding. This tendency can be attributed to the increase in pore pressure caused by the high injection rate. Furthermore, the same figure shows that high diffusivity impedes rapid sliding. This trend can be explained by the rapid dissipation of the injection fluid at high diffusivity, which impedes the buildup of pore pressure and inhibits the occurrence of rapid sliding. Notably, the linear boundary within the large diffusivity (⁠$log10⁡c>0$⁠) and horizontal boundary within the small diffusivity (⁠$log10⁡c<−2$⁠) were observed, which is consistent with the findings reported by Alghannam and Juanes (2020).

Fig. 3b  shows that a large fault stiffness impedes rapid sliding. A larger value of fault stiffness leads to sliding even with the short elongation of the spring, resulting in early sliding with a small discharge of elastic energy. Moreover, the figure indicates that a high normal stress acting on the fault leads to rapid sliding. According to Amonton’s law, a high normal stress results in a corresponding increase in the maximum static frictional force. Consequently, a larger shear force is required to initiate slip. However, the requirement of a larger shear force does not preclude rapid sliding in the spring-slider system. The accumulation of elastic energy within the system continues until the system experiences a slip, ultimately resulting in rapid sliding.

Fig. 3c  shows that a large value of a impedes rapid sliding. In this case, as the velocity of the slider increases, the increase in the frictional coefficient will be large, leading to the inhibition of the increase in slip velocity. On the other hand, the same figure shows that a large value of b leads to rapid sliding. In this case, as the velocity of the slider increases, the state variable decreases, resulting in a large decrease in the frictional coefficient. This, in turn, leads to an increase in slip velocity. The figure also shows that rapid sliding is observed only in the region where $a-b<0$⁠, while it is not observed in the region where $a-b>0$⁠. These preferences are consistent with velocity strengthening (⁠$a-b>0$⁠) and velocity weakening (⁠$a−b<0$⁠). In velocity strengthening, the frictional coefficient increases with slip velocity, leading to stable sliding. In velocity weakening, on the other hand, the frictional coefficient decreases with slip velocity, leading to rapid sliding.

Fig. 3d  shows the change in the vibration period within the range $10−8<ε<10−6$ is not a key factor in triggering rapid sliding. The range for the vibration period was chosen in accordance with the typical values from Alghannam and Juanes (2020). However, looking closer to Fig. 3d , a higher vibration period leads to a gradual decrease in the maximum velocity within the rapid sliding region. As shown in Table 3 , the vibration period is expressed as

where $m$ is the mass of the slider (kg). Assuming $ks$⁠, $dc$⁠, and $V*$ are constant, the vibration period depends solely on the mass of the slider. Consequently, this trend is interpreted as follows: A larger mass of the slider results in a longer vibration period, leading to a decrease in the sliding velocity. In addition, the figure shows that a large value of $α$ leads to rapid sliding. Basically, the high effective normal stress leads to a decrease in the state variable due to the initiation of newly formed asperities. The value of $α$ is proportional to the ratio of the decrease in the state variable to the increase in the effective normal stress. Therefore, a higher value of $α$ leads to a larger decrease in the state variable, leading to a larger decrease in the frictional coefficient, which, in turn, results in rapid sliding.

Fig. 3e  shows that a low initial pore pressure leads to rapid sliding. This trend is consistent with the finding of the effect of the normal stress acting on the fault on the velocity of the slider. A low initial pore pressure means a high effective normal stress, prompting rapid sliding with a larger release of elastic energy. Furthermore, the figure shows that a high loading velocity leads to rapid sliding. A high loading velocity accumulates elastic energy quickly, thereby resulting in rapid sliding.

The occurrence of induced seismicity associated with the multistage fracturing of horizontal wells in shale reservoirs is widely recognized. To validate our seismicity prediction model, we conducted a synthetic case study to analyze seismicity events during hydraulic fracturing in a shale reservoir. Field observations from various shale fields (GTI-NETL 2020; Rutledge et al. 2004; Maxwell and Cipolla 2011) indicate that the timing and magnitude of microseismicity events monitored during hydraulic fracturing can vary because of factors such as reservoir heterogeneity, mechanical properties of the rock formation being fractured, injection rate and pressure, completion design, and in-situ stress conditions. Typically, microseismic events occur in clusters, with multiple events occurring within a short period.

A review by Warpinski et al. (2012) summarized thousands of fracture treatments and associated microseismic measurements available from major shale basins in North America. The study found that microseismic events observed during hydraulic fracturing in shale formation range from −3.0 Mw to 1.0 Mw and are usually too small to be felt on the surface. Tan et al. (2021) analyzed field data obtained from a Hydraulic Fracturing Test Site-2 experiment and found that microseismic events occurred immediately after the injection started downhole and continued to occur during the injection.

In our study, we compared the timing and magnitude of microseismic events simulated using a dynamic seismicity prediction model with field observations.

Fig. 4  illustrates the setting of this case study. As shown in the figure, we consider a hydraulic fracturing treatment originating from a horizontal well within a shale formation. The horizontal fracture is propagating, because we assume that the direction of the minimum principal stress is vertical. In this case study, we examine the effect of fluid diffusion exclusively because the normal stress and loading velocity are assumed to be constant in the governing equations.

First, the model input parameters used in the governing equations were calculated based on typical values in shale reservoirs.  Appendix B includes an explanation of how to convert field data into dimensional physical quantities and Table B-1  summarizes the input data used in the calculations. Table B-1  enumerates nine values for the pressure diffusion length (⁠$L$⁠) and five values for the length of the sliding region (⁠$R$⁠), all of which will be used in the subsequent calculations. By using the embedded SIRK method to solve the governing equations, we computed the time required to induce a first rapid sliding (⁠$tslip$⁠) and the maximum velocity during the first rapid sliding (⁠$vmax$⁠) with the chosen values of $L$ and $R$⁠. It should be noted that the pressure diffusion length $L$ can be considered as the distance from the injection point, as illustrated in Fig. 4 .

The calculation results are shown in Figs. 5a and 5b . In Fig. 5a , $R$ is set at 8.0 m, while in Fig. 5b , $L$ is set at 100 m. According to Fig. 5a , a shorter distance from the injection point leads to a smaller value for $tslip$ and a larger value for $vmax$⁠. These results are consistent with the fact that a more rapid change in pore pressure is observed in the vicinity of the injection point. On the other hand, Fig. 5b  shows that a smaller length of the sliding region results in smaller values for $tslip$ and $vmax$⁠. Based on Eq. B-3, a smaller length of the sliding region means a larger value of $κ$⁠, leading to early sliding with a small discharge of elastic energy.

In addition, we calculated the displacement due to the first rapid sliding using the change in relative displacement (⁠$u$⁠) during this event. With this value, the moment magnitude of the “seismicity” due to the rapid sliding could be calculated using the following equation.

where $G$ (Pa) is the shear modulus, $D(m)$ is the slip distance due to rapid sliding, $S(m2)$ is the area of the sliding region, $E$ (Pa) is the Young’s modulus, $ν$ (–) is the Poisson’s ratio, and $R(m)$ is the length of the sliding region. The length of the displacement (⁠$D$⁠) can be calculated using the change in relative displacement during the first rapid sliding; the area of the sliding region (⁠$S$⁠) is assumed to be $R2$⁠. Based on Eq. 13, the moment magnitude of the “seismicity” due to the rapid sliding can be calculated within the range of zero to unity. Therefore, the rapid sliding in this case study can be referred to as “microseismicity.” The time required to induce this microseismicity (⁠$tslip$⁠) and the magnitude of the induced microseismicity (⁠$Mseis$⁠) are illustrated in Fig. 6 .

Comparing these values with the general characteristics of microseismicity triggered by hydraulic fracturing in shale reservoirs, as mentioned in the previous section, the magnitude of the induced microseismicity (⁠$Mseis$⁠) was within a reasonable range. However, the time required to induce microseismicity (⁠$tslip$⁠) exhibited significant error. Even in close proximity to the injection point with $L=0.1m$⁠, the calculated value of $tslip$ is 114.39 minutes. Therefore, the effect of the diffusion of pore fluid due to the injection is minimal and insufficient to account for the typical characteristics of microseismicity. This discrepancy may be attributed to three plausible factors—the dimensionless loading velocity $v0$⁠, the dimensionless pore pressure $p$⁠, and the dimensionless normal stress $σ$⁠. As we mentioned above, the poroelastic spring-slider model assumes a constant loading velocity and normal stress. However, these physical quantities vary with time in formations where hydraulic fracturing is conducted. Furthermore, considering the presence of natural fractures in the formation, the rate of increase in the pore pressure should be higher than that calculated using Eq. 8. In the following section, we describe how the model was modified to resolve these issues.

In the above case study, the calculation of $tslip$ and $Mseis$ was limited to 1D cases in which the pore pressure diffusion occurs along the L-axis (Fig. 4). To overcome such limitation, the analysis domain is discretized into a number of elements. Then, the spring-slider model is extended to allow us to calculate the 2D spatial distributions of the values of $tslip$ and $Mseis$ in the center of each element. In each element, the slip direction was predetermined, and we solved the modified governing equations (explained below) to calculate $tslip$ and $Mseis$ for every element. The value for $R$ is set at 8.0 m, while the other parameters that do not need modifications are the same as the above case study. It is assumed that these parameters remain uniform throughout the simulated region.

Herein, we explain the procedures of accommodating the required changes discussed in the above section. First, to accommodate the change in normal stress, the normal stress change was added to the second term on the right-hand side of Eq. 7 as follows:

Second, Eq. 8 (diffusion equation of pore fluid) was removed from the governing equations such that the temporal profile of the pore pressure could be assigned arbitrarily to the formation. This change allows us to evaluate the effect of the existence of natural fractures, even though simulating fluid penetration through these fractures is a challenging task. The effect of natural fractures was represented by the temporal profiles with the rapid increase in pore pressure. Finally, the changes in the loading velocity and normal stress due to hydraulic fracturing were determined based on the numerical simulation results obtained from a linear elastic FEM model. The FEM mesh was configured as shown in Fig. 7a , where the bottom blue elements correspond to the region where the horizontal fracture propagated. This FEM model solves the following elasticity equation:

where $km$ is the plane strain stiffness matrix, $u$ is the nodal displacements, and $f$ is the force components including in-situ stresses, gravity, pore pressure, and surface loads. Although fracture propagation is influenced by injection rate, heterogeneity of rock properties, temperature, and competition between other hydraulic fractures (Li et al. 2015), the fracture propagation rate in this case study was assumed to remain constant at 30 feet per 8 minutes, resulting in a complete 300-ft-long fracture after 80 minutes.

To examine the effect of fracture propagation, the FEM calculations were performed incrementally in 11 consecutive steps, as illustrated in Fig. 7b . The fracture propagation was modeled simply by assigning increased pore pressure in the bottom elements and surface loads caused by the net pressure on the fracture surface. At the initial condition, there is no increase in pore pressure or surface loads. At 8 minutes after the start of fluid injection, there are increases in pore pressure in the lower-left element and surface loads along the line between (0,0) and (30,0). Sixteen minutes later, there are increases in pore pressure in the lower-left element and the second lower element from the left, and surface loads along the line between (0,0) and (60,0). These incremental loading steps were repeated for 80 minutes.

The calculated values of normal stress acting on faults and strain rate in the fault direction at each time are connected by straight lines, as shown in Fig. 7c , and the calculation was performed using these kinds of time profiles. It should be noted that this FEM model calculates the normal stress and strain rate only considering the loads and the expansion due to increased pore pressure and does not solve the diffusion equation simultaneously. This is because pore pressure diffusion in formations with numerous natural fractures is too intricate to be feasibly incorporated in this work. To simply evaluate the effect of pore pressure diffusion through permeable natural fractures, two cases are investigated in this study. The first case, denoted as Case 1, assumes that the pore pressure in the unfractured formation does not change throughout the simulation. This case should arise when the permeability of the shale formation is extremely low, such that the pore pressure change occurs only in the fracture. The second case, referred to as Case 2, involved a linear rapid increase in the pore pressure across the entire region over time. This occurs in a formation where natural fractures exist and connect well with each other in the shale formation, thereby allowing the injection fluid to penetrate deeply into the formation. In either case, Eqs. 5, 6, and 14, with the simulated temporal profiles of the loading velocity and normal stress, and with predetermined and simplified pore pressure profiles rather than Eq. 8 are solved by the embedded SIRK method to evaluate rapid sliding in the prescribed fault direction in each element.

However, there are several limitations in this case study. First, Eq. 15 does not consider dynamic fracture growth. Second, the fluid flow equation is not coupled into the FEM model. In spite of these limitations, we can learn some insights into the spatiotemporal distribution of microseismicity induced by hydraulic fracturing.

We computed the time required to induce a first rapid sliding (⁠$tslip$⁠) in each element as illustrated in Fig. 8a . The maximum value of the color bar was set to 4,800 minutes. Thus, elements with a value of 4,800 minutes indicate that either the time required to rapid sliding exceeds 4,800 minutes or no rapid sliding occurs. Fig. 8b  shows the magnitude of the first rapid sliding (⁠$Mseis$⁠), which means that this rapid sliding can be called microseismicity. The results indicate that microseismic events occur only near the tip of growing hydraulic fractures, and it takes a long time to initiate microseismic events during the fracturing treatment.

Similarly, we computed the time required to induce a first rapid sliding (⁠$tslip$⁠) in each element as illustrated in Fig. 9a . Fig. 9b  shows the magnitude of the first rapid sliding (⁠$Mseis$⁠), which means that this rapid sliding can be called microseismicity. Compared to the simulation results for Case 1, microseismic events were observed in a wider region. In addition, microseismic events occurred shortly after the hydraulic fracture treatment began.

As illustrated in Fig. 8a , the simulation results for Case 1 suggested that microseismicity occurred only in a limited area near the growing hydraulic fractures. In addition, the slip time of the microseismic events was predicted to be very late, taking more than 49 minutes. These simulation results are not consistent with the field observations that microseismic events occurred immediately after the injection started downhole during the fracturing treatment. However, as shown in Fig. 9a , deeper pore pressure penetration through the natural fractures induced microseismicity in a wider area, taking fewer than 5 minutes across all elements.

The field observations showed the immediate occurrence of microseismicity after the start of the injection and the continued occurrence of the microseismic events during the injection. Despite the apparent inconsistency of the results in Case 2 with this temporal distribution, they can be reconciled by considering the following two assumptions. The first assumption is that the injected fluid starts flowing into all elements through natural fractures immediately after the start of fluid injection. The second assumption is that the rate of the pore pressure change was assumed to be constant in space and time. However, under the actual reservoir conditions, due to the reservoir heterogeneity, the timing when the injected fluid starts flowing through each element through natural fractures should vary. Furthermore, the rate of the pore pressure change varies spatially and temporally. In the case where the onset of pore pressure increase is early and the rate of pore pressure is large, the timing of seismicity is expected to be early, and vice versa. Considering these effects due to reservoir heterogeneity, the typical temporal distribution of microseismicity triggered by hydraulic fracturing can be explained.

From the above discussion, we cannot explain the typical temporal distribution of microseismicity by considering only the change in normal stress and strain rate due to fracture propagation. In addition, in Case Study 1, we found that the impacts of fluid diffusion without natural fractures are minimal. Therefore, to explain the temporal distribution of microseismicity, a rapid increase in pore pressure is required. Therefore, we can conclude that a significant number of natural fractures may exist with increased pore pressure owing to the presence of permeable natural fractures, which serve as the principal factor for the observed microseismicity.

Because our model computes the governing equations for each element, the large number of elements leads to prolonged computational time. However, the integration with FEM enables the consideration of heterogeneity, including spatiotemporal distributions of stress and pore pressure. This stands as a significant advantage over the original model proposed by Alghannam and Juanes (2020), which lacks the capability to account for heterogeneity.

However, it should be noted that we did not account for the effects of stress changes due to slip, including stress concentration or the stress shadow effect. In addition, concerning the invasion of pore fluid, we had to develop a more sophisticated model to fully describe the pore fluid invasion through natural fractures, which is also one of the limitations of this study. However, even the simple comparison in this study yields the finding that fluid invasion through natural fractures is one of the key factors triggering microseismicity during hydraulic fracturing in shale formations.

Fig. 10  shows the time required to induce microseismicity (⁠$tslip$⁠) at the elements in contact with the horizontal fracture in Case 2. The four elements in Fig. 10  are labeled as A, B, C, and D from left to right.

Fig. 11a  shows the temporal profiles of normal stress at Elements A, B, C, and D calculated by FEM, while Fig. 11b  shows the spatial distribution of normal stress acting on the faults at 8 minutes after the start of fluid injection.

Figs. 10 and 11a  show a consistent relationship between the time required to induce seismicity (⁠$tslip$⁠) and the rate of normal stress over time in the initial 8-minute period. Fig. 10  also shows that microseismicity occurred during the initial period; Elements A and C had the largest and smallest value, respectively, of $tslip$⁠. In this early stage, the load caused by the 30-ft-long fracture increases the normal stress of Element 30, resulting in the largest value of $tslip$⁠. In addition, a decrease in the normal stress was observed at the tip of the fracture at this early stage, specifically in Elements C and D. (Element B is the tip of the fracture, but no decrease in stress was observed owing to the effect of the load on the lower-left node of this element.) According to Alghannam and Juanes (2020), the likelihood of triggering earthquakes depends largely on the rate of increase in the pore pressure rather than its magnitude. Our findings suggest that both the rate of decrease in normal stress and the rate of increase in pore pressure are important factors in triggering microseismic events. Therefore, the stress decrease at the tip of the fracture is considered to be one of the key factors in tracking the evolution of microseismicity clouds. Considering this effect, time series seismic events can be used to identify the possible fracture geometry and extent of pore pressure penetration in field applications.

However, as reported from the Hydraulic Fracturing Test Site-2 experiment (Tan et al. 2021), microseismic events may represent a subset of the created fracture geometry or stimulated reservoir volume, especially in formations with less rock strength contrast. Our current model only captures microseismic events induced by tensile-opening fractures. Microseismic events caused by other mechanisms, such as slow-slips or shear-opening hydraulic fractures, need to be addressed in future study.

The conclusions drawn from this study are summarized as follows:

• Contrary to the popular understanding derived from Amonton’s law, the poroelastic spring-slider model revealed that a large normal stress on the fault leads to rapid sliding. The accumulation of elastic energy within the system continues until the system experiences a slip, ultimately resulting in rapid sliding.

• The governing equations of the poroelastic spring-slider model were solved numerically using the embedded SIRK method, and the results were compared with those obtained using an alternative solution method in MATLAB to confirm the accuracy of our computations. By using this numerical technique, we calculated the microseismicity triggered by hydraulic fracturing and demonstrated that our simulations produced magnitudes of microseismicity typically observed in real fields.

• The governing equations for the poroelastic spring-slider model were extended to accommodate changes in normal stress, which enabled us to conduct a more comprehensive case study by incorporating the dynamic slider model into a 2D FEM. The FEM analysis demonstrated that the timing of microseismicity was consistent with typical observations in shale fields when a significantly larger pore pressure diffusion into the formation was assumed, indicating the presence of natural fractures in the formation that enhanced its permeability.

• The decrease in stress at the tip of the fracture is considered as one of the factors that triggers seismicity. Therefore, time series seismic events can provide insights into the potential geometry of fractures and the extent of pore pressure penetration, which is useful for practical applications in the field.

 
• a

  dimensionless frictional parameter

 
• $a^$

  experimentally derived parameter relating frictional coefficient to the velocity of the slider

 
• Area

  cross-sectional area of fluid flow

 
• b

  dimensionless frictional parameter

 
• $b^$

  experimentally derived parameter relating frictional coefficient to the state variable

 
• c

  dimensionless diffusivity

 
• d_c

  characteristic slip distance

 
• D

  slip distance due to rapid sliding

 
• E

  young’s modulus

 
• {f}

  force components vector

 
• G

  shear modulus

 
• h

  step size

 
• k

  permeability

 
• k_i

  constant in the Runge-Kutta method for a single ordinary differential equation

 
• k_im

  constant in the Runge-Kutta method for simultaneous ordinary differential equations

 
• [k_m]

  plane strain stiffness matrix

 
• $kneff$

  effective normal stiffness

 
• k_s

  shear stiffness

 
• K

  bulk modulus

 
• $1/Kf$

  fluid compressibility

 
• $1/Kv$

  uniaxial compressibility

 
• L

  pressure diffusion length

 
• m

  mass of the slider

 
• M_seis

  magnitude of the induced seismicity

 
• p

  dimensionless pressure inside the slider

 
• p₀

  dimensionless ambient pressure

 
• P

  pressure inside the slider

 
• P₀

  ambient pressure

 
• Q

  volumetric injection rate per unit area

 
• Q_w

  volumetric injection rate

 
• rq

  dimensionless injection rate

 
• R

  length of the sliding region

 
• S

  area of the sliding region

 
• t_c

  characteristic time

 
• t_slip

  time required to induce a seismicity

 
• T

  vibration period

 
• u

  dimensionless relative displacement between the loading point and the slider

 
• {u}

  nodal displacements vector

 
• U

  relative displacement between the loading point and the slider

 
• v

  dimensionless velocity of the slider

 
• v₀

  dimensionless loading velocity

 
• v_max

  maximum velocity of the slider

 
• V

  velocity of the slider

 
• V₀

  loading velocity

 
• V_*

  velocity of the slider at steady state

 
• α

  dimensionless frictional parameter

 
• $α^$

  experimentally derived parameter relating the normal stress change to the frictional force change

 
• $ε$

  dimensionless vibration period

 
• η

  fluid dynamic viscosity

 
• θ

  dimensionless alternative state variable

 
• $Θ$

  alternative state variable

 
• $κ$

  dimensionless shear stiffness

 
• μ_*

  frictional coefficient when the slider is sliding at V∗

 
• v

  Pisson’s ratio

 
• σ

  dimensionless normal stress

 
• $Σ$

  normal stress

 
• τ_c

  characteristic shear stress

 
• $ϕ$

  porosity

This study was performed as part of the activities of the Research Institute of the Sustainable Future Society at the Waseda Research Institute for Science and Engineering, Waseda University. The authors thank INPEX Corporation for reviewing and approving the publication of this paper.

This study uses the embedded SIRK method to compute numerical solutions of Eqs. 5 through 8. The embedded SIRK method is a combination of the semi-implicit and embedded Runge-Kutta methods, which are briefly explained in this section. For ease of comprehension, the ordinary differential equation (Eq. A-1) is used to explain the solution technique rather than the governing equations in Eqs. 5 through 8.

To obtain a numerical solution, it is necessary to determine the estimated value $y1$ at $x=x1$ based on the known values $x0,y0$⁠.

The general form of the Runge-Kutta method is expressed as follows (Hairer et al. 1987):

where $h$ denotes the step size (⁠$h=x1-x0$⁠). The $s$th order square matrix $A$ and the $s$th order vectors $b$ and $c$ are introduced as follows:

The values in the matrix and the vectors are determined such that the numerical solution $y1$ obtained from the Runge-Kutta method matches the exact solution derived from Taylor series expansion. These arrays are organized using the Butcher array as shown in Fig. A-1a , which is a common representation of the Runge-Kutta method. Fig. A-1b  shows the Butcher array used in this study.

If the $s$th order square matrix $A$ is a strictly lower triangular matrix, where $aij=0j≥i$⁠, the calculation of $y1$ is performed as follows:

In this case, $k1,k2,…,ks$ are obtained successively. This method is known as the explicit Runge-Kutta method. If the $s$th order square matrix $A$ is a lower triangular matrix, where $aij=0j>i$⁠, the calculation of $y1$ is performed as follows:

In this case, $k1,k2,…,ks$ are obtained successively. However, in the process of obtaining $ki(i=1,…,s)$⁠, it is necessary to perform iterative calculations such as the Newton method. This method is called the SIRK method.

The dimensionless quantities $u,v,θ$⁠, and $p$ undergo rapid changes during rapid sliding. Differential equations characterized by intervals of such rapid changes are referred to as stiff equations (Hairer and Wanner 1991). Because of the stiffness in Eqs. 5 through 8, the SIRK method is used in obtaining numerical solutions for $u,v,θ$⁠, and $p$ rather than the explicit Runge-Kutta method.

The above explanation describes a numerical solution technique for a single ordinary differential equation. However, the governing equations in Eqs. 5 through 8 are four simultaneous ordinary differential equations. To address this issue, we consider the following simultaneous ordinary differential equations:

where $n$ denotes the number of ordinary differential equations. The general form of the Runge-Kutta method for solving simultaneous ordinary differential equations is as follows:

The Runge-Kutta method outlined above has one serious problem: It cannot automatically adjust step size $h$⁠. One possible approach to address this problem is the embedded Runge-Kutta method. In this method, the numerical solution $y^1$⁠, which is one order higher than $y1$⁠, is computed simultaneously, as follows (Celaya and Aguirrezabala 2015):

The step size is updated based on the error between $y1$ and $y^1$⁠. The larger the error, the smaller the next step size. For a more detailed explanation, please refer to Celaya and Aguirrezabala (2015). The $s$th order vector $b^$ is embedded in a Butcher array, as shown in Fig. A-2a . Fig. A-2b  shows the Butcher array used in this study.

In the case study, it is necessary to calculate or assign the nondimensional physical quantities (⁠$κ,c,rq,σ,p0,a,b,α,ε,$ and $v0$⁠) to solve the governing equations using the embedded SIRK method. This section provides a detailed explanation of the calculation and assignment procedures for these nondimensionalized physical quantities.

According to Table 3 , the dimensionless shear stiffness $κ$ is obtained from the following calculation:

Because $ks$ (Pa/m) is the shear stiffness, which corresponds to the stiffness of a fault, $ks$ can be calculated as follows (Scholz 2002):

where $E$ (Pa) is the Young’s modulus, $ν$ (–) is the Poisson’s ratio, and $R$ (m) is the length of the sliding region. By substituting Eq. B-2 and Eq. 10 into Eq. B-1, the dimensionless shear stiffness $κ$ can be expressed as

As shown in Table 3 , the dimensionless diffusivity $c$ is calculated using the following expression:

The effective normal stiffness $kneff$ (Pa/m) can be calculated as follows (Alghannam and Juanes 2020):

where $ϕ$ is the porosity, $1/Kv$ (1/Pa) is the uniaxial compressibility, and $1/Kf$ (1/Pa) is the fluid compressibility (Wang 2000). Uniaxial compressibility $1/Kv$ can be calculated using the following equation (Fjær et al. 2021):

where $K$ (Pa) is the bulk modulus and $G$ (Pa) is the shear modulus. The bulk modulus $K$ and the shear modulus $G$ are expressed as (Fjær et al. 2021):

By using Eqs. B-6 through B-8, the effective normal stiffness $kneff$ can be expressed as follows:

By substituting Eqs. B-9 and 9 into Eq. B-4, the dimensionless diffusivity $c$ can be expressed as

According to Table 3 , the dimensionless injection rate $rq$ is calculated as follows:

Because $Q$ (m/s) is the volumetric injection rate per unit area, $Q$ can be expressed as follows:

where $Qw$ (m³/s) is the volumetric injection rate and $Area$ (m²) is the cross-sectional area of the fluid flow. In this case study, the assigned value of $Area$ was based on the size of the horizontal fracture, which can be assumed by the size of the microseismicity cloud. For simplicity, we assumed a homogeneous formation in which the injected fluid permeated both upward and downward from the horizontal fracture. This assumption allows for a simplified representation of the fluid migration in a targeted formation. By substituting Eqs. B-12, B-9, and 9 into Eq. B-11, the dimensionless injection rate $rq$ can be expressed as

Furthermore, according to Table 3 , $σ,p0,a,b$⁠, and $α$ are calculated as follows:

In addition, the values of $ε$ and $v0$ are assigned as follows:

Table B-1  provides a summary of the assigned physical quantities necessary for calculating $κ,c,rq,σ,p0,a,b$⁠, and $α$⁠.

This paper (SPE 214891) was accepted for presentation at the SPE Annual Technical Conference and Exhibition, San Antonio, Texas, USA, 16–18 October 2023, and revised for publication. Original manuscript received for review 28 November 2023. Revised manuscript received for review 27 January 2024. Paper peer approved 7 March 2024.