3D seismic interpretation and fault slip potential analysis from hydraulic fracturing in the Bowland Shale, UK

The Bowland Shale Formation is one of the most promising targets for unconventional exploration in the United Kingdom, with estimated resources large enough to supply the country's entire natural gas consumption for 50 years. However, development of the Bowland Shale has stalled due to concerns over hydraulic-fracturing-induced seismicity. Only three wells have been drilled and hydraulic-fractured to date in the Bowland Shale, and all three have produced levels of seismicity of sufficient magnitude to be felt at the surface. Susceptibility to induced seismicity will be determined by the presence of critically stressed faults. However, such faults can go undetected in conventional interpretation of 2D or 3D seismic surveys if they are shorter than the resolution retrievable from a seismic survey, or if they have low (and in some cases even zero) vertical displacement. In such cases, the faults that cause induced seismicity may only be visible via microseismic observations once they are reactivated. To better identify fault planes from 3D seismic images, and their reactivation potential due to hydraulic fracturing, a high-resolution fault-detection attribute was tested in a 3D seismic survey that was acquired over the Preston New Road site, where two shale-gas wells were hydraulic-fractured in the Bowland Shale in 2018 and 2019, obtaining fault planes with lengths between 400 and 1500 m. Fault slip potential was then estimated by integrating the obtained faults with the formation's stress and pore pressure conditions (with the Bowland shale also being significantly overpressured), and several critically stressed faults were identified near the previously hydraulic fractured wells. However, the faults that induced the largest seismic events in the Preston New Road site, of c. 200 m in length for seismic events of magnitudes below 3.0 (as imaged with a multicomponent, downhole microseismic monitoring array deployed during the hydraulic-fracturing stimulations), could not be identified in the 3D seismic survey, which only mapped fault planes larger than 400 m in length.

The Bowland Shale Formation is one of the most promising targets for unconventional exploration in the UK because of its high total organic content (typically 1-3%, but can reach up to 8%), porosity (4-7%) and brittleness (>30% of silica content) (Smith et al. 2010;Andrews 2013). Located at depths from 2000 to 4750 m, stretching across central Great Britain (Fig. 1a), the Bowland Shale is composed of dark grey shale with interbedded limestone, siltstone, calcareous mudstone and sandstone (Andrews 2013). The study area is located in the NW England's Bowland Basin, which started to develop during the Late Devonian to Carboniferous by the rifting in the north-south to NW-SE direction (Anderson and Underhill 2020). After that, this area was influenced by various geological processes including the post-rift thermal subsidence in the Pennsylvanian and the Variscan Orogeny in the Early Permian, leaving the area heavily deformed (Anderson and Underhill 2020). The stratigraphic column of the Bowland Basin is shown in Figure 2, where the Bowland Shale is depicted sitting on top of the Worston Shale Group, and below the Millstone Grit Group. The Bowland Shale been estimated by the British Geological Survey (BGS) to potentially host P90, P50 and P10 resources of 822, 1329 and 2281 trillion cubic feet (23.3, 37.6 and 64.6 trillion cubic meters) of natural gas respectively (Andrews 2013). A 10% recovery rate of the P50 value translates to almost 50 years of UK gas consumption (Andrews 2013).
Despite the potential size of this resource, development of the Bowland Shale has stalled due to concerns over hydraulicfracturing-induced seismicity (HF-IS). The UK's Department of Energy and Climate Change (DECC) first imposed a moratorium in 2011 in response to a magnitude M L 2.3 earthquake at the Preese Hall well (Fig. 1c) (Clarke et al. 2014). This moratorium was later lifted by DECC in December 2012 (OGA 2019), with the implementation of a strict traffic light scheme (TLS) with a redlight threshold set at M L 0.5 to suspend injection (Kendall et al. 2019). However, further hydraulic-fracturing at the Preston New Road (PNR) wellsite ( Fig. 1c) produced an M L 1.5 earthquake in 2018 (Clarke et al. 2019a) and an M L 2.9 earthquake in 2019 . The fault planes which are believed to cause both M L 1.5-and M L 2.9-induced seismicity at the PNR site were identified by mapping the seismic event (M L >0) locations during an injection hiatus (Clarke et al. 2019a) and utilizing the locations of the aftershock events to delineate the fault plane . Also, the focal mechanisms of the largest induced earthquakes exhibit a strike-slip mechanism with nearly vertical nodal planes, making these fault planes difficult to detect in reflection seismic datasets as no strike-slip faults typically have small (and in some cases null) vertical offsets. In response to these anomalously large seismic events, some of which were also felt by local communities, the UK government has imposed a further indefinite moratorium on shale gas development. Hence, there is a clear need to better manage the risks posed by induced seismicity if further development of this resource is to take place.
Induced seismicity takes place when subsurface activities, such as hydraulic fracturing, create perturbations in the pore pressure and/ or stress field in the surrounding rocks. If a pre-existing fault is present that is near to its critical stress state, then these perturbations may be sufficient to trigger slip, resulting in induced seismicity. The magnitude of induced earthquakes will be determined by the dimensions of the resulting slip; as such, induced seismic magnitudes may be influenced by the scale of the perturbation (e.g. Shapiro et al. 2011), by the size of faults present in the vicinity of the stimulated well, or by the smoothness (or planarity) of the fault planes, as smoother ('simpler') faults produce bigger earthquakes (Wesnousky 1988). For HF-IS to occur, several necessary conditions must therefore be met: a fault of sufficient size must be present in the subsurface in the vicinity of the stimulated well, the in situ stress conditions must be such that the fault is close to its critical point, and the perturbation created by the stimulation must be of sufficient size to move the conditions on the fault to failure.
In a broad sense, there are two methods by which the risks of HF-IS could be mitigated: real-time modification of injection schemes, and pre-operation site selection based on mapping of subsurface faults. Real-time HF-IS mitigation involves decision-making based on observed microseismicity during operations. Perhaps the simplest form of real-time HF-IS decision-making is the TLS, where injection rates are reduced or stopped depending on the magnitudes of induced events (e.g. Verdon and Bommer 2021). Alternatively, statistical forecasting of expected magnitudes can also be used to guide decision-making during stimulation (e.g. Verdon and Budge 2018;Clarke et al. 2019a;Kwiatek et al. 2019;).
As described above, for HF-IS to occur, a pre-existing fault that is critically stressed must be present in the vicinity of the well. Hence, if we were able to obtain a complete understanding of the position and orientation of faults in the subsurface, then HF-IS could, in theory, be mitigated by avoiding sites where critically stressed faults are present (i.e. faults that are optimally oriented relative to the maximum principal stress, to have the maximum shear stress acting upon them, which could cause these faults to slip if the shear stress exceeds the faults' friction; discussed in more detail in the Fault Slip Potential section and in Fig. 11 below. It is this possibility that we investigate in this study, with a focus on the HF-IS observed at the PNR site in 2018-2019 (Clarke et al. 2019a;. This site offers several high-quality, publicly available datasets, including high-resolution microseismic data acquired during stimulation, 3D reflection seismic data and geophysical well logs. This study aims to compare the fault planes identified using microseismic observations (i.e. location of microseismic events, and focal planes of higher-magnitude events) as those responsible for the PNR HF-IS (Fig. 3), with structures identified in interpretation of 3D seismic data acquired at the site. Our interpretation of the 3D Fig. 1. Bowland Shale area in central Great Britain (a) with the 2D and 3D seismic exploration surveys (b) acquired for exploration of conventional oil and gas fields, coal and coalbed methane, and more recently for unconventional shale gas. The top view of the Bowland-12 3D reflection seismic data near Blackpool (c) shows the location of the wells hydraulic-fractured (to date) in the Bowland Shale (first Preese Hall 1 in 2011, followed by the wells in Preston New Road (PNR) in 2018 and 2019), and a time slice of the 3D seismic data at 1260 ms ( just below the two PNR wells, with their depths converted to two-way-travel time), with a distance of 25 m between in-lines and cross-lines. The microseismic events observed during the hydraulic stimulations of the PNR wells are shown in detail in Figures 3 and 5, and the structural interpretation of the vertical cross-sections are shown in Figure 4. seismic begins with well-to-seismic calibration, followed by manual horizon and fault picking. We follow this with the application of seismic attributes, including similarity, spectral decomposition, curvature, and an automated fault detection attribute. While these methods identify numerous faults around the PNR site, they struggle to identify the faults responsible for the HF-IS, highlighting the challenges of basing HF-IS mitigation on pre-operational site selection.

PNR site description and datasets
The Bowland Shale was deposited in a marine environment with relative sea level fluctuations, causing alternating changes in lithology between deltaic and deep marine facies (Andrews 2013). Its mineral content consists of 56-59% clay, 45% quartz and 10% carbonate (Smith et al. 2010). It was deposited in a very tectonically active environment: starting from the formation of the Bowland Basin during the Devonian and Carboniferous periods, thermal subsidence during the Pennsylvanian (Namurian and Westphalian stages), and exhumation from the Variscan Orogeny in the Late Carboniferous and Permian periods (see the stratigraphic column in Fig. 2) (Anderson and Underhill 2020). This tectonic complexity can be seen in 3D seismic data.
To date, activities in the Bowland Shale have been focused on the Fylde Peninsula, Lancashire, where stimulation has taken place at the Preese Hall and PNR sites near Blackpool (Fig. 1b). This area is covered by a 3D reflection seismic dataset, pre-stack depth migrated, of around 10 by 10 km (Fig. 1c). The dataset was acquired in 2012 after the first felt seismic events associated with hydraulic fracturing were reported in 2011 at the Preese Hall 1 well (Clarke et al. 2014). One fault near the Preese Hall 1 well was observed in this 3D seismic dataset, with very similar location and orientation as one nodal plane of one of the largest-magnitude events of M L 2.3 registered during the well's hydraulic stimulation (Clarke et al. 2014). This suggested that the 3D seismic data would be capable of detecting more critically stressed faults that could be reactivated during hydraulic stimulation of future shale-gas wells.
The location of the PNR wells, around 4 km south of the Preese Hall 1 well, was chosen in part to avoid major faults identified in the same 3D seismic dataset (Cuadrilla Bowland Ltd 2019). Smaller seismic discontinuities (SDs) observed closer to the PNR wells were also considered as potential faults, although the closest SD was 300 m away from the closest injection point and therefore considered to have a low risk of fault reactivation during hydraulic fracturing.
In the autumn of 2018, the hydraulic fracturing stimulation of the first PNR well, PNR-1z, produced multiple 'red-light' events (i.e. seismic events with magnitude above M L 0.5, which forced the operator to suspend the injection operations temporarily). Operations resumed with the operator skipping some stages to avoid a further reactivation of faults identified with microseismic observations (Clarke et al. 2019b). The maximum magnitude reported during the stimulation of the PNR-1z was M L 1.5, with the BGS reporting an intensity of 2 (i.e. scarcely felt, according to the European Macroseismic Scale -EMS). The event was felt by some members of the public living near to the well site, though the event was not strong enough to cause any disturbance or damage to the nearby communities or infrastructure.
The hydraulic fracturing stimulation of the second PNR well, PNR-2, in August 2019, also produced multiple red-light events, with a maximum magnitude 2.9 and an intensity level 6 on the same EMS scale (i.e. strong ground motions and possible minor damages to ordinary buildings). This event occurred after only seven hydraulic-fracturing stages (of 45 planned; Cuadrilla Bowland Ltd 2019), after which the regulator ordered a shutdown of injection operations, and later imposed a further moratorium on shale gas hydraulic fracturing in the UK that remains in place today.
The locations and geometries of the faults that were reactivated during stimulation of both PNR-1z and PNR-2 ( Fig. 2) have been constrained from observations of microseismic event locations and source mechanisms (Clarke et al. 2019a;. However, these structures were not identified in interpretations of the 3D seismic dataset that were produced prior to the hydraulic-fracturing operations. The local stress regime at the PNR site has been described by Clarke et al. (2019b). The stress gradients for the Bowland Shale, located at an average depth of 2200 m for the PNR wells, are S HMAX = 1.4 psi/ft, S V = 1.114 psi/ft, and S hmin = 0.725 psi/ft (31.7, 25, and 16.4 Mpa km -1 respectively), and assuming Andersonian stress conditions (i.e. the vertical stress, S V, is a principal stress, and the other two are horizontal; Anderson 1951). The vertical stress gradient (S V ) was calculated from density logs from the vertical wells in the area ( particularly from the well Preese Hall 1, shown in Fig. 1c), and the gradient of the minimum horizontal stress (S hmin ) was obtained from the fracture closure pressures measured from diagnostic fracture injection tests done at different depths within the Bowland shale (also from the Preese Hall 1 well). The gradient of the maximum horizontal stress (S HMAX ) has a higher uncertainty as it can only be constrained from the rock's internal strength, and it requires further measurements of the rock's internal friction (which can be highly variable in a formation as heterogeneous as a shale), and the implementation of a frictional faulting theory (as the Mohr-Coulomb faulting envelope, discussed in detail in the Fault Slip Potential section and in Fig. 11 below). The S HMAX orientation for the studied area, measured from the breakouts observed in image logs from the Preese Hall-1 well, is c. 173°. The pore pressure gradient is 0.58 psi/ft (13 Mpa km -1 ), which is significantly higher than a hydrostatic pressure of 0.45 psi/ft (10 Mpa km -1 ), meaning that the Bowland Shale is highly overpressured. Source mechanisms for the microseismic events showed that, for both PNR-1z and PNR-2, the largest induced events had strike-slip focal mechanisms.

Fault size and earthquake magnitude
The magnitudes of induced earthquakes will be limited by the size of the reactivated fault. Earthquake moment, M O is defined by where A is the rupture area, D is the average slip and G is the rock shear modulus. The resulting moment magnitude, M W , is given by: Equation (1) can be reformulated in terms of the earthquake stress drop, Δσ (Kanamori and Brodsky 2004), which is observed to be scale invariant, typically ranging between 0.01 and 100 MPa (1.45-14 500 psi) (e.g. Abercrombie and Leary 1993;Bohnhoff et al. 2016), such that  100 MPa, in Log 10 scale), and a fault aspect ratio of 1 (i.e. a square rupture plane).
The geophysical detectability of a fault will primarily depend on its length and on the amount of offset that has accumulated. Faults typically display fractal (i.e. self-similar) scaling behaviour, where fault dimensions scale with fault offset (e.g. Cowie and Scholz 1992). This occurs because both offset and length accumulate concurrently as a fault develops and grows. While this scaling between fault length and offset is observed to vary depending on tectonic settings (e.g. Anderson et al. 1996), Table 1 also lists the expected offset that would be present given a typical relationship between offset and fault length of 1% (e.g. Dawers et al. 1993). Note that these values represent minimum fault size for a given magnitude, since they represent faults rupturing in their entirety whereas, in reality, larger faults may host smaller earthquakes if only a portion of the fault is ruptured. Nevertheless, these values serve to provide an initial constraint on our ability to mitigate HF-IS through the geophysical mapping of pre-existing faults: if typical resolution of good-quality 3D reflection seismic data is of the order of tens of metres, then we might expect faults capable of hosting earthquakes of M 4 or lower to fall below the survey resolution.
However, fault identification in 3D seismic surveys is usually based on mapping vertical displacement of horizons. In the case of the Bowland shale area in central Britain, multiple faults with vertical offsets (i.e. normal and reverse faults) have been interpreted in regional 2D seismic images (Andrews 2013), which has been essential to reconstruct the area's complex tectonic history. However, if the present local stress regime is of a strike-slip nature, then any induced seismicity would be expected to occur on strike-slip faults (as those interpreted from the location and focal mechanisms of the microseismic events recorded during the hydraulic stimulation of the PNR wells; Fig. 3). Since these may not produce significant vertical displacement, even larger strike-slip faults may be missed when interpreting 3D seismic datasets.

3D seismic interpretation
We began our 3D seismic interpretation in this area by manually choosing key horizons from tied wells. Calibrated with the 3D seismic data, Preese Hall-1 (PH-1), Grange Hill-1 (GH-1), and Thistleton-1 (TH-1) wells provide locations of formation tops, which can be used as the starting points for horizon picking. From there, key horizons such as the Manchester Marl Formation, Coal Measure Group, Millstone Grit Group, Upper and Lower Bowland Shale Formation and Worston Shale Group have been chosen (see the stratigraphic column in Fig. 2). In addition, the Variscan Unconformity has also been tracked across the 3D seismic volume. Preliminary fault selection based on detecting discontinuities and offsets of horizons was then used to locate faults. We note that, due to the complicated structural and tectonic setting, manual selection of horizons and faults can be challenging even in these high-quality seismic data. Figure 4 shows the interpretation of the 3D seismic data including faults and horizons on a selection of cross-sections, and a close view of the 3D seismic and the microseismic events recorded during the hydraulic fracturing stimulations of the PNR wells (i.e., close view of the highlighted section in Fig. 4d) is shown in figure 5. The Variscan Unconformity is located below the Manchester Marl Formation and can be seen clearly cutting through the older formations below. The tectonic complexity of the strata underneath the unconformity is noticeably higher than those above, including faulting and folding in the older formations. The 3D seismic data reveal an anticline in the northwestern area of the data, with its folding axis lying in the NE-SW direction. The horizons below the unconformity generally incline towards the NE, terminating against the unconformity. The identified faults can be separated into two populations. The first group consists of thrust faults which strike NE-SW and dip in either NW or SE directions. Most of the thrust faults found in this region are located below the Variscan Unconformity. Larger normal faults are observed running through strata both above and below the unconformity. These faults can reach up to 7.5 km in length. In addition to these faults, we observe large fracture zones (dark, chaotic areas shown by the blue arrows in Fig. 6), where the reflection signals are much more chaotic. These zones are located mostly in the SE of the seismic data cube, with some in the north and NW.
However, the dimensions of the main faults reactivated during the hydraulic stimulations of the PNR wells, of less than 0.5 km in length (interpreted from the microseismicity observed during the stimulations, as shown in Fig. 3), are much smaller than the major faults of several kilometres in length, interpreted from vertical sections of the 3D reflection seismic data (Fig. 4). A close view of the same 3D seismic data near the PNR wells (Fig. 5) shows that the aforementioned reactivated faults at the PNR site are not clearly visible in the volume, at least from the raw seismic amplitude.
To complement this fault interpretation based on raw amplitudes, we tested multiple seismic attributes commonly implemented to enhance the detection of faults and seismic discontinuities in seismic images. In this study, similarity, spectral decomposition, and curvature have been applied to the seismic volume (Figs 6 and 7).
Introduced by Bahorich and Farmer (1995), similarity measures coherence by comparing signals between adjacent gathers using cross-correlation or semblance. A fault may create a difference in the signals between adjacent locations; thus, the similarity attribute would highlight faults as regions with low similarity. In this study, we used the similarity attribute provided by OpendTect (dGB Earth Sciences 2015) with a calculation time gate of [−28, 28] ms. The similarity values can be calculated by the following equation: where sim is the similarity value and X and Y are vectors of length N samples. Spectral decomposition is a method of converting seismic signals from time domain to frequency domain using the continuous wavelet transform (CWT) algorithm.
Given the wavelet c(t), where c is the complex conjugate of c, t is time, t is the time translation and s is the dilation of the wavelet (Sinha et al. 2005). The data in the frequency domain can then be tuned to a specific frequency value to manipulate the interference of signals for highlighting geological features. Different geological features have different tuning behaviours since the interference is influenced by We assume a stress drop of 1 MPa, a fault aspect ratio of 1 (i.e. square), and a fault offset/fault length ratio of 1%.
local distribution of the impedance contrasts and the wavelet. For example, faults with different sizes and locations would tune in or tune out at different frequencies. Therefore, it is commonly utilized for identifying lateral changes or discontinuities on a surface. In practice, the spectral decomposition attribute is commonly used with RGB-colour blending. The attribute volumes are computed using three different tuning frequencies: 15, 30 and 75 Hz (translating to low-, mid-and high-frequency ranges, respectively). These are represented by different colour schemes (red, green and blue in order). Then, they can be overlapped and displayed together as one attribute, which can improve the accuracy of the results and interpretations. Fig. 4. Structural interpretation of the cross-sections from the 3D seismic data shown in Figure 1c. A close view of the 3D reflection seismic data from the red dashed rectangle in (d) with the microseismic events from the hydraulic fracturing stimulations of the PNR wells shown in Figure 3, is shown in detail in Figure 5.  Figures 1 and 3, and the microseismic events observed during the hydraulic-fracturing stimulations of the PNR wells PNR-1z and PNR-2 (also shown in Fig. 3a and c), on top of the same 3D seismic data ((c) and (d) respectively). The yellow horizon shown in the front views ((a) and (c)) corresponds to the interpreted top of the Lower Bowland Shale.
Finally, the curvature attribute indicates the rate of change of dip of a surface. In other words, it is a measure of how deformed or bent a particular point of a surface is. A curvature value (K) of a surface can be calculated by using the following equation: where x and y are the variables in the quadratic equation that represents the surface or curve (Roberts 2001). Positive curvature values represent antiform features, while negative values represent synform features. Thus, the curvature attribute can be used to detect features that offset a surface, or to enhance the relief of geological features. In three dimensions, the curvature is calculated from two perpendicular vertical planes, referred to as normal curvatures, composed of maximum and minimum curvatures. The maximum curvature is a measure of the maximum bending of a surface at the given point, while the minimum curvature measures the curve perpendicular to the maximum curvature. The maximum curvature is typically utilized in fault detection. There are various ways of displaying curvature, including mean curvature, Gaussian curvature, most-positive curvature, most-negative curvature, dip curvature and strike curvature, though most-positive and mostnegative curvatures are the most convenient to use to detect geological features (Chopra and Marfurt 2007). Most-positive and most-negative curvature attributes calculate the most-positive and most-negative values from the normal curvature.
In addition to tracking these three seismic attributes (similarity, spectral decomposition, and curvature) across the horizons shown in Figure 4 (results of which are shown in Fig. 6), we choose additional sub-horizons sh-A and sh-B around the PNR wells (the results of which are shown in Fig. 7). These sub-horizons run through the centres of the microseismic clouds generated during stimulation of PNR-1z and PNR-2, and as such should be optimally positioned to identify the re-activated faults.
The similarity attribute is able to pick up both groups of faults (normal and thrust) described above, represented by black lines or dissimilarities in Figures 6a, 7a and b. Likewise, the spectral decomposition attribute (Figs 6b, c and 7c-f), especially at higher frequencies. The curvature attributes display some features which do not show up on the similarity and spectral decomposition attributes. These are represented by continuous lines of high-or  low-curvature values in Figure 7g and h. Our interpretation is that these geological features are offsets from faults and remnants of ancient channels.
However, none of attributes are able to confidently and unambiguously identify the reactivated faults near the PNR site (Fig. 7). Although the similarity and spectral decomposition attributes cannot directly pick up these faults, they can detect some evidence which could be used to infer their presence. The similarity and spectral decomposition attributes can pick up geological features (red dashed rectangles in Fig. 8) which terminate at the location of the PNR-1z fault plane identified by Clarke et al. (2018), perhaps indicating the offset of geological features. The high-frequency component of the spectral decomposition attribute can detect faint dark lines (white dashed rectangles in Fig. 8) near the location of the fault planes proposed by Clarke et al. (2019a) and . Nevertheless, without prior knowledge of the fault plane locations from microseismic observations, these faint features are not distinguishable above the typical background variability observed across the section. As such, we conclude that these seismic attributes are still not capable of reliably detecting the causative faults at PNR-1z and PNR-2 wells.

Automatic fault interpretation
To complement the manual, and attribute-assisted, 3D seismic interpretation, a thinned fault likelihood (TFL) attribute was calculated for the same dataset to detect fault surfaces automatically. This TFL is based on a semblance attribute, and therefore is also a coherence measurement between seismic traces (varying between 0 and 1). Then, a structure-oriented filter is applied to thin the fault likelihood, to generate fault images with higher resolution when compared with conventional discontinuity attributes as semblance or curvature (Hale 2013).
Hundreds of faults are clearly visible from the TFL attribute calculated from the 3D seismic (Fig. 9). To extract the geometry of the fault planes (i.e. centre location, length and azimuth), observed from the TFL attribute, we used the standard Hough transform (SHT) which is a feature extraction algorithm commonly implemented as part of image processing and computer vision workflows to extract analytically defined shapes, as circles or lines (the simplest case of the Hough transform), from images. This SHT method first calculates the edges from the input image to obtain a binary version (black and white background from Fig. 10a), and then uses the parametric equation of a line: where ρ is the perpendicular distance from the origin to the line, and θ is the angle of the perpendicular projection to the line, also from the origin. The SHT transforms the binary image into the ρ−θ domain (Fig. 10b), where each peak extracted from the ρ−θ image corresponds to one line (shown in green in Fig. 10a). Of the faults identified by the TFL attribute ( Fig. 9), some are compatible with the major faults identified in the manual geological interpretation described above (Fig. 4). Some NE-SW striking faults with high 'likelihood' (i.e. TFL closer to 1) are identified. Their location and orientation are close to that observed for the 2011  Table 1). These lengths confirm one of the issues of using 3D seismic data to guide site selection to mitigate induced seismicity: faults of sufficient size to cause M L 3.0 earthquakes may be at the very limit of survey resolution. While such events are unlikely to be of sufficient size to cause any significant damage (e.g. Nievas et al. 2020), since they occur at relatively shallow depths, they may be felt by people nearby, leading to significant public unease with the causative industries.

Fault slip potential
Once the fault lines are extracted from the TFL attribute from the 3D seismic survey, their fault slip potential (FSP) can be estimated by integrating them with the formation stress and pore pressure gradients. This FSP analysis is commonly done by integrating the effective normal and shear stress, the fault planes, and a failure boundary in a Mohr diagram (Walsh and Zoback 2016;Walsh et al. 2018), as shown in Figure 11b. To do so, first, the effective normal stresses (σ) must be calculated by subtracting the reservoir's pore pressure. At the depth of the horizontal section of the PNR wells Fig. 9. The thinned fault likelihood attribute (TFL) from the Bowland-12 3D seismic dataset (time slice at 1200 ms, just below the PNR wells as shown in Fig. 4d), shows a comprehensive understanding of the fracture network present in the Bowland Shale at 3D seismic-resolution scales (see the faultlength ranges obtained in Fig. 10c), which is compatible with the previous geological interpretations of major faults (black solid lines) shown in vertical sections of the same 3D seismic dataset (Fig. 4). The highlighted inset with the microseismic events observed during the hydraulic stimulations of the PNR wells is shown in detail in Figure 3, with the same vertical and horizontal slices of the 3D seismic dataset near the PNR wells shown in Figure 5. (shown in Fig. 3c), of c. 2200 m (7218 ft), and with stress gradients of S HMAX = 1.4 psi/ft, S V = 1.114 psi/ft, and S hmin = 0.725 psi/ft, and a pore pressure gradient (P p ) of 0.58 psf/ft, the effective normal stresses are ¼ 7218 ftÃ(0:725 psi=ft À 0:58 psi=ft) ¼ 1047 psi: From these normal stresses, the shear stress (τ) can be calculated for any fault orientation from a 3D Mohr diagram, by calculating three Mohr circles centred between the normal stresses and at a shear stress of zero (Fig. 11b). Then, each fault plane obtained from the 3D seismic dataset in the Bowland is plotted inside the same Mohr circles from their orientation relative to the orientation of the maximum horizontal stress, of 173°in this area. Then, we measure the distance between each fault plane (i.e. inside the Mohr circles) and a failure function (as the red line shown in Fig. 11b), which corresponds to the additional pore pressure required to fail (i.e. to slip). In this study we used a linear Coulomb failure function, where τ = μσ, and a friction coefficient (μ) of 0.75 that was first estimated for this area (Clarke et al. 2019b). However, this friction coefficient, which typically ranges between 0.6 and 1.0, can be highly variable in a formation as heterogeneous and highly fractured as the Bowland Shale, and the additional pore pressure required for each fault line to reach the failure envelope varies significantly depending on the friction coefficients. We find that several of the identified faults surpass the Mohr-Coulomb failure threshold (red line in Fig. 11b) meaning that they are critically stressed, due to their optimal orientation relative to the maximum horizontal stress (Fig. 11c, d), and the high pore pressure that significantly reduces the effective normal stresses. The abundance of critically stressed faults in the Bowland Shale likely explains why all wells that have been hydraulically fractured in this formation generated induced seismicity. Also, many smaller fractures in the same formation, which were undetected from the 3D seismic dataset in the Bowland shale, could also be critically stressed, which can cause wellbore stability issues while drilling the wells, as reported for the vertical wells drilled in the area (as in the Preese Hall 1 well shown in Fig. 1c; Clarke et al. 2019a).
Hydraulic fracturing for shale gas is currently subject to an ongoing moratorium in the UK, although we note that no such moratorium exists for similar activities such as geothermal stimulation (which has produced similar levels of induced seismicity, e.g. Holmgren and Werner 2021), and for smallerscale hydraulic fracturing in conventional reservoirs, which has been conducted in the UK for decades, with no cases of induced seismicity recorded (Mustanen et al. 2017). Similarly, hydraulicfracturing in offshore fields in the North Sea is not under moratorium. Based on our findings here, future work would be well-directed to investigating the relative abundance of faulting, and the in situ stress conditions, across the wider Bowland Shale formation, as areas with lower abundance of faulting and lower shear stresses may be less prone to induced seismicity. However, it should be recognized that faults of sufficient size to generate M L 3.0 earthquakes, which will be strongly felt by people nearby, may be  (Fig. 9) and extracted from a 2D standard Hough transform (Fig. 10), just below the PNR wells (PNR-1Z and PNR-2) hydraulic-fractured in 2018 and 2019 respectively. The colour code for each fault corresponds to the increase in pore pressure (ΔP P ) required in each fault to reduce their effective stress and reach the Mohr-Coulomb failure envelope (red line shown in (b)), and therefore to slip. In the linear failure envelope used in this study, where τ = μσ, a friction coefficient (μ) of 0.75 was first estimated for this area (Clarke et al. 2019a). However, this friction coefficient, which typically ranges between 0.6 and 1.0, can be highly variable in a formation as heterogeneous and highly fractured as the Bowland Shale, and the additional pore pressure required for each fault line to reach the failure envelope varies significantly depending on the friction coefficients. These pore pressures calculated from the same fault lines can also range from zero (as several faults are already reaching in some cases the Mohr-Coulomb failure envelope with μ = 0.75, meaning that they are critically stressed) to more than 1000 psi. The orientation of these critically stressed faults (shown in (c) and (d) in the normal composite and stereonet projection respectively) are also optimally oriented for a strike-slip faulting mechanism relative to the maximum horizontal stress (S HMAX ). below the limits of seismic resolution. As such, pro-active decision making based on real-time microseismic observations (e.g. Clarke et al. 2019b;Verdon and Bommer 2021), will be required.

Conclusions
A detailed structural interpretation of an unconventional reservoir, based on multiple exploration datasets including well logs and 2D and 3D seismic surveys, is an essential step required to determine an optimum well location and design of hydraulic-fracturing stimulations, especially in areas prone to induced seismic activity such as the Bowland Shale. The PNR wells near Blackpool, hydraulically fractured in 2018 and 2019, generated multiple 'red-light' seismic events, including a M L 2.9 event in August 2019. Analysis and interpretation of the 3D seismic dataset did not show any large faults (or seismic discontinuity) near these wells, and the faults were only visible from the microseismicity recorded with downhole arrays during the hydraulic stimulations. We tested multiple seismic attributes to complement the manual seismic interpretation of the main horizons and faults from the same 3D seismic dataset, observing a high structural complexity below the Variscan Unconformity.
To complement this manual seismic interpretation, we tested an automated fault detection method based on a high-resolution coherence attribute, and analysed the slip potential of faults that were identified. This method allowed the detection of numerous critically stressed faults in the Bowland Shale with a minimum fault length of 0.4 km ( potential for magnitude 3 earthquake, approximately), some of them near the hydraulic fractured wells in Preese Hall and PNR and with similar orientations to the nodal planes reported for the largest-magnitude seismic events detected near each well. This method could also be implemented to evaluate possible locations of future unconventional wells that require similar hydraulic-fracturing stimulations, by avoiding critically stressed faults that could trigger anomalously high seismic activity. However, it is also clear that multiple strike-slip faults could go undetected if they are shorter than the minimum fault length detected, or if they have a very low vertical displacement, some of which could also be critically stressed if they are optimally oriented relative to the maximum horizontal stress.
Even if the application of the TFL and the FSP in the Bowland Shale does not reliably pick up strike-slip faults, there is still a great potential that both methods can improve the drilling site selection and operation planning. The ability to better detect faults and calculate the probability of each fault to rupture is essential since avoiding critically stressed faults will tremendously reduce the chance of triggering HF-IS. Furthermore, these processes are performed before the drilling, meaning that operators can manage the problem before it actually occurs, unlike with TLS. We predict that, combining with the traditional, manual methods, the improved TFL and FSP will play a significant role in risk assessments of unconventional explorations in the future, not only for the Bowland Shale, but for other regions as well.