Coupled poroelastic modelling of hydraulic fracturing-induced seismicity: Implications for understanding the post shut-in ML 2.9 earthquake at the Preston New Road, UK

,


Introduction
Hydraulic fracturing has proved to be an effective technique to commercially exploit oil and gas resources from low-permeability reservoirs that are otherwise considered uneconomical.This technique involves the pumping of pressurised fluids into the subsurface to create a fracture network that acts as a permeable channel to increase the production of hydrocarbons from low-permeability formations.However, hydraulic fracturing operations in some formations have faced major challenges in terms of induced seismicity (Atkinson et al., 2020;Schultz et al., 2020).High pressurised fluids have the potential to activate pre-existing fractures/faults, either directly or indirectly, which results in induced seismicity, with moderate-size earthquakes that have been felt at the surface.Regulators have responded to induced seismicity concerns by imposing Traffic Light System (TLS) mitigation schemes (e.g., Verdon & Bommer, 2021) in some jurisdictions, which have resulted in the suspension or termination of operations at several different sites around the world (Schultz et al., 2020).
In the majority of field sites that have been affected by induced seismicity, the seismicity rate peaks during the hydraulic stimulation stage, followed by diminished level of seismicity after completion of the well (Schultz et al., 2020).However, it is not uncommon for seismicity to persist for days or even months during the shut-in phase, and this is thus referred to as the "trailing effect".For a number of hydraulic fracturing cases, in both shale gas development and enhanced geothermal systems (EGS), the largest earthquake has also occurred after cessation of fluid injection.Examples of hydraulic fracturing operations with significant trailing effects include the South Sichuan Basin, China (Lei et al., 2019) and Preston New Road, UK (Kettlety & Verdon, 2021); examples of EGS sites include Soultz-sous-Forêts, France (Evans et al., 2005), Basel, Switzerland (Häring et al., 2008), Paralana, Australia (Albaric et al., 2014), and Pohang, Korea (Grigoli et al., 2018).The trailing effect poses a challenge to the management of seismic risk using TLSs (Verdon & Bommer, 2021), since they represent a retroactive measure (operations are ceased after an event of a given magnitude).
The use of maximum seismic magnitude forecasting models provide an alternative to TLSs for induced seismicity mitigation (e.g., Cao et al., 2020;Clarke, Verdon, et al., 2019;Kwiatek et al., 2019;McGarr, 2014;Verdon & Budge, 2018).These models relate the cumulative seismic moment to the injected volume, and are thus highly dependent on the timely update of operations data.The cessation of fluid injection terminates the monitoring of operations data as model inputs, which impacts the applicability of such models.Regardless of whether TLSs or forecasting models are used, when trailing red light events are detected, few effective mitigation strategies have been identified to alleviate further seismicity, since, by definition, for trailing events the injection has already been stopped.
The driving mechanisms for the trailing effect have not been well understood.Some insights could be obtained by referring to three fundamental mechanisms for induced seismicity: (1) the direct increase in pore pressure (Elsworth et al., 2016;Talwani & Acree, 1985), (2) poroelastic stress perturbations (Segall, 1989;Segall & Lu, 2015), and (3) fault slippage induced stress transfer (Cao et al., 2021;Eyre et al., 2019;Guglielmi et al., 2015;Schoenball et al., 2012).The most common explanation for the trailing effect pertains to the delayed effect of pore pressure increase after fluid injection (Baisch et al., 2010;Hsieh & Bredehoeft, 1981;McClure & Horne, 2011;Parotidis et al., 2004).Pore pressure perturbations continue to propagate away from the injection point after the shut-in of the well, which can produce further seismic activity if additional faults are encountered.Another explanation is based on the poroelastic effect caused by fluid injection (Segall & Lu, 2015).Under certain circumstances where injection-induced poroelastic stresses inhibit slip, the abrupt shut-in would cause relaxation of poroelastic stresses, and in turn heightened seismic activity.Poroelastic coupled modelling results of injection-induced fault slippage have suggested that the poroelastic effect could lead to a surge in the post-injection seismicity rate, and that the permeability of faults and hydraulic connectivity of faults are crucial factors governing this process (Chang et al., 2018;Chang & Segall, 2016).As a secondorder triggering effect, the stress transfer through aseismic creep subjected to delayed pore pressure diffusion may also play a crucial role in triggering trailing events (Eyre et al., 2020).Elevated pore pressure results in stable sliding on a fault, and subsequently co-seismic slippage of unstable regions.
The persistent stable sliding of the fault and continuous loading on unstable regions account for the long-lived nature of post-injection seismic swarms.It was also argued that a large volume of seismic swarms previous linked to fluid diffusion can be alternatively explained by aseismic slip (Eyre et al., 2020).
In additions to explanations based on the three fundamental mechanisms for induced seismicity, a number of novel hypotheses that consider the shut-in conditions have also been proposed.One alternative explanation is that pressurised fluids in dead-end fractures backflow into larger fractures during the shut-in phase, potentially generating even larger events than those occurred during injection (McClure, 2015).Ucar et al. (2017) attributed the sustained post-injection seismicity to the normal closure of fractures after ceasing the injection, which acts as a fluid pressure support to advance the pressure front away from the injection point, increase apertures of fractures beyond the near-well region, and cause seismic events.It has also been recognised that the superposition of various mechanisms, such as the direct fluid pressure increase, stress transfer through fault slippage, and thermal effects (primarily in EGS) may have contributed to the persistent post-injection seismicity (De Simone et al., 2017).Some faults may tend to be destabilised by mechanical and thermal effects but held stable by the hydraulic effect during injection.The abrupt termination of injection resulting in sudden pore pressure decrease may trigger such faults to rupture.The various response times of reservoir formations to hydraulic, mechanical and thermal effects further complicate the variations of the superposed stress field and the fault stability.
The identification of causal mechanisms for post shut-in seismicity usually requires integrated interpretations of geophysical observations, hydrological properties and the geomechanical response.
However, insufficient field monitoring data and large uncertainties in hydrological and geomechanical properties often impede such efforts.It is also difficult to preclude the possibility that more than one mechanism is at play in complex geological settings and fault orientations.Nevertheless, certain characteristics of post-injection seismicity may help identify or at least constrain the plausible mechanisms involved in field observations.For example, the delayed occurrence of post shut-in seismicity favours the delayed pore pressure diffusion mechanism or the aseismic slippage stress transfer mechanism, whereas the immediate post-injection occurrence and long distance away from the injection region indicate the poroelastic stressing mechanism.Some features may be difficult to explain with the delayed pore pressure diffusion mechanism, such as prolonged duration (up to several months), steady seismicity rate over time, and lack of hypocentre migration, indicate a role for the aseismic slippage stress transfer mechanism (Eyre et al., 2020).Seismicity event locations that are beyond previously stimulated regions (e.g., in Basel and Paralana) may suggest one or more of the delayed pore pressure diffusion mechanism, aseismic slip stress transfer mechanism, and postinjection fracture normal closure mechanism.So far, the delayed pore pressure diffusion mechanism has been considered as the primary driving mechanism for post-injection seismicity in several hydraulic fracturing sites, e.g., the South Sichuan Basin, China (Lei et al., 2019), and the Red Deer region in Alberta, Canada (Wang et al., 2020).The aseismic slip mechanism was favoured to explain the persistent post-injection seismicity in the Fox Creek region in Alberta, Canada (Eyre et al., 2020).

Field observations of other recent novel hypotheses have not yet been reported.
Building upon different causal mechanisms of trailing events, various countermeasures have been proposed to mitigate against seismic risk.For example, in regards to the poroelastic stressing mechanism, tapering fluid injection rate instead of abrupt termination upon shut-in could lower or even eliminate the post shut-in spike in seismicity (Segall & Lu, 2015).Concerning the fluid backflow from dead-end fractures mechanism, reducing wellhead pressure (by initiating flowback immediately after injection) may be effective in alleviating post shut-in earthquakes (McClure, 2015).
In this work, the causal mechanism for a post shut-in M L 2.9 earthquake at the Preston New Road, UK, has been investigated through coupled poroelastic modelling that considers poroelastic solid deformation, fluid flow in both porous rocks and fault structures, and hydraulic fracture propagation.
A recent work by Kettlety and Verdon (2021), where elastostatic stress modelling was performed under a representative ambient stress field, has suggested that this earthquake was most likely driven by delayed pore pressure diffusion.The objective of the current work is to reveal the evolving stress, pore pressure and seismicity rate on the activated fault during and after the hydraulic fracturing operations, and to ascertain the respective contribution of delayed pore pressure diffusion and poroelastic stressing towards the post-injection fault slippage.In particular, we have examined three plausible mechanisms for the occurrence of the post shut-in earthquake, i.e., the post shut-in pore pressure diffusion, poroelastic stressing on a non-overpressurised fault, and poroelastic stressing on an overpressurised fault.We have also investigated the role of fault permeability and its connectivity to injection regions on the hydromechanical behaviour of reservoir formations and associated seismicity on the fault.

Post shut-in M L 2.9 earthquake at Preston New Road
In 2011, hydraulic fracturing operations commenced at the Preese Hall site in Lancashire, UK, marking the first onshore shale gas exploration in the UK (Clarke et al., 2014).In 2017, two horizontal wells (PNR-1z and PNR-2) were drilled in preparation for resumption of hydraulic fracturing tests at the Preston New Road (PNR) site, some 2.5 miles from Preese Hall.The PNR-1z well targeted the upper-most section of the Lower Bowland Shale at 2.3 km depth, while the PNR-2 well was drilled approximately 250 m to the north of the PNR-1z well through the lower-most section of the Upper Bowland Shale at 2.1 km depth.A surface monitoring array (broadband seismometers and geophones) and a downhole geophone array situated in the adjacent well were installed to monitor microseismicity associated with fracturing operations.All operations were regulated by a TLS, where operations would proceed with caution when the seismic magnitude reaches a M L 0 threshold, and be suspended for a minimum of 24 hours after reaching a M L 0.5 threshold.
Fracturing operations at the PNR-1z well commenced on 15 October until 17 December 2018 in 16 stages, with a maximum injected volume of 431 m 3 per stage.An M L 1.1 event occurred at the end of October, which triggered the TLS red light.Operations remained suspended throughout November.
Hydraulic fracturing resumed to complete 5 further injection stages at the heel of the well in December 2018.An M L 1.6 event occurred during this time, and operations were paused for around 48 hours.Both red light events were believed to be related to a seismogenic planar structure referred to as the PNR-1z fault (shown as the grey plane in Figure 1).The fault geometry was illuminated by microseismic event locations, and the trend of the fault is aligned with the focal mechanisms of the largest events.The microseismic monitoring, processing and interpretation at the PNR-1z well were detailed in Clarke et al. (2019) and Kettlety et al. (2020).
Hydraulic fracturing at the PNR-2 well took place during the period 15-23 August 2019, which sequentially stimulated 7 sleeves evenly spaced at 14.5 m from the toe of the well.The first 6 stages were operated to inject the full volumes of fluids and proppants as planned, with a maximum injected volume of 432 m 3 per stage.After stage 6, seismicity began to escalate and magnitude M L > 0.5 events occurred, resulting in a pause in injection.During stage 7, the well received a reduced volume of injection fluids with increased viscosity to alleviate the risk of seismicity.However, elevated levels of seismicity continued to occur after the end of injection, ultimately culminating in the occurrence of the M L 2.9 earthquake.The regions stimulated and fault structures activated from the PNR-1z and PNR-2 operations were believed to be hydraulically isolated, since almost no overlap exists between microseismic event locations from the two wells (Kettlety et al., 2021).Here, we focus on the spatiotemporal distribution of induced seismic events including the largest M L 2.9 earthquake during the PNR-2 operations.(2021) and Kettlety and Verdon (2021).We reprise key aspects on that analysis here.In the PNR-2 operations, microseismic events induced during stage 1 formed a cluster extending approximately 50 m above and below the well, and 150 m to the north and south, centred on the injection location.The cluster is closely aligned to the maximum principal stress direction, with a strike of 350° in the northward-propagating segment, and 155° in the southward-propagating one.This cluster was believed to represent the main hydraulic fracture emanating from the PNR-2 well, and was referred to by Kettlety et al. (2021) as the NS Zone (shown by the blue cluster in Figure 1).The NS Zone was driven by stage 2 operation to extend roughly 200 m northwards and 100 m southwards, followed by being maintained over stage 3 and 4 operations.Further to the west of the main NS zone, a cluster containing a smaller number of microseismic events also developed from stage 2 onwards.This cluster generally followed the maximum horizontal principal stress orientation, but manifested as a more diffusive feature.This microseismic cluster (shown by the green cluster in Figure 1) was also interpreted to result from hydraulic fracturing, and was referred to as the Western Cluster by Kettlety et al. (2021).
After the stage 4 operation had stopped, a new seismogenic zone emerged approximately 100 m to the east of the main NS Zone, and slightly deeper than the well.This cluster with a height of approximately 60 m, gradually propagated around 50 m southwards along the maximum horizontal principal stress orientation.This cluster (shown by the yellow cluster in Figure 1) was believed to represent another hydraulic fracture extending from the PNR-2 well, and was referred to as the Eastern Zone by Kettlety et al. (2021).During stages 5 and 6, microseismic events continued to occur within the NS Zone.In addition, the length of the Eastern Zone was extended to approximately 100 m and then to 300 m to both the north and south of the well.Most microseismic events induced during the stage 7 operation were restricted within both the NS and Eastern Zones.Roughly 5 hours after injection of stage 7 had ceased, a sequence of earthquakes of magnitude in excess of M L 1.0 occurred, including the largest M L 2.9 earthquake, which occurred 66 hours after the end of stage 7. Kettlety et al. (2021) used the aftershock locations determined from the downhole array to illuminate the fault structure as the source of the M L 2.9 earthquake.The aftershock cluster defined a near-vertical seismogenic planar fault measuring 330 m × 250 m (length × height), and extending to the southeast of the Eastern Zone (Kettlety et al., 2021).Integrated interpretations of the M L 2.9 earthquake focal mechanism and the fault plane fitting to the seismic cluster have suggested that the fault has a strike/dip/rake of 135°/80°/180°.This fault was denoted as the PNR-2 fault (shown as the red plane in Figure 1).
The fault activation mechanism responsible for the largest earthquakes at the PNR site is of particular interest.Kettlety and Verdon (2021) investigated the fault triggering mechanisms of both PNR-1z and PNR-2 faults through elastostatic stress modelling of the two hydraulic fracturing operations and the spatio-temporal evolution of microseismic event locations.To evaluate the impact of hydraulic fractures on stress conditions in the reservoir formation, they adopted a stochastic hydraulic fracture model, where a population of hydraulic fractures were generated following statistical distributions of fracture geometrical attributes.This modelling approach allowed them to examine median Coulomb stress changes on target faults and the variability from multiple model realisations.It was found that PNR-1z was activated by the compound effects of direct pore pressure increase and stress transfer caused by hydraulic fracture opening.The PNR-2 fault was most likely governed by the post shut-in diffusion of increased pore pressure through hydraulically stimulated regions.However, the stress transfer produced by hydraulic fracturing opening may also have contributed to destabilise the fault.
The difference between triggering behaviour of the two faults was attributed to the fault orientation in respect to the in-situ stress field: the PNR-1z fault is moderately well oriented to slip, whilst the PNR-2 fault is extremely well orientated to slip.
The elastostatic stress modelling of Kettlety and Verdon (2021) has qualitatively demonstrated the respective contribution of pore pressure change and poroelastic stress towards fault slippage by representing the most representative stress state that would exist during the hydraulic fracturing operations.However, this model did not capture the evolution of the pore pressure and stress state on target faults during the fracturing operations, and more importantly, after the end of injection.In addition, it is still unclear how the respective contribution of pore pressure change and poroelastic stress towards fault slippage varies, depending on the injection stages, stimulated regions, and hydraulic properties of the reservoir.This requires a more complex coupled hydromechanical model that considers both the hydraulic fracture propagation and injection pressure history in order to reveal the time-varying stress and pore pressure changes on target faults as the hydraulic fractures propagate and the reservoir is stimulated.

Governing equations
The theory of linear poroelasticity has been used to describe the hydromechanical behaviour of porous media such as subsurface rocks (Wang, 2017).The poroelastic constitutive equations consist of a set of six equations that describe the solid deformation as a function of the stress and pore pressure, and an equation that describes the pore fluid mass related to the pore pressure and mean stress.The first set of constitutive equations for an isotropic, linear elastic porous medium relate the strains ε ij to the stresses σ ij and pore pressure p: where ν, K and G are Poisson's ratio, bulk modulus and shear modulus of the porous medium, respectively,  is the Biot's coefficient, and δ ij is the Kronecker delta.
The other constitutive equation relates the increment of fluid content  to the pore pressure p and mean normal stress σ kk /3: where B is the Skempton's coefficient.
The geomechanical deformation of the poroelastic medium is based on stress equilibrium expressed in terms of the linear momentum balance equations: where f i is the body force.The effective stress σ' ij is defined by the stress σ ij and pore pressure p: By substituting Equation ( 1) and considering the compatibility relations between the strain and displacement ,, 1 () , the stress equilibrium equations can be expressed in terms of the displacement: The fluid flow in the porous medium is described by the mass conservation equation: where q i is the flux of fluid flow, and Q is a fluid mass source.Darcy's law for fluid flow in porous medium takes the form: where k is the permeability of rocks, and µ is the fluid viscosity.Substituting Equations (2)(7) to Equation ( 6) gives: The linear poroelasticity of the medium involves the two-way coupling between geomechanics and fluid flow.Both coupling terms are implemented in the constitutive equations: the fluid-to-solid coupling is reflected in the influence of the pore pressure on the strain (Equation ( 1)), and the solid-tofluid coupling is considered by the influence of the mean normal stress on the fluid mass (Equation ( 2)).The coupled poroelastic response could be simulated by solving the governing equations of geomechanics and fluid flow incorporating these constitutive relations.

Coupled poroelastic reservoir model
A 1,000 m-long cubic hydromechanical model was constructed to simulate the PNR-2 hydraulic fracturing operations and the associated hydromechanical response of shale formations (Figure 2).For simplification, the model was considered to be comprised of shale formations with uniform mechanical and hydrological properties.The fluid injection-induced geomechanical response was modelled using the linear elastic constitutive model, with Poisson's ratio υ = 0.29, and Young's modulus E = 25.7 GPa (Verdon et al., 2020).The matrix permeability of the Bowland Shale was estimated to be typically less than 1 × 10 -4 mD (Clarke et al., 2018), therefore, the shale formations at the PNR-2 site were assumed to have a permeability k = 1×10 -4 mD before stimulation.Generic values were used for other hydrological properties: porosity ϕ = 0.1, and Biot coefficient = 0.8.The tectonic stress at the PNR site is characterised by a strike-slip fault regime.Following the fracture growth trajectories delineated by microseismic clouds, the maximum principal stress was estimated to orient at 173 at the PNR site.The stress gradient of the maximum, intermediate (vertical) and minimum principal stresses are 0.032, 0.026 and 0.017 MPa/m, respectively (Clarke, Soroush, et al., 2019;Verdon et al., 2020).The pore pressure gradient is 0.012 MPa/m.The model was initialised with both the in-situ stress and the initial poroelastic strain caused by the in-situ pore pressure, so as to achieve a uniform initial stress distribution at the same depth before hydraulic fracturing operations begin.The boundary conditions were set up in such a manner that the model base is fixed, and normal and shear stress components calculated from the in-situ stresses were applied to top and lateral boundaries.The initial pore pressure was vertically distributed based on gravitational equilibrium of fluids using a density of 1,200 kg/m 3 .Fluid pressures computed from the pressure gradient was applied to all the outer boundaries to reach an initial pore pressure equilibrium.
The propagation of the two major hydraulic fractures was modelled in such a way that the fractures initiate upon the onset of hydraulic fracturing operations, followed by progressive extension to the maximum lengths over the fracturing period (Zeng et al., 2021).Considering that microseismicity began to appear along the designated hydraulic fractures during the main stages, hydraulic fractures mainly propagate during and shortly after main stages, when the bottomhole pressure exceeds the minimum principal stress at 2,100 m depth (around 35.7 MPa).For simplicity we assumed that hydraulic fractures only propagate within the main stages of hydraulic fracturing, when the maximum injection rate maintained for a period (Figure 3b).It is acknowledged that hydraulic fractures may propagate sublinearly over time as indicated by analytical solutions such as KGD and PKN models (Rahman & Rahman, 2010).Here we assumed a linear approximation and used a constant propagation velocity during each operation stage for simplicity.This model is not intended to accurately simulate the physical process of fracture propagation, but to represent the spatio-temporal overpressure distribution in designated hydraulic fracture paths and its influence on distant fault structures.As the hydraulic fractures propagate, the permeability of surrounding reservoir rocks (within a certain stimulated width given later) is elevated from 1×10 -4 mD to 100 mD.
Although the injection ports from the 7 injection stages are spaced by 14.5 m, it is believed that hydraulic fracture branches are well connected according to microseismic observations.Thus, fluid injection into the NS and Eastern Zones was modelled by applying overpressure on the two main hydraulic fractures planes, instead of on the respective injection ports along the PNR-2 wellpath.The bottomhole pressure history at the fracturing depth back calculated from the wellhead pressure history was used as injection pressure inputs (Figure 3a).The NS Zone initiates since stage 1, and the Eastern Zone is not activated until after stage 4. The field wellhead pressure history was recorded until 25 August, and the bottomhole pressure used after this date evolves with a leak-off type pressure decrease following an exponential function.The use of the injection pressure control in the numerical model allows for the accurate modelling of stress perturbations resulting from fluid injection, while the injection volume, which is not the focus of this work, is not explicitly represented.To understand the dominating factor of the fault slippage, the fault permeability and its hydraulic connection with injection regions were varied to assess the fault behaviour in various scenarios.Two end members of fault permeability, 100 mD and 110 -6 mD, were considered to represent conductive and sealing fault scenarios in the numerical model.The fault permeability was represented by assigning a uniform aperture over the PNR-2 fault plane, based on the relationship between the effective fracture permeability k and fracture geometry (width h and aperture e), given by k = e 3 /(12h).
Assuming a 10 m wide fault zone, apertures of 0.23 mm and 0.00049 mm provide effective permeabilities of 100 mD and 110 -6 mD, respectively.The hydraulic connection between hydraulic fractures and the PNR-2 fault was varied by controlling the width of hydraulically stimulated regions.
We used 0 m, 100 m and 200 m stimulated widths to represent hydraulic isolation of the PNR-2 fault, hydraulic connection only to the Eastern Zone, and hydraulic connection to both the NS and Eastern Zones, respectively.A total of six model scenarios were considered in the model, as listed in Table 1.
The finite element method-based solver COMSOL Multiphysics was used to solve the coupled poroelastic model.A maximum timestep of 10 mins was used due to accuracy considerations in the fluid flow modelling.

Coulomb failure stress evaluation
The potential for fracture slippage can be evaluated by the Coulomb failure stress change along fracture planes: where f is the friction coefficient, and σ

Seismicity rate model
Dieterich (1994) developed a model to quantify the rate of earthquake occurrence based on the assumption that the timing of a sequence of earthquake nucleation events is controlled by the initial conditions of nucleation sources and the stressing history.The in-situ effective normal stress at 2,100 m depth at the PNR site is 10.5 MPa.We assumed the constitutive parameter a = 0.005, and the background stressing rate 0  is 10 -3 MPa/yr, such that 1 MPa stress along the fracture plane accumulates in 10 3 years.As a result, a characteristic decay time is t a = 52.5 yr.In our model, we assumed that seismicity will only occur in regions with abundant fractures, such as within the NS and Eastern Zones and the PNR-2 fault.Figure 4 (d-f) presents the mean normal stress change σ kk /3 from the three models.In the two coupled models, the overpressure causes the expansion of the shale formation, which is resisted by surrounding rocks.This leads to more compression of rock matrix within stimulated regions but less compression outside, as shown in Figure 4 (e and f).In the two-way coupled model, the rock expansion outside the stimulated regions creates liquid sinks and thus results in pore pressure reduction immediately surrounding the stimulated regions (Figure 4 c).The pore pressure distribution in turn dictates the normal mean stress distribution, resulting in a much larger stress perturbation as compared to the one-way coupled model.In both coupled models, poroelastic stress generated is much smaller than the overpressure within the stimulated regions, but has a much larger extent of influence than the overpressure.In following sections, numerical results presented are from the twoway poroelastic coupled model.sealing or isolated fault scenarios, it begins to receive fluids as long as the hydraulic fractures impinge on the fault in conductive fault scenarios, such as at stage 7 for 100 m stimulated width (Figure 5f), and at stage 3 for 200 m stimulated width (Figure 5i).The pore pressure diffusion within the stimulated regions forms a pressure gradient, indicating less poroelastic stressing away from the hydraulic fractures.Consequently, the attenuation of pore pressure surrounding stimulated regions due to poroelastic stressing is less apparent for large stimulated widths, in particular a 200 m stimulated width.

Coulomb failure stress change
Figure 8 presents the Coulomb failure stress change due to poroelastic stressing τ s + fσ n after injection stages 3 and 7 for the six model scenarios.The resemblance between Figure 8 and Figure 6 suggests that the poroelastic stressing effect is dominated by the normal stress changes.In particular, the relief of the normal stress prevails over the elevation of shear stress within the stimulated regions.
Nevertheless, strong negative shear stress changes to both sides of the NS and Eastern zones outside the stimulated regions contribute to the inhibition of the potential for fault slippage away from the hydraulic fractures.In isolated or sealing fault scenarios, the northwest end of the PNR-2 fault falls within this seismic inhibited region.To examine the contribution of different stresses towards the Coulomb failure stress change τ, two measurement points B and C, spaced by 50 m apart, were set up on the measurement line A-A (Figure 9a).Point B is on the extension line of the Eastern Zone, thus would be subjected to the largest           In hydraulically connected scenarios, the seismicity rate within the PNR-2 fault jumps to high levels when the hydraulic connection is established at injection stage 6 for a stimulated width 100 m, and at stage 2 for a stimulated width 200 m, so as the seismicity rate computed based on pore pressure change rate fp  (Figure 16 c).In hydraulically isolated fault scenarios, the surge in the seismicity rate occurs when the poroelastic stress becomes prominent at injection stage 6, regardless of a conductive or sealing fault (Figure 16 c, d).The fault is subjected to minimal mean pore pressure disturbance in both scenarios, as shown in Figure 16 (a)(b).Interestingly, in the former scenario, the increased pore pressure spreads across the fault plane leading to a gentle overall seismicity rate increase (Figure 16 c), but in the latter, the increased pore pressure is localised surrounding fracture tips, causing a relatively large seismicity rate increase even after being averaged across the full fault plane (Figure 16 d).

Seismicity rate
The heightened seismicity rate in hydraulically isolated fault scenarios demonstrates that the seismicity rate is very sensitive to the stress change rate.This sensitivity can also be observed by comparing seismicity rate within the fault after stage 1 in all the scenarios, which is purely attributed to poroelastic stressing.Depending on the poroelastic stressing influenced by different stimulated widths, the post-stage 1 seismicity rate within the fault is only 10 0.33 = 2.1 for a 0 m stimulated width, and 10 0.55 = 3.5 for a 100 m stimulated width, but could reach up to 10 1.80 = 63 for a 200 m stimulated width (Figure 16 c, d).

Mechanism of the post shut-in M L 2.9 earthquake
Three plausible mechanisms examined for the occurrence of the post shut-in M L 2.9 earthquake at the PNR site include the post shut-in pore pressure diffusion, poroelastic stressing on a non- Therefore, none of three mechanisms could be ruled out in terms of the Coulomb stress change value.
Nevertheless, the focus lies in whether the fault slippage criterion is satisfied after injection at the PNR field conditions, equivalently, whether the maximum Coulomb failure stress change within the fault occurs after injection at the PNR field conditions.Under constant injection rate conditions, following shut-in an unfavourably oriented, hydraulically connected fault may experience rapidly increased normal and shear stress changes, both contributing to the destabilisation of the fault, before the pore pressure declines (Segall & Lu, 2015).This would result in a rapid increase of Coulomb stress, and thus a post-injection spike in seismicity rate.Even after the pore pressure begins to decline after shut-in, the combined action of rapid poroelastic stress changes and the delayed response of pore pressure may cause a post shut-in peak in Coulomb stress.In this case, the fault transmissivity plays a crucial role on the peak magnitude and duration of the post-injection increase in Coulomb stress (Wassing et al., 2021).A conductive fault, with a fast post-injection pore pressure decline, tends to Seismicity rate allows more straightforward comparison between recorded and modelled event counts.
Examination of seismicity rate indicates that all the three mechanisms could result in heightened seismicity rates over the majority of the PNR-2 fault plane (Figure 13).However, the seismicity rate estimated for hydraulically isolated fault scenarios is mostly below 10 4 , while that for hydraulically connected fault scenarios can reach above 10 6 .The field observation of the surge in event count surrounding the fault plane indicates that the poroelastic stressing on a non-overpressurised fault mechanism is less favourable.
The timing of post shut-in fault slippage differs for different mechanisms.For the poroelastic stressing on a non-overpressurised fault mechanism, the maximum Coulomb stress change occurs instantaneously when injection ceases, followed by a monotonical decline.For the post shut-in pore pressure diffusion mechanism, fault slippage usually happens sometime after the end of injection, depending on the permeability of the hydraulically connected fault.For the poroelastic stressing on an overpressurised fault mechanism, the delayed occurrence of fault slippage is also possible.The time after shut-in depends on the relative decline rate of pore pressure and normal stress, again modulated by the fault permeability (Wassing et al., 2021).Field records at the PNR site showed that the seismic magnitude began to increase around 5 hours after the end of injection.Before the M L 2.9 event occurring over 60 hours post shut-in, there was an M L 1.1 event around 9 hours post shut-in, an M L 0.5 event around 14 hours post shut-in, and an M L 2.1 event around 33 hours post shut-in.These observations indicate that the triggering mechanism was active over a prolonged period, thus at least the poroelastic stressing on a non-overpressurised fault mechanism could be ruled out.
Comparison between the field-derived rupture area and regions with positive modelled Coulomb stress change provides useful constraints on the underlying triggering mechanism.For the PNR site, of particular interest is the distribution of recorded seismicity around the poroelastic clamped fault section: recorded seismicity within this fault section would suggest that the PNR-2 fault is conductive and hydraulically connected, and seismic quiescence would suggest the opposite.The PNR-2 fault delineated by aftershocks is elliptically halo-shaped, within which seismicity is quiescent.The poroelastic clamped fault section in the models falls outside the halo and is not seismic active (see Figure 7 of Kettlety and Verdon, 2021).This suggests that the PNR-2 fault is likely to be partially sealing, which allows pore pressure to diffuse but at a slow rate, so that the fault slippage is promoted by gradually elevated pore pressure but the pore pressure is not sufficiently large to activate the poroelastic clamped fault section.
It is noteworthy that the stimulated width and fault permeability used in the models may not accurately represent the field conditions, but provide reasonable upper and lower bounds for extreme scenarios of hydromechanical behaviour.Building upon these modelling results and analyses, it is proposed that the occurrence of the post shut-in M L 2.9 earthquake was a three-staged process: hydraulic fracturing operations first stimulated surrounding reservoir formations and propagated fracture tips along the maximum horizontal principal stress orientation.Fracture tips then reached and established hydraulic connection with the partially sealing PNR-2 fault, followed by gradual pore pressure diffusion to the fault through stimulated regions.After the injection ceased, the pore pressure was significantly lowered, but it remained higher than the in-situ pressure and continued to drive the diffusion across the majority of the fault plane, eventually triggering the fault slippage and the earthquake.In view of this mechanism, continuous microseismic and hydrogeological monitoring is recommended over a prolonged post shut-in period.In case of continuously increasing seismic magnitude of post shut-in events, flowback of injected fluids might be performed to lower overpressure and prevent the delayed occurrence of large induced earthquakes.

Conclusions
A 3D fully coupled poroelastic model was developed to simulate the hydromechanical response of the shale reservoir formation embedded with the 330 m long, 250 m high PNR-2 fault during and after a one-week period of hydraulic fracturing operations in August 2019 at the PNR site, UK.Based on the stress and pore pressure modelled, the slippage potential of the PNR-2 fault responsible for the post shut-in M L 2.9 earthquake was evaluated in terms of the Coulomb failure stress change and seismicity rate.A total of six model scenarios, considering various fault permeabilities and hydraulic connection between injection regions and faults, were modelled to identify the causal mechanism amongst three hypotheses, i.e., the post shut-in pore pressure diffusion, poroelastic stressing on a nonoverpressurised fault, and poroelastic stressing on an overpressurised fault.
Coupled modelling results have shown that increased pore pressure plays a dominant role in triggering the fault slippage, although the poroelastic stress may have acted to promote the slippage.
Amongst the three plausible mechanisms, the post shut-in pore pressure diffusion is the most favoured in terms of Coulomb stress change, seismicity rate, timing of fault slippage and rupture area.
Comparison between various model scenarios has indicated that the occurrence of the post shut-in M L 2.9 earthquake was a three-staged process, where hydraulic fractures first stimulated surrounding reservoir formations, then hydraulically connected to the partially-sealing PNR-2 fault that allowed gradual pore pressure diffusion, and eventually the fault was activated primarily under the direct increase in pore pressure.Model results also highlighted the paramount role of the fault permeability and its connectivity to injection regions in promoting fault rupture, in addition to the fault orientation with respect to the ambient stress field.Co-seismic activation of faults of the same orientation may be attributed to different triggering mechanisms in different hydrogeological settings and stimulation conditions.

Figure 1 .
Figure 1.Map view of PNR microseismic event locations, focused on the PNR-2 events (Kettlety & Verdon, 2021).The two well paths are shown by black lines, and sleeve locations by yellow diamonds.Events are sized by magnitude and coloured by clusters.Coordinates used are based on the Ordnance Survey United Kingdom grid system.

A
3D fully coupled poroelastic model, considering poroelastic solid deformation, fluid flow in both porous rocks and fault structures, and hydraulic fracture propagation, was developed to model the hydromechanical behaviour of reservoir formations during and after the PNR-2 operations, and to further evaluate the potential for earthquakes on the PNR-2 fault.Section 3.1 introduces the mathematical formulation of the coupled poroelastic reservoir model.Section 3.2 presents the development of the coupled model used to simulate the fluid injection-induced hydromechanical behaviour of the shale reservoir.Based on the model results, the evaluation of potential for seismicity in terms of the Coulomb stress change and seismicity rate is described in Sections 3.3 and 3.4, respectively.

Figure 2 .
Figure 2. 3D model geometry for hydraulic fracturing operations at the Preston New Road, UK: (a) 3D view, (b) plan view, and (c) side view.

Figure 3 .
Figure 3. Injection pressure and hydraulic fracture growth histories used to simulate the hydraulic fracturing operations at the Preston New Road, UK: (a) bottomhole pressure history for both NS and Eastern Zones, and (b) hydraulic fracture length for northward-and southwardpropagating segments of the NS and Eastern Zones.

Table 1
Six model scenarios with various fault permeabilities k and stimulated widths w k n and  s are normal and shear stress changes resolved on the fracture plane, respectively.Here, negative normal stress changes σ n indicate rock compression.The potential for fracture slippage is elevated for a positive Coulomb failure stress change, and suppressed for a negative value.To isolate respective contributions of poroelastic stressing and pore pressure change, Equation (9) can be re-arranged in terms of poroelastic stress change and pore most vulnerable to rupture is 45°-φ/2 (φ is the internal friction angle of reservoir rocks given by f = tanφ) off the maximum principal stress direction, according to the Mohr-Coulomb failure criterion.Considering a friction coefficient f = 0.6(Verdon et al., 2020) and the maximum principal stress orientation of 173°N at the PNR site, the most vulnerable fracture plane is orientated at 144.5°N or 202.5°N.The Coulomb failure stress change  in response to hydraulic fracturing was resolved on the PNR-2 fault plane with the fault strike 130°N, which is well oriented to rupture.

4. 1
Uncoupled, one-way coupled and fully-coupled models Before investigating the full problem involving different model scenarios, we illustrate the effect of poroelastic coupling on the reservoir behaviour by comparing results of the baseline model scenario from an uncoupled model, a one-way fluid-to-solid coupled model, and a two-way poroelastic coupled model.

Figure 4
Figure4 (a-c) presents the pore pressure distribution immediately after injection stage 7 at 2,100 m depth of the PNR-2 well.The uncoupled and one-way coupled models yield the same overpressurised regions, constrained within the stimulated regions and the PNR-2 fault.The fully-coupled model presents a slightly smaller pore pressure increase within the stimulated regions.This is because the rock compression by the overpressure increases the volume fraction to host fluids and acts as a liquid source (Equation (2)), which causes less fluid injected under controlled injection pressure conditions, and thus lower overpressure within the stimulated regions.The uncoupled model does not represent the solid-to-fluid coupling effect, and thus slightly overestimates the overpressure.In contrast, under controlled injection rate conditions, the poroelastic effect would cause larger overpressure according to Equation (2)(Chang & Segall, 2016).

Figure 4 .
Figure 4. Pore pressure change p and mean normal stress change σ kk /3 immediately after injection stage 07 at 2,100 m depth of the PNR-2 well: (a, d) uncoupled, (b, e) one-way fluid-to-solid coupled, and (c, f) two-way poroelastic coupled models.The fault is conductive and hydraulically connected to the Eastern Zone (stimulated width 100 m, fault permeability 100 mD).

Figure 5
Figure 5 presents the pore pressure change p after injection stages 3 and 7 for the six model scenarios.The increased pore pressure only distributes within stimulated regions in all the scenarios.Whilst the PNR-2 fault is not overpressurised throughout the hydraulic fracturing operations in

Figure 5 .
Figure 5. Pore pressure change p immediately after injection stages 03 and 07 at 2,100 m depth of the PNR-2 well: (a, b) hydraulically isolated conductive fault, (c, d) hydraulically isolated sealing fault, (e, f) conductive fault hydraulically connected to the Eastern Zone, (g, h) sealing fault hydraulically connected to the Eastern Zone, (i, j) conductive fault hydraulically connected to both the NS and Eastern Zones, and (k, l) sealing fault hydraulically connected to both the NS and Eastern Zones.

Figure 6 and
Figure 6 and Figure 7 present the normal stress change σ n and shear stress change τ s after injection stages 03 and 07 for the six model scenarios.Displacement vectors are also indicated using the same length scale in the graphs.Fractures oriented in the fault direction are clamped within the stimulated regions, and relieved immediately surrounding these regions.Farther away from the stimulated

Figure 6 .
Figure 6.Normal stress change σ n immediately after injection stages 03 and 07 at 2,100 m depth of the PNR-2 well: (a, b) hydraulically isolated conductive fault, (c, d) hydraulically isolated sealing fault, (e, f) conductive fault hydraulically connected to the Eastern Zone, (g, h) sealing fault hydraulically connected to the Eastern Zone, (i, j) conductive fault hydraulically connected to both the NS and Eastern Zones, and (k, l) sealing fault hydraulically connected to both the NS and Eastern Zones.Grey arrows indicate displacement vectors.

Figure 7 .
Figure 7. Shear stress change τ s immediately after injection stages 03 and 07 at 2,100 m depth of the PNR-2 well: (a, b) hydraulically isolated conductive fault, (c, d) hydraulically isolated sealing fault, (e, f) conductive fault hydraulically connected to the Eastern Zone, (g, h) sealing fault hydraulically connected to the Eastern Zone, (i, j) conductive fault hydraulically connected to both the NS and Eastern Zones, and (k, l) sealing fault hydraulically connected to both the NS and Eastern Zones.Grey arrows indicate displacement vectors.

Figure 8 .
Figure 8. Coulomb failure stress change due to poroelastic stressing τ s + fσ n immediately after injection stages 03 and 07 at 2,100 m depth of the PNR-2 well: (a, b) hydraulically isolated conductive fault, (c, d) hydraulically isolated sealing fault, (e, f) conductive fault hydraulically connected to the Eastern Zone, (g, h) sealing fault hydraulically connected to the Eastern Zone, (i, j) conductive fault hydraulically connected to both the NS and Eastern Zones, and (k, l) sealing fault hydraulically connected to both the NS and Eastern Zones.

Figure 9
Figure 9 presents the Coulomb failure stress change τ after injection stages 3 and 7 for the six model scenarios.The Coulomb failure stress change is dominated by pore pressure change within stimulated regions and the hydraulically connected conductive fault, albeit being restricted by the normal stress change.Outside the stimulated regions, Coulomb failure stress change is primarily contributed by the shear stress change, where the potential for fault slippage is suppressed to the sides of the NS and Eastern zones but promoted at the propagation fronts.In isolated or sealing fault scenarios, the northwest end of the PNR-2 fault is suppressed to slip; but once hydraulic connection is established, the increased pore pressure overwhelms in favour of fault slippage (Figure 9 f, i, j).

Figure 9 .
Figure 9. Coulomb failure stress change τ immediately after injection stages 03 and 07 at 2,100 m depth of the PNR-2 well: (a, b) hydraulically isolated conductive fault, (c, d) hydraulically isolated sealing fault, (e, f) conductive fault hydraulically connected to the Eastern Zone, (g, h) sealing fault hydraulically connected to the Eastern Zone, (i, j) conductive fault hydraulically connected to both the NS and Eastern Zones, and (k, l) sealing fault hydraulically connected to both the NS and Eastern Zones.

Figure 10 Figure 10 .
Figure 10 presents the Coulomb failure stress change τ resolved on both hydraulic fracture zones and the PNR-2 fault immediately after each injection stage and at the end of modelling for the six model scenarios.The PNR-2 fault is not stimulated in the isolated fault scenarios (Figure 10 a and b), which provides a unique case to examine the effect of poroelastic stressing.Since the fluid injection into the NS Zone approaches the PNR-2 fault (at injection stage 02), the northwest end of the fault to the side of the NS Zone is clamped by the increased shear stress.As the Eastern Zone is stimulated (at injection stage 5), the northwest end of the fault in between the NS and Eastern Zones is further suppressed, whilst the front of Eastern Zone with elevated shear stress impinges the central part of the fault, promoting the potential for slippage.

Eastern
Zone has a larger influence on Coulomb failure stress changes than the NS Zone (Figure 11 a, b).If the PNR-2 fault is hydraulically connected to the hydraulic fracture zones, the Coulomb failure stress change τ dominated by the pore pressure change could be an order of magnitude larger than that due to poroelastic stressing, either stabilising or destabilising.In hydraulically connected sealing fault scenarios, the southeast end of the fault is less influenced by the increased pore pressure (Figure 11 d, f), as compared to conductive fault scenarios (Figure 11 c, e).Depending on the overpressure distribution, the extent of the seismicity suppressed fault section is around 150 m for hydraulically isolated fault scenarios (Figure 11 a and b), but could reach up to 200 m for hydraulically connected sealing fault scenarios (Figure 11 d and f).The magnitude of Coulomb failure stress change is also the minimum for the hydraulically isolated fault scenarios (Figure 11 a and b).

Figure 11 .
Figure 11.Evolution of Coulomb failure stress change τ along the PNR-2 fault plane (the dashed purple line A-A in Figure 9a): (a) hydraulically isolated conductive fault, (b) hydraulically isolated sealing fault, (c) conductive fault hydraulically connected to the Eastern Zone, (d) sealing fault hydraulically connected to the Eastern Zone, (e) conductive fault hydraulically connected to both the NS and Eastern Zones, and (f) sealing fault hydraulically connected to both the NS and Eastern Zones.
change before being hydraulically connected to the Eastern Zone.Point C is located within the clamped northwest end of the fault.Figure12presents the stress components at the two measurement points for three different hydraulic connection scenarios.Under injection pressurecontrolled conditions, no prominent distinction in hydromechanical behaviour was observed between a conductive fault and a sealing fault.In particular, shear stresses in conductive and sealing fault scenarios, respectively marked by yellow solid and dashed lines, are the same in all the graphs, because shear stress is independent of pore pressure change.Notably, when hydraulically connected, a sealing fault has sharper response to fluid injection than a conductive fault, in terms of both the injection-induced increase and the post-injection decrease in pore pressure change p and Coulomb failure stress change τ (Figure 12 c, d, e, f).

Figure 12 .
Figure 12.Stress changes at measurement points on the PNR-2 fault: (a) point B (hydraulically isolated conductive fault), (b) point B (fault hydraulically connected to the Eastern Zone), (c) point B (fault hydraulically connected to both the NS and Eastern Zones), (d) point C (hydraulically isolated conductive fault) , (e) point C (fault hydraulically connected to the Eastern Zone), and (f) point C (fault hydraulically connected to both the NS and Eastern Zones).Locations of measurement points B and C are shown in Figure 9 (a).

Figure 13
Figure 13 presents the seismicity rate R resolved on the PNR-2 fault as well as the hydraulic fracture zones immediately after each injection stage and at the end of modelling for the six model scenarios.Heightened seismic levels are observed in regions with a positive Coulomb stress change τ, as shown Figure 10.Pore pressure change results in much larger seismicity rates than the poroelastic stress change.In hydraulically isolated scenarios with the poroelastic stressing effect alone, the seismicity rate is mostly limited below 10 4 .In contrast, the seismicity rate can reach up to 10 7 in hydraulically connected scenarios where pore pressure change dominates.Due to the quadratic relation between the seismicity rate and its change rate in Equation (11), the Coulomb stress change rate   has a more pronounced effect on the seismicity rate R than on the Coulomb stress change τ, in particular at high seismicity rates.When the PNR-2 fault is hydraulically connected to the hydraulic fracture zones, the Coulomb stress change τ increases by an order of magnitude, but the seismicity rate dramatically increases by over 4 orders of magnitude following the surge in the Coulomb stress change rate   (Figure 10 and Figure 13 c, d, e, f).

Figure 14
Figure 14 presents the seismicity rate evolution along the full fault length over the hydraulic fracturing operation.The temporal evolution of seismicity rate is closely associated with that of the Coulomb stress change shown in Figure 11.Interestingly, although the seismicity rate R surges following a rapid increase in Coulomb stress change rate   , it does not fade off as fast following a rapid decline in Coulomb stress change rate   .Each injection stage represents a high Coulomb stress change rate, bringing the seismicity rate to a peak.The Coulomb stress change rate dramatically drops immediately after each injection stage, and the PNR-2 fault is characterised by a steady Coulomb stress only influenced by the pore pressure diffusion process.However, the post-injection

Figure 13 .
Figure 13.Seismicity rate R on hydraulic fracture planes and the PNR-2 fault plane after each injection stage of the PNR-2 well: (a) hydraulically isolated conductive fault, (b) hydraulically isolated sealing fault, (c) conductive fault hydraulically connected to the Eastern Zone, (d) sealing fault hydraulically connected to the Eastern Zone, (e) conductive fault hydraulically connected to both the NS and Eastern Zones, and (f) sealing fault hydraulically connected to both the NS and Eastern Zones.
pressure during fracturing would be much smoother than in wellhead pressure, so as the Coulomb stress change rate closely associated with the seismicity rate.

Figure 14 .
Figure 14.Evolution of seismicity rate R along the PNR-2 fault plane (the dashed purple line A-A in Figure 9a): (a) hydraulically isolated conductive fault, (b) hydraulically isolated sealing fault, (c) conductive fault hydraulically connected to the Eastern Zone, (d) sealing fault hydraulically connected to the Eastern Zone, (e) conductive fault hydraulically connected to both the NS and Eastern Zones, and (f) sealing fault hydraulically connected to both the NS and Eastern Zones.

Figure 15 .
Figure 15.The recorded and modelled cumulative seismic event counts over the hydraulic fracturing operation period.The PNR-2 fault is conductive and hydraulically connected to the Eastern Zone in the model.

1
Role of poroelastic stressing in triggering post shut-in earthquakesWe proceeded to isolate the contribution of pore pressure change on the Coulomb stress change on the PNR-2 fault.Figure16(a)(b) presents the mean Coulomb stress change  and the contribution from pore pressure change fp resolved on the PNR-2 fault for the six model scenarios.When both the pore pressure and poroelastic stressing are at play within the PNR-2 fault, such as in hydraulically connected conductive fault scenarios (Figure 16 a), the Coulomb stress change  is lower than the pore pressure change fp.When only the poroelastic stressing effect is active within the PNR-2 fault, such as in hydraulically isolated scenarios and sealing fault scenarios (Figure 16 b), the opposite is true.

Figure 16 .
Figure 16.The contribution of pore pressure change to mean Coulomb stress change  and seismicity rate R resolved on the PNR-2 fault.Mean Coulomb failure stress change  and the contribution from pore pressure change fp resolved on (a) a conductive fault, and (b) a sealing fault.Mean seismicity rate R calculated based on Coulomb failure stress change rate   and pore pressure change rate p  on (c) a conductive fault, and (d) a sealing fault.

Figure 16 (
Figure 16 (c)(d) presents the mean seismicity rate R computed based on Coulomb failure stress change rate   and pore pressure change rate p  on the PNR-2 fault for the six model scenarios.If the PNR-2 fault is conductive and hydraulically connected, the seismicity rate computed based on the Coulomb stress change rate   is slightly smaller than that based on pore pressure change rate fp  , overpressurised fault, and poroelastic stressing on an overpressurised fault.The first mechanism is represented by the model scenarios (c)(d), the second by the model scenarios (a)(b), and the last by the model scenarios (e)(f).We examined the possibility of the three mechanisms in terms of four factors based on modelling results: (1) Coulomb stress change, (2) seismicity rate, (3) timing of fault slippage, and (4) rupture area.A Coulomb failure stress change in excess of the generalised triggering threshold of 0.01 -0.1 MPa is considered to have high potential to trigger fault slippage(Kettlety & Verdon, 2021;Shapiro et al., 1997).Although Coulomb failure stress changes in the hydraulically isolated fault scenarios are much lower, they could reach well above 1 MPa ahead of fracture tips after stages 6 and 7 (Figure12 a).
have a small and narrow peak in Coulomb stress.In contrast, a sealing fault, with a slow postinjection pore pressure decline, could have a prominent and prolonged increase in Coulomb stress.It is noteworthy that the post-injection increase in Coulomb stress does not necessarily emerge across the full fault plane, but occur in a localised fault section.At the PNR hydraulic fracturing site, when the PNR-2 fault is hydraulically connected to the hydraulic fracture zones, its pore pressure change can be well constrained by the field-recorded wellhead pressure history.Under such conditions, the Coulomb failure stress change τ generated keeps decreasing after injection stage 7 at both points B and C (Figure12 c, d, e, f), suggesting that the release of the normal stress is not rapid enough to cause a post-injection peak in the Coulomb stress.This indicates that although the poroelastic stressing mechanism is active, it does not play a governing role in triggering the post shut-in M L 2.9 earthquake.
Segall and Lu (2015)re-formulated the seismicity rate framework by eliminating the state variable and expressing the equation in terms of the seismicity rate relative to the background rate R: Implementation of the model to the nucleation of accelerating slip on faults with the rate-and-state friction law yields a state-variable constitutive formulation of seismicity rate associated with the applied stressing history.