Determining reservoir intervals in the Bowland Shale using petrophysics and rock physics models

Finally, we consider the density of natural fractures in the shale by developing an anisotropic rock physics model to reﬂect high-angle fractures observed on micro-imagery logs. We invert crack density using shear wave splitting well log data and ﬁnd a crack density of up to 4 per cent which correlates well with micro-imagery observations. Our work further supports previous authors’ core-based studies in concluding that the Bowland Shale holds good reservoir characteristics, and we propose that there are multiple intervals within the shale that could be targeted with stacked horizontal wells, should those intervals’ mechanical properties also be suitable and there be adequate stress barriers between to restrict vertical hydraulic fracture growth. Finally, our rock physics templates may provide useful tools in interpreting pre-stack seismic data sets in prospective areas of the Bowland Shale and picking the best locations for drilling wells.

United Kingdom (Smith et al. 2010), though the size of the potential resource is disputed (Andrews 2013;Whitelaw et al. 2019). Exploration thus far has focused within the Craven Basin of NW England, where the Bowland Shale exhibits excellent shale reservoir properties (Clarke et al. 2014a(Clarke et al. , 2018 and has been subject to several sedimentological and microstructural studies (Fauchille et al. 2017;Newport et al. 2018;Emmings et al. 2019Emmings et al. , 2020a. However, the geological structure of the Craven Basin is complex, and the shale resource is compartmentalized by high-angle reverse faults into narrow blocks (Anderson & Underhill 2020). This compartmentalization suggests that long (>2 km) horizontal wells (as are traditionally used in shale gas plays) may not be possible in this area. Furthermore, all hydraulic fracturing attempts have reactivated either seismic-scale or subseismic-scale faults, leading to elevated levels of induced seismicity and premature halting of operations (Clarke et al. 2014b(Clarke et al. , 2019bKettley & Verdon 2021;.
This faulting may necessitate the use of shorter horizontal wells targeting different stratigraphic intervals ('stacked' lateral wells) to produce commercial volumes of gas from the Bowland Shale. Such a technique would ensure that as much shale could be stimulated as possible, whilst optimizing the horizontal well trajectories away from mapped faults. Production from the Bowland Shale using stacked lateral wells has been briefly mentioned in the literature (Clarke et al. 2014a(Clarke et al. , 2018 but without detailed analysis detailing if multiple prospective intervals exist in the shale and how they may be targeted. Selecting the best intervals to target with horizontal wells ultimately requires knowledge of the reservoir quality and geomechanical characteristics of the shale formation. In this work, we address the issue of selecting the best reservoir quality (RQ) intervals for stacked drilling at the well Preese Hall-1 (PH-1) by using petrophysical and rock physics analysis. We focus specifically on characterizing the Bowland Shale's mineralogy, organic content, porosity, permeability and water saturation from wireline logs. We then implement two rock physics models; first, an isotropic model to assess the link between rock properties and isotropic elastic properties, and secondly an anisotropic model to invert natural fracture density from shear splitting logs. We then conclude by presenting a classification of the intervals that hold the best RQ and consider how these might be predicted from seismic data.

PH-1 lithostratigraphy and wireline character
The PH-1 well was drilled by Cuadrilla Resources between 16th August and 3rd December 2010 in the Fylde region of Lancashire, NW England (Fig. 1). As it was designed with the purpose of shale gas evaluation, several specialized wireline logging runs (Table 1) and core studies (Table 2) were undertaken, and these data sets were accessed for this study via the British Geological Survey. It encountered the Bowland Shale at 1980 m TVD-MSL (True Vertical Depth below Mean Sea Level), and drilled through 725 m of the formation before reaching total depth at 2730 m TVD-MSL.
The shallowest wireline logging at PH-1 is within the Sherwood Sandstone (Fig. 2). The Manchester Marl is characterized by very low Gamma Ray (GR), low Bulk Density (RHOB) and high Deep Resistivity (DDLL) zones indicative of evaporites. Very low RHOB values and high Compressional (DTC) and Shear (DTS) Sonic Slowness readings within the Millstone Grit are accompanied by high Caliper values indicating washout within the associated shales (not shown in Fig. 2). The Upper Bowland Shale (UBS) displays clear heterogeneity shown by the varying GR response. An increase in Uranium Concentration (GRUR) near the top of the UBS suggests an enrichment in organic matter. Further washout sections evidenced by very low RHOB values exist but are confined to thin intervals. The Lower Bowland Shale (LBS) displays a different character; one of alternating thicker high GR units interbedded with low GR units. The low GR units belong to the Pendleside Sandstone Member, and while these are relatively thin in the PH-1 well, they can reach up to 300 m thickness in the Clitheroe area (Clarke et al. 2018). Three of the LBS high GR units (shales) were penetrated before the well reached total depth.

Lithological heterogeneity in the Bowland Shale
The wireline log signature of the Bowland Shale is one characterized by variability in rock properties and this is owed, in part, to the mineralogical heterogeneity of the formation. The heterogeneity observed can be related to the complex and dynamic relationship between tectonic and sedimentary processes active during deposition of the shale (the Visean and Namurian periods). The LBS is a hemipelagic mudstone that accumulated during the Visean within several isolated syn-rift basins separated by shallow platforms on which extensive shelf carbonates were deposited (Fraser & Gawthorpe 1990;Corfield et al. 1996). During this time, the carbonate platforms were flanked by reefs that effectively cut the supply of detrital carbonate into the basins, leading to the development of a dysaerobic, stratified water column (Kirby et al. 2000). As a result, LBS mudstones are typically calcareous, locally organic-rich and predominantly reflect low-energy turbidite deposition (Newport et al. 2018). The thick low GR units observed at the base of PH-1 are coarse sandstones of the Pendleside Sandstone Member (Fig. 2) (Clarke et al. 2018) which represent higher-energy turbidites, possibly sourced from coastal deltas to the north (Kirby et al. 2000).
The boundary between the LBS and UBS is taken at the Visean-Namurian boundary which also marks a change in regional tectonics from syn-rift to post-rift subsidence (Fraser & Gawthorpe 1990). Throughout the Namurian period, the deltas already established in the Midland Valley of Scotland and NE England prograded southwards and their influence was reflected in the progressive infilling of the previous syn-rift depocentres (including that of the Craven Basin) with increasingly coarse clastic material Fraser & Gawthorpe 2003). In PH-1, the Namurian UBS consists predominantly of mudstone with varying silt content, thin sandstones and thin dolomitized limestones (Clarke et al. 2018). Marine transgressions ('marine bands') are a prominent feature and are visible on the wireline logs as thin, high GR intervals. Emmings et al. (2020a) observed these features in adjacent boreholes and outcrops confirming these as pelagic/hemipelagic carbonate and fossil-rich mudstones.
This macroscale heterogeneity in lithology types described above also manifests at microscale. A range of microtextures (relating to the degree of lamination and/or predominant mineralogy), organic matter types and fracture styles have been characterized within samples taken from PH-1 cores (Fauchille et al. 2017(Fauchille et al. , 2018 and those authors suggest that depositional variability is a strong control on the heterogeneity observed. Taking these multiscale observations of rock character together, it is apparent that the Bowland Shale is a mineralogically and texturally complex formation. It is important to define how this heterogeneity influences reservoir and geophysical Downloaded from https://academic.oup.com/gji/article/228/1/39/6355446 by British Geological Survey Keyworth user on 27 September 2021 Figure 1. Two-way-time structure map for the Lower Bowland Shale on the Fylde region highlighting the Preese Hall-1 (PH-1) well, additional wells and interpreted faults (after Anderson & Underhill 2020). GH-1Z, Grange Hill-1Z; Th-1, Thistleton-1; PNR, Preston New Road-1; El-1, Elswick-1; RW, Roseacre Wood.  properties as this will in turn influence the placement of horizontal wells.

Seismic anisotropy in the Bowland Shale
An important phenomenon when considering the geophysical behaviour of shale is anisotropy in elastic properties. Seismic anisotropy can be most simply defined as the 'dependence of seismic velocity upon angle' (Thomsen 2002). The phenomenon is caused by the preferential alignment of minerals, pores, or fractures. It is particularly acute in mudstones as they are deposited in low-energy settings (giving rise to sedimentary processes that produce laminated textures), significantly compacted during burial (accentuating alignment of pores and minerals) and often extensively fractured at multiple scales (such as exhibiting microfractures due to hydrocarbon generation, or macrofractures due to burial diagenesis). A material with anisotropic elastic properties is theoretically represented in Hooke's Law by a complex, fourth-rank elastic stiffness tensor of 81 independent parameters (Mavko et al. 2009). However, in practice, a plane of symmetry is often imposed, allowing the tensor to be significantly reduced to a manageable five constants. Therefore, shale is commonly referred to as transversely isotropic with either a vertical plane of symmetry (VTI, i.e. fabric is bedding-parallel) or transversely isotropic with a horizontal plane of symmetry (HTI, i.e. fabric is bedding-normal). Voigt notation is commonly used to represent the elastic stiffness tensor (c ij ), which for a transversely isotropic (TI) medium is (Mavko et al. 2009): While there is limited data to confirm that seismic anisotropy may be a strong feature of the Bowland Shale, aligned fabrics within the shale have been documented in some studies. One such fabric is the presence of natural fractures, which are a common feature to shale plays  and are prevalent in the PH-1 well. Some of these are open, but many are mineral (usually calcite)-filled and are usually oriented perpendicular or at high angles to bedding (Fauchille et al. 2017;Clarke et al. 2018), suggesting HTI behaviour. Clarke et al. (2018) demonstrate how these can be observed on core samples (Fig. 3b) and thin sections, whereas Fauchille et al. (2017) identify them using SEM images (Fig. 3c). These fractures are also of sufficient aperture to be imaged on micro-imagery log data and can be identified as the presence of either bright (i.e. resistive, healed) or dim (i.e. conductive, open, fluid-filled) sinusoids usually with amplitudes noticeably greater than the bedding planes (Figs 3a and d).
Furthermore, there is some evidence that these sub-vertical features correspond to shear-wave anisotropy. This anisotropy is characterized by the cross-dipole sonic tool, which measures the shearwave splitting (SWS) effect; where the shear-wave is polarized into fast and slow shear waves upon entering an anisotropic medium (Liu & Martinez 2012). As this polarization can often be attributed to the presence of aligned fractures in the logged formation (Esmersoy et al. 1995), it has become a popular tool for unconventional reservoir characterization. SWS can also be used to determine fracture sets from reflection seismic (Davis et al. 2003;Shuck et al. 1996), vertical seismic profiles (Winterstein et al. 2001) and microseismic events (Verdon et al. 2009;Verdon & Kendall 2011).
While shear anisotropy is generally low at PH-1, locally high values can be observed in areas where subvertical fractures and/or wider fracture zones are present (Fig. 4). These manifest as high amplitude sinusoids on the micro-imagery log and chaotic zones without clear character (Figs 4b and c). However, the greatest shear wave anisotropy readings appear to correspond to zones of induced, rather than natural fractures, which produce symmetrical, boreholeparallel high conductivity features that cut across bedding (Fig. 4a). In each case, the presence of subvertical features approximately normal to bedding can be taken as evidence of HTI behaviour, which is then quantified by the shear wave anisotropy log. However, the impact of drilling-related features can clearly overprint the anisotropy response caused by natural fractures and thus these zones need to be treated with caution.
As discussed previously, VTI anisotropy is also known to be prevalent in shales due to the presence of bedding-parallel features (e.g. laminated clays or aligned kerogen and/or pores). However, there is no clear evidence to suggest that this is a significant phenomenon in the Bowland Shale. Clarke et al. (2019a) show how micro-imagery log data can be used to identify strong layering of shales and carbonates but also conclude that both the anisotropy wireline log data and the ultrasonic velocity experiments demonstrate negligible anisotropy in elastic properties. Fauchille et al. (2017) highlighted the presence of a 'strongly laminated' facies in the shale and identified bedding-parallel microcracks filled with bitumen which all may influence seismic anisotropy, but at present, there is no data to quantify the degree of influence.
Further work is needed to determine the effect of VTI fabrics on the anisotropy of the Bowland Shale. Ultrasonic velocity experiments conducted on core samples oriented at different angles to bedding/lamination would be required to fully quantify the degree of TI anisotropy. For an accurate characterization, measurements on three core plugs (oriented normal to, 45 • to and parallel to bedding) are required (Vernik & Nur 1992). However, the preparation of such samples (particularly those at 45 • to bedding) can be challenging due to shale's friability, which led Wang (2002) to propose a method whereby transducers could be placed at different orientations to bedding on a single sample. However, both these methods require access to core samples of sufficient size to obtain a large plug for the velocity experiments. The friable nature of shale, combined with the presence of natural fractures often leads to poor core recovery and makes the collection of such samples challenging.   The upscaling of such studies to wireline log and/or seismic data also poses an additional challenge. Accurately measuring the degree of VTI anisotropy and gaining an understanding of the inherent causes of it can better inform modelling work at the larger scale (i.e. wireline log or seismic scale), but upscaling these effects can be troublesome. Laboratory measurements can only be performed on either intact samples or samples with healed (i.e. mineralized) fractures that may not be representative of the fractured rocks at depth. Furthermore, the relative contributions of microscale features that affect VTI anisotropy to anisotropy measurements acquired at the broader scale (e.g. shear wave splitting) is uncertain.

Effective medium models
In this work, a rock physics models are built to characterize the elastic properties of the Bowland Shale, the foundations of which are built on effective medium theory (EMT). EMT allows for the mixing of multiple constituents (i.e. minerals, pores or fluids) considering their volume fraction, individual elastic properties and (generalized) shape. The background theory to many of the effective medium models used in this work can be found in Mavko et al. (2009) and recent reviews such as Zhang (2017), Zhang et al. (2017), Guo & Liu (2018) and Dvorkin et al. (2021). We outline some of the theory pertinent to this work within the Appendix.

Rock properties
Full characterization of a shale gas reservoir requires the integration of various scales of data and the degree of uncertainty associated with these data sets increases with increasing scale of coverage. Core studies provide accurate quantification of mineralogy, total organic content (TOC), porosity, permeability and saturation however sample sizes are often limited. Wireline logging provides continuous measurements of various physical properties every 0.1524 m (though actual bed resolution varies by tool type) however these often require transformation into relevant rock properties. Seismic data provides 3-D measurements of reflectivity at 10s-m resolution however well logs are required to produce a meaningful inversion and generate predictions of rock or reservoir properties. The first step in performing a petrophysical analysis of the RQ at PH-1 and generating data for the rock physics model is to calculate continuous logs of mineralogy, TOC, porosity, and saturation. Throughout the analysis, attempts were made to honour any information from PH-1 core reports, including mineralogical data from XRD tests, porosity, permeability and saturation data from Mercury immersion and Helium pycnometry tests and TOC data from Source Rock Analyzer (SRA) pyrolysis (Table 2).

Mineralogy logs-elemental log analysis (ELAN)
Determining complex mineralogy from wireline logs is a challenging task in the shale reservoirs. Spectroscopy logging provides quantitative lithology measurements derived from direct measurements of elemental yields within the formation (e.g. Elemental Capture Spectroscopy (Herron & Herron 1996) but may not be included as part of a standard wireline logging suite, and was not run at PH-1 (Table 1). Determining quantitative lithology using traditional logging suites is achieved more often using a multimineral inversion approach, such as that followed by Schlumberger's Elemental Log Analysis program (ELAN, Quirein et al. 1986). The ELAN program is packaged within Schlumberger's Techlog petrophysical software and was accessed for this work. Analysis of wireline logs is the inverse of the forward problem stated as (Quirein et al. 1986): where t is the tool vector consisting of the log readings at a certain depth, R is the response matrix of constant parameters and v is the volume component vector. The aim of multimineral inversion is determining v, given a series of tool measurements, t and a user-defined response matrix, R. The inversion approach (Mayer & Sibbit 1980) involves iteratively solving the above forward problem, minimizing the error between reconstructed tool response and actual tool response (Quirein et al. 1986). If a fixed response matrix is assumed, the volume of each component can be inverted for in this approach. Prior knowledge of the typical mineralogical assemblages within the formation is required, as should an understanding of the appropriate logs to be used as input. The program can invert for the same number of mineralogical components as input log measurements plus one (Quirein et al. 1986). The tool response for each mineral Table 3. Default response parameters for the three log inputs to the ELAN model for the three minerals to be inverted. The neutron porosity log is calibrated to limestone, which results in a zero-porosity value for pure limestone and a negative reading for quartz. must also be known (e.g. the RHOB log in a 100 per cent quartz formation will read 2.65 g cm -3 ). In practice, response parameters would be extracted from a database using adjacent wells and/or laboratory studies. While ideally, a bespoke set of response parameters would be determined for the Bowland Shale using wireline logs (including elemental logs 'Elemental Capture Spectroscopy') and core analysis from multiple wells, for this work, as only one well was available, with some limited core data, the response parameters were kept at program defaults (Schlumberger 2019, Table 3). Three wireline logs that are strongly influenced by lithology were selected for use in the ELAN model: NPOR, RHOB and DTC. Throughout the inversion, reconstructed tool measurements were compared with actual measurements to ensure errors were kept to a minimum. Log outputs were also compared with core XRD (Table 2) information to ensure reasonable results (Fig. 5).

TOC and kerogen concentration
Several approaches exist to determining TOC and Kerogen concentrations from wireline logs. The most popular methods include developing empirical correlations with laboratory TOC measurements (Fertl & Chilingar 1988) or the analysis of the separation between porosity-based and fluid-based logs ( logR, Passey et al. 1990). As this latter approach is largely based on the elevated resistivity response due to liquid hydrocarbons it is not always transferable to overmature gas shales and therefore is not followed in this work.
Of the different correlations proposed between wireline logs and TOC, one of the strongest involves Uranium concentration (GRUR; Fig. 2). As authigenic Uranium precipitates at the sediment-water interface under anoxic conditions and accumulates along with organic matter, Uranium and TOC have been shown to strongly correlate (Supernaw et al. 1978;Lüning & Kolonic 2003). In the Bowland Shale, we observe high Uranium concentrations, particularly within the UBS, up to 10 ppm and corresponding to total GR readings of ∼200 API (Fig. 2). Using this log data and TOC measurements available from core data (Table 2), a linear correlation is established: where TOC is measured as a fraction and Ur is in ppm. TOC was then converted to Kerogen using the approach of Herron & Le Tendre (1990) and Vernik & Nur (1992): where C is a TOC to Kerogen conversion factor. While the maturity of the Bowland Shale at PH-1 makes it difficult to infer original kerogen type (Clarke et al. 2018), it is generally believed to consist of a mixture of Types II (marine) and III (terrestrial) (Andrews 2013;Slowakiewicz et al. 2015;Hennissen et al. 2017;Clarke et al. 2018). For these kerogen types, situated in the gas window, Tissot & Welte (1978) suggest a conversion factor around 1.2 is suitable, which we follow in our calculations. ρ b is the RHOB log reading and ρ k is the density of kerogen. The density of kerogen will vary with maturity, which itself varies with depth of the Bowland Shale in the PH-1 well. Clarke et al. (2018) demonstrate vitrinite reflectance (R o ) to vary between 1.5 and 2.25 (corresponding to the early-wet gas to early-dry gas window). For this maturity range, kerogen densities of around 1.4-1.5 g cm -3 are reasonable (Alfred & Vernik 2013), and in our calculations we set the parameter to 1.4 g cm -3 . Kerogen concentration is plotted alongside measurements from core data (calculated from TOC using the same approach) for comparison in Fig. 5 (track 8).

Water saturation
Estimation of water saturation is a challenging task in shales with many of the well-adopted methods developed on clean or shaly reservoir sands. While this thesis does not address the calculation of water saturation in detail, a method considered most applicable in shaly settings [the Indonesia method proposed by Poupon & Leveaux (1971)] was picked to determine an approximation of saturation. The equation uses formation resistivity, the volume of shale/clay and porosity as follows: where V sh is shale volume (taken as clay volume from the ELAN model), ϕ is porosity and R t is formation resistivity. Picking appropriate values for tortuosity, cementation, and saturation exponents (a, m and n) can be a challenge. Values of 1, 2 and 2, are widely used in the conventional setting, however, these are appropriate only for clay-poor rocks with well-connected pores. Some authors have suggested that the cementation exponent is the fundamental parameter that should be adjusted in the case of shales, which led Malekimostaghim et al. (2019) to conduct a laboratory study to determine a reasonable value for shale. We use their determined value of 1.6 in our calculations. The resistivity of water (R w ) and of shale (R sh ) are fixed at 0.03 and 5, respectively. This method provides a reasonable fit to the lower values of core-measured water saturation ( Fig. 5; track 9) but does not accurately capture the higher core-derived values.

Permeability
The relationship of permeability and porosity is complex, and reservoir zones of similar porosity can exhibit markedly different permeabilities. This has led some authors to recognize the presence of 'flow units' using cross-plots of core-derived porosity and permeability data (Corbett & Potter 2004;Aguilera 2014). In the absence of substantial core test data, permeability models become more uncertain, though fitting an exponential trend is often considered a reasonable approximation (Halliburton 2019).  Some limited effective porosity and matrix permeability core test results are available for PH-1 (Table 2), and while they were not sufficient in number to carry out a detailed rock typing exercise, they did allow the building of an approximate permeability model (Fig. 6). Note that the core-derived matrix permeability values are very low, ranging from 2 to 50 nD. This exponential function was used to estimate a continuous permeability from the ELAN porosity. The permeability curve was used in the final classification but not in the rock physics model itself.

Rock physics model
In building a rock physics model for the Bowland Shale, two broad aims were pursued. First, to account for the lithological heterogeneity described in previous sections, while honouring observations made regarding the fabric of the shale, and to investigate its impact on elastic properties. Secondly, given the observations of natural fractures within the shale and the effects on shear wave anisotropy, to invert the Hudson crack density parameter (a metric that can be seen as a proxy for the density of natural fractures) and quantify the degree of anisotropy associated with it.

Isotropic model
A rock physics model essentially seeks to transform rock properties (such as lithology, porosity, and fluid saturation) to elastic properties (Grana 2014). In shales, the task is challenging due to the presence of complex microstructural fabrics, mineralogical variability, and the presence of microcracks and fractures which all, in turn, significantly affect the elastic properties of the rock. Therefore, rock physics models for shales usually involve multiple modelling steps and utilize several numerical algorithms to best approximate the rock's complex physical properties and account for anisotropy. The rock physics model implemented in this work sought to encompass the effects of porosity, kerogen and mineralogy on elastic properties and involves several steps outlined as follows: (1) Elastic constants for mixed clay and mixed carbonate are calculated using the Hashin-Strikman average. Clay is modelled as an illite/kaolinite mixture with the proportion of each mineral taken from average values obtained from core XRD data. Likewise, carbonate is modelled as a calcite/dolomite mixture. The proportions of each mineral within its mixed group (i.e. clay and carbonate) do not change within the model.
(2) The pore system within the Bowland Shale is mainly within or associated with organic matter (Ma et al. 2016). To represent this, a kerogen/gas mixture is created using the approach followed by Carcione (2000) whereby the mixture consists of hydrocarbon bubbles embedded within a kerogen matrix. Carcione (2000) first used Wood (1955)'s equation to create a water/oil mixture, which we follow but consider gas rather than oil. We take the water saturation log as the water fraction and assume the hydrocarbon saturation to be composed entirely of gas. Secondly, Carcione (2000) used Kuster & Toksöz (1974)'s model to embed spherical fluid (water and gas) pores into the kerogen matrix. This approach was also followed by Guo et al. (2013) in their model for the Barnett Shale. Following Guo et al. (2013), we assume that the percentage of the total porosity allocated to kerogen is fixed at 20 per cent throughout the shale. This will likely vary with maturity, but given their work focused on the Barnett Shale, we consider it a reasonable general value for gas mature shales. The remaining 80 per cent of total porosity is then added as isolated, fluid-filled (also water and gas) pores in step 3.
(3) Berryman (1980)'s SCA model is then used to create a mixture of the solid minerals (quartz, carbonate and clay), porous kerogen and remaining isolated pores. The fractions of each constituent are taken from the appropriate input log. We consider quartz and carbonate as spherical (aspect ratio of 1), and the kerogen, clays and pores as ellipsoidal in shape (aspect ratio of 0.1).

Anisotropic model
Anisotropy is considered within the rock physics model by adding vertical, aligned, water-filled cracks into the rock to simulate the effect of natural, sub-vertical fractures within the formation. The shear wave splitting log forms the basis for the inversion, with the aim being to quantify the degree of vertical/sub-vertical fractures observed on the micro-imagery log data by calibrating an anisotropic rock physics model to shear anisotropy observed from the cross-dipole sonic logging data. However, as was highlighted previously, there is evidence from the micro-imagery log that zones of drilling-induced tensile fractures are also contributing to the shear wave splitting behaviour within the shale. Attempt was not made to decouple the effects of natural and drilling-induced fractures on the shear wave anisotropy log, though effective medium techniques have proven efficient in distinguishing the two phenomena (Prioul et al. 2007). Instead, a prominent zone of intense drilling-induced fractures between 2127 and 2170 m (measured depth) was eliminated entirely as to avoid it dominating the inversion calculations. Fig. 4 (left) showed the micro-imagery log response and associate shear wave anisotropy character for a subsection of this zone.
To ensure that the inverted crack porosity and non-crack porosity (from the ELAN model) do not exceed the total porosity of the rock, our inversion scheme begins with initial estimates of crack density (per unit volume). We assume that crack density (ρ c ) is linearly related to shear wave anisotropy (γ ) as: where k is a factor that converts shear wave anisotropy (expressed as a per cent) to crack density (expressed as a fraction). According to Kachanov & Mishakin (2019, their eq. 11), crack porosity (ϕ c ) can then be estimated from crack density as: where α is the aspect ratio of the cracks. Once crack porosity is estimated, an isotropic rock is then constructed according to the method outlined in the previous section. The porosity added in the isotropic model is adjusted by subtracting the crack porosity from the total porosity to ensure that when cracks are added in subsequent modelling, the total porosity is not exceeded. Following this, the Hudson model is used to add vertical cracks, with the shape of the cracks fixed as very narrow ellipsoids (aspect ratio of 0.1).
The Hudson model provides an output of the full anisotropic stiffness tensor for that crack density. We use this to calculate equivalent slowness logs by assuming transverse isotropy with a horizontal plane of symmetry (HTI) whereby the fast shear wave (DTS fast ) propagates normal to bedding and the slow shear wave (DTS slow ) propagates parallel to bedding. Thus, the modelled shear wave anisotropy becomes: Therefore, a log of modelled shear wave anisotropy for a given k (eq. 8) is determined. The quality of fit between this modelled curve and the actual anisotropy log is assessed by calculating the mean squared error (MSE), and the process is repeated for different values of k (ranging from 0.002 and 0.011). Using this approach, the resulting MSE ranged between 6.2 and 0.03 (Fig. 7a), with the latter corresponding to a k of 0.08. At k values above this, the modelled overpredicts the degree of anisotropy and at k values below this, it underpredicts (Fig. 7b).

Reservoir properties
The results of the multimineral inversion using ELAN, reveal that quartz is the dominant mineral within the Bowland Shale (Fig. 8).
Most of the shale contains quartz concentrations between 50 and 80 per cent with means of 59.4 and 52.0 per cent for the UBS and LBS, respectively. Clay and carbonate concentrations are lower, with means around 24.5 and 20 per cent, respectively though clay concentration displays a normal distribution around the mean, whereas the carbonate distribution is skewed towards the lower values. The mean carbonate content in the LBS is greater (25.5 per cent) than in the UBS (18.0 per cent), and while carbonate content is generally low, there are some intervals with very high concentrations up to 80-90 per cent.
TOC content is largely between 1 and 4 per cent and is noticeably higher in the UBS (2.7 per cent mean) than the LBS (1.5 per cent mean). Total porosity is mainly less than 6 per cent with means of 4.2 and 3.8 per cent for the UBS and LBS, respectively. Lastly, our modelled water saturation is almost exclusively less than 30 per cent, which with means of 11.2 and 10.3 per cent for the UBS and LBS, respectively. These values appear quite low and fail to capture the higher water saturation core measurements as illustrated on Fig. 5  (track 9). The values of 100 per cent in this instance correspond to zero porosity sections of the shale.

Validation of isotropic model
To assess the durability of the isotropic rock physics model, compressional (P-wave) velocity (V P ) and shear (S-wave) velocity (V S ) curves were modelled for a section near the base of the UBS and checked against the actual V P and V S logs (Fig. 9). Both modelled curves show a good match with the actual data. The V P model displays the strongest match with a correlation coefficient (R) of 0.89 for this interval, an intercept close to zero and a gradient near one (Fig. 9c). A small number of low V P values are not captured in the model (e.g. at 2300 m, also evident as datapoints positioned far above the linear regression trend), but as these are few, the overall performance of the V P model was considered adequate. The V S model shows a slightly weaker match, with an R-value of 0.77 for this interval. A systematic, slight overprediction is observed as evidenced by the linear regression fit positioning further up the y-axis ('Modelled V S ') than the Y = X trend in Fig. 9(d).
The poorer correlation could be partially explained by differences between the log types. The shear velocity log is derived from the cross-dipole sonic tool which has a typical bed resolution of ∼1.4 m. However, the input logs to the ELAN model form part of a suite of tools (termed quad combo) most of which can achieve bed resolutions ∼30 cm. While all wireline logs are presented with a consistent 0.1524 m sampling rate, the bed resolution (i.e. the thinnest unit that can be identified) achieved by each log is variable. Coupling these differences in resolution with uncertainties regarding the depth-matching of separate logging runs could partly explain the poorer relationship observed in the V S model. Nevertheless, in the absence of additional information to assess this, the correlations were considered reasonable, and sufficient for use in the rock physics templates and the subsequent inversion.

Rock physics templates
Following assessment of rock physics model's performance, crossplots were produced to compare the measured elastic properties from wireline logs at PH-1 with the predicted elastic properties from the rock physics model and to provide some insight into the link between the Bowland Shale's rock and reservoir properties and geophysical data. Fig. 10 presents these cross-plots with a focus on the lithological character of the shale. The wireline log data is plotted as scattered data coloured by lithology classification (further outlined in Fig. 11). The rock physics model is presented as a pseudo-ternary diagram in rock physics space. A full range of clay, carbonate and quartz concentrations are iterated over; however, kerogen concentration and water saturation are fixed at their median values (3.8 and 9.8 per cent, respectively). Further argument around the median values is provided in the Discussion section. Two models of fixed porosity are considered; a zero-porosity model represented as black lines with colour-fill, and a 10 per cent porosity model represented as grey lines with no fill.
The colour-fill used in both the scattered data and the model is designed to highlight certain lithologies using a simplified version of the sCore scheme for shale lithology classification each data point using the Rickman et al. (2008) method where E min and ν max (ductile rock) were set to 26.1 GPa and 0.40, respectively and E max and ν min (brittle rock) were set to 67.6 GPa and 0.08, respectively. The results of this calculation were used to interpret zones of brittle and ductile zones elastic properties on each crossplot.
The most clay-rich sections of the Bowland Shale are easily differentiated from the rest of the shale (clay-dominated and argillaceous lithologies shown in green and grey, Fig. 10). These correspond to the ductile sections of the rock physics plots and are most easily identifiable in the λρ versus μρ (Fig. 10a) and ν versus E crossplots (Fig. 10d). They are characterized by the lowest E values (i.e. low strength; less stress is required to induce strain in the rock), highest ν values (i.e. low strength; high transverse strain upon axial loading) and lowest AI values (i.e. low resistance to a propagating compressional wave). V P to V S ratio alone does not appear a useful indicator for this lithology, as it spans a range of values between 1.7 and 2.1 (Fig. 10c) and displaying some overlap with the carbonate regions of the shale.
The predominant lithology is siliceous mudstone which displays greater rigidity (μρ > 40 GPa g cm -3 ) than the argillaceous mudstone and equal, to slightly lower incompressibility (λρ, Fig. 10a). It exhibits the lowest V P to V S ratio (<1.7, Fig. 10b), ν (<0.2, Fig. 10d) and moderate AI (between 11 and 14 km s -1 · g cm -3 , Fig. 10b). There is little change in E between this lithology and mixed or calcareous mudstones (Fig. 10d). Both lithologies also display high brittleness though there is a small preference for the highest brittleness values to be located towards the quartz section of the plot rather than the calcite section.
With increasing carbonate content (the transition from siliceous to mixed to carbonate mudstone), a clear increase in λρ is observed (up to 110 GPa g cm -3 ), accompanied by a subtle increase in μρ (Fig. 10a). An increase is observed for both AI (from 14 to 16 km s -1 · g cm -3 ) and V P to V S ratio (from 1.7 to 1.9, Fig. 10b), though as previously discussed, this trend in V P to V S ratio is also observed for the argillaceous mudstone. It is also characterized by high ν and E (Fig. 10d). This high Poisson's ratio translates to a slight decrease in brittleness, relative to the siliceous mudstone lithology.
We also observe errors between the template lithology and lithology of the data points (for example, some mixed mudstone data points fall within the calcareous mudstone template zone (Fig. 10a)). This could be explained by the fact that porosity, kerogen content, water saturation and some textural parameters (e.g. aspect ratios of minerals) are all fixed in the generation of the model template.
Fluctuations in these parameters will cause changes in elastic properties and affect the shape of the triangular mesh on the charts that represents the rock physics model. However, when producing such a generalized template this is difficult to avoid. Producing multiple charts with multiple distributions of these fixed parameters is possible, but in doing so, the predictive power of them is lost and it becomes harder to interpret the results in a general sense.
On the whole, the results suggest that lithology has an acute control on the elastic properties of the shale. The shale exhibits predominantly brittle characteristics as evidenced by the relationship between the well log data points and the generalized brittleness classification. The most brittle sections of the shale appear to lie within the mixed mudstone lithology. As quartz content is often considered the main contributor to shale brittleness, this is an unusual observation. However, these siliceous mudstone units also display moderate porosity, which would serve to reduce the rigidity of the rock. Brittleness appears to reach a maximum within low porosity, moderately calcareous mudstones, but then falls with further carbonate content. Carbonate minerals such as calcite and dolomite exhibit greater Poisson's ratio (∼0.30) than quartz (0.08, Mavko et al. 2009), which may explain the reduction in brittleness observed for carbonate-rich intervals when Rickman et al. (2008)'s approach is used. This behaviour was also observed by Perez Altamar & Marfurt (2014) in their Barnett Shale study and does suggest that interpreting siliceous from carbonate-rich intervals using elastic property cross-plots should ideally be checked against available mineralogical data.

Crack density inversion
The inverted crack density is predominantly below 0.04 but does exceed 0.06 in some intervals (e.g. 2200, 2260 and 2285 m, Fig. 12). It appears to show an inverse relationship with total porosity and clay fraction. In some intervals (e.g. those quoted above), high crack density appears to be associated with high carbonate fractions, however this relationship is less apparent at greater depths (e.g. 2540 m) where high carbonate fractions correspond to moderate crack density. On this plot, the zone of intense drilling-induced fractures (DIFs) that was omitted from the inversion is highlighted in red ( Fig. 12). On the basis of the above observations, given this interval has moderate-high porosity and relatively low carbonate concentrations, we would expect this interval to display low-moderate crack densities, however this cannot be confidently tested. Rock physics templates were again produced but in this instance with data points coloured by inverted crack density (Fig. 13). The relationship between elastic properties and crack density is not an obvious one, but they do clearly show a preference for plotting within the brittle sections of the cross-plot (see Fig. 10 for the full brittleness template). Areas of high crack density (yellows, greens) correlate to high μρ (<60 GPa g cm -3 ) but span a range of λρ (Fig. 13a); this range appearing to follow a trend similar to the boundary lines between sections of different brittleness characteristics (dashed grey lines). Similarly, they correlate with high AI (>12 km s -1 · g cm -3 , Fig. 13b) and Young's modulus (>50 GPa, Fig. 13d) but cover a range of V P to V S ratio (Fig. 13b) and Poisson's ratio (Fig. 13d). We previously described how high μρ, high Young's modulus and low AI correspond to both siliceous and calcareous mudstones, and that the clearest indicator of carbonate is high λρ. Taking this into account when interpreting Fig. 13, and in alignment with the observations of Fig. 12, it would appear that while natural fractures do show a preference for clay-poor zones, there is little clarity provided in these plots as to whether they are forming preferentially in siliceous or calcareous beds.

Reservoir quality evaluation
To evaluate the presence of high RQ intervals within the Bowland Shale section at PH-1, five parameters are used: clay concentration, porosity, permeability, TOC content and water saturation ( Table 4). The zones that pass all these cut-offs are shown in log view (Fig. 14) and rock physics space (Fig. 15). The rock physics model is shown in the same manner as Figs 10 and 13 but with data points coloured where a RQ interval exists. To better identify the portion of the plots where the density of RQ intervals is highest, a bivariate kernel density estimation is used, with shading strongest where density is greatest.
Clay concentration is a key factor in the evaluation of a shale gas reservoir. The best production within the Barnett Shale comes from zones with 27 per cent clay (Bowker 2003), and 50 per cent has been suggested as the maximum clay concentration for which a shale could be hydraulically fractured (Bowker 2007). Based on this, a 50 per cent maximum cut-off was taken for clay concentration. As the clay concentration in the majority of the Bowland Shale is beneath this cut-off (Fig. 14, track 3), this had minimal effect in reducing the prospective interval besides eliminating a clay-rich interval at the very top of the UBS and three clay-rich intervals in the LBS.
Shale porosity is an important parameter in assessing the storage potential of the formation. Jarvie (2012) suggests that porosity for an effective shale gas reservoir should be between 4 and 7 per cent. Taking 4 per cent as a minimum cut-off on its own has a marked effect on reducing the prospective interval as it is very close to the median value (3.8 per cent) thus eliminating near half of the shale. Notably, this eliminates the Pendleside Sandstone Member in the LBS [e.g. Fig. 14 (track 4) between 2500 and 2600 m].
TOC is a metric of the potential for a shale to generate hydrocarbons which is usually used in combination with Hydrogen Index, S2 (generation potential) and other geochemical parameters to assess source rock quality. It is purely a measure of the present-day organic richness of the formation which will have been reduced since deposition due to burial, maturation and expulsion of hydrocarbons. Jarvie (2012) suggests that 1 per cent present-day TOC is sufficient to be considered a shale gas reservoir. We take a slightly conservative cut-off similar to Andrews (2013) of 2 per cent which is very near the median calculated TOC for the Bowland Shale interval (1.9 per cent).
The remaining cut-offs were selected by taking the median value within the logged section. Minimum and maximum cut-offs were selected in this manner for permeability and water saturation, respectively. The final cut-offs are defined in Table 4 and Fig. 14. A total rock thickness of 125 m is found pass all five of these criteria (Fig. 14, track 9). We term this as RQ intervals, for future discussion. While RQ intervals are observed across the entire Bowland Shale interval, they are predominantly found (92 per cent of the total thickness) within the UBS.
The λρ versus μρ plot reveals a well-defined zone in which RQ intervals are clustered (Fig. 15a). The highest densities plot between 30 and 50 GPa · g cm -3 λρ and 50-60 GPa · g cm -3 μρ. In the remaining plots, the RQ interval is slightly less defined. On the AI versus V P to V S ratio plot (Fig. 15b), it manifests in a region ∼12 km s -1 · g cm -3 and with V P to V S ratios of 1.6 to 1.7. It is challenging to identify on the V P versus V S plot (Fig. 15c) alone due to the narrow bounds of the data but displays V P values between 4.5 and 5 km s -1 and V S values between 2.75 and 3 km s -1 . It forms a narrow range of E values (∼50 GPa) but a wide range of ν (0.15-0.3) (Fig. 15d).

Comparison with previous work
As mentioned earlier, previous studies have addressed the reservoir properties of the Bowland Shale, though focusing predominantly on the analysis of core samples. The high quartz concentrations and predominance of the siliceous mudstone lithology type reported herein is aligned with published ternary diagrams created following XRD analysis of core samples (Clarke et al. 2014a(Clarke et al. , 2018. The ternary diagram presented in Fig. 11 shows a similar character to those authors diagrams [e.g. 'Trend 2 in Clarke et al. (2018)'s Fig. 7(d)], but our results do include some argillaceous intervals not present in the other authors' analysis. While uncertainty regarding log response cannot be discarded, particularly in soft shales prone to washout, this may also reflect intervals where core recovery has either not been directed, or not possible due to the lithology type. Our mineralogical results are also aligned with Fauchille et al. (2017)'s analysis that the Bowland Shale typically exhibits quartz concentrations in excess of 50 per cent and carbonate concentrations around 20 per cent, though in contrast to our study, they found no samples with clay concentrations above 10 per cent.
Our calculated TOC content for the Bowland Shale is high and suggests good source rock potential. With an average value of 2.7 per cent (but up to >10 per cent), if this were analysed as a conventional source rock, the Bowland Shale would be considered a good-very good candidate according to the scheme of Peters & Cassa (1994). This range of TOC is consistent with Clarke et al. Andrews (2013) drew upon porosities typical for U.S. shale plays in their resource estimate for the Bowland Shale. Typical ranges include 1-5 per cent (Curtis 2002), 2.9-6 per cent (Jarvie 2012), and an average value of 3 per cent was chosen as applicable for the Bowland Shale (Andrews 2013). At PH-1, average total porosities are reported 3.70 and 3.31 per cent for UBS and LBS, respectively (Clarke et al. 2018). These values are comparable with those reported for the UBS and LBS (4.8 and 3.0 per cent, respectively) in the Gainsborough Trough (iGas Energy 2019), though the UBS demonstrates greater porosity in that case. Our calculated average total porosity is 4.2 per cent for the UBS and 3.8 per cent for the LBS which is greater than the values quoted by Clarke et al. (2018) but comparable on the whole.

Implications of natural fractures
The findings presented herein support the outcomes of previous studies concluding that the Bowland Shale contains natural fractures in the PH-1 well (Clarke et al. 2014a(Clarke et al. , 2018Fauchille et al. 2017). There are several zones in the logged section which have inverted crack density values around 0.06 (Figs 12 and 13). These magnitudes are lower than Guo et al. (2013)'s inverted crack density parameter for the Barnett Shale, however the Barnett Shale is known to be substantially fractured (Bowker 2003;Gale et al. 2007). The findings also suggest that the fractures develop in a wide range of lithology types, but in particular, those containing high amounts of either silica or carbonate. Direct comparison between microimagery data and crack density is challenging due to the differences in the vertical resolution of data (micro-imagery resolution is 0.5 cm whereas cross-dipole sonic resolution is ∼1.4 m). However, correlations can be observed (Fig. 16) whereby crack density increases in zones where natural fractures can be identified.
The impact of such macroscopic (i.e. can be identified visually using core samples and/or micro-imagery log data) natural fractures on hydraulic fracturing treatments and subsequent gas production from shales has been widely discussed in the literature (Gale et al. 2007Gale & Holder 2010;Gasparrini et al. 2014) but there remains dispute regarding whether they are a benefit or detriment to production. It has been suggested that sealed (carbonate-filled) fractures within the Barnett Shale do not contribute to gas production, and may in fact form barriers to fluid flow, impeding production (Bowker 2003;Montgomery et al. 2005), though Bowker (2003) does recognize that microfractures, likely created through the   (Vernik 1994;Jarvie et al. 2003) are likely an important factor on production. However, if such microfractures are assumed to be bedding-parallel, an additional fracture set of subvertical orientation would be required to connect the matrix to a horizontal production well to be drilled within a bed. Un-healed, open natural fractures, on the contrary, have been demonstrated to benefit production in other gas shales such as the Antrim (Curtis 2002) and Marcellus (Engelder et al. 2009) shales. However, laboratory and numerical studies support the concept that geological discontinuities in general, affect hydraulic fracture propagation (Warpinski & Teufel 1987;Zhang et al. 2007Zhang et al. , 2009). Reactivation of natural fracture networks within the Barnett Shale by hydraulic fractures has been demonstrated using microseismic data Downloaded from https://academic.oup.com/gji/article/228/1/39/6355446 by British Geological Survey Keyworth user on 27 September 2021 Figure 14. Log display illustrating the inputs to the RQ screening (tracks 3-7) with respective cut-offs annotated, the resulting RQ intervals (track 9) and the interpreted, broad zones where most RQ intervals are found (track 10). The inverted crack density is also shown for reference (track 8) (Contains OGA Well Log C Data accessed and published with permission of BGS). (Fisher et al. 2004;Warpinski et al. 2005), and it has been suggested this may ultimately benefit gas production as it would increase the number of hydraulic fracture 'strands', ultimately placing more host rock in contact with the hydraulic (Gale et al. 2007). However, as the energy imposed by fluid injection is dissipated over a larger network of fractures, their lengths would be expected to shorten and this led Gale & Holder (2010) to question if production would necessarily be enhanced. Nevertheless, numerical modelling has suggested that reactivated natural fractures can improve gas productivity by up to 25 per cent in the early stages (Ahmadi & Dahi Taleghani 2016).
Observations by other workers that natural fractures within the Bowland Shale are calcite-filled suggest there would not be a direct improvement on gas production as noted with open fractures of Antrim and Marcellus shales. The preferential failing of healed fractures may improve the development of a connected fracture network upon hydraulic fracturing, though thought needs to be given Figure 15. λρ versus μρ (a), AI versus V P to V S ratio (b), V P versus V S (c) and ν versus E (d) cross-plots for the Bowland Shale. The scattered data is well log elastic properties from PH-1, coloured where RQ interval is present. The triangular mesh represents the rock physics model iterating over concentrations of clay, quartz and carbonate with fixed kerogen concentration and water saturation. Porosity is fixed at 0 per cent (black lines with a colour fill) and 10 per cent (grey lines, no fill). A bivariate kernel density estimation is also shown to better illustrate the areas in which the RQ intervals are most densely situated. The as to the locations of naturally fractured intervals within the shale. The presence of naturally fractured thin carbonates interbedded with shale may form a detriment to gas delivery upon stimulation and production. A substantial question for the production potential of the Bowland Shale is that if a hydraulic fracturing treatment were targeted towards a certain shale interval, is there a risk of the hydraulic fractures reactivating naturally fractures within adjacent carbonate-rich beds rather than the desired shale? This phenomenon may lead to high initial well productivities as the fractured carbonates are produced from but followed by a sharp drop if the shale matrix itself is not sufficiently fractured.

Rock physics templates as a predictive tool
Seismic data has not been considered in this work; however, it is expected that as exploration within the Bowland Shale progresses, new 3-D data sets will be acquired and become available to researchers in due course. For example, Formby-16 and East Midlands-17, both acquired over Bowland Shale acreage are due for public release in 2021 and 2023, respectively.
Rock physics templates provided a helpful tool to geophysicists for the quantitative interpretation of seismic data. Templates for the Bowland Shale presented herein allow for lithology and porosity to be interpreted from elastic parameters but also provide a basis for Downloaded from https://academic.oup.com/gji/article/228/1/39/6355446 by British Geological Survey Keyworth user on 27 September 2021 Figure 16. Log display illustrating the relationship between dynamic Compact Micro-Imager log (CMI-Dynamic) (track 5) and inverted crack density (track 6). Fractures are interpreted as sinusoid shapes with amplitudes different from the bedding planes. The greater the amplitude, the greater the angle of the fracture relative to the well orientation. The ELAN mineralogy and lithotype logs are shown for reference (tracks 3 and 4). Note the discordance in vertical resolution between the traditional logging and the micro-imagery log data (Contains OGA Well Log C Data accessed and published with permission of BGS).
which RQ intervals in the shale may be defined from seismic data. Emphasis has been made on the use of λρ versus μρ plots as several works have presented that production from gas shales can correlate with specific regions of such plots, which we briefly discuss.
The concept of using λρ versus μρ cross-plots to interpret hydrocarbon zones and mineralogy was first proposed by Goodway et al. (1997) who identified gas zones within a conventional reservoir based on observing cross-over caused by a drop in λρ values below μρ. They also found that shaly gas sands plotted at λρ values of ∼30 GPa · g cm -3 and μρ values of ∼40 GPa · g cm -3 . We find our RQ intervals plot with similar λρ values to Goodway et al. (1997)'s shaly gas sand but with greater μρ values suggesting greater rigidity. Goodway et al. (2010) later extended this work to unconventional reservoirs; observing that brittle and ductile regions of the Barnett Shale could be separated reasonably well through the λρ versus μρ cross-plot and that lower λρ values could be correlated with increased Estimated Ultimate Recovery (EUR). They explained this may be caused by either the presence of gas-filled porosity, microfractures, overpressure or other geomechanical effects. Our proposed reservoir sections plot with greater λρ values than those considered in Goodway et al. (2010)'s analysis which may reflect a greater carbonate content within the Bowland Shale than the Barnett Shale.
This link between recovery and rock physics was further developed in a series of subsequent studies. Devonian shales of the Horn River Basin that hold low λρ values (20-25 GPa · g cm -3 ) have been demonstrated to deliver greater EUR per frack (Close et al. 2012). Furthermore, brittle zones on λρ versus μρ plots have been demonstrated to correlate well with locations of induced seismicity (Perez Altamar & Marfurt 2014, 2015. There does appear significant potential for such plots to assist with well planning in unconventional plays. However, currently, there is no production data available for the Bowland Shale, so it is impossible to analyse the link between geophysical parameters and recovery. The plots presented herein form a useful basis to begin defining the geophysical signature of the optimal reservoir properties of the Bowland Shale, but we acknowledge that with further drilling, production and therefore data, these will require re-examination.

Interpreting target zones
The areas in which most of the good RQ shale are found can be broadly divided into five zones [ Table 5 and Fig. 14 (track 10)]. These zones were identified through visually picking areas where the good RQ sections were clustered. Zones 1-3 are ∼80-100 m thick and are located throughout the UBS. Zones 1 and 2 are separated by a 100-m-thick interval whereas zones 2 and 3 are separated by a 45 m interval. Two thinner zones (4 and 5) are located at the base of the UBS and the LBS, respectively. These are thinner than the zones identified in the UBS. Note that zone 4 appears to show an unusual pattern of rock properties. It holds the highest clay fraction and porosity; however, water saturation is low. While the precise origin of this is uncertain, the zone exhibits elevated resistivity readings which would explain the low water saturation.
This analysis suggests that the UBS holds the better potential for shale gas reservoir potential than the LBS. Of the three distinct shales at the base of the PH-1 well and within the LBS, only one (zone 5) is interpreted as a prospective unit based on this classification. This zone is believed to be equivalent of the LB03 zone targeted by the PNR-1Z lateral well which was subject to the first multistage hydraulic fracturing operations within the Bowland Shale in 2018. High clay concentrations and variable TOC content limit the remaining sections of the LBS in being classified with the remaining zones of good RQ. As Fig. 14 demonstrates, a significant portion of these shales contain clay concentrations above the 50 per cent cut-off and TOC contents below the 2 per cent cut-off.

C O N C L U S I O N S
This paper presents a petrophysical examination using wireline logs at PH-1 to infer parameters describing the rock and reservoir properties of the Bowland Shale at PH-1. A rock physics model has been developed to study the effect of variations in some of these rock components (clay, quartz, carbonate and pores) and texture (the aspect ratio of each component); a necessary exercise due to the heterogeneity in those properties as observed in the shale. Seismic anisotropy was then accounted for in the model and used to invert a pseudo-fracture density parameter which can be related, in part, to natural fractures within the shale. From this work, the following conclusions can be drawn: (1) The Bowland Shale at PH-1 contains over 100 m of shale which is considered to have very good reservoir characteristics (such as high TOC, low clay and low water saturation), and porosities considered reasonable for tight mudstones. These intervals are located predominantly within the upper section of the shale.
(2) Crack density (considered a proxy for fracture density) mainly ranges from 0 to 0.04 but can exceed 0.06 in some areas. In some cases, crack density appears to correlate with natural, medium-high angle fractures from micro-imagery logs, but in other cases, it correlates to drilling-induced fractures. The fractures develop within clay-poor intervals, but our analysis fails to distinguish if they form preferentially within siliceous or calcareous lithology types.
(3) The variations in elastic properties observed from wireline logs can be broadly represented in a rock physics model by iterating over the clay, quartz and carbonate concentrations and porosity. Measures of rock strength (such as incompressibility, rigidity, Poisson's ratio, Young's modulus) show good correlations with lithology, and typical ranges for argillaceous, siliceous, and calcareous intervals of the shale are quoted herein.
(4) The intervals in which the best reservoir quality is found, plot in discrete zones of geophysical cross-plots suggesting the possibility of using such plots as a predictive tool in future seismic inversion studies. This may aid with well planning and ensuring the optimal target zones are drilled.

A C K N O W L E D G E M E N T S
David Griffiths is thanked for his input and discussion throughout the authors' PhD project. The British Geological Survey (BGS) are thanked for providing access to well data, and Schlumberger are thanked for the provision of Techlog software under academic license to Heriot-Watt University.

F U N D I N G
Iain Anderson acknowledges the support of a James Watt Scholarship from Heriot-Watt University (HWU) and the receipt of a British University Funding Initiative (BUFI) studentship award (grant number GA/16S/024) from the British Geological Survey (BGS), which provides the funding for the PhD project upon which this work is based. The PhD forms part of the Natural Environment Research Council (NERC) Centre for Doctoral Training (CDT) in Oil and Gas (grant number NE/M00578X/1). Jingsheng Ma acknowledges NERC grant number NE/R018022/1 for financial support.

DATA AVA I L A B I L I T Y
The data underlying this paper were provided by the British Geological Survey under license / by permission. Data may be shared on request to the corresponding author with permission of the British Geological Survey.
(a) (b) Figure A1. Theoretical plots of quartz/pore mixtures using Hashin-Strikman (HS) and Voigt-Reuss bounds. The Reuss and HS-lower bounds are equal, but the HS-upper forms a narrower bound than the Voigt bound. As water has no shear strength, the resulting shear modulus drops to zero upon addition of any pores (b).
(1966)'s tensor relating the far-field strain to the strain within the inclusion. Both approaches were presented by Berryman (1980). Fig. A2 illustrates the adding of fluid pores of different aspect ratios (ratio of short to long axis, α) into a quartz matrix. Note the elastic moduli fall beneath zero when the inclusion fraction is greater than the matrix fraction. A shortcoming of the KT approach is the assumption that the pore (or crack) concentration is dilute, that is the porosity must be substantially smaller than the aspect ratio of the pores (Xu & White 1995). This led Berryman (1980) to introduce the concept of self-consistency [self-consistent approximation-SCA, following the work of O'Connell & Budiansky (1974)] which removed the need to define a background material and allowed for the interactions between pores to be accounted for. Conceptually, the background host is instead replaced by an unknown effective medium and the solution developed iteratively (Fjaer et al. 2008): where K i and μ i are the bulk and shear moduli of the ith constituent of volume fraction χ i and K SC and μ SC are the bulk and shear moduli of the self-consistent approximate. Another approach to solving the porosity-aspect ratio problem is to incrementally add very small amounts of the inclusion into the host (differential effective medium-DEM), calculating the properties of the new effective medium at each step (Xu & White 1995. Both SCA and DEM allow for high concentrations of inclusions to be modelled without the resulting elastic moduli dropping beneath the elastic modulus of phase 2. This is exemplified in Fig. A3.

A3 Anisotropic models
Inclusion models such as SCA and DEM can be modified to consider the anisotropic effect of aligned inclusions. Hornby et al. (1994) extended the work of Berryman (1980) to consider the use of SCA and DEM to model anisotropic properties of aligned clay/fluid composites. They replaced the polarization factors, P and Q with the definition of a fourth-rank tensor G (after Eshelby (1957) and Mura (1987)), representing the strain felt by a single ellipsoidal inclusion in an unbounded isotropic medium. Jakobsen et al. (2000) and Bandyopadhay (2009) present the full equations for Hornby et al. (1994) vary the degree of alignment using an orientation distribution function (ODF) (Johansen et al. 2004;Sayers 1994). As previously mentioned, aligned fabrics observed in the Bowland Shale may cause anisotropy in elastic properties, however, there is little evidence available currently to support this and justify the need to incorporate anisotropy due to aligned clays, pores or kerogen. However, the anisotropic effect of aligned fractures is considered in the rock physics model. Crack models consider very thin ellipsoids added into a host matrix. The Hudson model (Hudson 1980(Hudson , 1981 for weak inclusions is based on a scattering theory analysis of the mean wavefield within an elastic solid with thin, aligned penny-shaped ellipsoidal cracks (Mavko et al. 2009). This modelling approach has been applied to investigate the anisotropic effect of bedding-parallel microcracks; known to develop in kerogen-rich, thermally mature shales (Vernik 1994) during the hydrocarbon generation (Vernik 1993;Vernik & Liu 1997). The Hudson model can consider bedding-parallel (VTI), bedding-normal (HTI) or multiple sets of cracks (orthorhombic) of different aspect ratios and crack densities so is, therefore, a very flexible model. In this work, it is used to add vertical, bedding-normal cracks to replicate the vertical or sub-vertical fractures discussed in later sections.
(a) (b) Figure A2. Theoretical plots of quartz/pore mixtures using KT theory. The different curves correspond to different aspect ratios (ratio of short to long axis, α) of pores. At high concentrations of low aspect ratio pores, the effective medium drops to less than zero strength.
(a) (b) Figure A3. Theoretical plots of quartz/pore mixtures using K-T, Berryman SCA and DEM theories. The aspect ratio of the pores is maintained at 0.1 in each case. The SCA and DEM techniques allow high concentrations of low aspect ratio pores to be modelled without the effective medium dropping to less than zero strength. Downloaded from https://academic.oup.com/gji/article/228/1/39/6355446 by British Geological Survey Keyworth user on 27 September 2021