Stress Transfer From Opening Hydraulic Fractures Controls the Distribution of Induced Seismicity

Understanding the dominant physical processes that cause fault reactivation due to fluid injection is vital to develop strategies to avoid and mitigate injection‐induced seismicity. Injection‐induced seismicity is a risk for several industries, including hydraulic fracturing, geothermal stimulation, oilfield waste disposal and carbon capture and storage, with hydraulic fracturing having been associated with some of the highest magnitude induced earthquakes ( M>5 ). As such, strict regulatory schemes have been implemented globally to limit the felt seismicity associated with operations. In the UK, a very strict “traffic light” system is currently in place. These procedures were employed several times during injection at the PNR‐1z well at Preston New Road, Lancashire, UK, from October to December 2018. As injection proceeded, it became apparent to the operator that stages were interacting with a seismogenic planar structure, interpreted as a fault zone, with several ML>0.5 events occurring. Microseismicity was clustered along this planar structure in a fashion that could not readily be explained through pore pressure diffusion or hydraulic fracture growth. Instead, we investigate the role of static elastic stress transfer created by the tensile opening of hydraulic fractures. We find that the spatial distributions of microseismicity are strongly correlated with areas that receive positive Mohr‐Coulomb stress changes from the tensile fracture opening, while areas that receive negative Mohr‐Coulomb stress change are quiescent. We conclude that the stressing due to tensile hydraulic fracture opening plays a significant role in controlling the spatiotemporal distribution of induced seismicity.


Introduction
Felt or damaging earthquakes have been induced or triggered by subsurface fluid injection related to a number of industrial activities. These include enhanced geothermal systems at Basel (Deichmann and Giardini, 2009) and Pohang (Grigoli et al., 2018;Kim et al., 2018), wastewater injection in the central United States (Keranen et al., 2013;Walsh and Zoback, 2015), carbon capture and storage at In Sala, Algeria (Stork et al., 2015), and hydraulic fracturing in central and western Canada (Bao and Eaton, 2016;Atkinson et al., 2016;Kao et al., 2018), the central United States (Holland, 2013;Skoumal et al., 2018), and the Sichuan Basin, China (Lei et al., 2017;Lei et al., 2019;Meng et al., 2019). However, while the links between fluid injection and seismicity are clear, the underlying physical processes by which injection causes fault reactivation are not yet well established. This matters because developing this understanding is crucial if we are to develop methods to prevent or mitigate injection-induced seismicity. In a broad sense, the mechanism of most injection-induced seismicity is well established: Fluid injection leads to an increase in pore pressure, decreasing the normal stress acting on critically stressed faults, and bringing them closer to failure (Raleigh et al., 1976). On large spatial scales in relatively permeable formations (as in the case of wastewater injection), pore pressure increases transmitted over large distances by diffusion would appear to be the dominant activation mechanism (Goebel et al., 2017;Goebel and Brodsky, 2018). In low permeability reservoirs and on smaller scales (on the order of hundreds of meters, within hours of injection), other mechanisms can dominate: the poroelastic expansion of the rock frame; direct pressure from the injected fluids; elastic stress changes from seismic events or fracture opening; and aseismic creep Bhattacharya and Viesca, 2019;Eyre et al., 2019).
Elastic stress change models have been used for decades to determine the triggering mechanism of tectonic earthquakes (Stein, 1999;Harris, 1998;Steacy et al., 2005;Meier et al., 2014;Wedmore et al., 2017), illuminating the sometimes unexpected spatiotemporal patterns which occur during seismic sequences. These models are regularly applied in physics-based earthquake hazard forecasts, using the observed slip on faults 10.1029/2019JB018794 to model the spatial distribution of subsequent, potentially damaging, earthquakes (Cattania et al., 2018;Mancini et al., 2019). Elastostatic modeling has also been applied with tensile sources, such as the analysis by Green et al. (2015) of a seismic sequence associated with dyke intrusion in Iceland. The areas receiving positive elastic Coulomb stress changes that resulted from the opening of the dyke were well correlated with the locations of seismic events throughout the sequence. As the sequence progressed and the dyke's orientation changed, earthquake rates were suppressed in areas experiencing negative Coulomb stress changes. In hydraulic fracturing, the tensile opening of hydraulic fractures produces perturbations to the stress state in a similar manner. Spatiotemporal observations in microseismicity that would be difficult to explain through any other mechanism could also be explained through the elastic stress changes that result from the tensile opening of fractures.
Such observations were made during hydraulic fracturing at the Preston New Road PNR-1z shale gas well in Lancashire, UK, in 2018 (described in Clarke, Verdon, et al., 2019). This was the first onshore well in the United Kingdom to be stimulated since a government review of this technique (The Royal Society, 2012). It was therefore the subject of extensive scrutiny by the public and by national media and was extensively monitored both by the operator and by independently funded organizations (Clarke, Soroush, & Wood, 2019).
Hydraulic fracturing at PNR-1z was subject to a traffic light scheme (TLS). This is a procedure developed to avoid felt seismicity (M L > 1.5) by taking mitigating actions (e.g., reducing injection rates, pausing injection, or skipping injection stages) when induced events of particular threshold magnitudes are observed. The "red-light" threshold in the United Kingdom is set at M L = 0.5, exceedance of which requires an 18-hr pause in operations. Microseismicity during injection at PNR-1z exceeded this limit on several occasions. During operations, the operator used a statistical model to forecast and manage induced seismicity (Clarke, Soroush, & Wood, 2019). One felt event did occur, with M L = 1.5 on 11 December 2018. Interestingly, the observed spatiotemporal distribution of microseismicity is not easily explained by the growth of hydraulic fractures or a diffusive pore pressure increase. Thus, in this study we examine the elastic stress changes in the vicinity of the well that occurred during the opening of hydraulic fractures and the potential impact these stress changes could have on the observed microseismicity. This is distinct from a poroelastic model, which would calculate the change to the stress state that results from increasing pore pressure deforming the rock mass itself, a continuously distributed inflation of the matrix due to increased pore fluid pressure. Here, we look at the propagation of elastic stress from discrete opening of finite model fractures.
Slip on faults, and tensile opening of fractures, will generate elastic stress changes in the surrounding rock. These changes can be resolved into changes in the normal stress n (defined here as positive extensive) and shear stress acting on nearby structures and combined to compute the Coulomb failure stress change ΔCFS: where ′ is the effective coefficient of friction.
Modeling of ΔCFS is a simple and effective tool for examining the effects of stress on surrounding faults or fractures-a positive value indicates that stress has changed in such a way as to promote failure, while a negative value means the stress change acts to inhibit failure. However, it is difficult to robustly model and interpret elastic stress changes. Defining a significance threshold for the effect on a population of events (Meier et al., 2014), quantifying model uncertainties (Catalli et al., 2013;Kettlety et al., 2019), and untangling the effects of other failure mechanisms, such as dynamic triggering or poroelasticity, all provide a significant challenge. Nonetheless, elastostatic stress modeling has repeatedly provided a robust explanation for the spatial distribution of earthquake sequences (Steacy et al., 2005;Meier et al., 2014;Wedmore et al., 2017;Cattania et al., 2018), and when applied carefully, can be an effective method of studying the triggering of induced seismicity (Schoenball et al., 2012;Catalli et al., 2013;Sumy et al., 2014;Pennington and Chen, 2017;Kettlety et al., 2019).
In this study, we examine the stress changes that result from the tensile opening of hydraulic fractures, modeled as displacement on finite patches within an elastic medium, and their effect on the distribution of microseismicity observed during the Preston New Road PNR-1z hydraulic fracturing operation in 2018 in the United Kingdom. We develop a stochastic, Monte Carlo procedure for generating model fractures as a set 10.1029/2019JB018794 of pure tensile opening discrete patches and calculate the resulting cumulative elastic stress changes from each fracturing stage. We compare the spatial patterns in ΔCFS with respect to the spatiotemporal evolution of the microseismicity. We show the areas of positive ΔCFS from prior and current stages correlate well with the hypocenters of the observed microseismicity and that areas where seismicity was unexpectedly quiescent received predominantly negative ΔCFS, suggesting areas are being clamped by the opening of fractures.

Hydraulic Fracturing at Preston New Road, UK
In October 2018, Cuadrilla Resources Ltd. began hydraulic fracturing operations at the Preston New Road PNR-1z well in Lancashire, UK. The operation targeted the upper section of the Bowland shale, a 1.2-km-thick Carboniferous natural gas-bearing formation (Andrews, 2013;Clarke et al., 2018). Hydraulic fracturing was monitored by a microseismic array of 24 three-component geophones housed in the adjacent well (PNR-2) (Clarke, Soroush, & Wood, 2019), shown in Figure 2. This was combined with a surface array, composed of the local UKArray (Baptie, 2018) broadband stations operated by the British Geological Survey (BGS), supplemented by a mix of eight broadband and three-component short-period instruments deployed by the operator as part of the monitoring program. The monitoring array, both surface and downhole, is detailed in Clarke, Soroush, and Wood (2019).
Over the course of 3 months, 17 stages were stimulated, with a planned injection program of 400 m 3 of slickwater fluid and 50 t of proppant per stage. Strict seismicity constraints-the TLS that is currently in place in the United Kingdom (Green et al., 2012)-restricted operations during many of the worked stages, with any event detected during pumping above M L 0.5 requiring a pause in injection for a minimum of 18 hr. More than 38,000 microseismic events were detected, with magnitudes ranging from −3.1 to 1.6 ( Figure 1). Data were processed in real-time by Schlumberger Ltd., providing event locations, M W magnitudes, and estimated source parameters. Estimates of location errors are around 10 to 50 m, typical of downhole microseismic monitoring. Focal mechanisms were independently calculated by both Schlumberger Ltd. and the BGS for 41 of the highest magnitude events using the surface station polarity data. These are also shown in Figure 1.
As successive stages were injected, it became apparent that the operations were interacting with preexisting seismogenic structures (Clarke, Soroush, & Wood, 2019). Seismicity was repeatedly occurring with magnitudes approaching or exceeding the red-light threshold. This resulted in the operator skipping stages, moving further toward the heel of the well to avoid repeatedly activating these features. In late October 2018, roughly 2 weeks after the start of operations, six events occurred that exceeded the TLS thresholds. After this, operations were paused for approximately 1 month, during which low levels of microseismicity continued to occur. The highest magnitude events, as well as the events during this hiatus, were predominantly located around a particular structure, a subvertical planar feature, striking to the NE of the injection well ( Figure 1). As detailed in Clarke, Soroush, and Wood (2019), we take a sample of events to calculate the orientation of this feature: the largest (M W > 0) events that took place after it was first encountered (from Stage 18); and all events that continued to occur in this zone during the month hiatus in operations. It was during this time that it became very clear that a more seismogenic planar feature was present, as the areas around each of the worked stages became quiescent except in vicinity of this feature. A least squares planar fit to the hypocenters of these events gives its orientation: a strike of 230 • and a dip of 70 • .
The majority of the focal mechanisms also have a similar orientation as this feature, showing left-lateral strike-slip motion (see Figure 1 and Figure 6a of Clarke, Soroush, and Wood (2019)). This feature appears to be relatively well oriented within the in situ stress state in the region. Given the S Hmax orientation H of approximately 170 • , and a strike-slip stress regime (Clarke et al., 2014;Fellgett et al., 2017), faults striking to the northeast will also produce left-lateral strike-slip motion (rake of 0 • ).
The location of this feature does not correlate with any discontinuities observed in the 3-D reflection seismic that was acquired at this site (Clarke, Verdon, et al., 2019). This may be because of its strike-slip nature, meaning there is little vertical offset to be imaged in the reflection seismic. This seismogenic feature could be described as a "fault," or potentially as a zone of preexisting fractures. Despite the feature being around 500 m in strike, and 200 m in dip, the largest event during the monitoring had a magnitude of M L = 1.5. The basic formulation of seismic moment release for a circular fault of radius r , shear modulus G, and slip d is given by equation (2) (Aki and Richard, 2002). Figure 1. Hypocenters of events recorded by the downhole monitoring array during hydraulic fracturing operations at the Preston New Road PNR-1z well with magnitudes greater than −0.5 and a signal-to-noise ratio greater than 5. Events are shown as circles, with marker sizes indicating the magnitude range, while color shows the injection stage with which the event time overlapped. Diamonds denote the center of the sleeve position on the well and are also colored by stage. The gray plane denotes the inferred seismogenic "fault zone," with a strike of 230 • and a dip of 70 • . This was found from the least squares fit to events with M W above 0 and the events which continued to occur during the month hiatus in operations (see Clarke, Soroush, & Wood, 2019, for a detailed discussion). Lower hemisphere focal mechanisms are shown as black and white beach balls, derived from the surface station polarity data (Clarke, Soroush, & Wood, 2019).
A M = 1.5 event roughly corresponds to a displacement of ∼1 mm over a rupture length of less than 100 m. Thus, seismic failure on this feature only ever occurred on a small section of the suspected fault's area. Despite many small events occurring along its length, there is no clear evidence distinguishing if this is a single contiguous fault or a dense zone of fractures. Clarke, Soroush, and Wood (2019) term this feature "northeast Fault 1" (NEF-1). Thorough out this paper, we will refer to it as the "fault zone" adjacent to the wells.
Location uncertainties are naturally a concern when interpreting structure from microseismic data. In this case, with a single, mostly vertical, downhole array (as shown in Figure 2), there is the potential for systematic bias or offsets, due to its limited azimuthal coverage. However, the 3-D hodogram analysis, as well as the beam-forming inversion used in the location calculation should provide more accurate back-azimuths and polarity data than simpler methods. The locations found were also relatively similar to those independently Figure 2. Hypocenters for 172 events located using data from both the surface and downhole arrays, and the same velocity model, allowing for comparison of the two locations. These surface-derived locations were calculated by the British Geological Survey (Baptie, 2019). Naturally, the lateral and depth resolution is far lower than that of the downhole locations. However, these surface locations generally mirror the spatial and temporal trends seen in the downhole locations, with a bias (74%) of events north of the PNR-1z well, and events trending further NE as the heel stages are injected.
calculated by the BGS using the surface stations. These locations are shown in Figure 2. Broad-scale structure is generally the same, though naturally the precision of the surface-derived locations is significantly lower than that from the downhole. The velocity model was calibrated from the extensive 3-D seismic data and was refined several times during the stimulation of the well-when operations on the sliding sleeves occurred, the known times and locations were used to check its calibration. As these more involved methods of location inversion and velocity model refinement were used, we feel the locations provided are adequate enough to interpret the spatial distribution of seismicity around the well.
The structures interpreted in the microseismic, for example, the northward propagation of events, would also be difficult to systematically shift given some velocity model or station orientation error. Rotating the event clusters around the axis of the monitoring well, in order to shift the events in the center of the well, would shift events at the toe of the well to be propagating only south of the well. The offset between the injection well and the heel stage events (Figure 3f), could be attributed to the velocity model being incorrect. However, any kind of systematic shift in the velocity model, which could counteract the separation of heel stage events far from the injection well, would shift the events at the middle and toe of the well even further from the injection well. Thus, it is difficult to envisage purely processing errors resulting in the structure interpreted above.
Some locations for stages greater than 38 are subject to a processing artifact produced by the fundamental 180 • ambiguity when locating events with a single downhole array (e.g., Jones et al., 2010). The P wave particle motion is used to determine the back-azimuth of the event from the monitoring array. Events could therefore be placed at mirrored positions either side of the monitoring array. Evidently, the processing contractor has placed all of the events to the south of the PNR-2 well, when in reality events will have occurred both to the north and the south. This artifact does not affect the observations presented above, as a gap between the injection well and the events will be present whether or not events are placed to the north of the monitoring well.

Microseismic Observations in Detail
In this section, we focus on some noteworthy aspects of the microseismic event locations. Event hypocenters from stages illustrating behavior of particular interest are shown in Figure 3. We will describe these observations sequentially, in the order the stages were injected. Full injection stages were effectively completed in ascending order; however, small-scale "minifracs" were conducted on Sleeves 35 through 40 just before the start of the month-long hiatus, prior to Stages 37 through 41. Only small numbers of events were generated during these minifracs, with no particularly noteworthy behavior.
For all stages conducted at PNR-1z, the microseismicity occurred asymmetrically, propagating to the north of the injection well. This is unlikely to be a detection effect, as the sensitivity of the array is such that it is capable of detecting events at least 1 km from the well. However, it does not detect events south of the well even for the heel-most stages, which are within 300 m of the array. This suggests that hydraulic fractures grew primarily asymmetrically in a northward direction. This could also be related to more seismically productive, shearing type events occurring in the inferred fault zone in the area approximately 250 m north of the well. Asymmetric fracture growth has been ascribed in previous work to a gradient in the geomechanical parameters, such as a laterally heterogeneous stress field, a change in the elastic properties of the rock, or the result of using sliding sleeve as opposed to plug-and-perf completions (e.g., Maxwell, 2011;Chorney et al., 2016).

10.1029/2019JB018794
As can be seen in Figure 3a, during Stages 2 and 3, an isolated cluster of microseismicity occurred around 200 m northeast of the injection, north of the location of sleeve 12. There is a clear gap between the events adjacent to the toe stages (1-3), and this anomalous cluster, with only a small number of low magnitude events sparsely connecting the two. Figure 3b shows the microseismicity that occurred when the operator skipped forwards to stimulate Stage 12, which was roughly adjacent to the anomalous microseismicity observed during Stages 1 through 3. Here we observe microseismicity to the north of the well, connecting into the same cluster of events that occurred to the northeast of Stages 1 to 3. However, we observe little microseismicity to the west back near these toe stages: What little microseismicity that is observed here is primarily the postinjection tailing of events from the earlier stimulation, not a reactivation of events. It is interesting, therefore, to consider why activity around Stages 1-3 was able to create a cluster of microseismicity adjacent to Stage 12, but activity near Stage 12 was not able to have the obverse effect on microseismicity near the toe stages.
During Stage 18, very little fluid was injected (around 8 m 3 ). However, this stage produced a significant microseismic response, with over 1,200 events occurring in a cluster extending over 150 m to the north of the injection point. This stage generated relatively high magnitude microseismicity, with eight events above M w 0, and a M L 0.5 trailing event around one hour after injection ceased. It is very unusual for an injection volume of around 8 m 3 to create a hydraulic fracture over 150 m in length, and to produce such significant amounts of microseismicity. Events that took place in the 6 hr after injection had a combined moment release of 3.10×10 10 Nm. This constituted a notably large increase in the ratio of seismic moment release to injection volume compared to the previous stages. This is also relatively close to the upper bound of moment release proposed by the McGarr et al. (2002) relation, which for this small injected volume and a shear modulus of 25 GPa, would be around 2 × 10 11 Nm. Previous stages had a far lower "seismic efficiency" (Shapiro et al., 2010;Hallo et al., 2014), with moment release less than 0.1% of this theoretical upper bound for each of their injected volumes.
During Stage 22 (Figure 3d), the full planned volume of just over 400 m 3 was injected, however, with only around a third of the planned proppant (∼17 t). This was conducted in two separate injection periods on 25 October 2018. This stage generated a large number of events, around 5,700, with 12 events with M w > 0. During the first period, events propagated perpendicular to the injection well, appearing to trace the hydraulic fracture growth northwards from the well. However, in the second period, events began to extend laterally, both east and west of the initial line of fracture growth, clustering along the seismogenic "fault zone" described above (Clarke, Soroush, & Wood, 2019). Events extended along ∼70% of the feature's length, tracing back toward Stages 12-14, and extending north of Stages 30-32.
Stages 30 through 41 continued to interact with this seismogenic zone, with large numbers of events clustering further north of the well. However, events rarely propagated westward, back along this structure, that is, toward the stages that had been previously stimulated. This is shown in Figure 3e, for Stage 32. If it is assumed that this planar feature is a preexisting fault or a zone of preexisting fractures, one would anticipate that when stages reconnect to this seismogenic area, events would again be stimulated along its length, especially as the pore pressure around these faults or fractures has been increased by the previous injection, so we might expect successive injection would continue to stimulate seismicity back westward along its length. Stress relaxation may contribute somewhat to the limited reactivation as subsequent stages reconnect along the fault's length. However, previous cases of fault reaction have observed repeated reactivation into the same fault as injection reconnects .
The clear clustering of events at a notable distance from the injection well is apparent in Figures 3e and 3f, for Stages 32 and 38, respectively: Clusters of microseismicity are not centered at the point of injection. If microseismicity were being driven directly by elevated fluid pressures, then we might expect more microseismicity to occur near to the well. These gaps between the well and the focus of the microseismicity are seen for stages all along the well, although they are particularly prominent for the latter stages at the heel of the well (Stages 37-41). This absence of microseismicity immediately adjacent to the well could be due to the tensile opening of fractures being a more aseismic process than shear slip on small faults or fractures that is occurring within the fault zone.  (4)). Shapiro et al. (1997) show that, where microseismicity is driven by diffusion of pore pressure, it should develop along a characteristic triggering front that extends a distance r from the injection point as a function of time t:

Spatiotemporal Evolution of Microseismicity
where D is the hydraulic diffusivity. It has also recently been shown that the hydraulic fracture growth can produce similar r-t behavior (Barthwal and van der Baan, 2019). In contrast, a simple model of hydraulic fracture growth can provide the upper bound for the seismicity distribution. Under constant flow conditions and assuming minimal leak-off of fracturing fluid, microseismicity driven directly by hydraulic fracture propagation might be expected to show a linear distance-time relationship, since the length of the hydraulic fracture L scales with the injection rate Q, the height of the fracture h , and its width w (Economides and Nolte, 2003;: (4) Figure 4 shows examples of the r versus t behavior for several stages: These plots are typical for the PNR-1z microseismicity. In Figure 4 we also show the expected r versus t produced by the diffusivity approach (equation (3)) using various values of D, and for the hydraulic fracture propagation approach with minimal leak off (equation (4)), using approximate values of h = 25 m and w = 2.5 mm.

10.1029/2019JB018794
We do not observe the r ∝ t 1∕2 behavior, characteristic of diffusion-controlled microseismicity. Realistic values of diffusivity for hydraulically fractured rock are considered to be 1.0 m 2 s −1 (∼ 1 D) or less, which Figure 4 shows is clearly not adequate to describe the observed spatiotemporal distribution (Gehne and Benson, 2017;Tan et al., 2018;Gehne and Benson, 2019). Instead, we observe microseismicity occurring near instantaneously across a range of distances from the injection point. This behavior is weakly consistent with the linear relationship between r and t posited by equation (4) for hydraulic fracture propagation with minimal leak-off, because in such circumstances, given a typical flow rate at PNR-1z of 0.07 m 3 s −1 , we might expect a hydraulic fracture to propagate a distance of 300 m in less than 10 min. Note, however, that this is an upper bound, because in reality we expect multiple hydraulic fractures to form, sharing the overall injection volume between the fractures, and because equation (4) assumes that no fluid is lost to the surrounding formation.
The near-instantaneous onset of microseismicity, regardless of hypocentral distance from the well, implies that pore pressure diffusion is not driving the microseismic activity, as this would produce microseismicity growing outward from the well with time. In contrast, stress transfer effects occur instantaneously and so might provide a mechanism for fault reactivation that is more consistent with these observations.

Stochastic Hydraulic Fracture Model
To produce the loading, or sources, for our stress transfer simulations, we require estimates of the number of hydraulic fractures, their orientation, length and height, and the amount of tensile fracture opening that takes place. This can be done using coupled hydromechanical fracture stimulation codes (e.g., Warpinski et al., 1994;Profit et al., 2016), as commonly used by industry. However, such models are highly dependent on poorly constrained geomechanical input parameters, which may be tuned based on observations made during operations (Profit et al., 2016). Detailed modeling of this kind is beyond the scope of this study, which aims primarily to evaluate not the hydraulic fractures themselves, but their impact on the stress conditions in the surrounding rock. Instead, we adopt a stochastic approach, generating hydraulic fracture populations by drawing their properties (positions, orientations, dimensions, etc.) from statistical distributions representing typical, expected hydraulic fracturing cases. The use of a stochastic approach allows us to create thousands of model instantiations, such that we can identify features in the resulting deformation that are consistent across a range of input hydraulic fracture models and so may be considered robust and not dependent on a single choice of model parameterization.
We assume that both the lateral (i.e., along well) and vertical locations of the fractures are normally distributed around the sleeve location, producing an ellipsoid which extends to match the observed microseismic clouds, as well as those observed from other hydraulic fracturing sites (Urbancic et al., 2003;Chorney et al., 2016;Kettlety et al., 2019). This truncated normal distribution has a mean of 0 m, a standard deviation of 25 m, and a limit of ±100 m. For the stages with an obvious gap in microseismicity between the well and the cluster (e.g., Stage 38 and onwards), this assumes that the initial propagation and opening of fractures is mostly aseismic, and then the seismicity observed is the result of changes in stress that occur during injection, promoting slip in a more seismogenic area. Fractures are modeled as uniformly opening rectangular patches, oriented in the direction of S Hmax (strike of 170 • and dip of 90 • ) with an on average 10 • von Mises random perturbation to the geometry. Fractures are randomly set to propagate either north or south from the well, with a bias of 80% extending north, to match the observations from the microseismic data.
We use the analytical solutions for the opening of a Griffith crack, commonly employed in fracture modeling, to approximate the fracture width (Perkins and Kern, 1961). For the injection rates at PNR (0.07 m 3 s −1 ), a shear modulus of 25 GPa, a Poisson's ratio of 0.25 (believed to be appropriate for this setting, as described in section 3.2), and a fracture aspect ratio of 0.2, the fracture width is around 2.1 mm. The total number of fractures is then calculated by dividing the total volume of fluid injected in the stage by the total volume within the average 75-m-long fracture. We set fractures to have a fixed aspect ratio AR of L dip ∕L str = 0.2. Fracture lengths L str are sampled from a truncated normal distribution, with a minimum value of 25 m, a maximum of 250 m, a mean of 50 m, and a deviation of 50 m, with at least one fracture above 100 m in length. L dip is then calculated from the L str and AR. These values were again chosen to approximate the expected stimulated zone for each stage, as well as being comparable to hydraulic fracture dimensions estimated at other sites (accounting for the smaller injection volumes used at PNR-1z (∼400 m 3 per stage), compared to many wells in North America (>1,000 m 3 per stage)). Fracture width for each of the model fractures is then defined as the total volume of fluid injected divided by the total area of all generated fractures (d = V tot ∕ ∑n i L str,i L dip,i ). This gives a width very similar to that found using the solutions of Perkins and Kern (1961) or Nordgren (1972), with normally distributed values of 2.6 ± 0.3 mm for each set of fractures.
The modeled fractures are then ordered, with the longest fractures located closer to the center of the sleeve, producing an ellipsoidal stimulated volume of tensile opening fractures around each stage. An example of a fracture set produced in this manner is shown in Figure 5.

Modeling Stress Change
These opening patches are treated as the sources in the elastic stress change model. We use PSCMP developed by Wang et al. (2006) to compute these changes in stress. This approach uses the analytical Okada solution (Okada, 1992) for the Green's function for a homogeneous elastic half-space to calculate the strain field, and Hooke's law to find the resulting change in the stress field.
The resulting elastostatic stress changes within the volume around the well are resolved onto the receiver geometry of the fault plane identified in Figure 1-a of 240 • , of 70 • , and of 0 • -in order to compute the ΔCFS using equation (1).
The effective coefficient of friction ′ in equation (1) is derived from by ′ = (1 − ) and is an attempt to account for the way in which a change in pore pressure p effects the change in the normal stress Δ n (Rice, 1992;Simpson and Reasenberg, 1994). This is achieved through the Skempton's coefficient (Skempton, 1954), where, through a series of assumptions concerning the material properties of faults, it can be found that = −p∕ n . The value of ′ can range from 0 to 0.8 and varies between tectonic settings and lithologies. Typical values of ′ are generally around ′ = 0.4 ( = 0.7 and = 0.4), which we adopt here (King et al., 1994;Harris, 1998;Stein, 1999). We assume a shear modulus of 25 GPa, and a Poissons ratio of 0.25. These values have been used in previous studies on induced seismicity (e.g., Schoenball et al., 2012;Catalli et al., 2013;Pennington and Chen, 2017) and are consistent with laboratory measurements of the frictional and mechanical properties of shales (Kohli and Zoback, 2013;Islam and Skalle, 2013). These values are also similar to those found from studies of the Bowland shale, the formation targeted by PNR-1z (Herrmann et al., 2018). Using the stochastic process described above, we model 1,000 fracture set realizations for each stage. We compute the ΔCFS for each case and compute the median ΔCFS value for each point in the subsurface for each stage. We also examine the variability of the ΔCFS change across the 1,000 model instances: ΔCFS values that do not change significantly across a wide population of models can be considered robust. Figure 6 shows an example of the median modeled ΔCFS changes for Stage 22, and the variability introduced by our stochastic modeling approach. Lobes of negative Coulomb stress change dominate to the east and west of the hydraulic fractures, while positive lobes extend north and south of the fracture tips, as well as above and below. The variability within the zone of hydraulic fracture propagation is high. This is because the ΔCFS values in close proximity to opening fractures can be very high, and so modeled stress changes within this zone will be strongly dependent on the particular stochastically generated fracture model used as the input. However, further from the fracture zone, the median absolute difference in ΔCFS values is low. In these areas, the stress change is not sensitive to the particular stochastic fracture model used, and so can be considered to be more robust. In other words, the general distribution and shape of the lobes of positive and negative ΔCFS seen in Figure 5 exist for all fracture models that have tensile fractures extending roughly 100 m from the well. Therefore, the use of the median value allows us to examine the typical effect of the fracture sets, without the perturbations produced by the generation of random fractures.
To assess the significance of stress transfer effects, we interpolate the median modeled ΔCFS changes onto the location of each microseismic event, assuming the left-lateral faulting mechanism on the inferred plane. From this we compute the Coulomb Index, CI, which gives the proportion of events within a population that received positive ΔCFS changes. If stress transfer effects are playing a significant role, then we would expect most microseismicity to occur within lobes of positive ΔCFS, and therefore the CI would be high-typically >70% (e.g., Harris, 1998;Steacy et al., 2005;Catalli et al., 2013).

Model Scenarios
For a given stage, we compute the median ΔCFS values for three points in time. We compute the stress change created by all of the preceding stages-this represents the stress conditions at the start of the selected stage. We refer to this as the prior ΔCFS. We compute the stress change created by hydraulic fracturing of the stage in question. This shows the ΔCFS produced by that stage. We refer to this as the "current" ΔCFS. Finally, we combine the stress change from all preceding stages and the stage in question. This represents the overall ΔCFS conditions that will be present at the end of a stage. We refer to this as the total ΔCFS. Obviously, the "total" stress conditions and the end of one stage will be the "prior" stress change for the following one. Included in the supporting information are the complete set of figures for each stage, showing the current, prior, and total ΔCFS maps in multiple orientations.   It can be seen that for stages from 18 (those that encountered the seismogenic fault zone), CI is largely well above 50%, and frequently in excess of 70%. The heel stage (37-41), while not appearing to be significantly effected by stress triggering during each of the stages, show strong signals for the prior stages. Stages 3 and 18, both of which showed anomalous seismicity, show significant correlation between positive stress change and event hypocenter location, with CI in excess of 70%. role in controlling where microseismicity occurs; this role can be further demonstrated by considering the CI values, shown on a stage-by-stage basis in Figure 8. We find that the majority of the stages have high values of CI, consistent with microseismicity that is triggered by stress transfer, especially when the cumulative impact of multiple stages is taken into account. This effect appears to be particularly strong for the latter stages where reactivation of the fault zone was taking place.

Results
In Figures 9-11 we examine some of these stress transfer effects in more detail, with particular focus on some of the observations presented in section 2.1. Figure 9a shows a map of ΔCFS produced by Stages 1 to 3 at the toe of the well. In Figure 3a we observed a cluster of events occurring roughly 100 m to the northeast of the main event cluster. In Figure 9a we see that this region is at the center of a large positive ΔCFS lobe created by the tensile fracture opening. In contrast, during Stages 12 and 13, we did not observe microseismicity back-propagating in the reciprocal direction. Figure 9b shows the ΔCFS produced by Stage 12. We note that this region is within a lobe of negative ΔCFS. This stress-shadowing effect (Green et al., 2015) as the ΔCFS shifts from positive to negative as the hydraulic fracturing moves from west to east might explain why microseismicity appears able to propagate to the northeast ahead of the fracturing but is suppressed in the region behind the active stage. What seismicity persists in that stress shadow may be continuing due to the large increase in pore pressures from the injection into Stages 1 to 3 at the toe of the well. In (c) we show a map of ΔCFS produced by Stage 12, with the microseismicity produced this stage. The region to the west of this stage is now in a lobe of negative ΔCFS, and microseismicity is suppressed here. Figure 9c shows a cross section of the median ΔCFS produced by Stage 3. Positive lobes extend above and below the well, with a plane of null ΔCFS dipping at about 45 • . The events around the well fall within this lobe, which results in a structure that appears to dip at the same angle. Our interpretation is that this angle does not represent dipping hydraulic fractures, since in this strike-slip environment the intermediate principal stress is oriented vertically, but instead is caused by microseismic events being limited to this lobe of positive ΔCFS. Figure 10 shows the ΔCFS produced by all of the previous stages prior to Stage 18, and the microseismicity that occurred during Stage 18. This stage produced a surprisingly large microseismic response from an injection volume of less than 10 m 3 , with eight events above M w > 0 and events extending over 150 m from the injection point. In Figure 10 we observe that the locations of these events are strongly portioned into the lobe of positive ΔCFS produced by these prior stages, with a CI = 80%. Our interpretation is that the earlier stages caused prestressing of fractures in this region, such that a small perturbation in the stress state caused by the small injection volume was able to produce such a large number and extent of events. Figure 11a shows the ΔCFS produced by Stage 22. As for Stages 1 through 3, we observe a lobe of positive ΔCFS extending both above and to the northeast of the modeled tensile fractures, within which most of the microseismicity falls, with CI = 74% for this stage. Figure 11b shows the cumulative ΔCFS from all previous stages and Stage 38, with microseismic events from Stage 38 overlain. Again, we observe a very high CI = 80% for this scenario. Whereas during Stage 22 we observed northeastward propagation of events along the Figure 10. Map of ΔCFS changes produced by all stages prior to Stage 18, with the Stage 18 microseismicity overlain. Stage 18 saw minimal injection, yet produced significant amounts of microseismicity. In this figure we see that the effect of the prior stages was to create positive ΔCFS in this region. Figure 11. Maps of ΔCFS in stages toward the heel of the well. In (a) we show the ΔCFS produced by Stage 22, overlain with the microseismicity from this stage: A lobe of positive ΔCFS extends to the northeast, in which microseismicity is observed. In (b) we show the ΔCFS produced by all stages up to 38 (inclusive), and the microseismicity produced by Stage 38: the area to the west, behind the active stage is now in a region of negative ΔCFS, and microseismicity in this region is suppressed. fault zone, in these latter stages we do not observe significant numbers of events propagating back to the southwest. Figure 11 shows that the cumulative impact of the latter stages is to place this portion of the fault zone within a lobe of negative ΔCFS, and therefore, seismicity is less prevalent. This significance of this effect can be seen in Figure 8b: for Stages 30 to 41, when considering the cumulative impact of prior stages, the CI values are consistently at approximately 80%, indicating event location is consistent with elastic stress transfer. As hydraulic fractures are created during each stage, a lobe of positive ΔCFS is pushed toward the northeast, while a lobe of negative ΔCFS is created behind (i.e. to the west) of the active stage. This geometry of positive and negative ΔCFS lobes appears to have a strong control on whether the fault zone is, and is not, reactivated.
For a number of stages, including the example of Stage 32 shown in Figure 7, a number of the largest events (M W > 0) occur in areas of consistently negative median elastic ΔCFS, mostly near the injection point and the injection well. Obviously, this stress transfer effect is occurring contemporaneously as injection of hundreds of cubic meters of fluid at over 50 MPa. Clearly, stress transfer from fracture opening will not be the sole driver for seismicity during this case of fault reactivation. The increase of pore pressure, and the associated poroelastic stress change, immediately adjacent to the well will naturally give rise to seismicity in areas that receive negative elastic stress change on the order of 1 MPa.
Using the derivations of Rudnicki (1986) for pore pressure and poroelastic stress change in a 3-D homogeneous poroelastic medium, we can estimate the approximate magnitude and extent of pore pressure change ΔP for a Q = 0.07 m 3 s −1 , 90-min injection (the rate and pump time of the largest stages during PNR-1z operations). For this estimate we use an average matrix permeability around the injection point of 5 mD, a Biot-Willis coefficient of 0.7, a shear modulus of 20 GPa, a drained Lame parameter of 20 GPa, an undrained Lame parameter of 25 GPa, and a dynamic viscosity of the fluid of 1 mPa s. At the end of pumping the stage, this simple model gives a ΔP of at least 0.5 MPa out to a radius of ∼50 m from the point of injection, and within 10 m, ΔP exceeds 10 MPa. The change to the stress tensor from increased pore pressure provides a poroelastic Coulomb stress change on the receiver fault geometry of at least 0.5 MPa around 70 m NNW-SSW from the injection point. Twelve hours after injection, a ΔP of at least 0.5 MPa will extend out ∼100 m from the point of injection. The poroelastic stress decays rapidly as elevated pore pressures diffuse into the surrounding medium and decrease in magnitude, so by 12 hr after injection, poroelastic ΔCFS is less than 0.1 MPa 50 m from of the injection. Thus, both during the stage and after, the magnitude of stress changes from both the diffusion of elevated pore pressures and poroelastic ΔCFS are comparable to the fracture opening elastic stress transfer. Without a complex model of the permeability structure around the well, providing conduits for increased ΔP, the spatiotemporal distribution of events does not clearly correlate with the areas of increased poroelastic stress or pore pressure.
Interevent static Coulomb stress increase is most likely another mechanism contributing to the failure of events within the fault zone that receive negative stress change from opening fractures. As the several M w > 0 events occur, failing in a left-lateral strike-slip fashion in the fault zone, positive stress changes will extend around 100 m from the tips of the fault, encouraging continued failure along its length. This effect will naturally be combined with the static stress change from opening fractures, however the magnitude of the interevent stress changes will be smaller in comparison due to the relatively small size of the events. There is also no clear aftershock-type sequences in the spatiotemporal distribution of events that occur after the M w > 0 events, which would be a clear indicator of interevent triggering.
The spatial distribution of seismicity will naturally reflect the multiple mechanisms at play, and thus, only the elastic model of fracture opening will not account for every event's location. What is notable, however, is that during most injection stages, the majority of events are located in areas that do receive positive stress from fracture opening and that this mechanism provides a possible explanation for the unexpected observations in the microseismic.

Discussion
Using a simplified model of distributed fracture opening around a hydraulic fracturing well, we have seen that microseismic event locations were predominantly distributed in regions of positive stress change when resolved onto the geometry of an inferred adjacent fault zone. Specifically, unexpected microseismic event locations during several stages, which would otherwise be difficult to explain, are located in regions of positive stress as generated by a simple model of tensile opening of hydraulic fractures.

Model Uncertainties
The input parameters used in this model, such as fracture dimensions and distribution, or elastic moduli, are not overly tuned to this specific location or site-they are broadly applicable to most hydraulic fracturing cases. Model fractures are centered on the injection point and their locations follow fairly generic distributions for stimulated ellipsoids around an injection point. Thus, it is noteworthy that, despite this generality, many of the observations are consistent with static stress transfer promoting failure on the inferred failure mechanism of the larger fault zone. Naturally, the extent of the ΔCFS lobes are dependent on the fracture modeling parameters, such as the average length of the fractures and could thus be varied in order to increase or decrease the significance of the results. For example, model fracture growth could be offset by small distances (tens of meters), within the uncertainty, to shift most events into the areas of positive stress change. However, we found that generic values gave a clear indication of stress triggering, through good agreement between areas of positive stress change and event location, and consistently high CI.
The magnitude of the ΔCFS change will be sensitive to model assumptions, such as the shear modulus, and the modeled fracture opening. We do not take into account the effects of leak off or proppant during injection, as in our model the total amount of fracture opening is sufficient to contain all of the injected fluid. In reality, some of this fluid will be lost to the formation, reducing the total volume of fluid available to cause fracture opening. Since our model fracture lengths are chosen from a fixed distribution, and the fracture widths are constrained by analytical solutions (Perkins and Kern, 1961), the net effect of a reduced injection volume would be to reduce the number of fractures in the stochastic model. The overall deformation is computed by adding the deformation produced by each hydraulic fracture, so a reduction in the number of fractures would reduce the magnitudes of the modeled stress change but would not change the polarity of the ΔCFS change. This magnitude is already sensitive to the elastic parameters used, as well as the simplistic uniform-slip source model, which can lead to unreliable stress changes within the near field of the source (Steacy et al., 2004;Meier et al., 2014;Kettlety et al., 2019). Thus, we deliberately choose not to interpret this magnitude. Instead, we focus on the sign of the modeled ΔCFS (i.e., if microseismic events occur in regions experiencing positive ΔCFS), since this is far more consistent and robust than the magnitude. Most events within the positive lobes do receive stress changes in excess of the triggering thresholds for critically stressed faults, which range from 0.001 to 0.5 MPa (Kilb et al., 2002;Freed, 2005;Shapiro, Kummerow, et al., 2006).
Accounting for the effects of leak-off and proppant in the fracturing fluid can also affect the calculation of fracture width. Reducing net flow into the fracture by accounting for leak-off would decrease the calculated width, while proppant increases the slurry viscosity and would act to increase the width (Nordgren, 1972). However, accounting for these effects would not significantly modify the overall stress change shape as we estimate that the width of each individual fracture would only change on the order of 0.1 mm. This would only have a small effect on the distance to which the lobes propagate, which is more sensitive to factors such as the spatial distribution of fractures and the shear modulus. Thus, the width parameter affects the magnitude of the stress, rather than the sign of ΔCFS.
When modeling the deformation produced by cumulative stages, we assume that the hydraulic fractures from each stage remain open, and we linearly sum maps for the previous stages. This situation is unlikely to be the case in reality, because as pressures reduce after each injection stage, fractures will begin to close. However, the flowback volumes between stages were small, typically less than 20-25% of the injected volume (over the course of weeks during the hiatus period specifically), and some (though not all, see Clarke, Soroush, & Wood, 2019) of the stages had proppant injected, which would serve to keep hydraulic fractures open after injection stops. Therefore, the extent to which fractures closed after injection, reducing the magnitude of stresses that are transferred to subsequent stages, is not well constrained. Naturally, adding the stress change from some earlier stages by a different factor would have the effect of altering the prior and total ΔCFS, shifting the positions of some of the positive and negative lobes somewhat. However, more complex fracture modeling would have to be conducted to determine the relative amount of fracture closing during each stage, and thus the scaling of the effect of each individual stage, with time.
Therefore, the magnitudes of the ΔCFS values could be higher or lower than those we describe here, depending on the assumptions concerning the factors described above. However, our study is primarily concerned with the polarity of the ΔCFS signal: whether events occur in regions that are experiencing positive or negative ΔCFS change, as described by the CI value. The shapes of the positive and negative ΔCFS lobes are primarily controlled by three factors: the orientations of the hydraulic fractures, the assumed length of the hydraulic fractures, and the orientations of the receiving fractures on which microseismicity occurs.
The orientation of the hydraulic fractures is determined from the in situ stress state, which has been well constrained from borehole measurements within the PNR-1z well (Clarke et al., 2014;Fellgett et al., 2017;Clarke, Verdon, et al., 2019). The orientations of the receiving fractures have been determined by consistent, well-constrained source mechanism observations (Clarke, Soroush, & Wood, 2019), as shown in Figure 1. The lengths of the hydraulic fractures that we have used in our model are based on generic assumptions about hydraulic fracture lengths given the injection volumes used. However, they are similar to the fracture lengths, between 100 to 300 m, which have been calculated by the operator based on their observed pumping parameters (Cuadrilla Resources Ltd., 2019). Therefore, while the magnitudes of the ΔCFS values may not be well constrained, the spatial distributions of positive and negative values, and therefore our results expressed in terms of the CI, can be considered to be robust.

Possible Impact on Fault Rupture Dimensions
Assuming the basic formulation of seismic moment given in equation (2) holds, maximum earthquake magnitude would be controlled purely by the dimensions of the fault on which induced seismicity is being triggered. For the feature identified in Figure 1, assuming a typical stress drop value of a rupture (∼ 1 MPa) along a 500-m by 200-m area, this corresponds roughly to a M 3 event. The largest event size during the operations had M L = 1.5, approximately 30 times smaller than this potential maximum magnitude, corresponding to a rupture radius of less than 100 m as discussed earlier. Our modeling shows that the ΔCFS values on the fault were positive in some places, but negative in others. This clamping at certain points along the fault, in particular the regions behind (i.e., to the west of) the active stage, could be seen as a mechanism for the limited rupture extent on this inferred fault plane. However, previous studies have shown that rupture extent is not limited to the portion of a fault zone receiving positive stress during failure along its length (Ripperger et al., 2007;Ampuero and Rubin, 2008). Dynamic stress changes during rupture can quickly overcome regional stress and local, smaller-scale stress changes (Meng et al., 2012;Preuss et al., 2019). Also, it is certainly not clear this zone is a well-connected fault surface or just a region of preexisting fractures that are oriented favorably in the present regional stress state. Thus, the likelihood of a M 3 event is not well constrained.
Many of the proposed mechanisms for constraining the maximum magnitude during an induced sequence (e.g. Shapiro et al., 2011) function under the assumption of a limited rock volume stimulated by injection. Shapiro et al. (2011) assume that seismicity is driven by pore pressure diffusion; however, an analogous argument could be made with respect to the dimensions of the portion of the fault that receives positive ΔCFS. Fracture opening does introduce significant changes to the stress state in hydraulic fracturing settings, and for well-oriented faults adjacent to the operations (i.e., where the magnitudes of stress transfer are significantly positive) this could modify the extent and shape of the "stimulated" rock volume greatly. While this clamping effect is a possibility for general cases of fracture opening stress transfer, the model proposed by Shapiro et al. (2011) produces a truncated Gutenberg and Richter (1944) distribution, which is not observed Journal of Geophysical Research: Solid Earth 10.1029/2019JB018794 at the PNR-1z site (Clarke, Soroush, & Wood, 2019). Thus, it is by no means clear that this is occurring in this case of injection-induced fault activation.

Conclusions
During hydraulic fracturing at PNR-1z, we observed the reactivation of a preexisting fault that produced tens of thousands of microseismic events, the largest of which was felt by nearby populations, and several of which required the operator to pause their activities under the conditions of the United Kingdom's TLS. Here, we have investigated the role of elastostatic stress transfer in triggering these events, as well as producing other microseismic observations that are not obviously driven solely by injection-induced pore pressure increases or the growth of hydraulic fractures.
To do this, we develop a stochastic approach to modeling hydraulic fractures as a loading source for the elastic stress transfer model. This allows us to assess the impact of expected, generic fracture sets, without being overly influenced by the results of a particular representation of the hydraulic fractures. We then look at the median ΔCFS of the 1,000 realizations that were conducted.
We find that the observed microseismicity occurs predominantly within volumes of rock that receive positive median ΔCFS. This indicates that stress transfer effects produced by the tensile opening of hydraulic fractures are in part driving the spatiotemporal distribution of induced seismicity at PNR-1z. These elastic effects, while often considered to be less significant than the increase in pore pressure, appear to play a role in prestressing nearby fractures or faults, as well as promoting failure near instantaneously at anomalously larger distances from the point of injection.
For the particular orientations of the hydraulic fractures and the preexisting fault at PNR-1z, the tensile fracture opening creates positive ΔCFS to the northeast of the active stage, with multiple stages adding cumulatively to this effect. Because stimulation progressed eastward along the well, each new stage was therefore injecting into a volume of rock that had been prestressed by the previous stage. This may have contributed to the repeated exceedance of the TLS threshold over multiple stages. In contrast, the regions to the west of the active stage were clamped by the tensile fracture opening, suppressing microseismic activity in these areas. This implies that if the wells were drilled in the opposite E-W direction, proceeding injection stages would have actively clamped the fault, rather than stimulating it further. The fault was not identified on any of the 3-D reflection seismic data that were acquired for the site, however, and thus, it was not possible to know its orientation prior to the fault being reactivated. These effects will be highly dependent on the specific orientations of both the hydraulic fractures and the receiving faults and so cannot easily be generalized to other sites. However, the stochastic modeling approach, combined with the PSCMP modeling code, is able to provide results at a speed that could plausibly be applied in near real time during injection operations. Doing so could enable operators to identify whether their planned stimulation program is likely to stress or to clamp any faults identified during injection, and potentially to make appropriate adjustments to their program to minimize induced seismicity.