Maximum Magnitude of Seismicity Induced by a Hydraulic Fracturing Stage in a Shale Reservoir: Insights from Numerical Simulations

A key unknown limiting assessment of risk posed by inducing anomalous seismicity during hydraulic fracturing is the potential maximum magnitude of an event. To provide insights into the variation in maximum magnitude that can be induced by a hydraulic fracturing stage, worst-case scenarios were simulated in 2D using coupled hydro-geomechanical models. The sensitivity of the magnitude to the hydro-geomechanical properties of the fault and matrix rock were quantitatively compared through parametric analysis. Our base model predicts a maximum event with moment magnitude (M w ) 4.31 and M w values range from 3.97 to 4.56 for the series of simulations. The highest magnitude is predicted for the model with a longer fault and the lowest magnitude for the model with a smaller Young’s modulus. For our models, the magnitude is most sensitive to changes in the Young’s modulus and length of the fault and least sensitive to changes in the initial reservoir pressure (i.e. pore pressure) and the Poisson’s ratio.


Introduction
The potential risk to the environment and public health and safety resulting from anomalous seismicity induced by hydraulic fracturing has become an increasing concern since the first reported occurrence in 2011 from Blackpool, England [1] [2] [3]. While the research effort has increased with the concern, various challenges still exist to accurately assess the true risk posed by anomal-ous seismicity induced by hydraulic fracturing. One of the key challenges, as well as a subject of considerable debate, is the maximum magnitude of an event that can be induced by a particular hydraulic fracturing operation. Some researchers have suggested that an upper limit exists for the maximum magnitude, which is controlled by the fluid injection operation [4] [5] [6] [7]. Other researchers argue that a deterministic upper bound is only applicable for self-arresting ruptures, whereas in high-stress environments, the rupture may be self-sustaining and the maximum magnitude is only bound by tectonics [8] [9].
In an attempt to better understand maximum magnitude, a parametric analysis was performed to quantify the sensitivity of the magnitude to variations in the hydro-geomechanical properties of the fault and matrix rock. In this study, the fluid injection from a hydraulic fracturing stage is simulated for a series of numerical models with a critically stressed fault within a fractured shale reservoir using the 2D discontinuum based distinct element model, UDEC [10]. The input parameters for the modelling were chosen to represent a worst-case scenario, providing insight into the maximum magnitude earthquake that can be induced by a hydraulic fracturing stage. The fluid injection rate and duration, length of the fault, distance from fault to injection source, maximum and minimum principal stresses, pore pressure, and the Poisson's ratio and Young's modulus of the shale, were systematically varied to end-member values for the parametric analysis. Understanding the sensitivity of the magnitude of induced events will aid in developing protocols for reducing the probability of inducing felt events and mitigation procedures to prevent larger magnitude induced events, as well as help to understand the variations in the magnitude of induced seismicity between regions. Following a brief explanation of the modelling procedure, the results of the simulations are presented, compared, and discussed.

Model
In order to simulate the interaction between natural and hydraulic fractures in combination with fault reactivation during a hydraulic fracturing stage, a coupled hydro-geomechanical approach is necessary. The two-dimensional (2D) distinct element model, UDEC version 5.0 (UDEC5, first presented by [11]), which uses the discontinuum, distinct-element method, was chosen for this study as it is specially designed to simulate fluid flow and the resulting behavior of fractured rock masses. The distinct-element method uses an explicit time-stepping algorithm, which allows for the deformation of the rock matrix blocks and large shear and closure/dilation displacements along the discontinuities between blocks.
In UDEC, Newton's equations of motion are directly solved using explicit time-stepping. The 2D model space is discretized into mutually interacting blocks (i.e. discrete elements) to simulate a discontinuous system with deformable boundaries representing planes for possible failure (i.e. fractures/joints/bedding planes/faults). Both shear and normal displacements are allowed along the shear slip and dilation along pre-specified failure planes caused by disturbance of the initial stress field resulting from fluid injection and the generation of hydraulic fractures. Therefore, the response of the modeled system to fluid injection includes the mechanical effects, due to the changing stress field, coupled with the flow effects, due to changes in fluid pressure. In this study, the pore pressure build-up and diffusion along a fault, near the injection source, within jointed shale, were simulated and the resulting slip along the fault was interpreted as a seismic event. The moment magnitude (M w ), calculated from the magnitude and extent of the slip, was then quantitatively compared for a series of simulations. The input parameters that were investigated during the parametric analysis include: the injection duration and flow rate; the length of the fault and its distance from the injection source; the magnitudes of the maximum and minimum principal stresses and the pore pressure; and the Young's modulus and Poisson's ratio of the shale matrix. Each parameter was systematically varied from the average value in the base model to the maximum and minimum end-member values shown in Table 1.
The input parameter values for our simulation were chosen to represent a possible worst-case scenario where the injected fluid is transported through the newly created hydraulic fractures and opened fracture network to a nearby, critically stressed fault. A 2D numerical model was built, 3 × 3 km in size, containing an orthogonal incipient fracture network, oriented to the principal stress directions, with a 50 m spacing and a fault oriented 30˚ from the maximum principal stress. The 2D model, which is shown in Figure 1, represents a horizontal plane with a strike-slip fault. The initial friction angle of the fractures was set as 30˚, the cohesion 1 MPa, and the tensile strength 0.5 MPa. Strain-softening was included in the models by setting the residual friction angle to 25˚, while the residual cohesion and tensile strength are 0. The fault was modeled having no cohesion, tensile strength, or residual friction and an initial friction angle of 20˚. The large difference between the residual and peak friction angles results in the brittle behavior of the fault. To allow sufficient fluid flow to model a worst-case scenario, the peak permeability of the fractures was input as 1 md, with a minimum value of 0.1 md and the peak permeability of the fault was input as 1 d, with a minimum value of 100 md. The density of the matrix blocks in the model was held constant at 2.5 g/cc throughout the simulations.     Table 1.
The magnitude earthquake induced by the fluid injection was calculated for each model at 5-minute intervals during the injection period and 10-minute intervals after shut-in by assuming all the shear slip along the fault, since the start of injection, was accommodated instantaneously at that time. At each time interval, the length of the fault that slipped > 1 cm and the average of the slips along the patches of that fault length (average slip) were determined to calculate a moment magnitude, M w , using: where M 0 , the scalar seismic moment in Nm, is the product of the shear modulus of the host rock, µ, the area of the fault that slip, A, and the slip along the fault, s: The fault area was calculated assuming a circular fault with a diameter equal to the modeled rupture length, L: The modeled slip and length of fault slip and the calculated M w were compared for the series of simulations providing insights into the relative importance of the hydro-geomechanical parameters on the maximum possible magni- While faults within shales are often argued to be impermeable (ex. [14] [15]), the goal of this study is to investigate the range of maximum possible magnitudes possible for varying hydro-geomechanical parameters. Therefore, instead of predicting multiple smaller events along the fault and/or stable slip, as commonly observed during hydraulic fracturing and in numerical modelling, this study assumes a permeable fault in order to predict the worst-case scenario.
The main uncertainty in the modelling is the 2D simplification, which results in difficulties estimating a representative injection rate and fault plane area. In addition, the fluid flow in UDEC5 is restricted to the discontinuities between blocks and the reservoir is assumed to be fully water saturated with no gas, which also provides limitations to the modelling. Therefore, the goal of this study is not to provide absolute magnitudes, but a comparison of the change in magnitude with hydro-geomechanical parameters.

Results and Discussion
The simulation of a hydraulic fracturing stage near a critically-stressed fault for the base case model predicts an average of 3.90 cm of slip (solid curve in Figure   3(a)) along 806 m of the fault (dotted curved in Figure 3(a)), after 5 minutes of injection, which produces a M w 3.59 event (Figure 3

Length of Fault
Three models were simulated with fault lengths of 500, 1500, and 5000 m for comparison with the base model with a fault length of 1000 m. The model with the 500 m long fault predicts the largest slip (solid blue curve in Figure 9), but the shortest fault rupture length (dotted blue curve on Figure 9), resulting in the lowest peak magnitude event of, M w 4.10 (dotted red curve in Figure 10). The model with the 1500 m fault predicts smaller slip (solid red curve in Figure 9), but longer fault rupture length (dotted red curve in Figure 9) than the models with shorter faults (i.e. 500 and 1000 m) resulting in a greater magnitude event of M w 4.56 (solid red curve in Figure 10). However, while the model with the longest fault length of 5000 m predicts the longest fault rupture length (dotted green curve in Figure 9), the average slip along the fault is much smaller (solid green curve in Figure 9), resulting in a lower magnitude peak event, M w 4.46 (dashed red curve in Figure 10), than for the model with a 1500 m fault. The results indicate that, for the tested model parameters, the magnitude of an induced event increases with fault length, up to a maximum length, after which the smaller slip along longer faults results in smaller magnitude events. Therefore, for the model in this study, the greatest hazard is posed by the moderate sized faults, ~1 -5 km in length.

Distance to Injection Source
Increasing the distance from the injection source to the center of the fault from

Duration and Rate of Fluid Injection
The  Figure 12) and 4.06 (dotted red curve in Figure 12). The injection rate; therefore, has a greater impact on the magnitude than the injection duration. The results also show that a lower magnitude event is predicted later in time if the same amount of fluid is injected with a slower rate over a longer duration. The entire length of the 1 km long fault is not predicted to significantly rupture until shut-in for the two models with 150 m 2 of injected fluid (dotted, blue and light blue curves in Figure 13).

Effective Principal Stresses
As  results show that the peak magnitude is least sensitive to the pore pressure of all parameters tested (Figure 4 and Figure 5). Increasing and decreasing the initial pore pressure by 5 MPa changes the peak magnitude from 4.31 to 4.30 (solid and dotted green curves in Figure 14). Higher initial pore pressures (i.e. 15 MPa) predict larger slip (solid red curve in Figure 15); however, shorter fault rupture lengths (dotted red curve in Figure 15), which trade-off resulting in a low impact   on the magnitude of induced event.
Increasing the differential stress, by increasing the maximum principal stress to 35 MPa from 30 MPa in the base model, predicts a lower peak magnitude event of M w 4.29 (solid blue curve in Figure 14), while decreasing the differential stress, by decreasing the maximum principal stress to 25 MPa, predicts a greater  Figure 14). Similarly, increasing the differential stress by, decreasing the minimum principal stress to 15 MPa from 20 MPa in the base model, predicts a slightly lower peak magnitude event (dotted red curve in Figure 14), while decreasing the differentail stress, by increasing the minimum principal stress to 25 MPa, predicts a greater peak magnitude event of M w 4.36 (solid red curve in F Figure 14). Such results may seem contrary to intuition and previous results from numerical studies, which found that higher differential stresses likely to result in large events. However, in this study, higher differential stresses result in failure of the incipient natural fractures in the models. A failure state for the incipient fractures is predicted before injection for the model with a minimum principal stress of 15 MPa (blue Mohr circle on inset plot in Figure 16) and with a small pore pressure increase for the model with a maximum principal stress of 35 MPa (blue Mohr circle on inset plot in Figure 17). The failure of the incipient fractures dissipates the stress and pore pressures resulting in smaller increases in stress along the fault plane and therefore lower magnitude events. In contrast, the incipient fractures remain intact for the models with lower differential stresses, resulting in a larger pore pressure increase along the fault and thus greater slip and magnitude events.
The models with critically stressed fractures resulting from the effective principal stresses (the model with a pore pressure of 15 MPa, the model with a minimum principal stress of 15 MPa, and the model with a maximum principal stress of 35 MPa) predict greater pore pressure diffusion following shut-in, resulting in a more gradual decrease in the slip (solid red curve in Figure 15, solid blue curve in Figure 16, and solid blue curve in Figure 17) and hence the predicted magnitude of an induced event with time (solid green curve, dotted red curve, and solid blue curve in Figure 14).

Conclusions
A series of 2D hydro-geomechanical models, representing potential worst-case scenarios where fluid is injected by a hydraulic fracturing stage nearby a permeable, brittle, critically-stressed fault, were simulated and compared to provide insights into the sensitivity of the magnitude of seismicity induced by hydraulic fracturing. The sensitivity of the magnitude to fluid injection rate and duration, fault length and its distance to the injector, effective stress state, and shale geomechanics were quantitatively compared through a parametric analysis.