Bowland Shale well placement strategy – Part 1: Determining landing intervals using geomechanical properties

The production performance of a shale reservoir is directly affected by the geomechanical characteristics of the formation. A target shale interval will ideally develop hydraulic fractures upon stimulation that stay open with the aid of injected proppant. However, shales are geomechanically complex due to heterogeneities in their rock properties such as mineralogy and porosity and the extent to which they may be naturally fractured. These characteristics can complicate the task of identifying the ideal target interval for placing horizontal wells. Whilst the Bowland Shale is the UK ’ s most prospective shale gas target, long horizontal wells are generally not feasible or practical in the Craven Basin, due to the existence of many, large-offset reverse faults and high bedding dips. An alternative to this approach could include drilling shorter, stacked horizontal wells targeting different stratigraphic intervals. However, it is unclear if there are enough intervals within the stratigraphic section with the desired geomechanical properties to target with stacked horizontal wells, nor if there are adequate intervals that can limit vertical hydraulic fracture growth between those wells. The absence of the latter may ultimately lead to well interference and reduced production. These issues were addressed by the creation of a series of wireline log-based geomechanical logs at well Preese Hall-1, calibrated to pressure test data. Aided by the results of a cluster analysis model, the upper section of the Bowland Shale was classified into geomechanical zones to identify the optimal intervals for hydraulic fracturing and barriers to vertical hydraulic fracture growth. Three intervals are highlighted with low effective stress, low fracture toughness and high brittleness which may form excellent landing zones. Importantly, these landing zones are also separated by intervals of high effective stress that may limit vertical hydraulic fracture growth and mitigate the risk of well interference.


Introduction
The exploitation of tight mudrocks, or shales, for hydrocarbon production, revolutionized the global energy industry during the early 21st Century and posed unique challenges not typically encountered when producing from conventional, higher permeability reservoirs. The USA became a net gas exporter for the first time in 2017 (EIA, 2020a), driven largely by the country's widespread adoption of multi-stage hydraulic fracturing, but this success has not yet been emulated to the same degree in other countries looking to harness their unconventional resources.
Production from shales is dependent on the creation of an interconnected hydraulic fracture network that serves to unlock gas from the very low-permeability matrix and deliver it to the production well. Thus, the key consideration when determining a production well placement strategy is stimulating as much of the reservoir as is commercially possible. In many of the prolific USA shale plays this involves drilling horizontal wells in the target formation and stimulating these wells at multiple stages. Such wells can be very long and, in 2019, the average horizontal well length in the USA reached 4.5 km (EIA, 2020b).
However, these shale developments benefit from relatively unstructured subsurface geology. In areas that are relatively highly structured, faulting poses a major challenge and can contribute to induced seismicity during hydraulic fracturing operations. In the UK, the British Geological Survey (BGS) produced an initial screening study of prospective candidates for domestic shale extraction, concluding that the Mississippian Bowland Shale Formation of north England posed the best opportunity (Smith et al., 2010). Many recent studies have addressed the formation's petrophysical, sedimentological, and microstructural properties within the Craven Basin (Clarke et al., 2014a(Clarke et al., , 2018Fauchille et al., 2017;Newport et al., 2018;Emmings et al., 2019Emmings et al., , 2020a and it is generally believed to possess many of the properties required for shale gas extraction (Clarke et al. 2014a(Clarke et al. , 2018. However, there is some substantial uncertainty regarding the total size of the resource with estimates ranging from 1329 TCF (median estimate and including the underlying Hodder Mudstone Formation) (Andrews, 2013) to 140 TCF (Whitelaw et al., 2019). Furthermore, in the Fylde area where hydraulic fracturing has been tested, the Bowland Shale has been compartmentalized by large-offset (Anderson and Underhill, 2020) and small-offset (Clarke et al., 2019b;Verdon et al., 2020) faults, due to the region having undergone multiple episodes of extension and uplift. This limits where horizontal production wells may be drilled ( Fig. 1) (Anderson and Underhill, 2020), and has led to elevated levels of induced seismicity during hydraulic fracture operations (Clarke et al, 2014b, Clarke et al., 2019bVerdon et al., 2019Verdon et al., , 2020.
These issues highlight the need for a well placement strategy that does not involve drilling long horizontal wells and maximizes stimulation of the shale resource in a spatially confined area. Stacked multilateral production wells may solve this issue and this concept has been briefly mentioned in the literature (Clarke et al, 2014a, Clarke et al., 2018, but within the context of minimizing the surface footprint of well pads, rather than avoiding geological faults. Stacked drilling allows for multiple stratigraphic intervals to be targeted from a single pad, but the strategy relies on there being several distinct intervals within the formation with geomechanical properties conducive to hydraulic fracturing. This paper is the first of a two-part series addressing the issue of developing a well placement strategy for the Bowland Shale. In this paper, the optimal landing zone locations for horizontal wells within the formation are determined. In another paper (in preparation), fracture simulations are deployed to study the effects of geomechanical properties on fracture geometries, the stress shadowing effect, and the presence of fracture barriers. The studies focus specifically on well Preese Hall-1 (PH-1) (Fig. 1); drilled by Cuadrilla Resources in 2010. Digital wireline log data and the results of core analyses were accessed through the BGS and were supplemented with published results of mini-frac testing for the Bowland Shale (de Pater and Baisch, 2011;de Pater and Pellicer, 2011). In this first paper, the results of wireline-derived geomechanical models are presented and cluster analysis is used to identify the optimal landing zones within well PH-1. Included in these models are profiles of effective stress, rock strength, and fracture toughness, which are considered some of the primary controls on hydraulic fracture geometry (Smith and Shlyapobersky, 2000), and, ultimately, commercial development of a shale resource. These characteristics can be grouped and assessed using a Completion Quality (CQ) index (Cipolla et al., 2011;Jochen et al., 2011;Miller et al., 2011;Suarez-Rivera et al., 2011). Fig. 1. Two-way time-structure map for the top Lower Bowland Shale on the Fylde region highlighting well Preese Hall-1 (PH-1), additional wells, and interpreted faults (after Anderson and Underhill, 2020). The inset map shows the location of the Becconsall-1Z well which lies outside the main study area. The map is projected to the British National Grid coordinate system with X and Y coordinates expressed in meters. Bec-1Z: Becconsall-1Z, El-1: Elswick-1, GH-1Z: Grange Hill-1Z, Th-1: Thistleton-1, PNR-1Z: Preston New Road-1Z, PNR-2: Preston New Road-2, RW: Roseacre Wood.

Geological setting
The Bowland Shale Formation is a Carboniferous (Visean-Serpukhovian), predominantly fine-grained unit that was deposited in a series of basins across northern England. These basins formed in response to Late Devonian to Lower Carboniferous extension and were separated by broad highs; a configuration often referred to as 'block and basin' topography (Waters and Davies, 2006). Well PH-1 was drilled in the Craven Basin, which is a SW-NE-striking graben defined by the Lake District Massif to the north-west, the Askrigg Block to the north-east, and the West Lancashire High to the south-east (Kirby et al., 2000). The basin was inverted both during the Late Carboniferous (forming the Ribblesdale Fold Belt (Arthurton, 1984)) and the Cenozoic (Kirby et al., 2000). Shale gas prospectivity exists in the west portion of the basin where Carboniferous sediments exposed at the Ribblesdale Fold Belt plunge beneath a Permo-Triassic cover, though inversion-related deformation persists in this area (Anderson and Underhill, 2020).
The formation is formally subdivided into lower (Lower Bowland Shale, LBS) and upper (Upper Bowland Shale, UBS) units at the Visean-Serpukhovian boundary, which marks a tectonic transition from extensional-driven subsidence to post-rift sag subsidence across the region (Fraser and Gawthorpe, 1990;Waters et al., 2009). Compositional changes can also be observed between the LBS and UBS with the former reflecting a carbonate-dominated system (Newport et al., 2018) and the latter reflecting an increasingly siliciclastic system Fraser and Gawthorpe, 2003;Waters et al., 2009;Emmings et al., 2020a). As a result, the formation exhibits acute compositional and textural heterogeneity on multiple scales (Ma et al., 2016;Fauchille et al., 2017;Clarke et al., 2018).

The Mechanical Earth Model concept
The concept of systematically deriving depth-profiles of the stress, elastic, and strength properties for a rock formation has been formally defined by Plumb et al. (2000) as a Mechanical Earth Model (MEM). The approach has been extensively applied within the unconventional setting including the use of MEMs to better predict stimulated rock volume in shale reservoirs (Li et al., 2014), determine landing points for horizontal wells (Alimahomed et al, 2017(Alimahomed et al, , 2018Qiu et al., 2013) and improve completions with coal-seam gas reservoirs (Johnson et al., 2010).
Elements of the MEM approach were adopted for this work to study the stratigraphic variability in geomechanical and stress properties at well PH-1 (Fig. 2). The model equations and relevant references are further outlined in Table 1 and background to the key concepts can be found within the supplementary material to this paper. The MEM can be separated into two categories: stress properties ( Fig. 2 -left) and mechanical properties (Fig. 2 -right). Within the stress properties category, four steps ultimately build towards determining the minimum horizontal stress. This is a key parameter to quantify, as stratigraphic variations in the least principal stress (assumed equal to minimum horizontal stress) have a direct impact on hydraulic fracture growth and unconventional reservoir stimulation (Zoback and Kohli, 2019). In calculating minimum horizontal stress, separate models for vertical stress (Eq. 1), pore pressure (Eq. 2), Biot coefficient (Eq. 3), and static elastic properties were first created ( Fig. 2 -left) and then deployed in Barree et al. (2009)'s derivation of the poroelastic equation for minimum horizontal stress (Eq. 4). These models mainly used wireline log data from well PH-1 but these were calibrated with the results of mini-frac tests or core measurements where available.
In addition to deriving stress profiles for the Bowland Shale, two key failure metrics were also determined. Brittleness is a parameter often considered in a shale gas evaluation as a brittle rock will fracture with greater ease and remain better propped open than a ductile rock will (Zhang et al., 2016), though its applicability in picking target zones for hydraulic fracturing has been challenged (Bai, 2016). Brittleness index was quantified in this study using log data and adopting Rickman et al. (2008)'s approach (Eq. 5). Fracability index (Jin et al., 2014) was also calculated, which accounts for both brittleness and fracture toughness (Eq. 7). Fracture toughness is an important parameter to define as it quantifies the ability of a rock to resist hydraulic fracture propagation (Yuan et al., 2017), and its calculation from wireline logs is usually achieved using empirical correlations derived from confining pressure and tensile strength (Eq. 6; Table 1) (Chen et al., 1997a(Chen et al., , 1997bJin et al., Fig. 2. Illustration of the workflow followed in the calculation of subsurface stress (left) and mechanical (right) properties using wireline log models calibrated to mini-frac test data. The workflow adopts elements of the Mechanical Earth Model (MEM) approach proposed by Plumb et al. (2000). 2011; Sui et al., 2019). Following the development of these models, a k-means cluster model (Pedregosa et al., 2011) was used to aid the picking of intervals with the best CQ characteristics in well PH-1.

The PH-1 hydraulic fracturing test
The models used to calculate pore pressure and minimum horizontal stress were calibrated with the results of a multi-stage hydraulic fracturing operation conducted at the vertical well PH-1 in 2011. The operation involved testing six stages of the Bowland Shale; a minifrac test was performed at all stages, and a full hydraulic fracturing test was performed at five of the six stages (de Pater and Baisch, 2011). The stages varied in length between 18 m and 73 m and each stage consisted of three, 2.7 m long perforated intervals, each with 27 perforations (de Pater and Pellicer, 2011, courtesy of Cuadrilla Resources). During these tests, two large seismic events of greater than 1.5 local magnitude were recorded and were later attributed to the reactivation of a critically stressed reverse fault near the base of the well (Clarke et al., 2014b)..
As a consequence of the induced seismicity at well PH-1, the operator commissioned several studies to determine the cause of the events (de Pater and Baisch, 2011;de Pater and Pellicer, 2011). de Pater and Pellicer (2011) determined fracture closure stress for five of the six mini-frac tests by analyses of pressure decline curves, however, the tests conducted in the upper stages (4 and 5) were affected by previous fracture tests, which led Clarke et al. (2019a) to question their reliability. Therefore, for this study, it was decided to focus on the results of the three lowermost stages (stages 1-3) for calibration. The closure stress readings taken from de Pater and Pellicer (2011) were combined with pore pressure measurements taken from Clarke et al. (2019a) to produce a small calibration dataset (Table 2).

Static and dynamic elastic properties
Young's modulus and Poisson's ratio are important components of the minimum horizontal stress equation (Eq. 4 (Table 1)). In producing a wireline log-based MEM, these are usually calculated from elastic-wave velocity and density wireline logs. However, elastic moduli derived in this manner (termed 'dynamic') and those derived more traditionally from deformation experiments ('static') are not equal, even for the same rock (Mavko et al., 2009). Eq. 4 (Table 1) handles static elastic properties, and therefore it was necessary to correct the dynamic elastic properties from logs before calculating minimum horizontal stress.
The results of static and dynamic (ultrasonic) tests on core samples at well Bec-1Z ( Fig. 1) were used to calculate dynamic to static corrections for Young's modulus (E) and Poisson's ratio (ν). In both cases, least squares regressions were calculated for the test results, which in the case of Young's modulus, resulted in a similar equation to Wang and Nur (2000)'s equation for soft rocks (Fig. 3). The two equations were then taken as follows: Eq. 8 Eq. 9

Vertical stress
The bulk density wireline log at well PH-1 was first extrapolated to the surface using a geometric fit to three points ( Fig. 4) and then integrated as per Eq. 1 (Table 1). This resulted in a vertical stress gradient of approximately 23.6 kPa/m, reaching a maximum stress of 64 MPa at total depth.

Variable
Brief Description Calculation Equation Key References

Vertical Stress
The stress normal to the earth's surface due to the weight of the overlying rock layers.
Integration of bulk density log with depth.
Jaeger et al. (2007) Pore Pressure The pressure exerted by fluids within a rock's pore space.
Sonic velocity deflection from a Normal Compaction Trendline. Eaton (1975); Zhang (2011) Biot Coefficient The effectiveness with which pore pressure counteracts confining pressure to produce volumetric strain.
Empirical correlation with the porosity log.

Horizontal Stress
The stress that is required to initiate a fracture.
Poroelastic theory combined with the tectonic microstrain offset parameter.

Brittleness Index
The ratio of compressive strength to tensile strength.
Young's modulus and Poisson's ratio wireline logs.
Eq The critical stress intensity factor at which a fracture becomes unstable and propagates with high velocity.
Empirical correlation with confining stress and tensile strength, for mode I cracks only. Metric that combines Brittleness Index with Fracture Toughness.
Brittleness index plus normalized fracture toughness.
Eq. 7 Jin et al. (2014) et al., 2004;Zhang, 2011): Where Δt m is the compressional transit time in the shale matrix, Δt ml is the mudline transit time, c is a constant and Z is true vertical depth. Shallow logging sections from well PH-1 and two additional nearby wells (wells El-1 and Th-1; the locations of which are shown in Fig. 1) were combined to derive a representative, normally-pressured compressional-slowness trendline that was then extrapolated to deeper (assumed overpressured) formations. Clarke et al. (2019a) state that overpressure is believed to start beneath the anhydrites of the Manchester Marl, therefore only the formations above the Manchester Marl were used. Shallow sections of the three wells were plotted together and a logarithmic curve was fitted to the data to produce the following relationship for the normal compaction trendline (Fig. 5): Eq. 11 Pore pressure was then calculated as per Eq. 2 (Table 1) with the hydrostatic pore pressure gradient fixed at 10.18 kPa/m. A GR cut-off of 75 API was used to filter the Bowland Shale interval to only include finegrained intervals for use in the pore pressure calculation. Such a cut-off is traditionally used for separating matrix from "shale" in conventional reservoirs, however, adopting such terminology for the Bowland Shale is problematic as many clay-poor intervals in the Bowland Shale also exhibit high GR responses. Synthetic GR responses were calculated from a Bowland Shale Th, U, and K dataset reported in Appendix 4 of Emmings et al. (2020a) which showed only intervals of Facies A (blocky limestones (Emmings et al. (2020b)) exhibited GR responses of <75 API. The remaining silt and carbonate-prone mudstone facies typically exhibited GR > 100 API. This suggests that a 75 API cut-off does not accurately define mudstone from other fine-grained facies, at least in the proximal part of the Craven Basin described by Emmings et al. (2020a). However, as this threshold is used to discriminate fine-grained (and prone to overpressure) intervals from coarse-grained intervals, it is considered adequate in filtering the coarse-grained carbonates and sandstones from the pore pressure calculation.
The Eaton exponent (β in Eq. 2; Table 1) within the pore pressure calculation is a parameter required to be fitted to any calibration data available. In Eaton (1975)'s paper, an exponent of 3 was proposed, however, for this work, a range of exponents between 0.5 and 3 were tested and compared with the three calibration data points (Fig. 6). This suggested that an exponent of 0.8 was suitable for the calibration data available and that value was subsequently used in the calculation of the entire pore pressure profile.
The start of overpressure is observed in the Manchester Marl, following which, pore pressure increases from ~10 MPa to ~40 MPa at a gradient of 15.2 kPa/m (Fig. 7). The UBS exhibits moderate overpressure of up to 10 MPa greater than hydrostatic pressure, but the greatest overpressure is observed within the clay-rich, thick shales of the LBS. These exhibit pore pressures up to 40 MPa; 13 MPa greater than hydrostatic pressure.
Several empirical models have been proposed to determine Biot coefficient from porosity (Krief et al., 1990;Laurent et al., 1993;Lee, 2002;Halliburton, 2019), and these were compared using a small sandstone dataset published by Detournay and Cheng (1993) (Fig. 8). None of the models accurately predicted α for low-porosity samples, besides the unpublished Halliburton (2019) equation, and therefore the linear regression model shown in Table 1 was selected. Applying this method resulted in a α range of between 0.63 and 0.85 for the Bowland Shale interval.

Minimum horizontal stress
Previous studies have determined that the subsurface stress state around well PH-1 is that of strike-slip (Clarke et al., 2014a(Clarke et al., , 2014b, and in such a setting, the least principal stress is equal to the minimum horizontal stress (Anderson, 1905;Hubbert and Willis, 1957). For this work, minimum horizontal stress was calculated using poroelastic theory combined with the microstrain parameter (Barree et al., 2009;Blanton and Olson, 1999) (Table 1; Eq. 4).
The microstrain parameter is a constant, applied to the entire formation that shifts the calculated minimum horizontal stress by a value proportional to its strength (by multiplying it by the layer's Young's modulus, E). This constant was determined by calibrating the modeled minimum horizontal stress log to data from minifrac pressure declines (de Pater and Pellicer, 2011) (Table 2, Fig. 9), which suggested that a microstrain of zero was appropriate.
Following calibration of the microstrain parameter, the modeled minimum horizontal stress curve (Fig. 10) was calculated according to Eq. 4 (Table 1). The general gradient was calculated as 17.4 kPa/m, however, the curve highlighted several high minimum horizontal stress zones (shown in red on Fig. 10), where the calculated stress approaches  Eissa and Kazi (1988) c: Wang and Nur (2000), d: Mese and Dvorkin (2000). The least-squares regression for well Bec-1Z samples is between the Wang and Nur (2000)'s soft rocks and Eissa and Kazi (1988)'s range of rocks relationships.

Table 2
Pore pressure and closure stress values from minifrac tests across three stages of well PH-1 and collated from de Pater and Pellicer (2011) and Clarke et al. (2019a). The depth quoted (expressed as Measured Depth (MD)) is the arithmetic average taken over the complete stage interval which can be up to 70 m in length. As the well is sub-vertical, MD is a reasonable approximation for True Vertical Depth (TVD).

Stage
Depth ( the 20 kPa/m gradient, and which corresponds to a maximum stress of ~50 MPa at the base of the well. These high minimum horizontal stress zones are predominantly located within the LBS, though three zones can also be observed within the UBS. For the majority of the Bowland Shale, the minimum horizontal stress gradient lies between 15 kPa/m and 20 kPa/m, though in many intervals, it drops to a gradient of 15 kPa/m. These lower minimum horizontal stress regions are highlighted in green in Fig. 10. Zones are also highlighted where pore pressure exceeds the minimum horizontal stress. While increasing the tectonic microstrain parameter would avoid this occurrence, it is hard to justify this, based on the fit to the minifrac data. Therefore, it is assumed these anomalies correspond to intervals where minimum horizontal stress effective stress is low and close to the pore pressure (i.e. very low effective stress).

Rock failure estimates
Brittleness Index (BI) and Fracability Index (FI) were calculated using Eqs. 5 and 7 in Table 1, respectively, where fracture toughness is considered for mode I cracks (K IC ) using Eq. 6. In Eq. 6, the confining stress was determined using empirical correlation (Sui et al., 2019) between estimated clay fraction and dynamic Young's modulus where the ratio of confining to tensile stresses was set to 12.26, following those   6. Crossplot of calculated pore pressure versus minifrac pore pressure for several values of Eaton exponent (β). Also included is a y = x trendline along which the cross-plotted values would be preferred to lie. Input values for the pore pressure model (Eq. 2) are averaged over the entire test interval to produce a single estimate for each minifrac data point. Reducing β to 0.8 produces the best fit for the available calibration data. Fig. 4. Log display illustrating the calculated vertical stress profile (track 4) derived from the integration of the bulk density log with depth (Eq. 1). The extrapolated bulk density log is also shown; calculated using a geometric fit to three data points (7.8 m (mudline): 1.65 g/cm 3 , 638 m: 2.4 g/cm 3, and 2704.5 m: 2.64 g/cm 3 ). Contains OGA Well Log © Data accessed and published with permission of BGS.
authors. Brittleness (track 4) is largely greater than 0.6 throughout the section except for some discrete intervals where it drops to ~0.2 (Fig. 11). Several low brittleness (i.e., ductile) intervals can be easily identified as plotting within the first quartile and are colored red. These are present throughout the stratigraphic section but are thickest within the LBS. The most brittle intervals were highlighted in green and are mostly within the Pendleside Sandstone Member of the LBS and at the base of the UBS.
Fracability index was also found to be high; exhibiting readings mainly between 0.3 and 0.5, but up to 0.9 in discrete intervals. It follows a very similar trend to brittleness index; however, it accentuates some intervals where fracture toughness (track 3) is notably low. The brittle section near the base of the UBS described above continues to exhibit excellent failure properties, exhibiting both low fracture toughness and high fracability index. While there is generally a strong correlation between brittleness and fracture toughness, there are some intervals (e.g, 2155 m) where they do not correlate; here, low fracture toughness is accompanied by low brittleness index, effectively canceling each parameter out in the calculation of fracability index, which shows a relatively flat response. This suggests that calculating brittleness or fracability in isolation is not sufficient and both need to be considered when picking intervals for hydraulic fracturing.
These rock failure metrics are challenging to assess when considered on their own and do not show a simple correlation. Therefore, a classification scheme may provide insight into the relationship between these that cannot be determined by simply analyzing the metrics in log view. To tackle this, brittleness and fracture toughness were combined with effective stress in a cluster model to attempt to pick the optimal Fig. 8. Plot illustrating three published empirical relationships for calculating the Biot coefficient from porosity (Krief et al., 1990;Lee, 2002;Laurent et al., 1993), the un-published equation implemented in GOHFER hydraulic fracturing simulator (Halliburton, 2019), and a linear least squares regression model fitted to a small dataset presented in Detournay and Cheng (1993). None of the published methods accurately model the low-porosity samples.   Eaton (1975) pore pressure model (shown in track 5). The logs are plotted in measured depth (MD); however, the well is sub-vertical and thus MD is approximately equal to TVD. Also shown are the GR (track 3) and the DTC log, together with the normal compaction trendline calculated using Eq. (11). Hydrostatic pressure (10.18 kPa/m) and a generic pressure gradient of 15 kPa/m are also shown for reference. High-pressure zones are highlighted in red from visual inspection of the pore pressure curve. Contains OGA Well Log © Data accessed and published with permission of BGS. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) geomechanical characteristics within the formation.

Classification
This paper has thus far presented the results of a series of geomechanical models to derive logs of in-situ stresses and rock failure indicators for the entire Bowland Shale interval at well PH-1. To focus on the identification of the best potential landing zones within the formation, modeled data were filtered to include only the best quality shale intervals (defined as Reservoir Quality (RQ) intervals; as outlined in Appendix A). Brittleness index, fracture toughness, and effective stress were selected as key parameters that categorize a rock's Completion Quality (CQ) and were input to an unsupervised classification algorithm to identify clusters that share similar values of input data.
Before deploying the algorithm, each input log was first normalized to a 0 to 1 scale to standardize the distributions of each input variable. One of the simplest clustering methods, the k-means algorithm, was selected which minimizes the sum of squared distances of samples to their nearest cluster center. It was implemented using the open-source python software package, scikit-learn (Pedregosa et al., 2011) and the code can be found within the supplementary material to this paper.
To select an appropriate k-value (the number of clusters used to classify the dataset), a range of k-values was tested, and the performance of the clustering was assessed using the inertia metric and silhouette score. Both metrics consider the distance of a classified sample within its cluster center and the cluster centers of neighboring clusters and are described in further detail within Pedregosa et al. (2011). The results of this analysis suggested that a k-value of 4 was optimal for the dataset, and thus the classification proceeded on this basis.
A series of two-dimensional cross-plots were then constructed to present the results of the cluster model (Fig. 12). As there were no Fig. 11. Profiles of Brittleness Index (track 4) and Fracability Index (track 7), together with the input logs involved in their calculation (tracks 3, 5, and 6). Red color shading is added between the minimum value and the first quartile (25%), and green shading is added between the third quartile (75%) and the maximum value. These assist in determining highly ductile/poorly fracable and highly brittle/highly fracable areas of the shale respectively. The logs are plotted in measured depth (MD); however, the well is sub-vertical and thus MD is approximately equal to TVD. Contains OGA Well Log © Data accessed and published with permission of BGS. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) Fig. 10. Well PH-1 wireline log display illustrating the results of the minimum horizontal stress calculation (shown in track 5). The logs are plotted in measured depth (MD); however, the well is sub-vertical and thus MD is approximately equal to TVD. Also shown are the static Young's modulus and Poisson's ratio logs (tracks 3 and 4) which are input to the poroelastic equation (Eq. 4). A tectonic microstrain parameter (ε h ) of 0 is used (see Fig. 12) and vertical Biot coefficient (α v ) is calculated from porosity using the linear function outlined in Fig. 5. The horizontal Biot coefficient (α h ) is set to unity. The pore pressure and vertical stress curves are also shown for reference in track 5, as are generic pressure gradients (10, 15, and 20 kPa/m). Track 6 shows the generalized gradients for each stress profile, determined through least-squares regression. Contains OGA Well Log © Data accessed and published with permission of BGS.
training data provided to the cluster model, a traffic-light labeling scheme was interpreted according to the typical geomechanical properties within that cluster, where red represents the poorest CQ and dark green represents the best CQ.
The degree to which the different clusters are separated is variable. There is good separation in the cross-plots where fracture toughness is an input parameter ( Fig. 12b and c). The cluster centers are wellseparated and there is minimal overlap in the classified points. In brittleness/effective stress space, however (Fig. 12a), there is overlap between the red, yellow, and pale green clusters suggesting a less effective classification on this basis.
The cluster shown in red ("poor CQ") exhibits the highest effective stress (>8 MPa), lowest brittleness index (<0.5), and highest fracture toughness (>1 MPa m 1/2 ). The yellow cluster also exhibits high, but slightly reduced effective stress (>7 MPa), greater brittleness (>0.55), and lower fracture toughness (<0.85 MPa m 1/2 ). This is interpreted as "moderate CQ" as while there is a clear drop in fracture toughness between this and the "poor CQ" cluster, the remaining properties exhibit only a slight improvement in CQ metrics.
The pale green cluster shows fracture toughness values almost as high as the "poor CQ" class but a further decrease in effective stress (<6 MPa) and an increase in brittleness index (>0.6). Given these latter two metrics show improvement in CQ properties, this is interpreted as a "good CQ" rock type. An increase in fracture toughness accompanied by an increase in brittleness in this manner may reflect a carbonate-rich section of the shale which requires greater energy to fracture due to its higher strength.
The dark green cluster shows the lowest effective stress (<5 MPa), the highest brittleness index (>0.6), and low fracture toughness (<0.85 MPa m 1/2 ). This cluster is interpreted as likely to hold the ideal mix of CQ properties and is subsequently labeled as "excellent CQ".
The "good" and "excellent" CQ intervals are located almost entirely in the UBS ( Fig. 13; track 7), and that figure shows the results for this unit only. The highest counts of "excellent CQ" intervals are within the lower section of the UBS (between ~2350 m and ~2420 m), where high brittleness index, low fracture toughness, and the lowest effective stress are observed. The highest counts of "good CQ" intervals are mainly within an interval in the middle of the UBS (between ~ 2220 m and ~2300 m) and at the top of the UBS (between ~ 2000 m and ~2075 m). Zones of poor RQ are observed between these intervals, which were not incorporated in the classification scheme but from visual inspection, can be seen to exhibit some poor CQ properties (e.g, low brittleness and high effective zones at 2150 m and 2320 m). Intervals that were flagged as good RQ but have "poor CQ" properties such as moderate-high effective stress and moderate-low brittleness index are shown as orange zones in track 9 and seem to separate the "good CQ" intervals from the poor RQ intervals shown in pink.

Picking landing points
The modeling results showed that the best geomechanical properties were located mainly within the UBS and within this unit, three landing points could be defined ( Fig. 13; track 9, Table 3). Landing 1 was picked near the top of the formation, immediately above the E 1b 2 marine band (Clarke et al., 2018), and within a "good CQ" zone. While two peaks in the rolling count of "excellent CQ" (Fig. 13; track 8) were observed, the landing zone was placed at the lower peak as this was farther from the top of the formation and thus poses a lower risk of hydraulic fractures Fig. 12. A series of cross-plots illustrating the results of the cluster model. The model runs only on the intervals identified as the best shale quality at well PH-1. The data points are colored by their allocated cluster, with the larger, black-outlined points representing the cluster centers. The color scheme is interpreted as different classes of CQ and is explained further in the text. Contains OGA Well Log © Data accessed and published with permission of BGS. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) propagating outside of the target interval.
Landing 2 was picked at 2290 m, immediately above a possible, but unverified, marine band at 2330 m. The landing zone was picked at a peak in the rolling CQ count ( Fig. 13; track 8), immediately below several "good CQ" intervals. Landings 1 and 2 are separated by an extensive, 105 m poor RQ interval, with a 30 m thick, high effective stress zone at 2138 m ( Fig. 13; track 5) which may form a hydraulic fracture barrier and prevent interference between the two landing zones if targeted.
Landing 3 was picked at 2420 m where a peak in the rolling CQ count ( Fig. 13; track 8) was again observed. The base of the UBS, marked by the E 1a 1 marine band, is known to exhibit intense quartz cementation (Emmings et al., 2020c), and this marine band has been identified within well PH-1 at 2500 m (Clarke et al., 2018). It appears likely that the cementation is contributing to the high CQ within the interval. Immediately overlying the interval is a 62 m section of "good CQ". The spacing between Landings 2 and 3 is less than between Landings 1 and 2; however, an interval of low brittleness index and high effective stress between 2299 m and 2348 m may form a hydraulic fracture barrier and limit interference between adjacent wells. This interval and other similar intervals in the UBS may reflect the weakly-cemented hemipelagic or density flow mudstones as described in Emmings et al. (2020b).

Overpressure in the Bowland Shale
The subject of pore pressure and overpressure has not been discussed in detail for the Bowland Shale, however, it is a key parameter in determining the amount of gas-filled porosity as part of a resource assessment. Some early studies suggested there was a lack of overpressure in Bowland Shale basins (Smith et al., 2010), leading to the assumption of hydrostatic pressure in the British Geological Survey's resource estimate for the shale (Andrews, 2013).
Observation of the pore pressures derived from minifrac testing ( Table 2) alone suggests that this is not the case; with these samples exhibiting overpressure of ~10-12 MPa (assuming a normal pore pressure trend of 10.18 kPa/m). Taking the onset of overpressure at the Manchester Marl, the modeling presented in this work suggested that the pore pressure increases at a gradient of ~15 kPa/m to a maximum of ~40 MPa at the base of the well, exhibiting overpressures up to 14.5 MPa within the LBS. In Clarke et al. (2018)'s synthesis of Bowland Shale exploration, a pore pressure curve from well PH-1 was presented which demonstrated the highest pore pressures within the base of the UBSperhaps steered by the calibration data-point 7895 ft (2406 m) depth which measured anomalously high pore pressure. That data point was considered poor quality in their subsequent study (Clarke et al., 2019a), and their amended pore pressure curve showed the greatest pore pressures in the LBS; matching the findings presented herein.
Pore pressure gradients for a variety of US shale plays were compiled and plotted for comparison with well PH-1 (Fig. 14; after Zhang, 2019). Fig. 13. Well log plot illustrating the geomechanical zonation scheme and proposed landing zones for the UBS. The plot shows the input logs to the cluster model (tracks 3-5), reservoir-quality intervals (track 6), CQ classifications output from the cluster model (track 7), a rolling count of "good" and "excellent" CQ sections over a 5 m interval (track 8) and the interpreted landing intervals and generalized geomechanical zones (track 9). In track 10, key marine bands (after Clarke et al. (2018)) and the approximate zone of abundant quartz cementation (after Emmings et al. (2020c)) are shown for reference. The logs are plotted in measured depth (MD); however, the well is sub-vertical and thus MD is approximately equal to TVD. Contains OGA Well Log © Data accessed and published with permission of BGS.

Table 3
Summary table of the interpreted, generalized geomechanical intervals, highlighting the three landing zones. These are also shown in Fig. 13. The Barnett Shale is an often-quoted analog to the Bowland Shale (Andrews, 2013) and is of similar mineralogy. It is considered a mildly over-pressured shale formation (11-12 kPa/m (Zhang, 2019);), suggesting it is generally less pressured than the Bowland Shale at this location. In areas where the Barnett Shale is mature, Bowker (2007) suggested the formation had an overpressure gradient of 11.8 kPa/m. Using an average reservoir depth of 2286 m (Bowker, 2007), this corresponds to 27 MPa pore pressure. After subtracting a hydrostatic pore pressure gradient of 10.18 kPa/m, this would correspond to 4 MPa of overpressure. The Eagle Ford shale is considered a more over-pressured shale play (9-18 kPa/m; (Zhang, 2019). Using Cander (2013)'s determination of a 16.7 kPa/m gradient and a reservoir depth of 3048 m this corresponds to 51 MPa pore pressure. After subtracting hydrostatic pore pressure, this corresponds to 20 MPa of overpressure. These values are greater than those modeled in this study and it would appear the overpressure exhibited by the Bowland Shale in this well is somewhere between that of the Barnett and Eagle Ford examples.

Top (M-MD) Base (M-MD) Thickness (M-MD) Interpretation
In many cases, overpressure contributes to low effective stress, which ultimately makes a shale easier to produce hydrocarbons (Cander, 2012;Wang and Gale, 2009). Wells within the highly overpressured Haynesville Shale can exhibit IP rates greater than 10 MMscf/day (approximately double that of less overpressured resource plays such as the Niobrara Shale (EIA, 2021)). But the link between pore pressure and effective stress is complex and will vary by shale resource play. The highest pore pressures in the well PH-1 correlated with high minimum horizontal stress (Fig. 10) and therefore moderate-high effective stress (Fig. 13), suggesting that overpressured intervals do not form the best targets for hydraulic fracturing. In this case, the zones of low effective stress shown in Fig. 13 appear to be caused by low minimum horizontal stress rather than high pore pressure.

Comparison with previous work
Following an induced seismic event near the base of well PH-1 (Clarke et al., 2014b), several technical reports were commissioned to study its cause, which were summarised in de Pater and Baisch (2011). A geomechanical model at well PH-1 was constructed as part of these studies and later expanded on by Clarke et al. (2018). This work used calibration data from these studies in independent geomechanical modeling to build a unique classification.
The minimum horizontal stress gradient of 17.4 kPa/m presented herein is aligned with these previous studies. de Pater and Baisch (2011) reported an average minimum horizontal stress gradient of 0.74-0.8 psi/ft (17.0-18.1 kPa/m) from minifrac pressure declines within the Bowland interval at well PH-1. Clarke et al. (2018) presented a minimum horizontal stress curve for which a gradient of 17.8 kPa/m (0.79 psi/ft) can be determined to calculate a maximum value of 49.3 MPa.

Uncertainty analysis
While care has been taken to minimize any uncertainty surrounding the results of this work, some key elements remain that need to be considered. The study relies heavily on wireline log datasets and the accuracy of the bulk density log in particular within enlarged (usually clay-rich) sections of the well is considered a key uncertainty. In such sections, the tool will record a lower response than expected, and this error will propagate into our calculated elastic properties (i.e., the interval would appear weaker than it should). As this will impact the most clay-rich intervals, it is unlikely that this jeopardizes the choice of landing zones, but it may result in an overestimation of the effective stress barriers between these.
Furthermore, the Eaton method for pore pressure prediction is known to perform best in young (Mesozoic-Cenozoic) sedimentary basins and does not consider the effects of unloading (Zhang, 2011). By contrast, the Craven Basin is of Paleozoic-fill and known to have undergone multiple uplift events, but the modeled curves herein for pore pressure did honor the calibration data available.
However, there is also some uncertainty regarding the accuracy of the stress calibration data itself. The results from de Pater and Pellicer (2011)'s analysis of minifrac pressure declines were taken for values of fracture closure pressure (see Table 2). However, those authors state that fracture closure pressure is "hard to determine uniquely from pressure decline" and that additional flow-back and step-rate tests would reduce this uncertainty (de Pater and Pellicer, 2011).

Conclusions
This paper presents new findings around the geomechanical stratigraphy of the Bowland Shale in northern England, where it forms a prospective shale gas target. Wireline logs from the Preese Hall-1 exploration well, located on the Fylde peninsula, Lancashire, were used in a series of geomechanical models. The results of mini-frac tests were used to calibrate some of these models where appropriate. Understanding the geomechanical characteristics of the Bowland Shale is important for assessing the shale's resource potential and this work addressed the implications for developing a well placement strategy, specifically.
The shale exhibits strong heterogeneity in geomechanical properties, together with a high degree of structuring, which can make it challenging to determine the ideal zones for horizontal drilling using methods commonly applied within USA shale operations. This was addressed using an unsupervised clustering algorithm that identified intervals of shared geomechanical characteristics. The intervals with the lowest effective stress, highest brittleness index, and lowest fracture toughness were interpreted as the intervals with the best qualities for hydraulic fracture stimulation and three potential landing zones were highlighted that exhibit high concentrations of such intervals. Importantly, these three landing zones are separated by intervals of high effective stress, low brittleness index, and high fracture toughness which could impede vertical fracture growth from one landing to another; something which could be tested with hydraulic fracture simulations.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. academic license to HWU. Jingsheng Ma acknowledges NERC grant number NE/R018022/1 for financial support.

Appendix A: Method for Defining Reservoir Zones
A series of petrophysical models were constructed to evaluate the presence of Reservoir Quality (RQ) intervals in well PH-1. The ELAN model (Quirein et al., 1986), packaged within Schlumberger's Techlog application, was used to invert mineralogy (clay, quartz, and carbonate) and porosity logs from the Neutron Porosity, Bulk Density, and Sonic Slowness wireline logs. Total Organic Carbon (TOC) was calculated by establishing a linear correlation between the Uranium Concentration wireline log and core-derived TOC measurements. Water saturation was estimated using the Indonesia method (Poupon and Leveaux, 1971).
RQ intervals were then computed using a series of cut-offs applied to this modeled dataset. The maximum cut-off for clay concentration is taken at 50 %, following Bowker (2007)'s suggestion that this forms the maximum concentration for prospective Barnett Shale intervals. A minimum porosity of 4 % was taken following Jarvie (2012) and a minimum TOC of 2 % was taken following Andrews (2013). The cut-off for water saturation was taken as its median value (Table A1).
A total thickness of 125 m was established following this cut-off scheme and their distribution is shown in track 6 of Fig. 13.