Modeling of fault activation and seismicity by injection directly into a fault zone associated with hydraulic fracturing of shale-gas reservoirs

Abstract We conducted three-dimensional coupled fluid-flow and geomechanical modeling of fault activation and seismicity associated with hydraulic fracturing stimulation of a shale-gas reservoir. We simulated a case in which a horizontal injection well intersects a steeply dipping fault, with hydraulic fracturing channeled within the fault, during a 3-h hydraulic fracturing stage. Consistent with field observations, the simulation results show that shale-gas hydraulic fracturing along faults does not likely induce seismic events that could be felt on the ground surface, but rather results in numerous small microseismic events, as well as aseismic deformations along with the fracture propagation. The calculated seismic moment magnitudes ranged from about −2.0 to 0.5, except for one case assuming a very brittle fault with low residual shear strength, for which the magnitude was 2.3, an event that would likely go unnoticed or might be barely felt by humans at its epicenter. The calculated moment magnitudes showed a dependency on injection depth and fault dip. We attribute such dependency to variation in shear stress on the fault plane and associated variation in stress drop upon reactivation. Our simulations showed that at the end of the 3-h injection, the rupture zone associated with tensile and shear failure extended to a maximum radius of about 200 m from the injection well. The results of this modeling study for steeply dipping faults at 1000 to 2500 m depth is in agreement with earlier studies and field observations showing that it is very unlikely that activation of a fault by shale-gas hydraulic fracturing at great depth (thousands of meters) could cause felt seismicity or create a new flow path (through fault rupture) that could reach shallow groundwater resources.


ABSTRACT
We conducted three-dimensional coupled fluid-flow and geomechanical modeling of fault activation and seismicity associated with hydraulic fracturing stimulation of a shale-gas reservoir. We simulated a case in which a horizontal injection well intersects a steeply dipping fault, with hydraulic fracturing channeled within the fault, during a 3-hour hydraulic fracturing stage. Consistent with field observations, the simulation results show that shale-gas hydraulic fracturing along faults does not likely induce seismic events that could be felt on the ground surface, but rather results in numerous small microseismic events, as well as aseismic deformations along with the fracture propagation. The calculated seismic moment magnitudes ranged from about -2.0 to 0.5, except for one case assuming a very brittle fault with low residual shear strength, for which the magnitude was 2.3, an event that would likely go unnoticed or might be barely felt by humans at its epicenter. The calculated moment magnitudes showed a dependency on injection depth and fault dip. We attribute such dependency to variation in shear stress on the fault plane and associated variation in stress drop upon reactivation. Our simulations showed that at the end of the 3-hour injection, the rupture zone associated with tensile and shear failure extended to a maximum radius of about 200 m from the injection well. The results of this modeling study for steeply dipping faults at 1000 to 2500 m depth is in agreement with earlier studies and field observations showing that it is very unlikely that activation of a fault by shale-gas hydraulic fracturing at great depth (thousands of meters) could cause felt seismicity or create a new flow path (through fault rupture) that could reach shallow groundwater resources.

Introduction
The rapid increase in North American shale-gas energy production has been made possible through new technology development, including extended-reach horizontal drilling and multistage hydraulic-fracture stimulation. But these new technologies have also raised concerns related to a range of local environmental problems (Arthur et al. 2008;Zoback et al. 2010)). One concern, investigated in this study, is whether shale-gas hydraulic fracturing could activate faults and thereby cause seismicity, opening up flow paths for upward fluid leakage and possible contamination of shallow potable groundwater resources (Arthur et al. 2008;Zoback et al. 2010;Davies et al., 2013;Rutqvist et al., 2013).
A first modeling study to investigate the potential consequences of fault reactivation during shale-gas hydraulic fracturing operations was presented in Rutqvist et al. (2013). Consistent with field observations, the study showed that a hydraulic fracturing operation to stimulate a deep shale-gas reservoir could only give rise to limited fault rupture, along with the possibility of (unfelt) microseismicity. In another study, Flewelling et al. (2013) used injection data and elastic fracture volume and length relationships to bound fracture-height data from 12,000 hydrofracturing stimulations conducted across North America. The hydraulic fracturing data showed that all microseismic events occurred less than 600 m above well perforation, although most were very much closer, and the farthest were usually associated with faults. These studies indicated that shale-gas hydraulic fracturing at great 2 depth (thousands of meters) could not create flow paths for leakage to reach shallow groundwater resources.
Studies have also concluded that the likelihood of inducing felt seismicity during shale-gas hydraulic fracturing operations, while not to be ruled out completely, is extremely small (National Research Council, 2012;Davies et al., 2013). Indeed, after hundreds of thousands of shale-gas fracturing stages conducted to date, only three examples of felt seismicity have been documented (Davies et al., 2013). In Lancashire County, UK, two seismic events of Richter scale magnitude M L = 2.3 and 1.5 were likely induced by direct injection into a fault zone that had not been previously mapped (De Pater and Baisch, 2011). In another case at the Eola Field of Garvin County, Oklahoma, in January 2011 (Holland, 2011), there was a clear temporal correlation between the time of stimulation and the occurrence of 43 earthquakes that ranged in magnitude from M D = 1.0 to 2.8 (M D is the duration magnitude).
Finally, the third case of felt seismicity occurred at Etsho and Kiwigan fields in Horn River, Canada, where 19 events between M L = 2 and 3 occurred having a clear temporal correlation with the shale-gas operation; the largest (and felt) event, occurring in May 2011, had a magnitude of M L = 3.8 (BC Oil and Gas Commission, 2012;Davies et al. 2013). Each of these three cases of felt seismicity have been associated with reactivation of faults.
The biggest modeling uncertainty in the previous fault-activation modeling by Rutqvist et al. (2013) was a 2D simplification of the full 3D field settings. In 2D plane-strain simulations, it is difficult to estimate a representative injection rate, and some assumptions have to be made about the shape of the rupture area (e.g., circular with diameter equal to 3 2D rupture length), which affects the calculated seismic magnitude. In this study we conduct, for the first time, a full 3D model simulation of fault activation associated with shale-gas fracturing. In such a 3D model simulation, the exact injection rate from the 3D field is a direct model input, and the seismic magnitude can be evaluated directly from the calculated rupture area and mean slip without the model uncertainties inherent in a 2D simplification. In this new 3D modeling study, we simulate the case in which a horizontal well intersects a subvertical fault, which then can be reactivated by injection directly into the fault. In addition, we investigate some issues not addressed in the previous 2D modeling in Rutqvist et al. (2013), including how the results correlate with fault and injection depth, fault dip, and fault frictional properties. We conclude with a discussion relating our modeling results to field observations and attempt to explain under which conditions a shale-gas fracturing stimulation could induce a felt seismic event.

Model Setup
We adopted the modeling approach that was applied in the previous 2D modeling study in Rutqvist et al. (2013). That is, we used the coupled multiphase fluid-flow and geomechanical simulator TOUGH-FLAC (Rutqvist, 2011) to model water-injection and fault responses, and we applied seismological theories to estimate the corresponding seismic magnitude. The fault was modeled as a discrete feature using finite thickness elements having anisotropic elasto-plastic properties. Shear failure was governed by a Mohr-Coulomb constitutive model with strain-softening frictional strength properties, consistent with a seismological slip-weakening fault model (Cappa and Rutqvist, 2011). 4 This allowed us to model sudden (seismic) slip events and to estimate their seismic magnitude. The adopted modeling approach has also been extensively applied for modeling fault activation associated with underground CO 2 injection (e.g. Cappa and Rutqvist, 2012;Mazzoldi et al., 2012;. The model domain and the material properties are presented in Fig. 1 and Table 1, respectively. We model a full 3D-geological system (x, y, z: 2 km×10 km×2 km) generally tuned towards conditions that could be encountered in the Marcellus shale-gas play in the Northeastern U.S. This includes model input of in situ stress, fluid pressure, temperature, material properties, and injection rates. In a base-case simulation, we adopt conditions consistent with areas where the Marcellus shale is located at a depth of about 2000 m (6562 ft). The model is representative of the Marcellus shale-gas play with a 30 m thick gasbearing shale, bounded at the top and bottom by other low-permeability formations (such as inorganic gray shale and limestone). This system is intersected by a steeply dipping fault, which in the base case has a dip of 80°. We simulate a case in which the horizontal injection well intersects the fault, and we directly inject the fluid volume related to a 3-hour hydraulic fracturing stage directly into the fault.
We set the initial conditions assuming linear pore pressure and temperature gradients (9.81 MPa/km and 25˚C/km, respectively), with constant hydraulic boundary conditions (i.e., open to fluid flow), except for the planes x = 0 and y = 0 where a no-flow condition is applied (Fig. 1). Mechanical boundary conditions are null displacement at x = 0 and y = 0 planes, and constant stress elsewhere. The initial stress field is selected to represent the 5 conditions at the Marcellus shale play as detailed and justified in Rutqvist et al. (2013). We first set the vertical stress gradient (maximum principal stress) to 26,487 Pa/m, corresponding to an overburden density of about 2700 kg/m 3 . We then consider the minimum principal stress to be horizontal and oriented parallel to the horizontal well, which would lead to vertical hydro-fractures perpendicular to the well, but which in this case follow the weak planes of the fault. This does also correspond to a normal faulting stress field, in which the minimum horizontal stress (and minimum principal stress) is directed normal to the strike of the fault. We set the magnitude of the initial minimum horizontal stress corresponding to a horizontal-over-vertical stress ratio of R =  h / V = 0.6.
There are uncertainties in the horizontal-over-vertical stress ratio and, as highlighted by Rutqvist et al. (2013), this ratio has an impact on the magnitude of fault shear activation.
However, several sources (e.g. Cipolla et al., 2010) indicate a fracture closure stress of about 0.7psi/ft and this corresponds to a horizontal-over-vertical stress ratio of R =  h / V = 0.6. In this study, we keep the horizontal-over-vertical stress ratio fixed at R =  h / V = 0.6, but vary the depth of the system, which also means a variation in stress magnitude at the depth of the injection. The magnitude of the intermediate stress, which in this case of a normal faulting stress regime would be oriented parallel to the fault strike does not affect the potential for shear failure along the fault.
Another important parameter in our analysis is the shear strength of the fault and how it evolves along with the reactivation. Here, we use the strain-softening Mohr-Coulomb model, in which the coefficient of friction and cohesion decreases with slip, i.e., once the peak shear strength is achieved and the fault slips, the cohesion drops to zero and the coefficient of friction drops to a residual value. In the numerical model, this is simulated by reducing the coefficient of friction and cohesion from peak to residual values over a plastic shear strain of 10 -3 (Cappa and Rutqvist, 2011). In the base case, we use a coefficient of friction of µ = 0.6, with a residual value (after slip) equal to µ R = 0.4, whereas the cohesion drops to zero from an initial value of 1 MPa. A larger difference between the peak and residual friction values represents a more brittle behavior that is expected to lead to a larger shear-stress drop and seismic event. The selection of the frictional coefficient parameters are also discussed and justified in Rutqvist et al. (2013), acknowledging that this is one possible set of reasonable values of the frictional coefficient. The fault shear strength and how it weakens with slip is defined by a set of parameters that are varied in this study Other fault properties as well as properties of the shale listed in Table 1 are equivalent to those used and justified in Rutqvist et al. (2013). In this study, we assume that the fault is nearly impermeable (hydraulically indistinguishable from the host rock), though the permeability and porosity can increase as a result of fracturing and shear. We consider a nearly impermeable fault a realistic assumption in this case. As pointed out by Flewelling et al. (2013), hydrocarbons cannot accumulate where there are permeable faults serving as pathways for buoyant oil and gas to leak upward. A relevant example of an impermeable fault in shale is a fault zone in Opalinus Clay exposed at the Mont Terri Rock Laboratory, Switzerland (Croisé et al., 2004). This zone is several meters thick, consists of intensively fractured rocks, has an inferred shear offset of 5 m, but is still hydraulically indistinguishable from the host rock, having an estimated permeability k = 210 -20 m 2 7 (Croisé et al., 2004). Consistent with such a conceptual model, we considered 3.3 m thick fault where fracturing can take place along weak planes such as healed fractures along the fault plane, meaning that the tensile strength is low and taken as 1.7 MPa, based on the extension of the Mohr-Coulomb failure envelope into tension. Such a low tensile strength is reasonable, considering published shale strength data, including laboratory measurements on bedding planes and tectonic shears in Opalinus Clay (Bock, 2009;Young et al., 2010), and past laboratory experiments on Marcellus black-shale samples (Heard and Lin, 1986).
We consider porosity and permeability changes in the fault with tensile and shear rupture, according to a conceptual model of fracturing and shear along a fault zone rather than opening of a single hydraulic fracture. In this model, as soon as tensile or shear failure propagates within the fault plane, the porosity increases with plastic tensile strain p t  and plastic shear strain a P S  according to: where  is the shear dilation angle. The approach of changing porosity with rupture is associated with the constitutive ubiquitous joint model, where plastic tensile and shear strain can be extracted and used for superimposing porosity changes by the rupture. The porosity increase due to the fracturing provides additional fluid storage that is important to consider in the fluid-flow simulation and significantly affects the fracturing propagation. 8 Similarly, the tensile and shear rupture also provides permeability that is superimposed on the initial fault permeability. The adopted permeability-change model described in Rutqvist et al. (2013), considers changes in equivalent fractured rock permeability along the fault as a function of plastic strain normal to the fault plane: where k 0 is the initial fault permeability, A is a constant, and t n  is a threshold strain related required crack opening displacement for onset of permeability changes. Here, we used t n  = 110 -4 and A = 110 -5 , meaning that permeability would increase to about 110 -14 to 110 -11 m 2 for a plastic strain normal to the fault on the order of 110 -3 to 110 -2 . This is a very substantial permeability change from an initial permeability of 110 -19 m 2 , one that provides rapid pressure diffusion along the fractured elements with the fracture propagation.
We simulated a typical hydraulic fracturing stimulation stage with injection at a rate of 200 kg/s (75 bbl/min) for 3 hours. According to Fig. 1, the horizontal section of the well is explicitly represented and connected to elements at the fault to simulate injection focused at a packed-off section of the well. We assumed that a fracture had already been initiated and extended within the fault to about 10 m from the well. Such fracture initiation would in the field be the result of an initial pressurization to formation breakdown pressure (e.g. Mayerhofer et al., 2011). We assume this initial fracture because the model is not refined enough around the well to properly simulate fracture initiation and formation breakdown. 9 In the field, the injection is typically shut down for a few minutes before beginning the main stimulation stage (e.g., Mayerhofer et al., 2011), which in our simulation is assumed to be 3 hours long.
During the stimulation, the fracture is extended within the fault plane, as illustrated in Fig.   1. We then evaluate magnitude and frequency of microseismic events as well as the total extent of the rupture zone. The seismic moment M 0 , and the moment magnitude M w are evaluated for newly ruptured patches on the fault (related to the propagation of the fracturing front). According to Kanamori and Anderson (1975) and Kanamori and Brodsky (2001)

Modeling results
In the following subsections, we first present the results of the base-case simulation, considering injection at 2000 m depth into a fault dipping 80 degrees, and having a coefficient of friction µ = 0.6 that upon reactivation drops to a residual value µ R = 0.4. We then present parameter studies to evaluate the sensitivity of the results to changes in injection depth (1000 to 2500 m), fault dip (75 to 90), and fault friction coefficients (0.2 to 0.8). In each case, we evaluate the results in terms of number and maximum magnitude of microseismic events, and the total vertical rupture length along the fault.  Fig. 2a shows how the injection pressure increases to reach a peak value of about 43 MPa and then slowly decreases along with the failure propagation. Fig. 2b shows that shear displacement is initiated after a few minutes of injection-when the pressure exceeds the stress normal to the fault-and then slips gradually along with the propagation of the failure zone. The biggest microseismic events occur during the first few minutes of injection, and then microseismicity continues during the 3-hour injection at reduced magnitudes and frequency-but also continues for at least 10 minutes after terminating the injection (after 180 minutes). The simulated magnitude of the microseismic events generally ranges between M w = -2 and 0. Fig. 2c to e show that at the end of the 3hour injection, the failure zone has extended to a radius of almost 200 m from the well upwards and laterally, while only about 100 m downwards. This can be explained by the vertical gradient in stress that makes it easier to propagate a fracture upwards. Moreover, at the end of the 3-hour injection, the fault slip is almost 5 cm at the center of the failure zone (near the injection well). The simulation in Fig. 2 captures results typical for hydraulic fracturing, including (1) small pressure fluctuations, especially during the first 20 minutes when failure and shear slip occur relatively close to the injection well, (2) a gradual reduction in the pressure required to propagate the fracture for a given injection rate, and 11 finally, (3) a typical pressure fall-off after terminating the injection. After terminating the injection and shutting in the well (at 180 min), the well pressure declines to below the stress normal to the fault plane (32.5 MPa), in a shut-in pressure behavior that could be used in the field to estimate the minimum principal stress magnitude.
Overall, the result from the base-case simulation shows that maximum moment magnitude is less than 0.5, and the total rupture is less than 200 m up from the injection. Field observations show that during typical stimulation operations in the North American shalegas basins (Marcellus, Barnett, Eagle Ford, Woodford, Haynesville, Horn River), the maximum moment magnitude in a stimulation stage ranges between M W -2.5 to (in rare cases) up to M W = 1.0, with the high-end values usually associated with faults (Warpinski, 2009;Warpinski et al., 2012). The simulated microseismic events are within this range, with the lower limit probably affected by the minimum size of element patches in the fault (10 by 10 m) that can rupture in one instance. This means that in our modeling, we model the relatively larger microseismic events occurring from shear slip along the fault plane, whereas in the field, there are numerous smaller-magnitude events perhaps occurring as a result of slip in small-scale fractures in the host rock surrounding the fault, and these are not resolved in our modeling.

Sensitivity to injection depth
In their review of thousands of hydraulic fracturing treatments in North American shale basins, Warpinski et al. (2012) noticed an increased microseismic magnitude and frequency with depth. In Fig. 3, we present an analysis in which we varied the injection depth from 12 1000 to 2500 m (3281 to 8202 ft). The trend is clear and shows an increasing maximum magnitude and frequency with depth (Fig. 3f). Note that, from Fig. 3a, b, and d, although the total shear slip magnitude increases with depth, the extent of the rupture zone is equivalent, i.e., always about 200 m up from the injection point. We attribute the increased total shear-slip magnitude with depth to a higher shear stress on the fault when going deeper. A higher shear results in a greater shear-stress drop upon reactivation, and the shear slip is proportional to the shear-stress drop for a given fracture extent and shear modulus.
The higher shear stress can also explain the increased moment magnitude of seismic events with depth. Overall, such depth dependency is consistent with field observations in Warpinski et al. (2012). longer for inclined faults, as well as both shear slip and microseismic magnitudes increase, because the initial shear stress is much higher (e.g., 5.3 MPa for the 75 fault compared to 0 for the 90 fault). However, even for an inclined fault, when the shear stress is the highest, the event magnitudes are less than about M W = 0.7.

Sensitivity to slip-weakening model parameters
The coefficient of friction and its evolution during failure may affect the amount of stress drop and thereby the amount of slip and seismic magnitude. In the base case, coefficient of frictions µ = 0.6 and µ R = 0.4 were used. A coefficient of friction of 0.6 is commonly applied in fault-stability analyses and has been identified as a lower limit value for the most common rocks (Zoback, 2007). However, clay-rich fault rock and fractures in shale could have a much lower coefficient of friction, especially under wet conditions (Zoback, 2007;Samuelson et al., 2012). Another data source are the investigations of Opalinus Clay at the Mont Terri underground research laboratory, which indicates a coefficient of friction of about 0.4 along bedding planes (Bock, 2009;Young et al., 2010). One additional complicating factor is that the coefficient of friction is rate dependent and is also dependent of the clay content (Kohli and Zoback, 2013). Kohli and Zoback (2013) determined coefficient of friction on fractures in Barnett, Haynesville, and Eagleford shale samples to be in the range of 0.4 to 0.8, for high to low clay content. Moreover, unstable slip behavior was noted for shale with lower clay content, i.e. when the initial coefficient of friction was in the range of 0.6 to 0.8. . 14 Because of the uncertainty in the selection of coefficient of friction, we conducted sensitivity studies involving both peak and residual friction values. Moreover, we varied the rate of frictional strength drop with plastic shear strain. In the base case it was assumed that the coefficient of friction drops from µ = 0.6 to µ R = 0.4 over a plastic shear strain of 10 -3 .
As discussed in Cappa and Rutqvist (2011), shear strength drop with plastic shear strain can be related to a fault slip weakening model, commonly applied in earthquake seismology, and the evolution of shear stress with slip (Kanamori and Bodsky, 2004). A plastic shear strain of 10 -3 across the fault width of 3.33 m corresponds to a fault shear slip of 3.33×10 -3 m (3.33 mm). This would correspond to the critical slip distance, Dc = 3.33×10 -3 m in a slip weakening model, meaning that the coefficient of friction reduces to its residual value (dynamic friction) over a slip displacement of 3.33 mm (Aochi et al., 2014). In the base case, the maximum shear strain after the 3-hour injection is on the order of 10 -2 and a shear maximum slip is up to 5 cm (Fig. 2). However, the incremental plastic shear strain and slip during each event is much smaller, with slip magnitudes on the order of 1x10 -5 to 5x10 -4 m for event magnitudes ranging from about M w = -2 to 0. This means that the full frictional strength drop from µ = 0.6 to µ R = 0.4 will not occur and this may in turn affect the shear stress drop during an event. Although a critical slip distance of 3.3 mm is within the range of laboratory measurements on rough fractures (Ohnaka, 2013), the value could vary by orders of magnitude depending on scale (Ohnaka 2003) and might be strongly heterogonous in a fault plane (Aochi et al., 2014). Moreover, recent laboratory data on shale by Kohli and Zoback (2013) indicated values on the order of tens of microns, i.e. very small values. Again, the values could vary by orders of magnitude when going from core sample scale to field scale (Ohnaka, 2003), leading to considerable uncertainty in this parameter.
Acknowledging the uncertainties we investigate the sensitivity to changes in peak and residual friction coefficient and the critical slip distance. Figure 5f shows that there is a trend of fewer but higher magnitude microseismic events if reducing the critical slip distance to as low as 10 -5 m. With such a low critical slip distance, the full shear strength drop corresponding to a drop in the coefficient of friction from µ = 0.6 to µ R = 0.4 would occur as soon as shear failure occurs. Moreover, quite surprisingly, the analysis showed the rather small influence of peak and residual friction coefficient, for variations in the initial coefficient of friction between µ = 0.5 to 0.8 and for residual coefficient of friction between µ R = 0.3 to 0.6. It is only when the residual coefficient of friction is reduced to 0.2 that we see more significant change in calculated moment magnitudes and rupture length, including a maximum moment magnitude of M W = 2.3 (Figure 6f). These relatively larger events can be observed in Fig. 5d as more significant slip displacements occur at 70, 100, and 160 minutes. The M W = 2.3 event is of similar magnitude as those felt at the aforementioned U.K. and Oklahoma incidents, though such small magnitude events may likely go unnoticed. Note that if the residual coefficient of friction is µ R = 0.2, then the shear stress on the fault from the tectonic stress field (about 3.6 MPa) is higher than the residual shear strength of the fault (about 2.5 MPa). This means that once fault activation is initiated, shear strength within the rupture area will drop to below the prevailing shear stress, and because of this weakening, the shear slip can self-propagate outside the zone of pressurization. This is the reason why in the case of µ R = 0.2, the shear rupture is somewhat 16 larger. Nevertheless, this result indicates that the only way to produce a microseismic event above M W = 2 in the current model setting is for the residual fault shear strength to be less than the prevailing shear stress.

Discussion
Overall, the results of this 3D simulation study are in agreement with the findings in the previous 2D study in Rutqvist et al. (2013), in terms of potential extent of the rupture and magnitude of seismic events. The moment magnitudes calculated in this study ranged from about M W = -2.5 to 0.5, except for the special case of very low residual coefficient of friction of 0.2, when a magnitude M W = 2.3 event was calculated. Over the duration of the 3-hour injection, repeated events and aseismic slip were shown to amount to up to 0.06 m, with the total radius of rupture extending up to 200 m. Consistent with field monitoring, the modeling shows that when faults are present, somewhat larger seismic events are possible-compared to those associated with regular hydraulic fracturing seismic events-because larger surface areas are available for rupturing (Fisher andWarpinski, 2011: Warpinski et al., 2012). Indeed, regular hydraulic fracturing seismic events average around M W = -2.5, whereas higher than usual events up to about M w = 0.5 are usually associated with fracturing along faults (Warpinski et al., 2012).
Moreover, the modeling results of repeated small microseismic events and aseismic slip and fracture opening is consistent with field observations-that the energy emitted from microseismic events represents only a small fraction of the energy input or the energy to 17 open the fracture (or fracture network) (Warpinski et al., 2012). Natural fractures present in shale (Gale et al., 2014) may also open aseismically by increased fluid pressure and associated reduction of effective stress (Rutqvist, 2014). In this context, weakly cemented fractures being common in shales (Gale et al, 2014) could have a permeability that is extremely stress dependent. Also consistent with our modeling, field observations of fracturing along faults have shown that microseismic events can travel several hundred meters upwards (Fisher and Warpinski, 2012;Flewelling et al., 2013;Yang and Zoback, 2014).
Our analysis showed that the frequency and magnitude of the microseismic events are dependent on the initial shear stress on the fault plane, while the total rupture length is limited by the volume of water injected. In the case of injection directly into the fault plane, the shear activation (shear failure) occurs simultaneously with the hydraulic fracturing (tensile failure), and these events takes place starting near the well and propagate away from the well bore in repeated microseismic events. At greater depth and for inclined faults that are more optimally oriented for higher shear stress, these microseismic events are of slightly higher magnitude. This depth dependency is also consistent with field observations at North American shale gas basins (Warpinski et al., 2012). An explanation from the modeling is that upon reactivation, shear stress is relieved, and this shear-stress drop is greater in the case of a higher initial shear stress, which in turn can results in larger slip for a given rupture area.
The total rupture length obtained in this study (about 200 m radius), will depend on the total injection volume and leak-off from the main hydraulic fracturing into surrounding 18 rock, as well as the initial gas saturation. In this case the conceptual model is a 30 m thick gas bearing formation surrounded by gray shale and limestone layers that does not contain gas, and in this analysis we simplified the system to be fully water saturated. If considering some gas saturation in the gas bearing formation (e.g. 50%), the hydraulic fracturing process could be affected by the fact that gas is much more compressible than water therefore delaying the pressurization for a given injection rate (e.g. . However, as shown in Kim et al. (2014), hydraulic fracturing of shale gas reservoirs may be affected by complex two-phase flow processes, including vertical gravity segregation within the created fractures. To investigate these effects Kim et al., (2014) conducted a hydraulic fracturing simulation of a shale-gas reservoir for the extreme cases of water saturated and gas saturated rock, and achieved similar radial extent of the stimulated fracture. Thus, we would not expect any major extension of the calculated rupture zone if considering some initial gas saturation in 30 m thick gas bearing formation.
Nevertheless, the study in Kim et al., (2014) shows that the created fracture volume may not be equal the injected volume due to two-phase flow effects and fluid leak-off to the surrounding rock. Here, we injected 2160 m 3 (571,611 gallons) of water over the 3-hour injection, which can be considered a typical average injection volume per fracturing stage (e.g., compared with data shown in Mayerhofer et al. (2011) andBC Oil andGas Commission (2012)). If injection were to continue or were at a higher rate, the rupture would propagate slightly longer. However, this also means that the total length of the rupture zone can be controlled by limiting the total injection volume and guided by seismic monitoring.
One important question is, under what conditions could it be possible to induce felt seismic events? We would emphasize that only three such cases have been documented among hundreds of thousands of hydraulic fracturing treatments, meaning that these are very rare events, and hence it might be difficult to produce such events even in a numerical model.
Our analysis for the conditions considered in this study indicates that the only way to produce events on the order of magnitude M W = 2 to 3, is a very brittle fault with a residual coefficient of friction lower than the prevailing shear stress on the fault (see results for a residual coefficient µ R = 0.2 in Fig. 5). In this case, sufficient new surface area can be ruptured in one instance to produce a seismic event on the order of M W = 2 to 3. Moreover, consistent with field observations at the three known cases of felt events, these simulated events occurred hours into the hydraulic fracturing stage, i.e., after 70, 100, and 160 minutes, in Fig. 5d. Each simulated event in Fig. 5d resulted in up to 1 cm fault slip, and the total slip is greater than 6 cm, which would certainly put some strain on the horizontal well.
Deformations of the horizontal wells were also observed at the incidents at Lancashire County, U.K. and Horn River, Canada. However, some of the shale-gas-fracturing-related events reported in the literature, e.g., the events at the Eola Field, Oklahoma, occurred many hours after shut-in of the injection and farther away from the injection wells. This may indicate channeling of injection fluid that, after some time, then enters a critically stressed fault that then can be reactivated. Such a scenario is more in line with the conceptual model in previous 2D modeling in Rutqvist et al. (2013), where we considered a hydraulic fracture propagating to intersect a fault, not injection directly into the fault. In that study, such a fault could be reactivated many hours after continuous injection to produce larger than usual microseismic events (e.g., magnitude 2). That is, if that fault had 20 some permeability, were near critically stressed, and had brittle slip properties with a low residual strength. Thus, it could be that if fluid pressure migrates out of the shale layer it might encounter more brittle, faulted rock types that could be a source of larger than usual events.

Conclusions
We have conducted 3D model simulations of fault activation during a shale-gas hydraulic fracturing operation that eliminates some important model uncertainties associated with previous 2D modeling studies. We simulated a case tuned toward the conditions at the Marcellus shale simulation, a hydraulic fracturing stage of 2000 m, but with direct injection into a steeply dipping fault. Thus, the entire water volume from one hydraulic fracturing stage was injected to propagate a fracture along the fault. The results of the simulations are consistent with field observations regarding the possible microseismic magnitude and rupture length, indicating the occurrence of small (unfelt) microseismic events and fault ruptures that could propagate upwards a few hundred meters from the injection interval.
The following main results were achieved in this modeling study:  During the stimulation injection, the shear activation (shear failure) occurred simultaneously with the hydraulic fracturing (tensile failure), starting near the well and propagating away from the well in repeated microseismic events.
 The moment magnitudes ranged from about M W = -2.5 to 0.5, except for the special case of a very brittle fault with a residual shear strength below the prevailing shear stress, when a magnitude M W = 2.3 event was calculated. 21  The microseismic magnitude increased with depth and for inclined faults that are more optimally oriented for higher shear stress, indicating the dependency of shear-stress drop on the seismic magnitude.
 Over the course of the 3-hour injection, repeated events and aseismic slip amounted to up to 0.06 m, with the total radius of the shear rupture extending up to 200 m.
Overall, this 3D modeling study in agreement with previous findings showing that it is very unlikely that activation of a fault by shale-gas hydraulic fracturing at great depth (thousands of meters) could cause felt seismicity or create a new flow path (through fault rupture) that could reach shallow groundwater resources.
The analysis indicates that the very rare observations of felt seismicity associated with shale-gas fracturing (three documented out of several hundred-thousand fracturing stages) might have been caused by locally very unfavorable conditions, wherein fluid-pressure changes could reach a large section of a fault that was close to critically stressed for shear and having very brittle slip weakening properties in which the residual strength was lower than the prevailing shear stress.
Though our results in terms of seismic magnitudes are in agreement with field observations at various U.S. shale-gas plays exposed to various stress regimes, our analysis was for the case of steeply dipping faults under a normal-faulting stress regime at 1000 to 2500 m depth. It might be different when operating a lot deeper such as could be the case in future shale gas production in China. In such case, the stress field will be much higher, the shale 22 even tighter and potentially more brittle, and therefore potentially leading to larger magnitude events.
Nevertheless, adequate site characterization for identifying and avoiding faults should be a priority in any shale gas development. Moreover, continuous monitoring of induced seismicity from the start of the injection can be used to detect any runaway fracturing along faults. That is, to detect fracturing that propagates far above and below the gas-bearing shale layer, indicating potential reactivation of a fault.     Figure 5. Sensitivity of simulation results to changes in the critical slip distance (or equivalent plastic shear strain threshold) for shear strength drop : Vertical extent of (a) plastic tensile strain (signifying tensile failure) and (b) plastic shear strain (signifying shear slip), and time evolution of (c) well pressure and (d) slip magnitude, and (e) magnitude distribution for each value of the plastic shear strain threshold and (f) frequency and maximum magnitude as a function of plastic shear strain threshold. Figure 6. Sensitivity of simulation results to changes in residual coefficient of friction: Vertical extent of (a) plastic tensile strain (signifying tensile failure) and (b) plastic shear strain (signifying shear slip), and time evolution of (c) well pressure and (d) slip magnitude, and (e) magnitude distribution for each value of the residual friction and (f) frequency and maximum magnitude as a function of residual coefficient of friction.