Investigating the role of elastostatic stress transfer during hydraulic fracturing-induced fault activation.

SUMMARY We investigate the physical processes that generate seismicity during hydraulic fracturing. Fluid processes (increases in pore pressure and poroelastic stress) are often considered to be the primary drivers. However, some recent studies have suggested that elastic stress interactions may signiﬁcantly contribute to further seismicity. In this work we use a microseismic data set acquired during hydraulic fracturing to calculate elastic stress transfer during a period of fault activation and induced seismicity. We ﬁnd that elastic stress changes may have weakly promoted initial failure, but at later times stress changes generally acted to inhibit further slip. Sources from within tight clusters are found to be the most signiﬁcant contributor to the cumulative elastic stress changes. Given the estimated in situ stress ﬁeld, relatively large increases in pore pressure are required to reach the failure envelope for these faults—on the order of 10 MPa. This threshold is far greater than the reliable cumulative elastic stress changes found in this study, with the vast majority of events receiving no more than 0.1 MPa of positive (cid:2) CFS , further indicating that elastic stress changes were not a signiﬁcant driver, and that interaction with the pressurized ﬂuid was required to initiate failure. Thus, cumulative stress transfer from small events near the injection well does not appear to play a signiﬁcant role in the reactivation of nearby faults.


I N T RO D U C T I O N
Hydraulic fracturing (HF; Clarke et al. 2014;Bao & Eaton 2016), waste water injection (Keranen et al. 2013), carbon capture and storage (CCS; Verdon et al. 2013a) and the enhancement of geothermal systems (EGS; Schoenball et al. 2012) have all shown the ability to stimulate seismicity. Numerous felt earthquakes have been linked to injection across North America and Europe, with magnitudes reaching as high as M W 5.8 (Keranen & Weingarten 2018). Increases in pore pressure and poroelastic stresses are considered to be the dominant physical mechanism responsible in areas of anthropogenic seismicity, reducing the effective normal stresses acting on nearby faults and fractures and thus promoting failure (Raleigh et al. 1976). However, the extent to which faulting and fracturing near the well can promote more distant seismicity is not well established. The uncertainty surrounding the continued propagation of seismicity once initial faulting has occurred has significance for the mitigation strategies currently being used in certain jurisdictions. 'Traffic-light' schemes function under the assumption that if a smaller event (e.g. M 0.5 in the UK) is detected during injection, halting further operations can act to mitigate a larger, possibly damaging, earthquake (Green et al. 2012;Clarke et al. 2014). This assumption hinges on a transient change to the stress state produced from the fluid permeating the rock mass, as well as stress being relieved from the small faulting and fracturing that immediately accompanies fluid injection. However, if nearby faults and fractures receive enough elastic Coulomb stress change from the preceding events, bringing them closer to failure, by the time injection is halted due to the traffic light being exceeded, it may already be too late. In tectonic settings the transfer of static stress is widely used in the estimation of seismic hazard, where it is applied to the distribution of aftershock sequences (Stein 1999;Steacy et al. 2004). These Coulomb models have also been applied to fluid related triggering processes such as the seismicity associated with magmatism (Toda et al. 2002). Recent studies (e.g. Catalli et al. 2013;Sumy et al. 2014;Pennington & Chen 2017) have concluded that elastic stress transfer between induced events can also play a significant role during a sequence of induced seismicity.
Displacement on a fault changes the surrounding stress field. This stress change will act on nearby faults and fractures, and can promote or inhibit further slip. The Coulomb failure stress change, CFS, acting on a fault is an effective tool for examining the stress accumulated during a sequence of events: with change in shear stress τ and change in normal stress σ n (positive extensive) resolved onto the fault plane, and effective co-Stress transfer during hydraulic fracturing 1201 efficient of friction μ (Harris 1998). A positive CFS means the stress conditions have changed in such a way as to promote slip, whereas negative acts to inhibit slip. Here, where μ is used to summarize the complex assumptions relating pore pressure change, the frictional properties of the fault and the effective normal stress. Skempton's ratio β (Skempton 1954) is used to relate the change in pore pressure change to the change in the applied stress: P = −β σ kk /3. Assuming that the fault zone is homogeneous and isotropic, one can relate the diagonal components of the stress tensor σ kk to the normal stress σ n (Rice 1992;Simpson & Reasenberg 1994). This gives σ kk /3 = σ n , and eqs (1) and (2) are satisfied (Harris 1998). This β parameter can be thought of as simulating the frictional stability of a fault in a porous medium.
With β > 0, pore pressure changes act against the compressive stress, reducing the effective friction on the fault. The elastic Coulomb stress changes acting on a plane due to any previous slip on surrounding faults can then be calculated using the Okada equations (Okada 1992), identifying if slip was encouraged or discouraged by the preceding fault movements. The overall effect of stress changes acting on an earthquake population is measured using the Coulomb index (CI)-the percentage of events with positive CFS (Hardebeck et al. 1998). A CI significantly above 50 per cent is consistent with a population where elastic stress changes actively promoted slip.
This appears to be a relatively simple procedure, but determining the significance of a positive CFS signal for a large number of events is non-trivial. One should appropriately consider the uncertainties in fault orientation, material properties and slip behaviour. A study by Meier et al. (2014) highlighted the importance of a rigorous statistical treatment of static stress modelling by conducting a thorough analysis of the stress changes associated with the 1992 Landers earthquake sequence, assessing the significance of positive stress changes with regard to the modelling assumptions and the uncertainties in the input parameters. They explored techniques to determine significance and the sensitivity of elastic stress models, especially to focal mechanism uncertainty. Further studies have highlighted the effect of elastic stress transfer in some injection-induced earthquake sequences, with relatively high magnitude events correlated with positive Coulomb stress changes. Catalli et al. (2013) calculated CI in excess of 75 per cent for induced events associated with the Basel EGS, as well as giving a treatment of uncertainties in focal mechanism (FM), location, and frictional parameters. During the 2011 Prague, OK sequence, Sumy et al. (2014) found that the M W 5.0 foreshock, which is considered to have been induced by the nearby waste water injection, promoted the M W 5.7 main shock by elastic Coulomb stress transfer. Pennington & Chen (2017) also studied the events associated with waste water injection in Oklahoma, during the 2016 Pawnee sequence. They found that the three foreshocks promoted failure of the M W 5.8 main shock, which went on to promote failure on 67 per cent of the events that occurred after it. Coulomb stress changes resulting from microseismic events recorded during HF were calculated by Vasudevan & Eaton (2011), showing that non-double-couple source mechanisms have a significant effect on the resulting stress maps. Schoenball et al. (2012) investigated the static Coulomb stress changes for the seismicity associated with an EGS, neglecting any pore pressure change. They found that elastic Coulomb stresses may have only played a small role in promoting failure for the majority of the seismicity associated with the injection, but were more significant in the triggering of specific event clusters after shut-in occurred. Maghsoudi et al. (2018) used statistical techniques to show that spatiotemporal clusters of events, associated with an HF stimulation in the Horn River Basin (HRB), exhibited evidence for interevent triggering. They showed that the occurrence of events deviated from a purely Poisson process, and temporal clustering was present. In their study, while the microseismic event clusters exhibit this aftershock sequence-like pattern, the events deviate from the typical mainshock-aftershock behaviour observed in tectonic sequences. Large numbers of events also exhibit swarm-like behaviour, where the triggered events have very similar magnitudes to those that triggered them and, within a given cluster, the largest event is not typically the first.
In our study we use a microseismic data set acquired during an HF operation in the HRB. The HRB is one of the largest natural gas plays in British Columbia, Canada, and is part of the Western Canadian Sedimentary Basin (WCSB). Most of the wells in the HRB are horizontal, multiwell operations deployed from single pads, employing multistage HF to stimulate production from three Devonian shale layers: the Muskwa, Otter Park and Evie formations (Barker 2014). In recent years, multiple felt earthquakes associated with HF have been recorded in the WCSB (Schultz et al. 2015;Bao & Eaton 2016), with magnitudes reaching as high as M W 4.6 (Kao et al. 2018). Microseismic monitoring during operations has also repeatedly indicated the presence of smaller scale injection-induced fault activation (BC Oil and Gas Commission 2012).
We use the PSCMP code of Wang et al. (2006) to model the cumulative elastic stress changes during a period of HF-induced fault reactivation. We use the full moment tensor solutions, including the non-double couple fault motion. Sensitivity analyses are then conducted to constrain the effects of focal mechanism uncertainty on CFS through Monte Carlo methods, as developed in Catalli et al. (2013) and Meier et al. (2014). We examine the effect of the fluids on the resulting CFS by varying Skempton's ratio β. Due to the observed tight clustering of events, we also compare the elastic stress changes within clusters to the total cumulative CFS to investigate the significance of intracluster triggering.
The variability of induced seismicity between different sites and regions (e.g. Verdon et al. 2016) has been linked to differences in the background stress field (King et al. 2014;Göbel 2015). Fault orientation has been shown to be connected to the occurrence of induced seismicity, with the M W > 5 events associated with wastewater injection in Oklahoma occurring on planes well oriented to the background stresses (Keranen et al. 2013;Alt & Zoback 2017). Thus, the calculated elastic stress changes need to be considered in the context of the in situ stress field acting on the individual fracture or fault planes. The conditional probability of failure given the in situ stress field can be determined through a quantitative risk assessment (QRA) methodology developed in Chiaramonte et al. (2007) and Walsh & Zoback (2016). QRA methods have been applied to various geomechanical issues, such as wellbore stability (Moos et al. 2003) and CCS caprock integrity (Chae & Lee 2015). Following the approach of Walsh & Zoback (2016), the uncertainties are included in the calculation of the required pore pressure change to initiate failure (given the Mohr-Coulomb slip criterion). A 'probability' of failure for a given P can be determined by a Monte Carlo method of randomly sampling from distributions of the input parameters. This fault slip potential (FSP) can then be used as a proxy for the likelihood of failure for a given fluid pressure increase.
We conclude our study by interpreting the estimated elastic stress changes in the context of the estimated FSP for the activated faults.

DATA
We use a microseismic data set that was acquired during a multiwell, multistage hydraulic fracture treatment conducted in the HRB (described in Baird et al. 2017;Verdon & Budge 2018). 237 stages were completed in the 10 wells drilled, using a toe-heel, zipper-frac injection pattern. Continuous monitoring was provided by three downhole microseismic arrays. For the 119 stages for which we have data, around 92 700 events were recorded, processed and catalogued. These stages of injection constitute the 10 to 12 most proximate stages to the monitoring arrays for each well, and thus are expected to give the best quality data (shown in Fig. 1). Magnitudes (ranging from −2.4 < M W < 0.5), fault radii, stress drops, slips and locations were determined for all events in the catalogue, and around 35 per cent have full moment tensor solutions (MTs). The processing was conducted by a service provider, ESG Solutions. The data is proprietary and is not available for public release.
Event hypocentres were first calculated by inverting P-and Swave traveltimes through a layered anisotropic velocity model. This model was determined from sonic logs, and further calibrated using the perforation shots. Once preliminary locations were found, the velocity model was then linearly tapered at the transitions between the shale layers to more accurately reflect the smooth transitions observed in the sonic logs. Events were then relocated using this refined velocity model. Double difference corrections were then applied, refining the relative locations of the events. Seismic moment, source radius and slip were found for each event by fitting a source model to the event displacement spectra (see Stork et al. 2014). The source radii for the microseismic events used in this study are distributed around 15 m, with the majority (96 per cent) under 20 m. A small number (36 of the 923) of higher magnitude events have source radii up to 110 m. Moment tensors were determined by inversion of polarity data and P-to S-wave amplitudes for an estimated radiation pattern (as in Vavryčuk 2014Vavryčuk , 2015b. The fault plane and auxiliary planes were discriminated using the in situ stresses, finding the preferentially aligned nodal plane (NP) in the local stress field. This can be used as an estimate of the fault plane, over the auxiliary plane. However, the choice of nodal plane may not be significant for this elastic modelling, as discussed in Section 3.1. Throughout this work, we refer to a nominal 'primary', preferred NP chosen during the inversion, and a 'secondary' nodal plane.
Some of the determined MTs include a volumetric component. The vast majority (90 per cent) of events are effectively doublecouple (DC) sources (with tensile angles less than 3 • ), though some events, especially smaller events closer to the well, have larger non-DC components. Events with significant tensile angles (α > 10 • ), most likely associated with the opening of fractures though hydraulic stimulation, are generally far smaller in magnitude than those associated with the fault reactivation, with 84 per cent having magnitudes less than M W −0.5. Thus, while the non-DC component of slip is included in the analysis, it will have less of an effect on the overall pattern of stress change than the larger, DC sources.
No uncertainties were provided by the operating contractor. Thus, we estimate uncertainties throughout, informed by the methods through which the parameters were determined. As the location determination went through multiple stages of refinement, the absolute location uncertainty would be approximately ±50 m, while the relative error may be as low as ±10 m. We go on to estimate the location uncertainty further, by quantifying how diffuse each of the clusters are that map out the presumed fault planes (Fig. 4b). The combined effect of a low relative location error and the modelling method is discussed in Section 3.1, and more in Section 4. The uncertainty in the fault geometries, as inferred from the MTs, is the subject of the sensitivity analysis conducted later in this work (Section 3.3). Magnitudes, and thus source radii and slip, are considered to be well constrained, given the standard method by which they were calculated, and the low noise for these events detected on approximately 80 stations in deep borehole arrays.

Identifying fault reactivation
Initial examination of the microseismic event hypocentres indicated the re-activation of pre-existing faults, with clusters of events mapping out planar features that extend approximately 400 m below the wells, into the underlying Keg River limestone. The largest events were found in these clusters, with 27 events having magnitudes greater than M W = 0.0, whereas the overwhelming majority of events, those associated with hydraulic fracture propagation, have magnitudes around −2 < M W < −1. The largest events, and the most obvious downward growth of microseismicity, were seen in proximity to Wells 1 to 4, associated with the stimulation stages in the central portions of these wells (see Fig. 1). Because these events represent the clearest example of fault reactivation within this data set, we focus our analysis on the events associated with these wells and stages.
To further discriminate fault reactivation from the operationally induced microseismicity associated with fracturing, we examined variations in the b-value of the event magnitude-frequency distribution (Gutenberg & Richter 1944). For each stage, we clustered events based on their spatial positions using a density-based clustering algorithm DBSCAN (Ester et al. 1996). For each cluster we computed the b-value using the maximum likelihood approach (Aki 1965). To determine the minimum magnitude of completeness (M MIN ) for each cluster, we use a Kolmogorov-Smirnov test, with M MIN being the smallest magnitude at which the observed population can be fit with an exponential distribution with a confidence level of 95 per cent.
The resulting magnitude-frequency distributions are shown in Fig. 2. It is immediately apparent that the distribution of b-values is bimodal. The majority of event clusters have high b-values (b > 1.5), implying that they are associated with normal hydraulic fracture propagation, while a smaller number of clusters have b ≈ 1.0, a value similar to tectonic earthquake sequences, indicating that they are associated with fault reactivation (Verdon et al. 2013b;Verdon 2013;Eaton et al. 2014;Eaton & Maghsoudi 2015). We therefore use b < 1.5 as our criteria for identifying microseismic event clusters that represent fault reactivation.
These clusters occurred mostly during two periods of injection, as can be seen in Fig. 3. These are denoted by reactivation periods 1 and 2 (or, RP1 and RP2) and are separated by around 42 hours. Little seismicity occurred in the vicinity of these clusters except during these two periods. The events are shown in Fig. 3 with the same colour as the well with which they are associated. Each reactivation period is associated with an injection stage in Well 2, with fault reactivation clusters located beneath these injection stages. Seismicity continued in these clusters after the Well 2 stages had stopped, at which point injection switched to Well 4. However, the injection stages from Well 4 are over 500 m NE from Well 2, and there is clear lateral separation between the Well 4 events and the Well 2 fault reactivation, with the Well 4 events all appearing to show operationally induced HF events. We conclude that the fault reactivation events are associated solely with the Well 2 stages, and  not those of Well 4. Their apparent contemporaneity with Well 4 injection simply represents the fact that seismicity continued after injection had ceased in Well 2, as is commonly observed during injection-induced fault reactivation (Deichmann & Giardini 2009).
An additional cluster of low b-value seismicity occurred during a later injection stage in Well 2, to the NW of the clusters shown in Fig. 3, forming another cluster extending downwards from the well. This appears to connect to the diffuse cluster visible at around 2200 m depth in Fig. 3. Considering all of the events that were located in this area, two main structures become apparent: a planar feature, and a more diffuse cluster, both extending between 1900 and 2400 m depth. Fig. 4 shows the locations of the low b-value, fault related events, showing that these events map out two pre-existing fault structures both of which strike NE-SW, and dip towards each other.
The occurrence of events appears to be correlated tightly with the injection. We do not observe any clear aftershock-type sequences in either RP, where one large event triggers subsequent smaller ones in a spatiotemporal cluster. This type of behaviour would be indicative of a sequence with significant event-event triggering. Event clusters appear more swarm-like, where similar sized events are happening throughout the periods of activity. This swarm-like behaviour is contrasted with the results found by previous studies Maghsoudi et al. 2018), where there were clear patterns of triggered temporal clusters during injection-induced seismicity.
There is also no clear sign of the characteristic r ∝ √ t triggering wave front emanating from the point of injection, as exploited in many previous studies to estimate the diffusivity of the medium (Shapiro et al. 1997;Goertz-Allmann et al. 2011). This would be indicative of events triggered by the expanding wave front of pore pressure that results during injection. The elastic Coulomb modelling that follows will attempt to probe the ambiguous mechanism behind this sequence of events.
We use a least squares minimization algorithm to find the planes that best fit the seismicity. The RP1 plane has a strike φ of 250±10 • and a dip δ of 80±5 • . The RP2 plane has a φ of 60±15 • and a δ of 75±10 • . Uncertainties for these geometries are estimated given the spread of the events. These two suspected fault structures are denoted fault 1 and 2 (F1 and F2) respectively, named for the period of seismicity in which they were activated: events on F1 were seen mostly during RP1, whereas events of F2 occurred primarily during RP2, although some RP1 events also took place on F2. These planes can be seen projected onto their associated events in Fig. 4.
If events are thought to have occurred purely on these conjugate planes, the spread of events around the faults may provide an estimate for location uncertainty. This estimate would be ±25 m for F1, and ±46 m for F2, as shown by Fig. 4b. However, given the estimated rupture dimensions (∼20 m), it is uncertain as to whether the events observed are actually slip on pre-existing faults, or rather smaller fractures failing in a zone around these structures. Despite their size, all of the largest events (M W > 0) in the catalogue took place along these underlying planar features, we conclude that this region must be fundamentally different from that of a diffuse cloud of fractures, such as those near the wells. To compare these two interpretations, the following analysis will consider both the individual microseismic event geometries (Section 3.1), as well as the fault planes themselves (Section 3.2). This will also help to constrain the effect of the uncertainty in the geometries determined from the event moment tensors. The effect of Coulomb stress on the fault planes and the microseismic events should be similar if the event geometries reflect that of the larger structures.

Stress model for Horn River Basin
Assuming that the vertical stress is a principal stress, which is usually the case in sedimentary basins unless significant deformation or stress rotation is present, an in situ stress model is defined by the orientation and magnitudes of the principal stresses-the vertical stress S V , the maximum horizontal stress S Hmax , and the minimum horizontal stress S Hmin -as well as the pore pressure P, and the S Hmax orientation φ H . Borehole breakouts, density log measurements, leak-off tests, shut-in pressures, mud-weights, and the Downloaded from https://academic.oup.com/gji/article-abstract/217/2/1200/5315764 by University of Bristol Library user on 14 March 2019 Figure 3. (a) Hypocentre locations of the studied seismicity, grouped into 'reactivation periods' 1 and 2 (RP1 and RP2). Perforation shots are shown as squares, coloured in the same manner as their associated RP. Circles show nearby events considered not to be associated with the fault activation, and are coloured by their associated injection well. The map view of the seismicity is inset. (b) Temporal distribution of the seismicity. Injection durations are shown above, with colours denoting the well from which injection was occurring.
modelling of breakout rotations can be used to constrain estimates in these parameters (Haimson & Cornet 2003;Zoback 2010).
Previous studies have estimated the stress regime in the HRB for particular operations (Snelling et al. 2013;Sayers et al. 2016), as well as the frictional properties of the shale and underlying limestone (Chou et al. 2011;Hurd & Zoback 2012). There are a number of at-depth stress measurements in the region, which arise from the large number of shale gas operations in the HRB, as well as the adjacent Laird Basin and Montney shale play.
The stress model used in this study is calculated using data from the World Stress Map (Heidbach et al. 2016). At-depth stress measurements are taken from operator reports catalogued by Bell (2015), with each well located less than 30 km laterally from the site considered here. These measurements are shown in Fig. 5. These data provide a φ H of 55±10 • . The pore pressures given appear hydrostatic. This may be the case only for the formations from which the measurements were taken, or could be the product of a modelling assumption on behalf of the operator. Eaton & Schultz (2018) found that, for two other reservoirs in the WCSB, regions of overpressure were common, and correlated with an increased occurrence rate of earthquakes. We were not provided with drilling data for the field, and thus do not have a measure of the in situ pore pressure for the shales. We thus use hydrostatic pore pressure as a mean value, and give it an appropriate uncertainty to account for the possibility for overpressure. Uncertainties are also estimated for the other principal stresses (see Walsh & Zoback 2016). These will then be used in the Monte Carlo FSP analysis.
One factor not directly accounted for is the increase in pore pressure that accompanies HF. The magnitude of the P change, and its spatial extent, is highly dependent on the permeability structure of the reservoir Brown & Ge (2018a). More importantly, it is very sensitive the fracture permeability and extent. If matrix permeability is very low (∼10 nd; Chalmers et al. 2012;Dong et al. 2017), and fracture permeability is ignored, the P perturbation will be in far more than 10 MPa over a scale of tens of metres (Brown & Ge 2018b). Fractures act to increase bulk permeability, increasing the distance over which the fluid pressure is applied, and reducing the magnitude of the resulting P. Constraining these exact parameters requires complex hydromechanical modelling, and is beyond the scope of this work.
A least squares fit is applied to find the gradients of S Hmax , S Hmin , S V , and P, and their respective uncertainties (York et al. 2004). The stresses acting at the depth of the reservoir are then found by extrapolation. This gave an S Hmax of 77±12 MPa, S Hmin of 51±6 MPa, S V of 66±5 MPa, and P of 27±7 MPa. These values broadly agree with those found by Chou et al. (2011), who constructed a stress model for the HRB shales using borehole observations provided by industry, and Sayers et al. (2016), who used anisotropic seismic attributes to estimate stress conditions. In the Anderson classification scheme this is a strike-slip regime, with S Hmax > S V > S Hmin . Within the uncertainties, the regime could be considered normal faulting, however, that is taken into account in the later Monte-Carlo FSP analysis.

Elastic stress modelling
Using the estimated slips, rupture lengths, and fault plane orientations, the deformation resulting from each event can be modelled. The code PSCMP by Wang et al. (2006) is used to calculate the displacement field u due to slip on a square patch. This uses the analytical Okada solution of the Green's functions for a homogeneous elastic half-space (Okada 1992). The model then derives the stain tensor ij from u, which is then related to the stress sensor σ ij using Hooke's law for a linear elastic medium. This stress tensor is then resolved onto a receiver plane and rake, and the resulting Coulomb stress change CFS is calculated (eq. 1). All events during RP1 and RP2 are defined as 'receivers', and each event that occurred prior to each receiver is treated as 'source' of deformation. This gave a total of 1119 receivers (of which 923 have MTs), and potentially 6602 sources (of which 1876 have MTs). Proceeding iteratively through each receiver with a MT, the cumulative elastic stress change is calculated at that receiver's hypocentre.
The elastic parameters (μ = 0.7 and β = 0.4) used are based on measurements made by (Chou et al. 2011) in the HRB, as well as values obtained for formations of similar type, depth, pressures and porosities (Kohli & Zoback 2013). Sources are modelled as square patches, with dimensions calculated to match the rupture area given by the source spectra. The average patch length is approximately 25 m, with the largest around 190 m. Slip is resolved into the strike and dip components of the plane using the rake and tensile angles (see Vavryčuk 2011). Sources are excluded with sourcereceiver distances less than one source length because CFS values computed using a uniform slip model become unreliable in the very near-field (see Meier et al. 2014). This step will have a significant effect on the resulting CFS magnitudes and polarities, but will also aid in removing artefacts which result from uncertainties in the slip behaviour and source mechanism (Steacy et al. 2004), and stress drop (Schoenball et al. 2012). The ramifications of this are discussed in detail in Section 4.1. Fig. 6 shows resulting CFS values found for RP1 and RP2 when the preferred nodal plane, as identified during the MT inversion, is used. The Coulomb index for all of the events in RP1 is just below However, there is no correlation between the event magnitude and CFS magnitude. For both reactivation periods, no CFS received from the microseismicity exceeds ±1 MPa. There are a number of events (36 per cent) however that received positive CFS greater than the commonly used triggering threshold of 10 kPa (see King et al. 1994;Stein 1999).
Previous work in this field (e.g. Catalli et al. 2013;Pennington & Chen 2017), observe CIs generally in excess of 70 per cent. Given the above results, elastic deformation may have only modestly promoted slip in this case, during the initial period of activity. However, during the latter period, the events were generally inhibited by the preceding elastic deformation. This implies that fluid processes such as aseismic tensile slip, pore pressure changes and poroelastic effects were most likely the dominant cause of fault reactivation.
These results are only slightly altered when the secondary nodal plane is chosen, with a CI of 61 per cent positive for the first period and 36 per cent for the second. The windowed CI values also show very similar temporal behaviour as for the preferred NP, with peaks around 75 per cent early in RP1, and CI not exceeding 45 per cent in RP2. The similarity should not be surprising however due to the inherent symmetry in the stress change pattern when using a uniform slip model (Meier et al. 2014). This results in generally similar stress change values. The NP ambiguity is further explored in Section 3.3.

Resolving stresses onto fault structures
In addition to computing the Coulomb stress changes for each event, we also evaluate the stress interactions between the microseismic events and the two fault planes (see Section 2.1 for description), as in the latter section of Schoenball et al. (2012). The same geomechanical parameters as in Section 3.1 are used, along with the preferred NP geometries. As events on F1 only occurred during the first reactivation period, only stress changes up to the end of RP1 are calculated for this plane. Events on F2 occurred during both reactivation periods, so stress changes are calculated for this plane for the entire duration of both reactivation periods. In the model, the receiver geometries were those of the individual fault planes, and receiver rake was given by the average rake of events used to define the two fault planes (see Section 2.1). Cumulative CFS was calculated in half-hour increments, accumulating sources through successive time windows. The resulting stress changes acting at a given time on the plane are then calculated at sample points gridded at 20 m intervals along the surface. This grid was defined by examining the extent of the cloud of seismicity. F1 extended for 180 m along its strike and 500 m in the dip direction. F2 extended 270 m along its strike, and 550 m in the dip. The fault dimensions can be seen in Fig. 4. These gridded stresses can then be analysed in a similar manner to that of an event population, by examining the proportion of positive CFS received (Coulomb index), or the sum, or average of the stress changes at the sample points.
The results of this analysis show a similar pattern to the elastic stress changes as calculated in Section 3.1. In RP1, both the sum and the average stresses resolved across the plane are positive, and the effective CI (the proportion of sample points on the plane with CFS > 0) reaches as high as 78 per cent around the middle of RP1 (Fig. 7). This suggests that the elastic stress changes modestly promoted failure on F1. Similar temporal variations in CI are observed as for Section 3.1 for the individual events. However, there is a short period during RP1 where the opposite signal is seen. At around 13 : 00, when stresses are resolved across F1, the effective CI acting is at its lowest value: ∼35 per cent. But when the individual events are considered, this early part of RP1 however had the relatively high CI of ∼70 per cent.
The similarities between the two methods of stress change determinations continue for F2 through RP1 and RP2. The CI across the second plane never exceeded 50 per cent and rarely went above 30 per cent, indicating that across the surface of this plane, the majority of the elastic stress changes were inhibiting failure throughout both reactivation periods. The average and total stress change magnitudes acting on the two fault planes mirror these trends, with consistently positive CFS for F1 and negative for F2.
Thus, we reach the same conclusion -deformation may have modestly contributed to the failure on F1, but fluid and/or pore pressure effects must have dominated for F2, counteracting the effect of inhibiting elastic CFS. This also highlights that, even when the uncertainty of resolving stresses onto the individual receiver geometries is removed, roughly the same behaviour holds. However, geometry uncertainty is obviously not removed entirely here, as the sources of deformation are still the catalogued events.

Sensitivity analysis
We use a Monte Carlo approach to estimate the sensitivity of CFS to focal mechanism uncertainties, as in Catalli et al. (2013) and Meier et al. (2014). Using the first 500 events with MTs that occurred during the first reactivation period, 3000 test catalogues are produced. Each catalogue contains the same 500 events, but with the geometries permuted in three different fashions, giving 1000 catalogues for each type. In the first set of catalogues, the NP is randomly selected from the original FM in the catalogue. In the second case, Von Mises distributed rotations are applied to the preferred NP, representing uncertainties in the focal mechanism inversion. As no uncertainties were provided from the processing, a single FM uncertainty of σ = 30 • is given, representing a reasonable, perhaps pessimistic, value for this form of MT inversion using borehole-acquired microseismic data (Vavryčuk 2014). A final set of test catalogues are generated with entirely random strikes, dips and rakes (Kagan 2005). This allows us to compare the perturbed FMs stress changes to an extreme uncertainty scenario, as well as a potential null hypothesis case, where there is no preferential alignment of failure planes. The CFS modelling is then conducted in the same manner as in Section 3.1 for each catalogue. The CI of each of the catalogues is then calculated, and the distributions in CI can be compared for each of the three FM perturbation methods.
The results of the sensitivity analysis are shown in Fig. 8. For this subset of events, the CI when using the preferred nodal plane is 64 per cent, and when using the secondary NP, are 63 per cent. When the nodal plane are randomly selected for each event, the CI values  are approximately normally distributed around a mean of 63±1 per cent. This indicates that the choice of NP does not have a significant impact on the resulting Coulomb index when a uniform slip model is used. When the FMs are rotated by Von Mises distributed angles, the mean CI is decreased to 59±3 per cent. Simulating entirely random event geometries gave CI values distributed around 50±3 per cent, as expected for entirely randomly aligned fracture sets.
As this perturbation in fault geometry is applied, the CI distribution shifts toward that of the random focal mechanism scenario. A FM uncertainty of σ = 30 • would appear to correspond to an approximately ±4 per cent uncertainty in the observed CI values. However, even when the focal mechanisms are rotated by around 30 • , the positive signal observed during RP1 is still present. As shown in Catalli et al. (2013), any perturbation to the FMs appears to systematically decrease the resulting average CI. Meier et al. (2014) interpreted this result as the outcome of applying further randomness to an already uncertain, noisy FM catalogue. We see the same result here. This challenge of quantifying the effect of focal mechanism uncertainty, and the significance of the signals observed, is further discussed in Section 4.
With injection fluids interacting with the proposed faults in this system, the effect of changing pore pressure cannot be neglected. Thus, further sensitivity analysis is conducted to examine the effect of varying Skempton's ratio β. Using the same subset of events as in the above FM analysis, stress changes are repeatedly modelled with a varying β. Values of μ from 0 and 0.7 are used, in 0.1 increments, which corresponds to Skempton's ratio varying from 0 ≤ β ≤ 1.
Varying Skempton's ratio appears to have a similar scale of effect on the CI as applying rotations to the event geometries (Fig. 9). CI changes from 64 per cent to 56 per cent, as β varies from 0 to 1, simulating entirely unsaturated to saturated faults. Note that when β > 0.4, CI decreases at a more rapid rate than when 0 < β < 0.4. When looking at the values of the individual event CFS, changing β generally acted to change the magnitude, but infrequently changed the sign -only 14 per cent of receivers changed sign across the whole range of β. This minority are responsible for the change in CI as measured for the population. The sign of all receivers in this sample with | CFS| > 0.05 MPa remained constant from one extreme β to the other. This indicates that if injection fluids are interacting with the faults near the injection points (increasing β from 0.4) the value of the CI may decrease up to approximately 4 per cent. This would weaken the positive signal observed during the first reactivation period.

Cluster stress budget
To investigate the effect of tight clustering on the elastic stress change budget, the DBSCAN algorithm (Ester et al. 1996) is used again to identify clusters within the event population. Whereas in Section 3.1 we used only spatial dimensions to initially identify fault reactivation, we now incorporate time as an additional dimension for the clustering analysis. A scaling factor of c = 250 ms −1 is applied to the time dimension. This value is found to best separate the clusters in the time domain. Following the method of Daszykowski et al. (2001), an optimal neighbourhood distance is calculated using the sorted k-distance approach, giving a value of 35 m. The resulting clusters are shown in Fig. 10.
We calculate the stress changes in the same manner as Section 3.1, however, for each receiver, only sources from within the same cluster are used. This will allow for the full cumulative stress change to be compared to the change when only the deformation from within each cluster is considered. These resulting stress change value is termed the 'intracluster' stress change. Fig. 11 shows that, for the majority of events (89 per cent), the intracluster CFS makes up effectively all of the stress change received. This indicates the importance of tight clustering for microseismic stress transfer, as even when sources are excluded within one source length, the stress change from within the cluster contributes the majority of the CFS. The relationship in Fig. 11 also shows that when CFS is significant (i.e. greater than ±0.1 MPa), effectively all of the stress change results from within the cluster.
One would expect this to be true for clusters 1 and 2, where few events outside of the cluster preceded the seismicity, but this linear relationship is also present for cluster 5 (Fig. 12). This cluster will have received all of the CFS from all of those proceeding it, but still, the majority of the stress changes comes from within that cluster. This agrees with Vasudevan & Eaton (2011), in that for microseismic events, where the spacial propagation of stress is not laterally extensive, the majority of the stress change comes from nearby sources.

Fault slip potential
The code FaultSlipPotential (Walsh et al. 2017) is used to conduct a probabilistic geomechanical analysis, using the methodology described in Walsh & Zoback (2016). The stress gradients, S Hmax azimuth, fault geometries, and the friction coefficient are used as inputs, as well as their associated uncertainties. The principal stress uncertainties used are those given in Section 2.2, along with a μ of 0.7±0.1. The uncertainty in the orientation of the faults is given in Section 2.1. These input parameters are randomly sampled 1000 times within their uncertainties, and a resulting critical pore pressure required to initiate failure, P C , is then calculated using Coulomb's law: For the range of P C perturbations, a cumulative probability of exceeding the Mohr-Coulomb failure criterion is then found. Fault-SlipPotential uses uniform distributions from which to sample the input parameters. This will give a narrower distribution in the resulting P C values, than would be found when using standard errors. We use eq. (3) to examine the in situ stress field with respect to the fault orientations (shown in Fig. 4), as delineated by the microseismicity. Fig. 13 shows the computed σ n , τ , and P C as a function of fault strike and dip. Optimally oriented faults, where P C is at its minimum value of 14.5 MPa, have strikes of 28 • , 81 • , 208 • and 261 • . These are relatively close to the strikes of the planar structures: 250 • and 60 • . Fig. 14, shows the result of the probabilistic analysis using Fault-SlipPotential, showing the cumulative density function for the probability of slip given some pore pressure perturbation P for F1 and F2. The probability of failure does not exceed 50 per cent until the pore pressure change is in excess of 15 MPa for both of the fault planes. F1 is marginally more optimally aligned for failure, with F2 requiring an additional ∼1 MPa of pore pressure change to give the same probability of failure. This scale of pore pressure change (∼10 MPa) is similar to that found by Chiaramonte et al. (2007), for a trap bounding fault in north-central North America. Across the range of input values, the S Hmin and in situ pore pressure uncertainties had the largest effect on the P to slip for both faults. Significantly reducing S Hmin will have the effect of shifting the Mohr circle towards the failure criterion. Increasing the in situ pore pressure will act to decrease the effective stresses acting to clamp the fault.
A P required to fail of approximately 15 MPa would indicate that the two hypothesized faults here are not critically stressed within their stress field. A large pore pressure increase would be required for the Mohr-Coulomb failure criterion to be reached. With the average observed elastic stress increase on the fault planes being approximately 0.1 MPa, these planes would require two orders of magnitude larger pore pressure changes to initiate failure. Even at the extremes of the stress model and fault parameter uncertainty, the required stress change for failure is an order of magnitude Downloaded from https://academic.oup.com/gji/article-abstract/217/2/1200/5315764 by University of Bristol Library user on 14 March 2019  above the average elastic CFS. This suggests that the interaction with high pressure injection fluids was likely to have been responsible for the activity in both reactivation periods, and elastic stress transfer did not play a major role in the fault activation. As discussed in Section 2.2, this analysis does not take into account the potential for overpressure in the formation, other than in the uncertainty applied to the in situ P (±7 MPa). During injection downhole pressures were measured to be consistently greater than 50 MPa, giving a differential pressure between that and the in situ pore pressure of over 20 MPa. Given the top of F1 is ∼50 m from the nearest perforation, and the fact that fractures are likely to extend at least that distance, it is plausible that the pressurized fluid interacted with the pre-existing structure and exceeded this failure criterion.
This analysis however uses a simple model of fault stability, and does not take into account the possibility of frictional asperities on the fault surface. These areas of weakness would allow the failure threshold to be reached on a small fraction of the fault surface and enable failure to propagate outwards. This type of simple modelling also does not explain the spatiotemporal clustering observed heresmall events near the well start the sequence, but immediately after and 400 m below, larger magnitude events occur at the base of the faults, with very little intermediate seismicity. This would also suggest a rapid transport of fluid along the faults, and then subsequent failure back up along the plane of weakness. This type of behaviour has been observed in EGS-induced microseismic sequences (Majer et al. 2007), and is described as a potential factor in the maximum observed magnitude of the injection induced seismicity ). More complex fluid and fracture modelling is required to fully examine the pressure changes occurring within the fault system, however this slip potential model does indicate that these faults are not optimally oriented for failure in the regional stress state.

Elastic stress modelling
The positive signal we observe in the first period of activity is a CI of around 60 per cent, as well as small periods wherein CI reached around 70 per cent. Ascribing significance to this signal, a slightly higher proportion of events receiving positive CFS, is difficult, especially given the uncertainties in the FMs and the assumptions inherent in modelling elastic stress change. We see similar   magnitudes of stress changes as previous studies (e.g. Schoenball et al. 2012), however our results diverge from studies which find CIs in excess of 75 per cent, implying that elastic stress transfer played a more significant role during injection induced seismicity (e.g. Catalli et al. 2013;Pennington & Chen 2017). We have also tested the reliability of this signal by conducting a separate analysis, wherein individual receiver geometry is bypassed, and stress changes are instead mapped onto the visible fault planes. Whilst the measured quantities are different, the same pattern emergesa positive signal in the first reactivation period, and a generally negative signal in the second. Testing the sensitivity to particular uncertainties conventionally, using the method employed in Catalli et al. (2013), we observe a decrease in CI by both the focal mechanism (∼4 per cent) and frictional uncertainty (∼4 per cent), which weaken the positive signal observed. However, this approach may not be a reflection of the effect of FM uncertainty on the resulting CI. As described in Meier et al. (2014), adding random noise to the FMs by applying a perturbation just adds a smearing effect, always resulting in a CI closer to 50 per cent, no matter how significant the original signal. This method however, does show that even when FMs are rotated by, on average, a very significant amount (30 • ), the signal is still present. This indicates that despite FM uncertainty, the events in the first reactivation period are spatially distributed in such a way that elastic stress changes still modestly promote failure.
The inherent assumptions used and lack of data must also be considered, and their effect on the accuracy of our results. Both observational precision, and lack of data, constrain the accuracy of our modelling. Uncertainties in both the FMs and locations were not provided by the processing company, and MTs were not computed for a number of events in the original catalogues in this study due to low signal-to-noise ratios. That said, these 'missing' events generally have considerably lower magnitudes (M W < −1.5) than the majority of events used here (those with MTs), so the effect on the stress change from these sources is expected to be small. Their contribution to the stress changes will therefore be minimal and is assumed to not have a significant effect on the results. The uniform slip model used, which gives rise to a number of the limitations described above (NP ambiguity, source exclusion within 1 source length), has a dramatic effect on the accuracy of our results, as it does for all Coulomb stress modelling. A fuller treatment of modelling slip distributions, however, would have been impractical in this case, with up to ∼1400 sources being modelled for later receivers. Also, the nature of slip distribution on the fault plane is highly uncertain in itself. Ascribing some assumed rupture behaviour in the hope of more accurately reflecting the near-field stress change could in fact produce significantly more uncertain CFS values within one source length.
Thus, we followed the methodology of previous works (such as Meier et al. 2014), where slip is treated as uniform, and near-field sources are ignored. This method obviously results in an lower estimate for the resulting CFS magnitude, but also removes the significant effect of the uncertainty in the stress drop, slip model, and to a certain extent source location (Steacy et al. 2004;Schoenball et al. 2012). These parameters, each of which have significant uncertainties, mostly effect the near-field stress change, due to the r −3 power law describing the decay of CFS. The resulting values can thus be considered a more accurate reflection of the triggering effect of elastic stress change. One could imagine that this exclusion step is amplified in this study in particular, where events are relatively tightly clustered. However, as shown in Fig. 12, the CFS from within clusters still makes up a significant proportion of the stress change received, and thus effect from sources within a cluster is not negated using this method.

Focal mechanism sensitivity analysis
Comparing our permuted focal mechanism catalogues to those with entirely random FMs (Fig. 8) also may not provide an accurate measure of significance when considering a possible null hypothesis, as described by Meier et al. (2014). A more suitable method could involve comparing to the CFS computed for events not thought to be triggered by elastic stress changes, and thus not associated with the fault activation. This could include populations near to the wells, assumed to be caused by tensile fracturing. For these events, fluid effects are expected to dominate and no considerable elastic triggering, that is, the CI for this population would be around 50 per cent. However, the accuracy of the focal mechanisms for these events would be even more questionable, with event size generally around one to two magnitudes smaller than the fault reactivation events. We already have considered a relatively pessimistic uncertainty for the FMs (30 • ), however, this would most likely increase dramatically for those events where the signal-to-noise ratio is even smaller. Another consideration is the alignment of failure planes for tensile fracturing events. These events will be generally well aligned parallel to the maximum stress direction, and so their relative alignments will deviate from that a zone of random oriented fractures. Thus, due to the spatial pattern of the stress change, they will predominantly receive positive stress change, and static stress transfer will appear to play a significant role.

Fault orientations
The geometry of the mapped faults and an estimate of the in situ stress regime indicate that the activity seen here was initiated by the high-pressure injection near the failure plane, as faults require in excess of 10 MPa of pore pressure change to reach the Mohr failure criterion. This is in contrast to the stress environment and fault system found in central the US, where waste-water injection has stimulated high magnitude (>3 M W ) events regularly in recent years (Alt & Zoback 2017). Both Walsh & Zoback (2016) and Schoenball et al. (2018) found, using the QRA methods employed here, that multiple fault systems mapped in north-central Oklahoma are preferentially aligned for failure, with some faults requiring only ∼2 MPa of P to reach a failure potential of 0.5. In our case, the required pore pressure changes are approximately six times greater.
This could provide a possible explanation as to the difference in significance of elastic CFS seen in this study compared to other sites (e.g. Catalli et al. 2013;Sumy et al. 2014;Pennington & Chen 2017). If the faults are already close to being critically stressed, then the small amounts of stress change associated with static elastic stress transfer can play a significant role. If the faults are not close to failure, then larger stress changes are required. In such cases the stress changes associated with static elastic stress transfer are not large enough, and instead the very high pressures associated with fluids dominate the activation of faults, meaning that static stress transfer does not play a significant role.

C O N C L U S I O N S
For this case of injection induced fault activation, elastic stress transfer modelling suggests that during the first period of seismicity, deformation on faults and fractures may have weakly promoted further failure. However during the second period, these stress changes inhibited failure. This result is insensitive to the choice of event geometry (from the two nodal plane solutions of the MT), due to the symmetrical nature of the stress change pattern that results from a uniform slip model. Applying Von Mises distributed uncertainties to the focal mechanisms of the events appears to diminish, but does not completely eliminate, the observed CI signals. However, this is due to the nature of measuring the FM sensitivity on a population of events with already uncertain FMs, which acts to consistently shift the CI towards 50 per cent. Variability in the Skempton's ratio is also found to diminish, but not completely eliminate, the positive CI signal. The scale to which the resulting CI distributions shift while varying the model parameters highlights the importance of properly incorporating the various uncertainties that affect CFS calculations.
Intracluster stress transfer is found to be significant in this case, meaning that the elastic stress change budget is chiefly derived from events in the same cluster. This is expected for tightly clustered microseismic events because of the small area and (size of) slip on the fault plane, resulting in spatially contained elastic stresses. Thus, events only receive significant CFS when they are within ∼100 m from the source.
The weak contribution of elastic stress transfer to the failure can be interpreted with respect to the in situ stresses acting on the hypothesized fault structures. Using at-depth measurements of the regional principal stresses, the probability of slip for a given pore pressure change on the two planar structures is examined using a quantitative risk assessment technique. We find that for both faults, the pore pressure changes required to reach the failure criterion are around 15 MPa, far in excess of the reliable values of elastic stress change. This supports that, in this case, the failure was stimulated by the presence of highly pressurized fluids in the system. These faults are not critically stressed, unlike those found in the central US, which have been found to require only a few MPa of pore pressure to reach the slip criterion (Walsh & Zoback 2016;Schoenball et al. 2018). This could provide an explanation for the difference between the CFS behaviour found here and studies which concluded that this mechanism played a more significant role (Catalli et al. 2013;Pennington & Chen 2017).

A C K N O W L E D G E M E N T S
The authors would like to thank the operator for supplying the data, and ESG Solutions Ltd., the service provider who collected and processed the data set. TK was supported by the Natural Environment Research Council (NERC) GW4+ Doctoral Training Partnership (grant number NE/L002434/1). MJW was supported by NERC (grant number NE/R017956/1) and the Southern California Earthquake Center (SCEC) (Contribution No. 8961). SCEC is funded by NSF Cooperative Agreement EAR-1600087 and USGS Cooperative Agreement G17AC00047. JPV was supported by NERC (grant number NE/R018162/1). JMK was supported by NERC (grant number NE/R018006/1). This project was supported in part by Bristol University Microseismicity Projects (BUMPs). This work used data from the World Stress Map Database Release 2016 of Heidbach et al. (2016) . The code MTDecompostion of Vavryčuk (2015a) and MATLAB scripts produced by Meier et al. (2014) were both modified for use in this work. The DBSCAN algorithm of Ester et al. (1996) was used to determine spatiotemporal clusters of events. The code FaultSlipPotentialv1.07 of Walsh et al. (2017) was used in the latter sections. We would like to thank the reviewers for their very constructive feedback, which improved the quality of the manuscript.