Linking multi-scale 3D microstructure to potential enhanced natural gas recovery and subsurface CO 2 storage for Bowland shale, UK

Injection of CO 2 into shale reservoirs to enhance gas recovery and simultaneously sequester greenhouse gases is a potential contributor towards the carbon-neutral target. It oﬀers a low-carbon, low-cost, low-waste and large-scale solution during the energy transition period. A precondition to eﬃcient gas storage and flow is a sound understanding of how the shale’s micro-scale impacts on these phenomena. However, the heterogeneous and complex nature of shales limits the understanding of microstructure and pore systems, making feasibility analysis challenging. This study qualitatively and quantitatively investigates the Bowland shale microstructure in 3D at five length scales: artificial fractures at 10–100 m m scale, matrix fabric at 1–10 m m-scale, individual mineral grains and organic matter particles at 100 nm–1 m m scale, macropores and micro-cracks at 10–100 nm scale and organic matter and mineral pores at 1–10 nm-scale. For each feature, the volume fraction variations along the bedding normal orientation, the fractal dimensions and the degrees of anisotropy were analysed at all corresponding scales for a multi-scale heterogeneity analysis. The results are combined with other bulk laboratory measurements, including supercritical CO 2 and CH 4 adsorption at reservoir conditions, pressure-dependent permeability and nitrogen adsorption pore size distribution, to perform a comprehensive analysis on the storage space and flow pathways. A cross-scale pore size distribution, ranging from 2 nm to 3 m m, was calculated with quantified microstructure. The cumulative porosity is calculated to be 8%. The cumulative surface area is 17.6 m 2 g (cid:2) 1 . A model of CH 4 and CO 2 flow pathways


Introduction
Climate change has resulted in many countries aiming to reach net-zero greenhouse gas emissions by 2050. In the energy transition from fossil-based fuel to zero-carbon, several strategies need to be considered. Natural gas produced from shales is considered to be a relatively 'clean' energy source, which could increase the energy supply and reduce CO 2 emissions significantly, when compared with oil. 1 Additionally, it is generally accepted that to meet net-zero, carbon capture, utilisation and storage is needed.
The concept of CO 2 injection into shale reservoirs has the potential to achieve these two goals of enhanced natural gas extraction and carbon sequestration simultaneously e.g. ref. 2 and 3. It could potentially protect fresh water, enhance gas recovery, and sequester CO 2 into the subsurface. 3 The process of preferentially adsorbing CO 2 while desorbing CH 4 has been suggested in shale reservoirs previously. 4,5 Globally there has been interest from government, industry and the academic community, but high costs and outstanding research questions restrict wide-scale implementation of this technology e.g. 6. A core issue for commercialisation is the limited understanding of pore networks, gas transport pathways and gas behaviours in shales. 7 The heterogeneous nature, fine-grain size and small pores within shales makes characterization of microstructure, pore networks and gas transport challenging. ref. 8 Shale microstructure classically consists of pores and fractures distributed within a matrix in which non porous and porous grains are embedded at different scales. All can contribute differently to transport properties and storage properties of CO 2 and CH 4 . Furthermore, recent experimental results by  suggest that using CO 2 as a fracturing fluid rather than water may have a significant impact driven by a lower fracturing pressure, while leading to a denser network of largely shearmode fractures.
A range of techniques have been applied to quantify microstructure and pore networks in shales, for example, gas adsorption, mercury intrusion, small angle neutron or X-ray scattering e.g. ref. [9][10][11]. They can provide data on pore size, porosity, adsorption, and permeability, but the pore structure, spatial distributions and connectivity cannot be observed or quantified directly. The complexity of shale makes it necessary to additionally quantify the microstructure and their pore system in 3D by means of advanced imaging techniques. High-resolution SEM imaging techniques, such as broad ion beam-scanning electron microscopy (BIB-SEM), can characterise the microstructure over a range of scales in 2D. 12 X-Ray tomography has been applied to fractures, 13 organic matter, 14 large mineral grains 8 and clay minerals 15 at micro-scale in 3D. Ga + focused ion beam-scanning electron microscopy (FIB-SEM) 16 and Xe + plasma FIB (PFIB, or P-FIB) 17 can be applied to resolve macro-and meso-pores at nm-scale. 16 Nanopores within organic matter or between clay minerals can be imaged by Transmission Electronic Microscopy (TEM) tomography 17 and Helium Ion Microscopy (HIM). 18 Owing to the complex microstructure and high heterogeneous nature of shales, the assessment of potential of CO 2 storage and enhanced gas recovery in shales requires a multi-scale, multi-approach investigation of the open pore spaces for free molecules, the pore surfaces for adsorbed molecules and the flow pathways between them. Particularly, a detailed quantification is needed of the nm-scale pores and the mm-scale phases associated with these pores (e.g. minerals and organic matter) alongside with the mm-scale sample fabric and stimulated fractures. To this end, in this study we combine imaging across multiple scales with other laboratory measurements, for example, porosity, permeability and adsorption, to characterize shale microstructure and pore system in 3D and to link these with the potential for enhanced gas recovery and subsurface CO 2 storage.
As the largest shale gas resource in the UK, the Bowland shale is estimated to contain a gas resource of 1300 trillion cubic feet (at surface temperature and pressure). 19 The gas resource is comparable to those of major producing shales in North America, such as Barnett, Woodford and Fayetteville Shales. 20 However, previous studies mostly focused on the shale gas potential, while the carbon sequestration potential has not been considered. The aims of the research presented here are: (1) to characterise the shale microstructure across a range of scales (nm-to mm-) and to assess its heterogeneity; (2) to quantify the pore structure and network to investigate the possible gas storage and pathways for methane and CO 2 flow; and (3) to link the microstructure and reactive transport to further assess the potential for simultaneous enhanced gas recovery and subsurface CO 2 storage.

Sample selection
Recent studies of the Carboniferous Bowland Shale Formation (herein referred to as the Bowland Shale) e.g. ref. 20 and 21, have demonstrated that the formation contains organic-rich horizons. Samples extracted from these horizons contain greater than 1% TOC, are quartz-rich (450 wt%) and clay-poor (o20 wt%). Based on detailed core observation and microfacies analysis, three typical samples were selected in this study. These consist of (1) one typical organic-rich Bowland Shale (BOR; TOC 6.1 wt%) sample, refer to B6 in ref. 22 from the upper Bowland Shale Formation, which was selected for cross-scale 3D microstructure characterisation and adsorption experiments (Table 1). This sample was taken from a depth of 2344.7 m in the Preese Hall-1 borehole which was the UK's first shale gas exploration well (Clarke et al., 2018). (2) One typical organic-lean (BOL; TOC 1.0 wt%) sample refer to B8 in ref. 22 at a depth of 2495.3 in Preese Hall-1, which was used in adsorption experiments for comparison. (3) An adjacent sample (BF) to BOR with similar compositions in the same core (2354.5 m) was selected for fracturing characterisation owing to the limited volume of available BOR material.

Imaging techniques
The selected typical organic-rich sample was firstly examined in 2D for microstructure and pore network characterization and then in 3D for cross-scale quantification. (2) Synchrotron source micro-XCT. BOR sample was cut into a 2 Â 2 Â 10 mm 3 stick and scanned in I13 Beamline (MT15506-1) at DLS using PCO Edge 4Â camera and Pink beam. The energy was 20 keV and exposure time is 0.04 s. 4000 projections were taken for each scan with voxel size 1.6 mm. Four tomograms were taken at different vertical locations with 10% overlap to access the larger range of microstructure variation.
(3) Lab source nano-XCT. BOR sample was cut into a 60 mm diameter cylinder using a laser cutter and scanned in a Zeiss Ultra 810 nano-CT scanner in Manchester X-ray Imaging Facility (MXIF). Two scans using phase contrast mode and adsorption mode were performed for the same volume. These two images were registered together for segmentation. 8 kV monochromatic beam was used for the scan with 160 s exposition time. The voxel sizes of both images were 64 nm.

FIB-SEM.
To characterize the nanoscale features in shales, advanced Focused Ion Beam Scanning Electron Microscopy (FIB-SEM) is used to perform high-magnification 3D imaging. Specifically, Xe + plasma FIB (PFIB) and Ga + FIB are both used in this study at diffferent scales. 17 (1) Plasma FIB-SEM. The BOR sample scanned in nano-CT was used for plasma FIB. Pore and organic matter cannot be distinguished in the nano-CT images, but are identifiable in the PFIB data. Secondary electron images were taken at 1 kV and the organic matter and pores were quantified after image processing. The SEM images were taken at 3449 Â 3735 voxels with 33 Â 42 nm 2 resolution. 815 slices were collected during rock milling with a milling distance 50 nm.
(2) Ga + FIB-SEM. Two sites were selected to be imaged with Ga + FIB-SEM, including one for mineral matrix and one for a single organic matter particle. The organic matter particle site was slightly smaller than the mineral matrix site (Table 2), and both of the data set were collected at 2 kV.
2.2.4 Image processing. X-Ray tomography images required specific image processing to align the images into the same locations. Four micro-CT tomograms were registered into a whole image using the overlap volumes of the top and bottom of adjacent tomograms in Avizo (Thermo Fisher Scientific, United States). Two nano-CT tomograms of adsorption and phasecontrast scans were registered into the same location for an enhanced quantification of minerals and organic matter. Specifically, granular minerals and clay minerals were segmented in the phase-contrast image while the heavy minerals and organic matter were segmented at the adsorption image. PFIB and both Ga + images were aligned (rotation) and sheared to correct the drift during the imaging.
Subsequently, all images were filtered and segmented using varied filtering (e.g. non-local means, FFT, median filters) and segmentation methods (e.g. thresholding, tophat, watershed) in different series of images to present the identified features at  BOR  6  67  6  2  1  9  7  2  BOL  1  51  3  11  2  5  18  9  BF  1  69  7  5  1  6  8 3 the corresponding scale, as a standard workflow. 24 Features below 27 voxels (3 Â 3 Â 3 for XYZ) for XCT data and 9 voxels (3 Â 3 for XY plane) for FIB data were not considered. The size, volume, surface area, volume fractures of each feature were measured in voxels using Avizo (Standard and Fire versions, FEI, Hillsboro, United States). The sizes of pores or particles are measured in sphere equivalent diameter in this paper. This is defined as the diameter of a sphere with the equivalent volume. 2.2.5 Heterogeneity analysis. The heterogeneity of the microstructure and pore system in the 3D images was quantified with respect to vertical variation, fractal dimensions and degree of anisotropy.
(1) Vertical variation. Horizontal bedding was recognized in the Bowland shale samples (e.g. 1-2 mm), and therefore vertical variations of different phases were measured for heterogeneity analysis of microstructure. In each image dataset, the volume fractions of identified phases at this scale were quantified in all horizontal slices and were then plotted with the corresponding vertical positions. The maximum, minimum, median values and standard deviations for each phase at each scale were also identified. The variations of phase distribution at different scales were then quantified.
(2) Fractal dimensions. The fractal dimension was determined by a box-counting algorithm to analyse the spatial arrangement of porous phases or pores within the 3D binary images using 'BoneJ' plug-in in ImageJ. 25 Certain sizes of boxes with the identified feature scan the whole image and the number of the boxes were counted. The box counting starts from the size equalling to 80% of the side length of the original image and ends at the 27 voxels (3 Â 3 Â 3). Features smaller than 27 voxels are not considered owing to the small sizes. The scaling factor is set to be 1.2 for each of the image to get enough steps for the analysis, which means the box size was divided by 1.2 after each iteration.
(3) Degree of anisotropy. The degree of anisotropy is calculated to quantify the directionality of the shale components and pores using the 'BoneJ' plugin in ImageJ. 26 It consists of four steps to determine the degree of anisotropy including (1) find mean intercept length (MIL) of vectors in random directions; (2) plot MIL vectors from background to foreground into a point cloud; (3) solve the equation of an ellipsoid that best fits the point cloud and (4) calculate the degree of anisotropy from the radii of the ellipsoid (longest/shortest). The degree of anisotropy ranges from 0-1 and the closer to 1 means higher anisotropy of the identified feature.

Laboratory bulk property measurement
Bulk laboratory measurements were conducted to provide the essential information for the assessment of gas recovery and CO 2 sequestration along with the image results. These include (i) high-pressure adsorption measurements (CO 2 and CH 4 ); (ii) low-pressure N 2 and CO 2 adsorption measurements for pore size distribution analysis; and (iii) porosity and permeability.
(1) High-pressure adsorption. CO 2 and CH 4 adsorption measurements were made BOR and BOL at 7 MPa and 80 1C (reservoir conditions). The measurements have been carried out on powdered samples (2 g) in a Rubotherm Magnetic Suspension Balance following the procedure described in ref. 27. Samples were dried at 120 1C under vacuum for at least 12 hours in situ prior to the measurements. The obtained excess amount adsorbed, n ex , was converted to absolute amount adsorbed, n a , using the following equation: n a = n ex /(1 À r/r a ), where r is the bulk density of the given gas at the selected pressure and temperature, while r a is the adsorbed-phase density, presumed to be constant at 21.1 mol L À1 for CO 2 and 26.3 mol L À1 for CH 4 .
(2) Pore-size distribution. Cryogenic adsorption of N 2 at 77 K and CO 2 at 273 K was conducted on about 1 g of both samples using a Micromeritics 3Flex Surface Characterisation Analyzer. Samples were first degassed at 120 1C under vacuum ex situ for at least 12 hours and for a further 3 hours in situ at the same temperature. Gases were probed in the pressure range 2 Â 10 À6 -1 Â 10 À1 MPa. The distribution of pore sizes above 2 nm (mesopores and macropores) was obtained by the application of the BJH theory to the adsorption branch of the N 2 isotherms 28 while the micropore (pores less than 2 nm in size) volume was obtained by the application of the linearised Dubinin-Radushkevich (DR) equation to the CO 2 adsorption isotherm. 29 The Harkins-Jura equation was used to describe the thickness of the adsorbed layer. 30 A liquid density of 21.1 mol L À1 was assumed for CO 2 .
(3) Porosity and permeability measurement. Porosity was measured by helium porosimeter and gas permeability was measured by pore pressure oscillation 31 using Argon gas, and in two orientations, parallel and normal to bedding. Six permeability measurements were recorded in each orientation to determine a mean value and its uncertainty with at pore pressure 5 MPa and confining pressure 10.5 MPa.

Porosity and surface area upscaling from 1-100 nm to 10-100 lm scale
The porosity, pore size distribution and surface area measured at 1-100 nm from images and lab measurements were upscaled to pore associated phases (organic matter and minerals) at mm-scale; and, then to the sample fabric at 10-100 mm scale using an image-based upscaling method that has been shown to be applicable to shales. 24 The correlation coefficients (i.e. volume and surface area ratios) of pore and pore associated phases were quantified at 1-100 nm. Those were then assigned into the correlation coefficients of pore associated phases and sample matrix at mm-scale, and finally to the whole sample at 10-100 mm scale. Subsequently, pore size distribution was normalised as number of pores per mm 3 , and as surface areas per mm 3 for the comparisons between scales.  (3) and (4), which are modified based on the US-DOE-NETL method [eqn (1) and (2)]; 32 where G CO 2 is the mass of CO 2 that is stored in the available volume; A t and H g are the area and thickness of gas shale; r CO 2 is the density of liquid or supercritical CO 2 at downhole conditions; F is the porosity for free gas; 1 À F is the solid fraction; and r sCO 2 is the mass of CO 2 adsorbed per unit volume of solid rock (corresponding to 1 À F). E A and E h are the fieldscale efficiency factors (ratios of efficient and total values) for the area and thickness. E F and E s are pore-scale efficiency factors for pore volume and sorbed volume, which are used to correct values measured in the laboratory, including effects arising from inaccessible porosity, variable mineralogy and the presence of moisture. Furthermore, E s can be calculated by the equation below: where E s(d) refers to the adsorbed CO 2 at dry conditions and accounts for the variability in sample mineralogy. E s(m) is the moisture equilibrated adsorption efficiency, and refers to the ratio of the adsorbed CO 2 at moisture-equilibrated and dry conditions. It additionally accounts for a reduction in gas uptake due to the presence of moisture. Therefore eqn (1) and (2) can be further re-written as where V is the efficient volume for shale gas, E v is the fraction of shale volume used for CO 2 storage, which also refers to the practical efficiency compared with the theoretical maximum value. 2.5.2 Parameters used in the equations. Chen et al. 2019 33 report that pores below 2 nm aperture are predominantly filled by adsorbed gas and the majority of free gas resides in pores 45 nm aperture (which make up 77.7% of total pore volume). In this study, F is therefore computed as the cumulative pore space found in pores larger than 2 nm and 5 nm and these two thresholds are used to estimate to find the range of the stored gas volume. V, r CO2 , and E v were determined as the average values reported in the previous studies of Bowland shale. 19,20 E v , the practical efficiency, is a reservoir scale efficiency factor not measured in this study and therefore set to be 1 for the theoretical maximum value estimation. Different values of E v were used for practical estimation in the next step. r sCO 2 was estimated from the adsorption experiments by the gas uptake measured on a dry shale. E s(d) was estimated from measurements carried out on BOR and BOL samples. It is expected that the primary formation water reduces the adsorption capacity, and this effect is included in E s(m) . The efficiency factor E s was considered at moisture equilibrated status in the subsurface. The adsorption was measured on dry samples in the laboratory for E s(d) . In reality, the exiting of primary formation water in the surface reduces the adsorption ability, which affects E s(m) .
The water saturation in Bowland shale is relatively low with an average value 32.4% across all formations. 20 Therefore, the final adsorbed values of Bowland shale in the subsurface are assumed based on moisture equilibrated status according to other studies with a similar water content, for example, 22% less CO 2 adsorbed in a 32.10% moisture sample than the dried samples e.g. ref. 34. Therefore E s(m) in eqn (3) is set to be 0.78 (= 1-0.22) in this study.
Then, two simple calculations were performed for practical volume estimations. The first calculation was based on an estimation that technically accessible CO 2 storage was around 32% in some commercial shale plays 35 and therefore E v was set to be 0.32. The second estimation was performed using P10 (10 percent probability that an estimated value is less than the P10) and P90 (90 percent probability that the estimated value is less than the P90 value) of E v for free gas (0.15; 0.36) and adsorbed gas (0.11; 0.24), which were extracted from the modelling results of different shale reservoirs 36 for practical storage potentials estimation of stored CO 2 volume.

Results
Artificial fractures, sample fabric, mineral matrix and organic matter, macro-pores, and nano-pores were imaged and quantified at different scales. Specifically, artificial fractures were quantified at 10-100 mm scale; sample fabric, including the distribution of granular minerals, clay minerals, organic matter particles and heavy minerals, were characterised at low-resolution micro-scale (1-10 mm scale); the grain sizes and higher-resolution distribution of these porous phases were further characterised and quantified at high-resolution micron-scale (100 nm-1 mm scale); macropores (450 nm) within these phases were quantified at low-resolution nano-scale (10-100 nm scale) and mesopores (2-50 nm) were characterised at high-resolution nano-scale (1-10 nm scale). Additionally, micropores (o2 nm) were measured by sorption experiments.

Identification of microstructures
The major compositions are identified as granular minerals (GM), clay minerals (CM), organic matter (OM) and heavy minerals (HM). The distribution of these compositions and pores were recognized as sample microstructure. The overall microstructure of the sample is relatively similar horizontally with clusters of large grains and matrix (Fig. 1a). Granular minerals and heavy minerals scatter in the sample, while clay minerals and organic matter distribute partially along the beddings in the matrix (Fig. 1b).
3.1.1 Bowland shale components: (1) Granular minerals. Quartz, calcite, ankerite and albite are observed as granular minerals in these Bowland Shale samples (Fig. 1a-d). Small pores exist in the crystal gaps, with spherical or irregular polygon shapes, which are unlikely to have an interconnected pore network (Fig. 1e).
(2) Clay minerals. Kaolinite and muscovite dominate in clay minerals in these samples, partly orientated to the bedding  Fig. 1a-d). Inter-mineral pores are present between clay mineral grains or between clay mineral and other mineral grains with elongated wedge shapes (Fig. 1e). The latter of the two are generally slightly larger, especially when fine-grain granular minerals inlay the clay mineral matrix.
(3) Organic matter. Organic matter particles are present scattered in the matrix (Fig. 1a-d). Some are porous while others are non-porous. Many of the porous organic matter particles are large and spherical in shape. Some are lamellar, interbedded with chlorite or other clay minerals. Non-porous organic matter is dominated by elongate shape or irregular polygons. Some inter particle pores are present in the interface of non-porous organic matter and adjacent minerals. They are commonly present as curved flakes or possessing a crack-like geometry, aligned along the edge of non-porous organic matter particles (Fig. 1f).
(4) Heavy minerals. The major heavy mineral detected in the Bowland Shale is pyrite, present in the forms of framboids and single crystals (Fig. 1a-d). Framboids host some polygon shaped pores while single crystals do not host pores.
3.1.2 Microcracks and pores. Microcracks and pores were observed in and between minerals and organic matter ( Fig. 1e and f). Microcracks in the matrix follows the mineral boundaries primarily, while those in organic matter have different distributions. Some of cracks present at the interface of organic matter and clay minerals, and some cracks open in the middle or organic matter particles. Pores are observed as mineral pores and organic matter pores. The minerals pores are generally irregular in shape with wedge-shape corners while the organic matter pores are predominately spherical. Pores and microcracks are not distinguished in this study owing to the ambiguous definitions and comparable roles they are playing for carbon storage and sequestration.

Multi-scale 3D quantifications
The spatial distribution of the four major components was quantified at 1-10 mm scale (1.6 mm voxel size) and the individual grains of these components were characterized at 100 nm-1 mm scale (0.1 mm voxel size). Macropores were imaged at 10-100 nm scale, and then mesopores in minerals and organic matter then quantified respectively at 1-10 nm scale. To assess the heterogeneity, the volume fractions of each phase were plotted with the vertical distance from the top of the images at LR-mm, HR-mm and 10-100 nm scales. This was not present either for the 10-100 mm scale as the matrix is not able to be identified or for 1-10 nm scale owing to the small sample volume below the representativity elementary volume (REV).
3.2.1 Fractures at 10-100 lm scale. Stimulated fractures opened at 53 MPa injection pressure from a borehole oriented normal to bedding, and a confining pressure of 51 MPa. Two vertical (cross-bedding; along the injection direction) fractures and one horizontal (along bedding; vertical to the injection direction) were identified in the image (Fig. 2). The fractures are relatively planar across the sample section. This image was captured immediately after the fracture opening and the width and volume of the fractures may vary over time, therefore the quantification of fracture apertures and volumes is not presented here.
3.2.2 Sample fabric at 1-10 lm scale. The four phases in the selected volume of the sample, (organic matter, granular minerals, clay minerals, and heavy minerals), were imaged using synchrotron-based X-ray tomography at 1-10 mm scale (Fig. 3). Micro-cracks or pores are not distinguishable from the above four phases at this scale. The quantification uncertainty of volume percentages is AE4 vol% at this scale granular minerals occupy 71 vol% while clay minerals account for 26 vol%. Organic matter and heavy minerals only contribute 2 and 2 vol% in volume fractions respectively. Granular minerals show a decreasing trend upwards which ranges from 85% to 48% with a median value 72%, while the clay mineral increase from 12% to 45% with a median value 25%. Organic matter slightly increases from 1% to 3%, correspondingly with clay minerals, with a median value 14%. Heavy minerals distribute relatively evenly, and the volume fraction ranges from 1% to 3%, with a median value 2%. In this study, note here that the organic matter content is able to be higher than 1% in local intervals of 500 mm to 1.5 mm length.

3.2.3
Porous phases at 100 nm-1 lm scale. The four major components, organic matter, granular minerals, clay minerals, and heavy minerals, were imaged at high-resolution micro-scale using nano-CT (Fig. 4) for higher resolution quantification at 100 nm-1 mm scale. The quantification uncertainty of volume percentages is AE2 vol% at this scale. Granular minerals, clay minerals, organic matter and heavy minerals occupy 71, 17, 10 and 2 vol% respectively. Pores and micro-cracks are not distinguished at this scale, but are separated in the same sample at the next scale.
The vertical variation for granular minerals ranges from 46% to 52% with a median value 46% and for clay minerals ranges from 1% to 2% with a median value 1%. Organic matter varies from 43% to 8% with a median value 7%. Heavy minerals max at 15% and min at 0% with a median 4%.
3.2.4 Pore size, structure, and networks at nm-scale 3.2.4.1 Macro-pores at 10-100 nm scale. The sample scanned in PFIB was part of the one scanned in nano-CT, and macropores (4100 nm) are distinguishable from the matrix (Fig. 6a), which occupies 0.64% (Fig. 5b). The quantification uncertainty of volume percentages is AE0.5 vol% at this scale. 16 056 matter particles and 5832 pores were identified in the PFIB images at 10-100 nm scale. The vertical variation of the porosity ranges from 0.0 to 3.3% at this scale with the step of Z axis 64 nm (Fig. 5c). A good correlation between organic matter and pore volume fraction was found with a correlation coefficient 37 value 0.84. The largest pore sizes were quantified to be 2.8 mm with bimodal distribution peaking at 800 nm and 2 mm (Fig. 8b).
3.2.4.2 Meso-pores at 1-10 nm scale. Based on the observations under high-resolution SEM, organic matter and minerals are two typical phases hosting pores (Fig. 1). Imaging and quantification of the pore system was performed in an organic matter dominated site and a mineral dominated site ( Fig. 6a and d). Due to the image resolution limitation (equivalent pixel size B12 nm), only pores above 30 nm aperture are quantified and described in this section. 5574 mineral pores and 3353 organic matter pores were identified in the two FIB images respectively at HR-nanoscale. The quantification uncertainty of volume percentages is AE0.3 vol% at this scale. The normalized cumulative porosity and surface area of the whole range of pores is shown in Fig. 8.
(1) Mineral pores. Pores inside minerals (intra-mineral pores) or between minerals (inter-mineral pores) are both observed in these images (Fig. 6a). Intra-mineral pores are commonly found in granular minerals such as quartz, calcite, ankerite and albite, with wedge shapes. Small pores exist in the crystal gaps, with polygon shapes. They only occupy 8.5% of the total mineral pores. Inter-mineral pores are observed between many mineral types, and account for 91.5% of total mineral pores. Pores between granular minerals display irregular shapes, sometimes featuring sharp corners, while pores between clay mineral grains or between clay minerals and other mineral grains with wedge shapes (Fig. 6b).
The pore network shows a low connectivity at the resolution of the study and a relatively homogenous texture (Fig. 6c). The mineral pores occupy 1.4% of the volume within the mineral phase, while accounting for 65.8% of total pores.
(2) Organic matter pores. Organic matter pores are mainly spherical, with variable sizes (Fig. 6d and e). The distribution of OM pores is more heterogeneous compared with the mineral pores, and the size ranges from 30 nm to 3 mm in this image, with a major peak at 40 nm and 200 nm. They form a highdensity map of pores although not a globally connected network was observed in the FIB image (Fig. 6f). Organic matter pores occupy 9.0% within the organic matter phase.
Although there is a higher density and larger total volume of organic matter pores than mineral pores in the single image dataset, they tend to be lower in total number of pores per mm 3 when the whole sample is taken into account (Fig. 6g) owing to the lower volume fraction of organic matter overall. Number of pores show a bimodal distribution, with peaks at 30 nm and 200 nm (Fig. 6g). The incremental porosity peaks at 2000 nm, which is greater in sizes compared with the mineral pores (Fig. 6h). Cumulative organic matter pores in FIB image is 1.2% (Fig. 6i) and the cumulative surface area is 21 mm À1 (E54.6 Â 10 3 mm 2 g À1 ) (Fig. 6j).

Heterogeneity across scales.
The fractal dimensions and the degree of anisotropy of each feature were obtained at multiple scales for heterogeneity studies. Generally, fractal dimensions decrease with higher resolution and the degree of anisotropy increases for each phase due to the sample's limited size (Table 3) Organic matter is the only feature which can be characterized crossing three scales owing to the relatively large particles sizes, the fractal dimensions of which slightly decrease from 2.546 at 1-10 mm scale to 2.510 at 100 nm-1 mm scale and then to 2.458 at 10-100 nm scale. The degree of anisotropy of organic matter range between 0.435 to 0.342 and the standard deviations are 0.008 and 0.035. Fractal dimension of total pores is 2.709 and it is slightly higher for OM pores (2.315) than mineral pores (2.195). The degree of anisotropy increases from 0.429 at 10-100 nm scale to 0.736 for OM pores and 0.504 for mineral pores. As a consequence, particle and pore anisotropy (a) significantly and quantitatively changes across the scales from the mm to nm, and (b) it is heterogeneous for the different phases at the same scale.

Laboratory bulk property measurements
3.3.1 CO 2 and CH 4 adsorption. CO 2 and CH 4 adsorption of BOR and BOL was measured at reservoir conditions (80 1C and 7 MPa). The CO 2 is present as supercritical phase ( sc CO 2 ) under this condition. The amount of sc CO 2 adsorbed in BOR sample was 4.5 kg per t while the CH 4 is only 0.6 kg per t. In contrast,

View Article Online
BOL was found to only adsorbed 2.9 kg per t sc CO 2 and 0.3 kg per t CH 4 (Fig. 7A).

3.3.2
Permeability. The permeability result shows two orders magnitude heterogeneity between flow orientations. The mean value is 3.9 Â 10 À21 m 2 in the direction of normal to beddings and 3.3 Â 10 À19 m 2 in the direction of parallel to beddings (Fig. 7B).
3.3.3 Pore size distribution (PSD) extracted from adsorption data. The PSD measured by adsorption shows that the major pores distribute between 2-5 nm with a small hump around 25 nm, which corresponds to the peaks shown in the image results (Fig. 8a). When this is converted into incremental pore volume, the peaks are found at 2-3 nm and 10-50 nm (Fig. 8b).

Upscaled pore size, volume and surface area
The porosity and surface area measured in the images at 1-100 nm were upscaled to 10-100 mm scale using an image based upscaling method. 24 The number of pores shows an increasing trend when pore sizes decrease, and the peak appears at 2 nm which is the limit of the techniques used in this study (Fig. 8a). The incremental porosity has the major

View Article Online
peak at 200-500 nm and the second peak at 20-40 nm (Fig. 8b). The incremental surface area shows a continuously increasing trend with decreased pore size (Fig. 8c). The cumulative porosity is calculated to be B7.8% and the cumulative surface area is 1.79 m 2 g À1 (Fig. 8d). Among them, pore above 100 nm contribute 47.5% in total porosity (3.7% of sample volume) and less than 0.1% for cumulative surface area (1.94 Â 10 À6 m 2 g À1 ); pores between 10-100 nm contribute 35.9% in total porosity (2.8% of the sample volume) and less than 0.1% for cumulative surface area for cumulative surface area (7.67 Â 10 À4 m 2 g À1 ); pore below 10 nm contribute 16.6% in total porosity (1.3% of the sample volume) and over 99.9% for cumulative surface area (1.79 m 2 g À1 ) (Fig. 8d).

CO 2 storage calculation
The potential CO 2 storage volumes in shales are calculated through eqn (4), first for theoretical maximum values, and then the practical values using P10 and P90 values for free gas and adsorbed gas simulated in other studies. 36 3.5.1 Theoretical maximum values for CO 2 storage. V, the effective CH 4 shale volume, for the upper Bowland unit is 9.31 Â 10 11 m 3 and for the lower Bowland unit is 3.45 Â 10 12 m 3 . 19 E v , the practical efficiency is set to be 1 for theoretical maximum value. This assumption implies that the potential storage volumes estimated here should be thought of as upper bounds. The mean porosity, F, is 4.0%. 20 r CO 2 is the density of the liquid or sc CO 2 which is 134.1 kg m À3 at reservoir conditions (80 1C and 7 MPa). The r sCO 2 values measured in the laboratory on BOR (TOC 6%; 4.5 kg per t) and BOL (TOC 1%, 2.9 kt per t) at 7 MPa and 80 1C (Fig. 8a) were combined to estimate r sCO 2 E s . Specifically, the TOC content of the Bowland Shale varies mainly between 2 and 2.5 wt%. We assume here a shale composition of 75% BOL (TOC 1%) and 25% BOR (TOC 6%), yielding a mean TOC 2.25%, as observed in the field. The mean density of Bowland shale is 2.6 t per m 3 . 19 E s(m) , the moisture  91.5% of total mineral pores and all organic matter pores are CO 2 and CH 4 accessible, so E F , the pore connectivity efficiency, are calculated to be 0.95 (= 1 À 0.082 (non-assessable pores in mineral pores) Â 0.658 (mineral porosity in total porosity)).
Total porosity of BOR sample is calculated to be 7.8% of total sample volume (Fig. 8d). Specifically, pores above 2 nm and 5 nm occupied 94% and 85% of the total porosity respectively, corresponding to the min and max values of the free gas porosity (Section 2.5). When converted to mean porosity of Bowland shale (4%), F (above 2 nm) is 3.8% (= 4% Â 94%) and F (above 5 nm) is 3.4% (= 4% Â 85%). These values were used to define the boundaries to calculate the total  View Article Online amount of stored CO 2 using the eqn (5) and (6) which were derived from eqn (4).
In each 1 m 3 shale with 32.4% water saturation, 4.3-4.8 kg CO 2 can be stored as free molecules and B6.4 kg CO 2 can be adsorbed, assuming full efficiency of CO 2 shale and sorbed volume. This can be converted into 2.0-2.2 kg per t and 3.1 kg per t respectively as the average density for shale is 2.6 kg m À3 . Therefore, 19.0-21.2 Gt CO 2 is estimated to be stored as free gas and 28.3-28.5 Gt CO 2 is estimated to be adsorbed on pore surface, consisting of 47.5-49.5 Gt carbon storage as a theoretical maximum in Bowland shale. It is noted that these values were calculated based on the assumption of all gas shales in Bowland shale are used for CO 2 (i.e. E v = 1) owing to the lack of reservoir-scale measurements. 3.5.2 A simple practical estimation for CO 2 storage. The practical storage potential will be smaller than the theoretical maximum values shown above. Here we provide two simple estimations for practical stored volumes. Estimate (1) was calculated on basis of E v = 0.32. 35 Estimate (2) was calculated using the simulated P10 (0.15 for E F and 0.11 for E s ) and P90 (0.36 for E F and 0.24 for E s ) values of E v after 60 years' injection based on the studies from Myshakin et al., 2018 ( Table 4).
The practical estimation of total gas storage is 6.0-15.8 Gt with 2.9-7.6 Gt free gas and 3.1-9.1 Gt adsorbed gas. It is noted that the exact practical values will need to be refined when the reservoirs-scale efficiency factors are measured in future studies and the ethically and economically stored volumes will consequently be left for those studies.

Heterogeneity at multiple scales
The heterogeneity analysis relies on the targeting features and the corresponding scales. The variation of formation depths and thickness at basin-scale, 19 lithofacies and microfacies at core-scale, 21 and porosity, permeability and mineral compositions of bulk samples 20 can be found in many other studies. Based on these understandings, the heterogeneity is further extended to sub-cm-scales in this study through the distribution of mineral  clusters and matrix at 1-10 mm scale, individual grains of minerals and organic matter at 100 nm-1 mm scale, macropores and microcracks at 10-100 nm scale, and mesopores at 1-10 nm scale. These features are crucial to understand the pore-structure, flow pathways, and gas behaviors for accurate prediction of carbon sequestration and enhanced gas recovery. For the same component, increased resolution of images leads to a lower fractal dimension owing to more details of complexity and roughness quantified. At 1-10 mm scale, each component is present as assemblies without clear boundaries of individual grains, therefore the fractal dimensions do not vary dramatically between components (B2.5 to 2.6). Individual grains of organic matter and minerals can be distinguished at 100 nm-1 mm scale, and the fractal dimensions, representing for the surface roughness of these grains, decrease to B2.3-2.5. Pores can be only observed at LR-and HR-nanoscale, and the fractal dimension of macropores and fractures is relatively high when the voxel size is 50 nm and it decreases to 2.2-2.3 and when the voxel sizes reach 20 nm and the mesopores are taken into account.
The degree of anisotropy varies among different components. Clay mineral shows higher degree of anisotropy at both LR-and 100 nm-1 mm scales compared with other components, owing to the elongate geometry and the preferred orientation to the bedding planes. It is as high as 0.758 at high-resolution microscale, much higher value than 0.181 for granular minerals and 0.171 heavy minerals which have relatively regular crystal geometries. Organic matter pores show slightly higher degree of anisotropy than minerals pores, owing to the different mechanism of origin.

Potential flow pathways and storage of CO 2 and CH 4
A full understanding of the pore system is the basis of understanding flow pathways and storage volumes. The mechanism of CH 4 and CO 2 flow transport in shale reservoirs can take several forms (e.g. Darcy flow, surface diffusion and molecular diffusion), and each depends on the pore sizes and porenetwork size. 38,39 CO 2 is present as supercritical state at the reservoir conditions in Bowland shale, which has low viscosity and high diffusivity like gas but density similar to liquid CO 2 . 40 The variations of subsurface condition factors such as temperature, pressure and moisture status will impact the storage volume, flow pathways and fluid behaviours in the subsurface. The adsorption and permeability were both measured at reservoir pressures to simulate the subsurface environment, in order to provide better estimates. The pore systems of Bowland Shale samples are primarily associated with mineral matrix (60.2%) with wedge shapes or within organic matter (34.2%) with spherical shapes in terms of total pore numbers. Three size ranges of pores were detected in this study playing different roles for the fluid pathways. The micro-cracks and pores above 100 nm, shown in the PFIB images, contribute 47.5% in total pore volume and many of them exist at the interface of organic matter and minerals with relatively elongate shape (Fig. 5). These pores have extremely low connectivity themselves but they can be connected to smaller pores. They may act as flow porosity, in the form of Darcy flow, allowing the CH 4 and CO2 exchange between different compositions. The CH 4 /CO 2 displacement efficiency in these pores changes relatively rapid. 41 Therefore once the CO 2 reaches these pores, the CH 4 adsorbed on the pore surface can be extracted more easily. The injected CO 2 flow travels between mineral matrix and organic matter through these pores and the extracted CH 4 gas (both free and adsorbed gas) is transported from the organic matter to mineral matrix. The pores between 10 nm to 100 nm, majority of which were shown in the FIB images, contribute 35.9% for the pore volume and therefore provide the storage space for the free molecules of CH 4 and CO 2 and the pathways to the smaller pores (less than 10 nm size). Pores below 10 nm, identified by the nitrogen adsorption, SEM and HIM images, are associated with organic matter and clay minerals. They contribute over 99.9% of cumulative surface area and may provide the potential pathways in the alternative forms of molecular diffusion cross pores and adsorption-desorption at the surface of pores. [42][43][44] The adsorption capacity is also the highest in these pores. The capillary condensation phenomenon may be able to be characterised in these pores using relevant techniques. Additionally, 5.6% of pores (0.4% of the sample volume) are identified as intra-minerals (Fig. 1e), which are bounded in the mineral crystal, are likely to contribute little to the overall storage and flow transport for either CH 4 or CO 2 due to the extremely low connectivity relationship to other pores. In addition, the artificial fractures may connect many pores with different sizes and provide Darcy flow pathways for CH 4 and CO 2 .
The CO 2 and CH 4 flow transport in shale reservoirs are complex. Based on the previous understanding e.g. ref. 2, 45 and 46-48 and the results in this study, a flow pathways and storage mechanism is proposed below (Fig. 9). At the beginning of the CO 2 injection, high concentration of CO 2 will start the transport from the induced fractures (20-80 mm) to the preexisting microcracks (1-20 mm) and macropores (20 nm-1 mm), as a form of free molecules. Then CO 2 transports into the mesopores and micropores (o20 nm) in the forms of free molecules and surface diffusion, stored in the open pore space and adsorbed to the pore surface. As the concentration of CO 2 increases, CH 4 desorbs and CO 2 adsorbs, and CO 2 penetrates to smaller pores gradually. The surface diffusion takes more weight in the finer pores. When the transport reaches the equilibrium, the CO 2 may be stored in shales as the following forms: hydrodynamic trapping of buoyant CO 2 remains as a mobile fluid, adsorption trapping to the surface of the host rock, solution trapping of dissolving CO 2 in the formation water, mineral trapping due to mineral precipitation; and the residual or capillary trapping into an immobile fraction.

Feasibility of simultaneous CH 4 extraction and CO 2 storage in Bowland shale
The Bowland Shale meets the criteria for CO 2 storage in shale proposed by the US-DOE: 32 (1) hydrocarbons can be commercially produced from the shale through stimulated fracture networks, which has been established by numerous studies e.g. ref. 19 and 20. The gas in place is estimated to be 822 (P90) -1329 (P50) -2281(P10) tcf. 19 (2) The Bowland shale reservoir is found at a depth 41500 m, which can maintain CO 2 in a supercritical state at the reservoir's pressure and temperature conditions. These conditions were attained in this study for the adsorption experiments with CO 2 and CH 4 . (3) A suitable seal system is in place to prevent the vertical migration of CO 2 to the surface. The extremely low vertical permeability of Bowland Shale measured in this study, 3.9 Â 10 À21 m 2 (B3.9 nD, two orders of magnitudes smaller than the horizonal permeability), indicate that CO 2 will be retained in the subsurface once injected. Sufficient pore volume for free molecules and pore surface areas for adsorption are two important parameters to assess the potential of CO 2 storage in shales. Based on the estimation in this study, B40% of the CO 2 will be stored as free molecules and B60% will be stored in as adsorbed molecules in the Bowland shale. The adsorption of these gases largely depends on the organic and clay content in shale, which at the nanoscale is reflected in distinct surface functionality. 49 Oxygencontaining groups including -COOH, -H and -OH bonds, have been shown to have a severe impact on the adsorption of CO 2 and CH 4 on organic matter and mineral surface at reservoir conditions. 50 Inorganic vibrational bands of quartz, illite and kaolinite (Si-O, Al-O-H and Si-O-Al bonds) as well as carbonates can exhibit higher CO 2 adsorption than hydrocarbons (C-H) on mineral surfaces. 51-53 CO 2 injection causes a reduction in the intensity of C-O and C-H bonds within organic matter in shales, 54 leading to an increase in hydrocarbon recovery from organic matter. 55 In practice, clay minerals have a significant impact on gas adsorption on shale, 56 although in the presence of water the organic content reveals a stronger correlation with gas adsorption than clay minerals. 57 In the Bowland shale, the clay mineral content is relatively low with an average value 20% and the TOC value normally varies between 1% and 6%. The selected organic-rich (silicate and carbonate 75%, kaolinite 7%, muscovite 10% and TOC 6%) and lean (silicate and carbonate 65%, kaolinite 18%, muscovite 5% and TOC 1%) samples are typical for the whole Bowland shale formation. 20,21 In our study, we observed that CO 2 adsorption is approximately 8 and 10 times larger than CH 4 adsorption in BOR and BOL, respectively. This is in line with other studies reporting 3-10 times greater adsorption of CO 2 than CH 4 in shales. 58,59 Importantly for this application, enhanced CH 4 gas recovery could result from CO 2 injection in Bowland shale. Furthermore, the low clay mineral content in Bowland Shale is expected to minimise the impact of clay swelling following water injection, meaning that the differences between water and CO 2 injection observed in the Longmaxi shale formation (Song et al., 2019) are likely to be significantly smaller here.
The cumulative porosity (including micro-cracks) is calculated to be 7.8% and the cumulative surface area is 1.79 m 2 g À1 in the selected samples (Fig. 8c). We have estimated a theoretical maximum storage capacity in Bowland shale of 47.5-49.5 Gt, including 19.0-21.2 Gt CO 2 (38-45%) stored in the open pore space and 28.3-28.5 Gt CO 2 (55-62%) adsorbed onto the pore surfaces. Additional porosity can be produced by artificial fracturing (4.8% in the 20 mm diameter sample considered in this study), leading to an increase of storage space for free molecules. The enhancement of the storage capacity will therefore also depend on the implemented fracturing job.

Overall assessment of simultaneous CH 4 extraction and CO 2 storage in different shale plays
Compared to many successfully commercial gas plays in North America, Bowland shale is similar in composition (low clay mineral content), depth (41000 m), and water saturation (o35%). While the porosity of Bowland shale is lower than Marcellus and Haynesville shale, it is similar to values reported for Ohio shale, New Albany shale and Barnett shale. This may lead to a reduced gas recovery using conventional extraction methods (e.g. hydraulic fracturing), but the overall CO 2 storage volume is still high ( Table 5).

View Article Online
The theoretical CO 2 storage value of Bowland shale (47.5-49.5 Gt CO 2 ) is comparable to many commercial shale plays in the United States. For Marcellus shale, the CO 2 storage capacity was estimated to be 171 Gt (99 Gt of adsorbed gas and 72 Gt of free gas), while 340 Gt, 72 Gt and 21 Gt CO 2 were estimated fro Utica shale, Antrim Shale and Devonian Ohio shale, respectively. 3 A total of approximately 28 Gt of CO 2 was estimated to be potentially stored in Ohio shale and New Albany Shale. 58 The average CO 2 storage capacity of Bowland shale is 5.4 CO 2 kg per t, which represent the lower limit of the range reported for many shale plays (5 to 10 kg per t). 60 Additionally, the average recovery rate of shale gas is 25% with a CO 2 enhanced rate 7% which makes the final CO 2 recovery gas 32%. Meanwhile, the technically accessible CO 2 storage is around 32% in Marcellus shale. 35 The estimated total gas inplace estimate (P50) in Bowland shale is 1329 tcf (B37.6 tcm). Based on technical values estimated above, 438.8 tcf (B12.4 tcm) methane may be recovered by CO 2 injection and 6.5 (P10)-15.8 (P90) Gt CO 2 may be stored practically. This is equivalent to over 10 years' CO 2 emission in the UK according to the statistical figure (0.364 Gt = 80% Â 0.455 Gt) in 2019 published by UK government. 76 Some inherent limitations to the extraction operation include (i) the low viscosity of sc CO 2 (1/4 to 1/5 of water) leading to a poor sand-carrying ability; 40 and (ii) the low matrix permeability (3.9 Â 10 À21 m 2 normal to beddings and 3.3 Â 10 À19 m 2 parallel to beddings) which introduces severe mass transfer limitations. 61 Technical developments such as cyclic injection of carbon dioxide, and huff and puff operations, may be considered to further increase the efficiency in field operations. Compared with water, CO 2 injection require a lower pressure for fracturing and generates more complex fracturing networks, 62 while variations between induced microseismicity during water or CO 2 injection requires careful studies.

Conclusions
(1) This paper qualitatively and quantitatively investigates the 3D microstructure and pore system of the Bowland shale and characterise the possible flow pathways and storage capacity of CH 4 and CO 2 in shales. It demonstrates that the multi-scale 3D characterization is a powerful tool to assess the feasibility of CO 2 storage in shale and enhanced gas recovery and this can be used for further technical development to enhance the energy efficiency.
(2) The representative features at corresponding scales were identified, including fractures induced by fluid injection, sample fabric, individual grains, macropores, mesopores and micropores at 10-100 mm scale, 1-10 mm, 100 nm-1 mm, 10-100 nm, and 1-10 nm respectively. The particle and pore anisotropy is heterogeneous for the different phases and pore sizes at the same scale, and significantly changes across the scales from the mm to nm.
(3) The cumulative porosity is calculated to be 7.8% while the cumulative surface area is 1.79 m 2 g À1 . Pores associated with mineral matrix within organic matter occupy 60.2% and 34.2% of the total porosity respectively with another 5.6% intra-mineral pores which are not connected to the whole pore network. A pore size distribution spanning the range of 2 nm to 3 mm shows pore above 100 nm, between 10 to 100 nm and below 10 nm contribute 3.7%, 2.8% and 1.3% of the sample volume respectively, and pore below 10 nm contribute over 99.9% of the cumulative surface area.
(4) At the study conditions, a theoretical maximum of 49.5-49.5 Gt carbon storage can be stored in Bowland shale, consisting of B40% as free molecules in the open pore space and B60% as adsorbed molecules to the pore surfaces. The technical and economical storage of CO 2 in Bowland shales may be less (e.g. 6.0-15.8 Gt based on simple estimations) and the exact numbers can be estimated in future studies when reservoirs-scale factors are measured. The practical value might be equivalent to over 10 years' CO 2 emission in the UK.

Conflicts of interest
There are no conflicts to declare.