Normal fault reactivation induced by hydraulic fracturing: Poroelastic effects

Seismic fault slip has been attributed to engineered subsurface fluid disturbances as early as 1929 when an event with a magnitude of 5.3 $ML$ occurred during solution mining in Attica-Dale, New York (Street and Turcotte, 1977). Since then, induced earthquakes have also been reported in possible association with mining operations, oil and gas production, enhanced oil recovery fluid injection, saltwater disposal, reservoir impoundment, groundwater extraction, and hydraulic fracturing (Davies et al., 2013; Foulger et al., 2018; Haddad and Eichhubl, 2020a). The recent hydraulic-fracturing cases that have reportedly caused surface-felt earthquakes and their interpretations have been reviewed by Schultz et al. (2020). Seismic survey arrays installed to monitor hydraulic-fracture (HF) propagation have clearly revealed earthquake clusters, including events of a wide magnitude range that can be attributed to fault reactivation, with examples reported from Jonah Field, Wyoming (Wolhart et al., 2005; Maxwell et al., 2008); Lancashire, UK (Green et al., 2012); Montney Formation, Canada (Maxwell et al., 2011; Ghofrani and Atkinson, 2016); Duvernay Formation, Canada (Bao and Eaton, 2016; Eyre et al., 2019); Barnett Shale, Texas (Kratz et al., 2012); Appalachian Basin, Ohio (Brudzinski and Kozlowska, 2019); and southern Delaware Basin, Texas (Savvaidis et al., 2020). Seismicity has been linked to only a small fraction of the hundreds of thousands of hydraulic stimulations that have been conducted worldwide. For instance, earthquakes with magnitudes larger than 3.0 in Western Canada from 2009 to 2019 are spatiotemporally correlated with only 0.8% of the stimulation wells (Ghofrani and Atkinson, 2020).

Because of the relatively recent widespread application of hydraulic fracturing in unconventional plays, there are limited studies on the causative mechanisms of the seismic fault reactivation associated with hydraulic fracturing. These mechanisms can be influenced by a wide range of geologic parameters, including the mechanical properties of the target layer relative to the underlying or overlying layers, as well as operational ones such as the hydraulic stimulation design, the number of fracture clusters, and the spacing between the fault and HF (Schultz et al., 2017, 2020; Fasola et al., 2019). It is expected that this spacing markedly contributes to the likelihood of induced seismicity, as demonstrated by Hui et al. (2021a) in a hydraulic-fracturing case study in Fox Creek, Alberta, Canada. This spacing is considered a controlling parameter for regulatory guidelines put in place in some U.S. states (e.g., Ohio) for managing induced seismicity risk (Westwood et al., 2017).

Seismicity induced by hydraulic fracturing is differentiated from that by other subsurface fluid injection processes through its proximity to well-stimulation activity and by relatively small fluid volumes compared with wastewater injection, thus affecting smaller rock volumes. In contrast to the typically larger delays of fault reactivation between wastewater injection and seismic events in the order of months or years (Haddad and Eichhubl, 2020b, 2023; Zhai et al., 2020), hydraulic-fracturing-induced seismicity usually occurs within hours or days of hydraulic stimulation (de Pater and Baisch, 2011; Holland, 2011; BC Oil and Gas Commission, 2012; Lacazette and Geiser, 2013; Lomax and Savvaidis, 2019). The typically closer spatiotemporal link between hydraulic-fracture stimulation and seismic activity facilitates the establishment of a firm temporal correlation between seismicity and hydraulic fracturing in contrast to wastewater injection, for which such correlation is often more difficult to establish because of local variations in pressure diffusion rates. Compilations of microseismic data from thousands of stimulations reveal that seismic events generally cluster within 600 m of the well perforation in the depth direction (Davies et al., 2012; Flewelling et al., 2013), in contrast to disposal-induced earthquakes that may occur tens of kilometers away from the disposal wells.

Various numerical techniques have been used to study hydraulic-fracturing-induced seismicity. The Haskell propagator method offers a mesh-free solution for the coupled elastic subsurface deformation and Darcy flow through porous media based on surface displacement (Wang and Kumpel, 2003). The elastic deformation and Darcy flow equations can also be solved iteratively using a numerical finite-element mesh in the COMSOL Multiphysics® Software Program (2020), as conducted by Hui et al. (2021b, 2021c) for an induced seismicity problem. Another approach is based on the distinct-element method (Damjanac and Cundall, 2016), used by Maxwell et al. (2015) in reservoirs with a discrete fracture network. Rutqvist (2011) applies an iterative coupling scheme between the fluid flow and geomechanical simulators, TOUGH and FLAC, and a continuum hydraulic-fracturing-induced fault reactivation criterion based on models to increase porosity and permeability as functions of the plastic tensile and shear strains. Rutqvist et al. (2015) extend their 2D numerical model (Rutqvist et al., 2013) to a 3D case to resolve the intrinsic uncertainties of simplifying an actual 3D geometry to a 2D model. However, their model considers only hydraulic-fracturing injection into a subvertical fault plane. Andres et al. (2019) also use 3D simulations in TOUGH coupled with FLAC linking hydraulic stimulation at the Basel geothermal site with induced seismicity. However, the models based on TOUGH and FLAC cannot simulate large-magnitude seismic events due to hydraulic fracturing (Rutqvist et al., 2013, 2015); in these models, events up to $Mw$ 2–3 do not occur unless the fault is extremely brittle, with residual friction coefficients of less than 0.2. This may be attributed to the inability of these models to rigorously simulate the convective fluid transport through the fracture and fault conduits. In other words, these model choices result in substantial pore-pressure diffusion in large elements within the fault, which unduly reduces the excess pore pressure within the fault gaps and, hence, the fault reactivation potential. In this study, we do not simulate seismic events; however, to overcome the limitations of TOUGH and FLAC in simulating fluid convection, we developed a new numerical approach based on the finite-element method and the use of a cohesive zone model, which provides a fully coupled solution to the fluid flow and poroelastic stresses in the porous matrix as well as the fluid flow and mechanical deformation of the HFs. We then postprocess this solution to assess a fault reactivation criterion (e.g., the Coulomb failure stress [CFS]) over a fault core.

In this work, we pursue possible mechanisms for the reactivation of a fault due to nearby single- or multiple-stage hydraulic fracturing using numerical conceptual simulations. Our models leverage the simulation and in-situ stress estimation tools already developed for the Barnett Shale in the Fort Worth Basin, northeastern Texas (Fan et al., 2019; Haddad and Eichhubl, 2020b, 2023) for HF propagation. Through these numerical simulations, we intend to identify the most influential parameters for the increase in the fault reactivation potential due to hydraulic fracturing in the vicinity of an almost critically stressed normal fault. We focus on the effect of the following hydrogeologic, geomechanical, and operational parameters on hydraulic-fracturing-induced seismicity: fault damage zone conductivity, toughness- versus stress-controlling height growth, the number of HFs, and the fracturing sequence. We assume an optimally oriented fault case where $SH,max$ azimuth is parallel with the fault azimuth or $Sh,min$ azimuth is normal to the fault azimuth. Subsequently, we assume that the horizontal well, in our case, is placed normal to the fault azimuth, which is consistent with the drilling of horizontal wells in unconventional reservoirs in the $Sh,min$ azimuth to initiate HFs perpendicular to these horizontal wells.

In contrast to Rutqvist et al. (2013, 2015), who simulate fracturing-induced fault reactivation due to the direct injection of the fracturing fluid into the fault zone, we test if fault reactivation can occur by poroelastic stress interaction ahead of the HF tip, with and without the intersection of the HF with the fault, and without direct injection of the fracturing fluid into the fault zone. To avoid direct injection into the fault zone, we constructed the models by assuming a relatively large offset of 50–200 m between the fracturing perforation and the fault zone. In this parametric study, we include the effect of a fault damage zone of varying width to account for a geologically more realistic fault permeability structure that properly accounts for higher fault-zone permeability on the pressurized fluid distribution (Flodin et al., 2001; Eichhubl et al., 2009). The fully coupled poroelastic finite-element analysis enhanced by a fracture-fault intersection model allows us to evaluate whether seismic fault reactivation requires pore-pressure changes within the fault zone due to fracturing fluid invasion or if a poroelastic stress disturbance with only a small direct pore-pressure effect is sufficient to lead to fault reactivation, without fluid invasion.

For all of the model scenarios tested here, we consider a conductive normal fault striking parallel to the maximum compressive horizontal stress (⁠$SH,max$⁠) and dipping 60°. Here, $Sh,min$⁠, $SH,max$⁠, and $Svert$ gradients through depth are 15.37, 19.38, and 26.2 kPa/m, equivalent to 0.68, 0.86, and 1.16 psi/ft, respectively. We assume that these stress gradients are uniform across the model domain, meaning that the presence of the fault does not perturb in-situ stresses. The fault is initially stable but close to failure, with a friction coefficient of 0.6. The initial CFS expressed as the shear stress subtracted by frictional strength with the formulation in the “Method” section over the fault core varies from −0.0034 to −0.0037 MPa from the top to the bottom of the fault. We consider three stacked layers consisting of a middle target layer (where the horizontal well is drilled and perforated), an underlying layer, and an overlying layer, which are all crossed by the normal fault, as shown in Figure 1. These layers are assumed to have hydromechanical properties, as provided in Table 1. Fluid viscosity and bulk modulus values for the different layers were selected, assuming that the target layer was to be saturated with natural gas and the overlying and underlying layers were to be saturated with saltwater (Table 1). Although this work is a conceptual study, for the bulk modulus of natural gas in the target layer and water in the other layers, we considered the values suggested by Haddad and Eichhubl (2023) for the Barnett Shale and its adjacent layers within the Fort Worth Basin. In addition, for simplicity, we considered 1 cp for the viscosity of water and 0.1 cp for the viscosity of the mixture of natural gas and resident water in the target layer.

The target layer is a 60 m thick shale layer, and the underlying and overlying sedimentary layers are each 170 m thick. The overlying and underlying layers perform as mechanical barriers for HF height growth in some simulation scenarios through the higher than the target layer opening-mode fracture toughness of these layers. The fault damage zone permeability in all these conceptual models is assumed to be twice the hosting layer permeability in all the scenarios; however, we evaluate the effect of the damage zone conductivity as the product of the damage zone permeability and thickness by varying the fault damage zone thickness between 1 and 20 m.

We use the solution outcome within each time increment to calculate the CFS over the fault core to assess the reactivation potential of the fault through time. CFS is defined as $τ−μs(σn−Pp)$⁠, where $τ$ and $σn$ denote the shear and total normal stresses over the fault core, respectively; $Pp$ is the pore pressure; and $μs$ represents the static friction coefficient, which is assumed equal to 0.6 in this work. The term $(σn−Pp)$ is Terzaghi’s effective normal stress, which is one of the governing factors of fault stability and brittle failure (Brace and Martin, 1968; Paterson and Wong, 2005). Stress components $τ$ and $σn$ are calculated through the inner product of the Cauchy stress tensor and a unit vector normal to the fault plane. CFS values smaller than zero suggest that the fault is fully stable, whereas CFS values equal to or larger than zero show the onset of fault reactivation.

Fracture growth is simulated with the use of a pore-pressure cohesive zone model, as extensively adopted for this purpose in numerous works (Haddad and Sepehrnoori, 2015, 2016; Haddad et al., 2017; Sobhani Aragh et al., 2018). The plane over which an HF grows is parallel to the $σ1−σ2$ plane and normal to $σ3$ azimuth, where $σ1$⁠, $σ2$⁠, and $σ3$ are the maximum, intermediate, and minimum compressive principal stresses, respectively (Economides and Nolte, 2000). This assumption is supported by the observation of the linear HFs in a consistent overall orientation normal to $σ3$ azimuth in the cores through the stimulated reservoir volumes, for instance, in the Eagle Ford Shale (Raterman et al., 2018) and the Upper and Middle Wolfcamp (Gale et al., 2018). In our models, $σ1$ is assumed completely along the depth direction and, therefore, $σ2$ and $σ3$ in the horizontal directions. The plane on which the HF is expected to grow is modeled by a layer of the pore-pressure cohesive elements normal to the $σ3$ or $Sh,min$ azimuth. In addition, the plane representing the fault core is modeled by a layer of the pore-pressure cohesive elements, making a nonorthogonal angle with respect to the principal stresses. These enable the simulation of the hydraulic-fracture opening and fault-core opening as well as the shear- and tearing-mode failures. In general, the intersection of a weak fault plane by an HF may divert the fracture growth direction toward the fault plane (Haddad et al., 2019; Zhang et al., 2020). The cohesive elements within the fracture and fault-core domains follow a macroscopic cohesive law representing microscopic progressive damage within the fracture tip process zone. This law is described as a triangular traction-separation response composed of two parts: a linear elastic part prior to the attainment of the fracture initiation stress, where the stress across an interface monotonically increases, and a progressive damage part posterior to the fracture initiation, where the damage nucleation and coalescence lead to a gradual reduction in the interface load carrying capacity (Haddad and Sepehrnoori, 2015).

To improve the numerical convergence in the presence of the cohesive layer, we increased the load-carrying capacity of the cohesive layer by increasing the stiffness of the cohesive elements $Km$ in all the modes of failure by 120 times relative to the host rock for the cohesive layers that model hydraulic fracturing within the host rock and 60 times relative to the host rock for the cohesive layers that model the fault core. In other words, we weakened the fault-core fracture initiation stress and stiffness by half when compared with the fracture initiation stress and stiffness for the HF in the host rock. Notably, the effect of this elastic response on the numerical solution is likely negligible, as it is effective only prior to the rock failure in a thin rock volume that is modeled by cohesive elements or within the small fracture process zone at the fracture tip.

Fluid flow within the fracture or the fault core is coupled with the flow with the porous media through the continuum-based leakoff model that is governed by the fluid pressure gradient normal to the fracture wall multiplied by a constant equal to 5.879 × 10⁻¹⁰ m³/s/kPa in this work. We selected this small leakoff coefficient to simulate negligible fracturing fluid leakoff into the formation in an unconventional reservoir. Hence, we can assume that the hydraulic-fracture propagation negligibly wets the reservoir through leakoff. To honor the fluid flow continuity at the intersection of the fracture and the fault core, we developed a search plugin to search for the closest middle-edge nodes of the cohesive elements at the intersection and to tie the pore-pressure degrees of freedom over these middle-edge nodes together. Haddad et al. (2017) further elaborate on the details of this pore-pressure coupling at the intersection.

Hydraulic fracturing is commonly conducted through multiple stages with multiple perforation clusters within each stage. However, for the simplicity of the conceptual simulations, we assume one perforation cluster within each stage. In addition, we assume that one perforation cluster yields one major HF, which grows over a predefined plane of cohesive elements. To model the propagation of multiple HFs in various stages, an arbitrary number of these planes are defined in the finite-element framework by partitioning the model into several fracture domains, identified as planes with negligible initial thickness and the remaining poroelastic domain.

We assume that the pore space is occupied by a single-phase fluid: gas in the target layer and brine in the overlying and underlying layers. The fracturing fluid is water and comes in contact with the gas-saturated rock in the target layer. Due to the low permeability of shale formations and relatively short stimulation durations of 1 h per stage (hence, insufficient time for the fracturing fluid to substantially permeate into the matrix), we can assume that increasing the number of pore fluid phases (gas versus water and gas) very close to the fracture wall does not affect the HF growth and hence does not change the potential for hydraulic-fracturing-induced fault reactivation. Therefore, our choice to assume a single-phase pore fluid in our model does not substantially influence the outcomes for induced fault reactivation. In addition, for simplicity, we assume that the fracturing fluid is proppant-free and, therefore, the injection-induced dilation of the HF can be followed by fracture closure during shut-in.

The selected anisotropic in-situ stress state assumes an overburden stress gradient of 26.2 kPa/m (1.16 psi/ft) and a pore-pressure gradient of 10.26 kPa/m (0.45 psi/ft) (Table 2). The estimated minimum horizontal stress of 15.37 kPa/m (0.68 psi/ft) assumes fault stress criticality and a widely accepted friction coefficient of 0.6 (Zoback et al., 2002). The maximum horizontal stress is 19.38 kPa/m (0.86 psi/ft), assuming the generalized Angelier’s shape parameter $Aφ$ of 0.37, which we adopted from the Venus earthquake cluster in a normal faulting stress regime (Quinones et al., 2018; Haddad and Eichhubl, 2023). An extensive discussion about $Aφ$ is provided by Simpson (1997) based on the original introduction of the Angelier (1979) shape parameter $φ$⁠. Although the considered $Aφ$ is slightly lower than the suggested value of 0.5 for normal faulting in the Barnett Shale (Lund Snee and Zoback, 2022), the effect of $Aφ$ on CFS seems negligible because the assumed fault strikes parallel to the maximum horizontal stress.

A cubic model geometry is assumed, consisting of vertically stacked layers that are each homogeneous and bounding the shale layer in between (Figure 2). All of these layers are partitioned vertically for the assignment of exclusive HF and obliquely for the fault domains (highlighted by the blue and yellow lines in Figure 2). We assume that the minimum principal stress is horizontal and that the horizontal well is parallel to this principal stress and located within the XZ symmetry plane of the model. This would lead to vertical HFs perpendicular to the horizontal well and the modeled normal fault being crossed by this horizontal well.

The target layer is assumed to be saturated with gas with 0.1 cp viscosity, and the underlying and overlying layers are saturated with saltwater with 1 cp viscosity. We assume orthotropic matrix permeability with $kzz/kxx(yy)$ of 0.1, where $kii$ denotes the permeability in the ith direction, z represents the depth direction, and x and y are the horizontal directions. We also assume that permeability follows an exponential function of porosity $k/k0=(φ/φ0)8$⁠, where k and $k0$ denote the current and initial permeability, respectively, and $φ$ and $φ0$ are the current and initial porosity, respectively. The exponent 8 is chosen after Bourbie and Zinszner (1985) based on laboratory experiments on clastic rocks. The dry-rock bulk modulus $Kdry$ used in the material property assignment is calculated as $E/3(1−2ν)$⁠, where E and $ν$ denote Young’s modulus and Poisson’s ratio, respectively. In addition, the grain bulk modulus is obtained as $Kdry/(1−α)$⁠, where $α$ denotes the Biot-Willis coefficient (Table 1). Water density $ρw$ is 1046 kg/m³, and the solid density is expressed as $∂zSvert/g−φρw$⁠, where $∂zSvert$ is the assumed vertical stress gradient of 26.2 kPa/m (1.16 psi/ft) (Tables 1 and 2). We inject water for 3650 s (approximately 1 h) with the volumetric rate of 0.42 m³/s (160 U.S. barrels per minute) into one or multiple fracture clusters, meaning that the entire stimulation time is almost 4 h in scenario 6 (sequential fracturing) and 1 h in the other scenarios. Each fracture cluster comprises only one growing biwing fracture. The injected fluid in a cluster flows into a 3D fracture plane in the fracture wings and the height direction. The injection rate ramps linearly up to the maximum rate in the first 200 s of injection. This volumetric injection profile leads to 1505 m³ of total injected volumes during almost 1 h of injection time. Although hydraulic stimulations are conducted using multiple fracture clusters within a stage, we assume only one growing fracture within a stage. Simulating the allocation of the injected water into multiple fractures or clusters per stage is outside the scope of this study and its focus on fault reactivation. This assumption is justified by models predicting that dominant fractures receive most of the fluid injected into a fracture stage at small cluster spacings (Yang et al., 2022).

All the boundaries are under a hydrostatic pore-pressure boundary condition (i.e., open-flow boundary condition), the three lateral boundaries are under horizontal stresses increasing with depth, the upper boundary is under the overburden stress at a depth of 3600 m, and the lower and right-lateral boundaries are under a zero normal displacement boundary condition (Figure 2). To eliminate the free body motion in the y-direction, three nodes close to the lower-right corner of the XZ symmetry plane are fixed against y-displacement.

We consider nine scenarios that differ in geometric and operational aspects (Figure 1; Table 3):

1. Fault damage zone thickness: 1 or 20 m.

2. Intersection location of the HF and fault: within the layer overlying the shale reservoir or the underlying layer.

3. Upward or downward growth of the HF controlled by adjusting the opening-mode fracture initiation stress $tI$ and toughness $KIc$ of the overlying and underlying layers (Table 4).

4. Number of fractures: 1 or 4.

5. The sequence of fracturing clusters: simultaneous or sequential.

6. Spacing between the fault and HF, $dPerf-F$⁠, measured at the middle depth of the target layer: 50, 100, 150, and 200 m.

Notably, the parameter space that is generated by variations in the preceding model aspects does not cover many other factors that may influence fault reactivation, such as the hydromechanical properties of the host rock, except for $tI$ and $KIc$ of the overlying or underlying layers or the target layer thickness. However, the defined simulation models account for a range of parameters expected for different geologic settings and stimulation designs. Furthermore, these parameters are chosen for this study because they seemingly contribute significantly to fracturing-induced fault reactivation due to two mechanisms in our point of view: (1) poroelastic stress disturbance as a result of hydraulic-fracture tensile opening and (2) pore-pressure increase within the fault core as a result of fracturing fluid flow along an HF intersecting the fault zone. In general, pore-pressure diffusion through the host rock is considered the primary mechanism for injection-induced seismicity (e.g., in saltwater disposal cases), and host rock permeability seems the dominant controlling factor for this pore fluid diffusion (Chang and Yoon, 2022; Ge et al., 2022; Haddad and Eichhubl, 2023; Stokes et al., 2023). However, in unconventional reservoirs, pore-pressure diffusion solely is too slow to explain fracturing-induced seismicity in these reservoirs (Kettlety and Verdon, 2021); hence, we did not include host rock permeability in the parameter space. In addition, we did not consider the change in Young’s modulus in our parameter space because of its absence in the stress shadow analytical solutions developed by Sneddon and Elliott (1946) and Pollard and Segall (1987). Nevertheless, our conclusions are derived based on the difference of a reactivation criterion on a fault core across the considered simulation scenarios.

The effects of hydraulic fracturing on the stability of adjacent faults as a function of different geometric and operational parameters are quantified through changes in CFS over the fault core (Figures 3, 4, and 5). In Figures 3–5, the green contours designate the stable regions of the fault core with $CFS<0$⁠, and the red contours designate the unstable regions with $CFS>0.1MPa$⁠. Although there is a complex transition between these two states, for simplicity we focus on these two extreme states for comparison between scenarios. We selected $CFS>0.1MPa$ as the criterion for fault reactivation after Kilb et al. (2002). Notably, although this criterion determines the onset of fault reactivation, our model does not simulate the actual seismic events and, therefore, cannot predict the occurrence of specific events even if this criterion is satisfied. The unstable region of the fault core in Figures 3–5 is centralized in the y-direction because the injection is conducted into a central perforation within the target layer over the intersection line of an HF plane and XZ symmetry plane of the model (Figure 2). In addition, the larger model size in the y-direction compared to the fracture propagation in this direction, as shown in Figure 4 for three simulation scenarios, leads to the restriction of the unstable region to a central area of the fault core.

To show the effect of the fault damage zone thickness and the hydraulic-fracture-fault intersection location on the fault stability, we plot CFS across the fault core in normal view to the fault core in scenarios 1–4 after 1 h of injection (Figure 3). Increasing the fault damage zone thickness from 1 to 20 m (scenarios 1 and 3 versus scenarios 2 and 4) substantially extends the unstable area of the fault core. In addition, changing the hydraulic-fracture-fault intersection location from the overlying layer to the underlying layer reduces the unstable area of the fault core and moves this unstable zone from above to below the target layer (scenarios 1–3 and scenarios 2–4).

Varying the number of HFs in scenarios 1 and 5 (Figure 4) demonstrates that, given our specific model parameters and geometry, simultaneous fracturing of more than one cluster may slightly decrease the unstable area of the fault core (Figure 4b and 4d). In addition, changing the fracturing sequence of four clusters from simultaneous to sequential (scenarios 5 and 6) enlarges the unstable area of the fault core due to the staggered height growth of the HFs at various clusters (Figure 4d–4f). The slightly unsymmetrical growth of the HF in scenario 1, as shown in Figure 4a, is likely due to asymmetrical mesh with respect to the perforation. Considering that the boundary condition and material properties are symmetrical with respect to the perforation in the y-direction, the effect of the mesh asymmetricity on the fracture growth in the y-direction could be minimized by a larger viscous regularization parameter. Nevertheless, this asymmetrical fracture growth does not influence the overall comparison of CFS over the fault plane in this scenario with that in the other scenarios.

Increasing the spacing between the single-stage HF and the fault leads to a significant decrease in the unstable area of the fault core, as shown by comparing CFS in scenarios 1, 7, 8, and 9 (Figure 5). Moreover, increasing this spacing leads to the expansion of the fault-core stable area, shrinkage of the unstable area, and shift toward a shallower depth. The expansion of the stable area with the increase in the spacing between the HF and the fault is due to the increase in the rock volume between the HF and the fault under the compressional stresses induced by the HF tensile opening.

Results presented so far reflect the changes in poroelastic stresses over a smooth or planar fault with a well-quantified dip angle. However, faults are corrugated, and the deviation of a fault from a smooth plane can lead to irregular, unstable areas, considering that seismic event hypocenters are usually concentrated at fault seismogenic asperities as shown for the Tohoku-Oki and Loma-Prieta earthquake clusters (Wald et al., 1991; Xie et al., 2019). To test the robustness of our fault stability predictions as a function of fault corrugation, without running additional simulations, we accounted for the effect of fault corrugation on CFS explicitly. We assumed that the principal stresses at the fault location change only due to the stress shadow effect of a nearby dilating HF. The stress shadow effect is defined as poroelastic stress changes in a relatively large rock volume around an HF, induced by the opening-, shearing-, or tearing-mode displacement of this HF under large fracturing pressures. In addition, we assumed that the presence of the fault does not influence these stresses and that normal and shear stresses at a specific depth are only functions of the fault dip at that depth. We generated five fault dip realizations using 10 control points through depth along the dip section of the fault (Figure 6a) by calculating the fault dip at each control point using Monte Carlo sampling from a Gaussian cumulative distribution function with a mean of 1.047 radians (60°) and standard deviation of 0.1047 (6°), respectively. From these fault dip realizations, we obtained fault cross sections over the XZ symmetry plane, as shown in Figure 6b. Comparison of the CFS over the vertical symmetry line of the fault at the different dip realizations and various scenarios in Figure 7 shows that the CFS values are between −1 and +1 MPa for most of the depth interval and locally spike up to large values (e.g., 11 MPa in scenario 3, Figure 7a). Considering that these spikes could initiate fault reactivation by themselves, we can rank the scenarios that lead to fault reactivation from earliest to latest as follows: scenarios 3, 7, 2, 4, 5, 1, 6, 8, and 9. In this ranking, we assume a monotonous stress disturbance of each fault point through time due to the HF growth, which would lead to the earliest fault reactivation (or $CFS>0$⁠) in the scenario with the maximum CFS peak at an arbitrary injection time. As a reminder, scenario 3 was defined with $dPerf-F=50m$⁠, $tFDZ=1m$⁠, and the intersection of the fault and a single HF plane within the underlying layer, and scenario 9 was defined with $dPerf-F=200m$⁠, $tFDZ=1m$⁠, and no intersection of the fault and a single HF plane. Figure 7b–7j confirms this ranking regardless of the fault dip realization. Figure 7b–7j also shows that changing the fault dip realization leads to the observable variation in CFS between −1 and +1 MPa, and the CFS value and depth of the spike remain almost unchanged.

By changing only one model parameter in each scenario, we isolated the effect of each parameter on the fault reactivation induced by hydraulic fracturing. Following this approach, we find that thicker fault damage zones lead to larger unstable fault areas (Figure 3). Due to the large permeability contrast between the fault damage zone and the host rock, the pore fluid that is trapped within the fault damage zone and pressurized by the hydraulic-fracturing-induced compression can permeate more effectively within a damage zone with higher conductivity. For a fixed hydraulic-fracturing-induced compression, this mechanism would lead to more spread pore-pressure elevation within a more conductive damage zone in contrast to a more concentrated pore pressure increase within a restricted area of a less conductive damage zone. We also find that hydraulic-fracture height growth is closely associated with fault reactivation in the overlying layers (Figures 3 and 4a). In addition, for multiple-stage fracturing near a normal fault in a normal faulting stress regime, as we study here, we show that sequential fracturing leads to a larger area of the fault being under enhanced shear stresses (or higher CFS) because of nonuniform staggered fracture growth compared with the almost uniform multifracture growth associated with simultaneous fracturing. This variation in multifracture growth can be explained through the different stress shadow effects of the fractures in these stimulation scenarios. The nonuniform staggered fracture growth in sequential fracturing is due to the preexisting stress shadow effect of previous fully stimulated fracture(s) on a currently stimulated fracture, in contrast to an almost uniform multifracture growth in simultaneous fracturing due to gradually increasing the reciprocal stress shadow effect of all growing fractures on each other. Simultaneous fracturing may thus be preferred over sequential fracturing because it leads to more uniform fracture dimensions and less height growth driven by the stress shadow effect (Figure 4). In other words, sequential fracturing leads to a larger area of the fault under enhanced shear stresses due to nonuniform staggered fracture growth compared with the enhanced shear-stress area over this fault due to almost uniform multifracture growth by simultaneous fracturing. Furthermore, for single-stage hydraulic fracturing near a critically stressed fault, a spacing between the fault and HF of approximately 200 m is required to avoid fault reactivation during stimulation (Figure 5). We ranked the simulated scenarios in their potential to lead to fault reactivation based on the size of the fault area that exceeds the CFS of 0.1 MPa (Table 3). However, in all these scenarios, we simultaneously observed a large stabilizing fault area with a declining CFS on the other parts of the fault (the green area in Figures 3–5). The stabilizing forces over this fault area may resist the occurrence of a large earthquake and rather restrict fault reactivation to the surrounding areas. These large variations in the fault reactivation potential are in contrast to the relatively uniform increase of the reactivation potential of a fault associated with saltwater disposal (Haddad and Eichhubl, 2020b, 2023). However, these large variations are consistent with the inhibited near-field seismicity and promoted remote seismicity by poroelastic effects in the seismological modeling framework for fractured porous media proposed by Jin (2022).

To account for local variations in fault orientation, we evaluated the effect of fault corrugation on the change in the reactivated fault area. For the fault area that undergoes slightly positive CFS values in the case of a smooth fault with a favorable dip, a small deviation of the fault dip from the favorable dip can lead to the marked drop of CFS and, hence, the stabilization of the fault in the associated area (Figure 7). Because of the large variation in CFS between −1 and +1 MPa at different realizations of the fault corrugation, the stimulation-related fault reactivation potential should be assessed based on stress disturbances resulting in CFS exceeding 1 MPa. This assessment criterion is also supported by the relatively fast changes in the induced poroelastic stresses during hydraulic fracturing and the uncertainty in fault corrugation characterization.

The variation in rock properties, the geometric configuration of a fault with respect to HFs, and operational parameters may lead to fault reactivation in relation to hydraulic fracturing via three different mechanisms: (1) pore-pressure diffusion across the fault zone from a nearby HF or from an injection point within the fault zone at adequately low injection pressures not to induce the fault tensile opening; (2) tensile opening of a fault patch, gradual reduction of the load carrying capacity over the fault core, and fault slip due to the direct injection of fracturing fluid into the fault zone at adequately high injection rates and pressures for the tensile opening of the fault patch, and subsequent shear reactivation; and (3) poroelastic stress disturbance associated with an HF growth near a fault, without the fracturing fluid invasion into the fault zone (i.e., the direct pore-pressure effect).

The first mechanism, also known as the direct pore-pressure effect, involves the transmission of elevated fluid pressure from the stimulation wells to the fault planes through the permeable bedding planes between the laminated shale layers, a natural fracture system, or extensive HFs (de Pater and Baisch, 2011; Westwood et al., 2017; Zbinden et al., 2020; Hui et al., 2021c). This can be inferred by the correlation of the seismicity with the pump rate, the delay of the onset of the seismicity compared with the start or completion of the fracturing, or the fluid pressure wave detection, the latter of which uses a color-coded tomographic fracture imaging (TFI) map (Geiser et al., 2012). In TFI, seismic emission tomography is used to monitor fluid-pressure waves (i.e., the pressure pulse) through the reservoir fairway system. These fluid pressure waves hold soliton-like properties and can travel through the porous media for kilometers as fast as 10 m/s (Geiser et al., 2006). Soliton behavior is explained as small, rapid fluctuations of fluid volumes in conjunction with injection and extraction. This fluid pressure wave velocity is affected by the fairway system transmissibility, which is governed by the rheological characteristics of the fracturing fluid and the hydromechanical properties of the flow pathways (e.g., the bedding plane and hydraulic/natural-fracture conductivities and the bulk modulus of rock). In the Barnett Shale, the maximum number of low-frequency (i.e., slow shear slip) events were recorded during the hydraulic-fracturing stages with the largest pump pressure (correlated with the instantaneous shut-in pressure) and natural fracture density (Das and Zoback, 2011). The onset of fault reactivation induced by this mechanism may be delayed by several hours with respect to the start of the hydraulic-fracturing operation, such as has been documented in the Eola Field, southern Garvin County, Oklahoma (Holland, 2011); the Etsho and Kiwigana Fields in Horn River, Canada (BC Oil and Gas Commission, 2012); the Preese Hall 1 Well in Lancashire, UK (de Pater and Baisch, 2011), southwest of Pecos, Texas; the Delaware Basin (Lomax and Savvaidis, 2019); and the Eagle Ford Shale in Texas (Geiser et al., 2012). This delay can be partially due to (1) the large storage and transmissibility of the fault core and damage zone and (2) the required time for pressure transmission through porous media or a fairway system (Geiser et al., 2012; Lacazette and Geiser, 2013). The same reasons can contribute to the persistence of seismic events after pump shut-in (Bao and Eaton, 2016; Kozlowska et al., 2018). The delay in fault reactivation with respect to the stimulation time is likely controlled by this pressure wave velocity and the spacing between the fault and the active stimulation stage. This spacing could be as large as 400 and 433 m, as demonstrated so far in the Bowland Shale near Blackpool in the UK (Clarke et al., 2014; Westwood et al., 2017). As the CFS threshold increases (i.e., as the fault becomes less critically stressed), the lateral spacing required for the fault reactivation decreases and also becomes exponentially a weaker function of this threshold (Westwood et al., 2017). Evaluating the combined effect of the mentioned parameters on the delay of the fault reactivation with respect to the start of the stimulation would require a joint assessment of the hydraulic-fracturing-induced fault reactivation, using not only TFI, microseismic monitoring, and spatiotemporal data analyses but also fully coupled poromechanical modeling to quantify the hydromechanical properties of the media by history matching of the fault reactivation time with respect to the stimulation time.

The second mechanism involves direct injection at large injection rates and pressures to nucleate and gradually enlarge a fault patch area under the tensile opening, leading to the reduced load-carrying capacity of the fault and the ultimate shear reactivation at failure. This mechanism seems to be effective only in low-permeability fault zones and host rocks to contain injections within a confined volume of the fault damage zone, leading to adequately large pressure buildups for the nucleation and expansion of the fault area under the tensile opening. Unlike the first mechanism, in the case of an HF intercepting the fault, the reduction in the Terzaghi effective normal stress is limited to the fault patch under the tensile opening. This mechanism likely dominated during the experimental fault reactivation in the low-permeability Mont Terri Rock Laboratory, as reproduced using hydromechanical models (Park et al., 2020).

The third mechanism involves the poroelastic stress disturbance (i.e., the stress shadow effect) associated with the HF growth in the close vicinity of the fault but without the physical intersection of the HF with the fault. This poroelastic stress disturbance may increase shear stress over the fault core and, consequently, the CFS. Although this mechanism does not involve a direct pore-pressure effect on fault reactivation, the stress shadow effect due to the HF tensile opening may induce an undrained indirect pore pressure increase within the fault core. An example of this mechanism is presented by Kettlety et al. (2020) for the PNR-1z shale gas well in the UK. As discussed in this work, evaluating this mechanism is critical for the assessment of fault reactivation in the vicinity of hydraulic-fracturing operations, especially where multiple-stage fracturing is planned.

In permeable rock units, the first mechanism has the potential to reactivate a large fault area and, thus, large seismic events. For the second mechanism, the area of fault reactivation is limited by the fluid volume injected into the fault, the magnitude of fault-opening displacement, and the porosity of the fault core. The third mechanism is limited to the vicinity of the HF. As the HF approaches the fault, the zone of poroelastic stress disturbance would increasingly engulf the fault until a sufficiently large area of the fault zone is affected for slip nucleation to occur.

In our simulations, we observed fault reactivation due to the third mechanism: poroelastic stress disturbance associated with the HF tensile opening. The second mechanism was negligible in the simulations because we did not directly inject into a fault in any of our simulation cases, and the underlying fault tensile opening close to the injection point in the second mechanism would become effective in a relatively isotropic in-situ stress state, whereas in all the simulation scenarios, we considered large differential stresses. For an optimally oriented normal fault, tensile opening acts against the fault-normal stress controlled by $Svert$ and $Sh,min$⁠. The second mechanism would be effective where $Svert$ converges to $Sh,min$⁠, leading to the intermediate stress (⁠$SH,max$⁠) converging to $Sh,min$ as well. This mechanism is favored under low differential stresses. Although our model is capable of simulating this fault reactivation mechanism, as demonstrated by Haddad et al. (2019), we did not observe this mechanism in this work due to the large differential stresses chosen for these specific simulations. We emphasize that these simulations do not simulate the process of the slip nucleation or slip instability, only the stress evolution leading up to the slip instability. As such, we can address the size of the fault area affected by the induced stress changes, which may affect the onset of a seismic or aseismic fault slip but not the size of the actual rupture area and, thus, of the slip magnitude.

Our models are capable of simulating HF height growth due to the horizontal stress gradient and the contrast in the mechanical properties between stacked layers (Figure 4a, 4c, and 4e). Therefore, we can evaluate if induced seismicity is promoted on a steeply dipping normal fault either above or below the simulated reservoir layer. Related to the third mechanism for fault reactivation, as discussed previously, this extension of the induced seismicity beyond the stimulation interval is the poroelastic effect caused by the stress disturbance in large rock volumes due to the HF tensile opening. Updip fault reactivation has been documented by Eyre et al. (2019) and numerically simulated by Rutqvist et al. (2015) and Haddad et al. (2017). Relating to the second fault-reactivation mechanism that was discussed previously, the effect of the intersection of an HF and a fault in the overlying or underlying layers on fault reactivation was possible to evaluate in this work. Our results demonstrate that fault reactivation by intersection with an HF is more likely to occur in the overlying layer (scenarios 1 and 2) than in the underlying layer (scenarios 3 and 4). This is because CFS exceeded zero across a larger area of the fault core within the overlying layer in scenarios 1 and 2 than within the underlying layer in scenarios 3 and 4 (Figure 3). The unstable red areas in scenarios 1 and 3 may seem comparable; however, this area in scenario 3 engulfs a relatively large, stable, green area below the HF-fault intersection, leading to a cumulatively more stable fault in scenario 3 than in scenario 1. In addition, an HF intersecting a fault damage zone with high permeability can substantially grow in length within the fault damage zone, leading to fracture-driven interactions (i.e., frac-hit) with the offset wells. This fracture interaction and the associated channeling of fracture fluid along a fault may lead to fault reactivation not only at shallower depths but also at distant locations from the injection sites. This effect can be compounded by a fault receiving hydraulic-fracturing fluid from multiple fracture clusters (e.g., in a zipper frac scheme), thus reactivating a larger fault area and resulting in potentially higher seismic magnitudes.

In our simulations, we observed that CFS remains above zero for at least 1 h after the injection shut-in with contours close to those in Figure 3, which shows the potential for post-shut-in fault reactivation. This agrees with the hour-long recorded duration of seismic activity after hydraulic fracturing in the Preese Hall 1 Well in Lancashire, UK (Clarke et al., 2014). This hour-scale duration of positive CFS over the fault after the stimulation shut-in, however, is very short compared with the persisting fault reactivation after the wastewater injection rates drop at the disposal fields, usually lasting for several months (Haddad and Eichhubl, 2020b, 2023). This large difference may originate from the local effect of injecting small fluid volumes in hydraulic stimulations versus the spatially more extended effect of elevated pore pressure affecting large rock volumes in disposal sites.

In the context of the current work, a progressive poroelastic stress disturbance over the fault plane can be assumed by the fracture height growth toward the adopted normal fault. Nevertheless, this stress effect on fault stability could be projected into a case with a vertically dipping strike-slip fault provided that the HF preferentially grows in length, rather than height, to intersect this vertically dipping strike-slip fault. Considering this rotation of the dominant fracture growth direction requires a revisit of the fracture-growth limiting mechanisms, as discussed extensively by Fisher and Warpinski (2012). They listed the following favorable conditions for hydraulic-fracture length growth rather than height: (1) a higher fracture toughness or minimum horizontal stress of the bounding layers compared with the target layer, (2) weak interfaces between the bounding layers and the target layer or (3) a favorable fluid-pressure gradient. Considering a change in each of these conditions in addition to fault orientation and stress regime in a poroelastic model would generate a combined effect on the change of the fault reactivation criterion. Deconvoluting this combined effect could be the objective of a future study.

Using a nonwater-based fracturing fluid (e.g., supercritical CO₂) may seem a robust remedy to reduce hydraulic-fracturing-induced seismicity by replacing water as a nearly incompressible fluid (leading to sharp stress disturbances during hydraulic fracturing) with a more compressible and less dense fluid (i.e., soft fluid). However, Verdon et al. (2010) show that CO₂ and water fracturing under very similar conditions can result in similar induced microseismic event magnitudes, and larger event magnitudes are correlated with higher injection pressures, not the fracturing fluid type. In addition, to the best of our knowledge, nonwater-based fracturing is not commercialized yet on a large scale, mostly due to the complex fluid viscosity effect on proppant transport and HF geometry. Using a less viscous fracturing fluid may lead to larger fractures, which disturb the stresses over a larger fault area and result in more fault reactivation; however, these fractures may not be effective because of less fracturing fluid capability for proppant transport. In contrast, using more viscous fracturing fluid leads to shorter and wider fractures (Settari and Cleary, 1986), which leads to larger stress disturbance over a smaller fault area (Sneddon, 1946; Sneddon and Elliott, 1946; Warpinski and Branagan, 1989), likely resulting in a smaller magnitude event. In the absence of water-saturated intervals, which would impose a risk of producing abundant water through out-of-zone fractures, large HFs (in height and length) are likely favorable from a production point of view. This is because large fractures would open an extended drainage rock volume to the wellbore. Generating these large fractures promotes the use of slickwater instead of highly viscous fracturing fluids. Using slickwater not only leads to a large fracture area (favorable for production) but may also lead to fault interception and direct fracturing fluid discharge into the fault core (unfavorable for induced seismicity). However, due to the operational complications in the use of nonwater-based fracturing fluids (e.g., their complex rheological and proppant transport characteristics), slickwater would remain the most common fracturing fluid.

We provided a physical understanding of the fundamental processes of fault reactivation associated with reservoir stimulation for various operational parameters (the number of the stimulated fractures and the stimulation sequence), the hydromechanical rock characteristics (fault damage zone conductivity and toughness of the target layer, overlying layer, and underlying layer), and the spacing between the fault and HF. This study can guide operators in the informed collection of the subsurface data to be used for the derivation of these suggested factors that influence fault reactivation. Hydromechanical properties of the host rock and the fault zone may be obtained from the core and, in some cases, from the surface exposures of similar rock systems, recognizing the frequent dissimilarity between surface-exposed and reservoir-depth rock units of the same stratigraphic unit. Of greater benefit than comparing surface-to-outcrop systems of stratigraphic equivalence may be the comparison of systems of comparable depositional, diagenetic, and structural history, even in geographically and stratigraphically disparate units. Determining the spacing between the HF stimulation and a seismogenic fault may be equally challenging for smaller “subseismic-scale” faults. A combination of high-resolution seismic imaging and stratigraphic well correlation (Horne et al., 2021) is a prudent step prior to well planning in minimizing the risk of seismic fault reactivation during well completion.

In this work, we developed fully coupled poroelastic models to simulate the reactivation of a normal fault by hydraulic fracturing. We used poroelastic continuum solid elements for the porous media and previously validated pore-pressure cohesive elements to simulate the tensile failure of one or multiple HFs and the fault core during the hydraulic stimulation. The shear failure of the fault plane was assessed spatiotemporally by the CFS criterion, obtained based on the Cauchy stress tensor. We considered three stacked layers at different fracture mechanical properties being crossed by a critically stressed fault under a normal faulting stress regime. For sensitivity analysis, we varied the tensile-mode strength and the toughness of the fractures within the overlying and underlying layers, the spacing between the fault and the HF, the number of HFs, the fracturing strategy (simultaneous versus sequential), and the fault damage zone thickness.

We find that poroelastic stress changes induced by hydraulic fracturing exert a dominant control on fault reactivation, whereas fault dilation by fluid injection into the fault after the fault-fracture intersection may not occur in formations under an anisotropic in-situ stress state. In all nine simulated conceptual scenarios, a relatively large, stabilized region, identified by $CFS<0$⁠, was observed over the fault core in the same depth range of the HF due to the compression generated by the HF. In contrast, an extensive unstable area, identified by $CFS>0$⁠, was detected surrounding this stabilized region resulting from the poroelastic stress disturbance. This poroelastic effect allows for fault slip ahead of the tip of the propagating HF and prior to, or in the absence of, a physical fault-fracture intersection.

We conclude that hydraulic-fracturing-induced fault reactivation is not limited to HFs intersecting a fault thereby resulting in fracturing fluid leaking into the fault. Instead, for a single HF, fault reactivation can occur at a spacing of 200 m from the fault without any fracture-fault intersection. The inclusion of poroelastic effects into seismic hazard assessments thus increases the potential for seismic fault reactivation compared with established procedures that do not include poroelasticity. Increasing the fracture-to-fault spacing lowers the fault reactivation potential due to more uniform compressional stresses over the fault core induced by the hydraulic-fracture dilation. In addition, increasing the number of fracture clusters within the 200 m spacing to the fault leads to a larger reactivated fault zone, thus increasing the potential for fault reactivation, regardless of whether hydraulic fracturing is simultaneous or sequential. Compared with simultaneous fracturing, sequential fracturing leads to (1) a larger stabilized fault zone due to the alternate upward-downward HF growth under stress shadows and (2) a larger area of $CFS>0$ over the fault core due to the poroelastic stress disturbance over a larger area of the fault core. Furthermore, a larger fault-damage-zone thickness and the upward growth of the HF (due to the stratigraphic fracture toughness heterogeneity) lead to a larger reactivated fault zone.

Financial support was provided through the Center for Injection and Seismicity Research (CISR) at the Bureau of Economic Geology, The University of Texas at Austin. The directions and outcomes of this research were not influenced by the industrial associates of this Center. We thank P. Hennings for technical and logistic support. We extend our gratitude to J.-E. Lundstern and the other anonymous reviewers and editors of this paper for their constructive comments. The Texas Advanced Computing Center provided access to high-performance computing. We acknowledge Dassault Systèmes for providing a research license for the Abaqus software program. Publication was authorized by the director, Bureau of Economic Geology.

No data have been required for this paper.