Comparing organic-hosted and intergranular pore networks: topography and topology in grains, gaps and bubbles

Abstract The relationship between pore structures was examined using a combination of normalized topographical and topological measurements in two qualitatively different pore systems: organic-hosted porosity, common in unconventional shale reservoirs; and intergranular porosity, common in conventional siliciclastic reservoirs. The organic-hosted pore network was found to be less well connected than the intergranular pore network, with volume-weighted coordination numbers of 1.16 and 8.14 for organic-hosted and intergranular pore systems, respectively. This disparity in coordination number was explained by differences in the pore shapes that are caused by variations in the geological processes associated with the generation of the pore network. Measurements of pore shape showed that the pores in the organic-hosted network were both significantly more spherical and had a more positive curvature distribution than the pores present within the intergranular network. The impact of such changes in pore shape on pore-network connectivity was examined by creating a suite of synthetic pore geometries using both erosion/dilation of the existing network and image-guided object-based methods. Coordination number, Euler characteristic and aggregate porosity analyses performed on these synthetic networks showed that organic-type pore networks become connected at much higher aggregate porosities (35–50%) than intergranular-type pore networks (5–10%).

The last five years has seen a dramatic change in the nature and genesis of the petroleum systems driving the worldwide hydrocarbon market, with an increasing emphasis on unconventional resources, particularly artificially stimulated shale reservoirs. Over the last ten years production from onshore shale formations has expanded greatly, increasing from around 1.3 trillion cubic feet (tcf ) per day in 2007 to around 15.8 tcf per day in 2015, contributing approximately 60% of US day-to-day natural gas production (EIA 2017). The recent success of shale gas has caused a dramatic shift in global energy prices by greatly increasing global hydrocarbon supply. In this new low hydrocarbon price environment, technology must be made more efficient and reservoir management more targeted for unconventional petroleum systems to remain competitive relative to mature conventional resources. The flow processes governing recovery are governed at the scale of the tiny tortuous porous pathways through which the flow occurs, and so an understanding at this scale is critical for understanding upscaled reservoir behaviour.
Pore-scale imaging has developed over the last 20 years from a primarily academic technique used to visualize pore structures for fundamental research into an increasingly crucial industrial tool.
High-resolution imaging is today used for the characterization of a range of geological and petrophysical properties, including porosity and pore connectivity (e.g. Fredrich et al. 1995), absolute and relative permeability (e.g. Andrä et al. 2013), mineralogy (e.g. Lai et al. 2015), geomechanical response (e.g. Fonseca et al. 2012), geological facies type, and diagenetic history (Vos et al. 2014), and while much progress has been made, many challenges remain (Blunt et al. 2013).
One particular area of interest is how differences in rock type, geological origin and diagenetic history impact the topology of the pore network, and, thus, single and multiphase flow and transport phenomena. It is these properties which ultimately govern and control variations in hydrocarbon production rates and recovery factors. The goal of this study is to compare two qualitatively different pore systems: those frequently found in conventional rocks, and those frequently found in unconventional reservoirs. Specifically, we examine the intergranular pore structures found in simple conventional sandstones and compare them to the organic-hosted 'bubble' pores frequently described in unconventional shale reservoirs (e.g. Loucks et al. 2012;Milliken et al. 2013;Chen et al. 2016). To perform this analysis we use classic pore-connectivity analysis (in the form of aggregate connectivity, coordination number and Euler characteristic distribution) and pore-shape analysis (in the form of pore sphericity and a novel geometrical curvature analysis), both on the real (imaged) geometries and synthetic geometries created using geometrical and object-based techniques.
An understanding of the contrast between these two pore structures is critical as intergranular pore structures dominate recovery from many classical clean reservoirs (e.g. Øren et al. 1998;Øren & Bakke 2002), whereas organic-matter pores can indicate strong reservoir storage capacities in shale reservoirs (e.g. Jiao et al. 2014;Wang et al. 2016;Ju et al. 2017) and contribute the most space for methane accumulation (Milliken et al. 2013).

Materials and methods
The two rock types compared in this study were Bentheimer Sandstone, an early Cretaceous highporosity, high-permeability quartz arenite (Andrew et al. 2014c), and the Vaca Muerta Shale, a late Triassic-Early Jurassic economic low-intermediate maturity onshore Argentinian shale.
For the sandstone sample, the central portion of a micro-plug 10 mm in diameter and 50 mm in length was imaged using a ZEISS Versa XRM 520 X-ray microscope (Zeiss X-Ray Microscopy, Pleasanton, CA, USA). 3D volumes with a voxel size of 8.96 µm were reconstructed using a proprietary implementation of the Feldkamp-Davis-Kress (FDK) algorithm (Feldkamp et al. 1984) from a series of 1600 projections, sequentially acquired at equally spaced angular increments during a 360°rotation. Each projection consisted of 1024 × 1024 pixels and was reconstructed into a 3D volume of size 1024 × 1024 × 1024 voxels, representing a physical volume of 9.14 × 9.14 × 9.14 mm.
3D imaging of the shale pore network poses a significantly larger challenge than the conventional systems (Saif et al. 2017) due to its geological heterogeneity (Ringrose et al. 2008) at multiple scales. Frequently, the resolution required to resolve fundamental pore structures comes at the expense of a field of view representative of the heterogeneity inherent in real geological systems. This issue is particularly problematic in organic-rich shale reservoirs where pore structures range can from the millimetre to the sub-nanometre in scale. What is more, the organic-hosted porosity frequently governing hydrocarbon flow and recovery occurs at the smallest scales, only accessible using 3D imaging techniques such as focused ion beam-scanning electron microscopy (FIB-SEM). The impact of scale in subsurface shale pore structures is complex (Ma et al. 2017), and thus the scope of this study was limited to examining the structure of organic-hosted 'bubble' porosity (Loucks et al. 2012;Milliken et al. 2013). All multiscale imaging was performed using the Atlas 5.2 correlative imaging platform.
A core-plug end trim 25 mm in diameter and 10 mm in thickness was first mechanically polished before being argon-ion milled for 20 min to produce a sample surface that was flat on the micrometre scale. It has been noted that such milling can change perceived sample maturity by imparting thermal energy to the surface of the sample (Mastalerz & Schieber 2017). The impact of such artificial maturation was reduced by minimizing both the end-trim thickness (10 mm: maximizing conductive heat loss through the ion beam stage) and the milling time (20 min: decreasing the total amount of thermal energy imparted to the sample).
The sample then was imaged in 2D across its entire surface with a pixel size of around 150 nm using a backscattered electron detector (BSD) on the ZEISS Crossbeam 540 (Fig. 1). This enabled the macroscopic heterogeneity of the sample to be fully described, and an organic-rich region to be identified. This sample was then imaged in 3D with a voxel size of 2.5 × 2.5 nm in the X and Y directions, and 5 nm in the Z direction.
One challenge when imaging porous materials using FIB-SEM is that of 'pore backs'an artefact where out-of-plane information is visible, interfering with segmentation. The impact of this artefact is minimized when using a secondary electron (SE) detector. However, this has a much lower sensitivity to material contrasts (such as those between different minerals, or between the organic regions and the inorganic matrix) than a BSE detector. To resolve this problem, each image consisted of two image channels, acquired in parallel from an in-lens SE and energy selective backscatter (ESB) detector. These images were then blended to produce an image which both maximizes contrast and minimizes artefacts (Fig. 2).
This stack was approximately aligned using multiple chevron fiducials, following the technique outlined in Narayan et al. (2014), implemented within Atlas V 3D. A final automated stack alignment was performed in the ORS Dragonfly software package, creating a volume of size 2093 × 1760 × 606 voxels, representing a real physical volume of around 5.2 × 4.4 × 3.0 µm.
Machine-learning image segmentation is a new state-of-the-art technique that partitions challenging image datasets into segments (labels), representing different groups of features of the rock microstructure (e.g. pores or minerals), previously too difficult to separate using traditional techniques (Kan 2017). In this study we use ZEISS Zen Intellesis, an advanced tool which utilizes machine learning to construct a classifier from a range of different features extracted from the image, including local and non-local greyscale, gradient, and texture information. Once the classifier is trained, it can be used for segmenting the same or similar imaging data.
When compared with traditional techniques such as Otsu universal thresholding (Otsu 1979) or watershed-based region-growing algorithms (Jones et al. 2007) frequently used in digital rock-image processing (Schlüter et al. 2014), machine-learning techniques are more artefact and noise tolerant, frequently resulting in fewer misclassified voxels (Chauhan et al. 2016;Andrew et al. 2017). Both the XRM and FIB-SEM images were segmented using this machine-learning software into pore and grain in the case of the XRM dataset, and into pore, organic and mineral grain in the case of the FIB-SEM dataset. The image-processing workflow is shown below in Figure 3.
The resulting connected pore networks of the sandstone and the shale were then separated using a watershed algorithm to create a network of discreet, separated pores and throats (Beucher & Lanteujoul 1979;Rabbani et al. 2014;Andrew et al. 2015). While such algorithms have their limitations for complex pore networks (Raeini et al. 2017), those of relative geometrical simplicity can be separated using watershed relatively effectively (Baldwin et al. 1996;Lindquist & Venkatarangan 1999;Dong & Blunt 2009). These separated pore networks could then be examined for various geometrical, topological and topographical characteristics. First, the size of each pore was measured, both in terms of the volume of each individual pore and the radius of the largest possible sphere within the pore. This was the equivalent to the maximum of a Euclidian distance map evaluated through each pore. Network connectivity was then measured in two ways. A simple evaluation of network connectivity was calculated by measuring the relative contribution to the total porosity made by the largest connected cluster-this is  termed 'aggregate connectivity'. A more detailed analysis of connectivity was then computed by calculating both the distribution of pore coordination numbers (the number of pores connected to each pore) and the volume-weighted coordination number. Finally, pore shape was measured using pore sphericity and curvature.
Pore sphericity (Ψ) is shown by equation (1), where V p is the object's volume and A p is its surface area. Pore sphericity gives a volume-normalized, dimensionless measure of how close a particular component of the pore space is to an ideal sphere (with a sphere having a value of 1.0): Interface curvature has been used extensively to measure the shape of fluid-fluid and fluid-solid interfaces for the measurement of local capillary pressure (Armstrong et al. 2012;Andrew et al. 2014aAndrew et al. , 2015 and grain shape (Andrew et al. 2014b). To measure curvature, a mesh of the segmented pore network was first created using a marching cubes technique. This mesh was smoothed using a Gaussian filter, before the mean curvature was measured at every point across the meshed surface.
One challenge is how to compare interface curvature from pore structures which exist at different scales. To normalize across the different length scales, a dimensionless curvature measurement (C n ) (equation 1) was created by multiplying the mean curvature C of each element on the surface (units m −1 ) by the volume-weighted average pore radius r of the network, where V p is a pore's volume and r p is its maximal inscribed radius. All these parameters were measured for all of the separated pores within each pore network: To investigate the intrinsic differences between different pore networks constructed from elements of different shapes, a suite of different synthetic pore structures was constructed to simulate the sort of pore networks expected in both organic regions (common in shales) and intergranular spaces (common in conventional reservoirs). These networks were constructed in two ways. The first technique involved either eroding (to decrease total porosity) or dilating (to increase total porosity) the existing pore networks until the desired porosity was reached. The second technique used to generate pore networks was an object-based approach where both types of pore network were generated in a way guided by the rock's geological origin. 'Objectbased' approaches to pore-network generation refer to a group of techniques revolving around the treatment of geological bodies (e.g. mineral grains and pores) and processes (e.g. compaction and diagenesis) explicitly. Such approaches have a long and well-established history in digital rock technology (Bakke & Øren 1997;Lerdahl et al. 2000;Øren & Bakke 2002); however, as imaging technologies improved their predominance as the primary tool for the identification of a pore network has declined (relative to direct imaging). They still retain a valuable role, however, in extrapolating and extending understanding generated using direct imaging through a wider range of conditions, for which they are used in this study.
The authigenic 'bubble' shapes present within the organic-hosted pore network were simulated by randomly placing differently sized spheres within a volume. The spheres' volume distribution was given by the measured pore-volume distribution measured from the 3D FIB-SEM image. New pores were placed until the desired porosity was reached and allowed to overlap with existing pores within the network. The intergranular pore network was modelled as the spaces between convex polyhedra, randomly placed into the pore space.
The size distribution of these convex polyhedra was defined by the grain-size distribution of the original geometry, defined using the same watershed technique used for pore-network analysis, only applied to the grain network (Fonseca et al. 2012;Andrew 2015). Multiple total porosity networks were created by randomly placing new grains until the desired porosity was reached. These different networks were then analysed to see the relationship between connectivity and aggregate porosity.
The networks created using object-based techniques were then further analysed, first to assess their statistical similarity to the original images, then to assess the trends present in the suites of synthetic networks. To do this analysis, coordinationnumber distributions were assessed on each of the synthetic networks. To examine trends within the network series, each of normalized frequency distributions was modelled using a Poisson distribution: where F n is the normalized frequency of pores occurring at a coordination number N c . λ, the sole free parameter in the distribution, corresponds to the average coordination number. The Poisson distribution is a useful tool for modelling the distribution of coordination numbers as it matches the observed pore-network coordination-number distributions well, it is asymmetrical (like the observed distributions) and it is simple, aiding an intuitive interpretation.
To mitigate uncertainty in any conclusions drawn on the relationship between network structure and connectivity introduced by the network separation algorithm, the Euler characteristic was also used to assess the evolution of connectivity over the network series. The Euler characteristic is a topological invariant, classically defined for polyhedra, given by: where χ is the Euler number, V is the number of vertices, E is the number of edges and F is the number of Fig. 4. Connectivity in intergranular porosity (a) and organic-hosted 'bubble' porosity (b). The total pore network is shown in green, and the largest connected cluster shown in red. A mostly red rendering corresponds to a well-connected pore network, whereas a mostly green rendering corresponds to a poorly connected network.
faces. As it is topologically invariant, it can be extended to examine the topology of any 3D object by polygonizing the surface of the (segmented) voxelized image. Generally, a more negative Euler characteristic corresponds to a geometry with more connecting strands within it (it is better connected), whereas a more positive Euler number corresponds to a geometry with fewer connecting strands (it is more poorly connected). A major disadvantage of using the Euler number is that it is extremely challenging to normalize between individual instances of differing geometries, and as such it has been largely used when looking at trends within 3D geometries evolving due to some processes (e.g. the evolution of non-wetting-phase topology during multiphase fluid flow: Herring et al. 2013Herring et al. , 2014. This makes it unsuitable for a comparison between the original images on their own but perfect for a comparison between the trends observed in two suites of (synthetically created) geometries. In order to compare the Euler characteristic distribution between the two network types, the distributions were normalized to a range of 0-1 using their maximum and minimum values.
Multiple connectivity metrics (aggregate connectivity, coordination number and Euler characteristic) were used in the network analyses due to their complementary advantages and disadvantages. As such, their agreement gives a greater confidence in any resulting conclusions.
All geometry creation was performed in GeoDict software (Math2Market GmbH).

Results and discussion
The sandstone had a measured image porosity of 22.6%, of which the largest connected cluster contributes 99.98% of the volume. The shale had a measured porosity of only 5.37%; however, this porosity was exclusively hosted within the organic matrix, which made up 25.01% of the rock, meaning the pore space and organics together make up 30.38%. The porosity therefore makes up 17.68% of the volume of the organic-hosted region; however, its largest connected cluster only makes up 6.52% of the overall pore network. This is in contrast to the intergranular network, which was extremely well connected, with the largest connected organic cluster contributing 99.64% to the total pore volume (Fig. 4).
This difference in connectivity is even more pronounced when looking at the distribution in pore coordination numbers (the number of pores that each pore is connected to) (Fig. 5). Organic-hosted pores have a volume-weighted coordinated number of just 1.16, so each voxel within the network      9. Synthetic pore networks created by modelling (a) intergranular pore networks, as the spaces between overlapping convex polyhedral, and (b) organic-hosted pore networks, as overlapping spheres. Total porosity in each case is shown as a percentage. The total pore network is shown in green, with the largest connected cluster shown in red. A mostly red rendering represents a well-connected network, whereas a mostly green rendering represents a poorly connected network. belongs to a pore which is only, on average, connected to just over one pore. In contrast, the intergranular pore network has a volume-weighted coordination number of 8.14, consistent with existing measurements of pore networks in siliciclastic quartz arenites (Andrew et al. 2014c).
The relatively poor connectivity of the organichosted network can be explained by the differences in pore topography (caused by different geological processes associated with pore-network genesis), and the impact that this has on network connectivity. The differences in pore topography (shape) can be investigated by computing two properties from the imagessphericity and pore curvature. The pores within the organic-hosted pore system are significantly more spherical than those within the intergranular pore network, with modal sphericities of 0.6 and 0.3, respectively (Fig. 6). This trend can be understood by examining the normalized curvature distribution of the two pore networks (Fig. 7).
While both intergranular and organic-hosted pore structures exhibit convex and concave regions of the pore space, the organic-hosted normalized porecurvature distribution is significantly more positive than the intergranular curvature distribution, with mean curvatures of 0.95 and 2.14, respectively. In the intergranular pore structure, the concave regions of the pore space correspond to the regions where large polyhedral grains bulge into the pore network, whereas the convex regions of the network correspond to pore corners, where multiple grains are touching, causing the direction of the pore wall to rapidly change direction. This shape, caused by the nature of deposition and diagenesis within the intergranular pore network, leads to the creation of abundant long, thin regions of the pore structure which connect one pore body to the nextpore throats (Fig. 8a). They are formed by the spaces left between matrix grains bulging into the pore space. The organic-hosted porosity, in contrast, has a much more convex (positive curvature) shape, with its surface bulging out into the surrounding matrix. In contrast to intergranular porosity (formed as the negative to the result of a deposition and diagenesis of the grain matrix), organic-hosted porosity is authigenic (forming in place during the maturation of the organic matrix). As the organic matrix is deformable, the hosted fluid tends to form a quasi-spherical shape (explaining the observed trends in pore sphericity). While the formation of this structure is complex, and the pores may then be deformed by other diagenetic processes, they retain their characteristic positive overall surface curvatures. In this case, connectivity only occurs when these pores overlap during growth, leading to a relatively poorly connected matrix (Fig. 8b).
To explore the impact of these differences in pore-body shape, suites of synthetic pore systems were created for both the intergranular and organichosted pore systems using both the dilation and erosion of the existing network and object-based techniques (Fig. 9).
The pore coordination-number distribution was used to assess the statistical similarity between the networks created using object-based techniques and the real (imaged) networks. To do this, a network was generated to a target porosity given by the original (imaged) pore networks, and pores separated and identified, with good agreement between the coordination numbers being found (Fig. 10). Connectivity of the simulated networks was first assessed by measuring the contribution made to the overall porosity by the largest connected cluster (Fig. 11a). The transition between a poorly connected pore network and a well-connected pore network occurs at porosities of between 5 and 10% in the intergranular pore network, and between 35 and 50% in the organic pore network, with the trends created by both network generation algorithms (object based and erosion/dilation) being consistent. The trends in average coordination number and Euler characteristic distribution show very similar results. The characteristic coordination number (Fig. 12a) of the intergranular pore system reaches a value greater than 3 at 5% aggregate porosity, whereas the same number is only reached at a porosity of 50% in the organic-hosted pore system. Similarly, when examining the distribution within the Euler characteristic (Fig. 12b), values transition from their maximum values to their minimum values between 5 and 15% total porosity in the intergranular network, and between 35 and 50% in the organic-hosted network. All metrics of network connectivity are consistent with the originally imaged networks and with each other, showing that organic-hosted pore networks are significantly less well connected at a given porosity than their equivalent intergranular network.
One potential source of additional hydraulic connectivity in real systems is sub-resolution porosity. This is of particular concern in the organic-hosted pore system, where connectivity is poor and pore sizes in such systems have been reported to the nanometre (or even sub-nanometre) scale (Nelson 2009). To assess the impact of such porosity on these systems, the 'resolvedness' of the imaged network was examined by looking at contributions to the total network volume of pores (separated using watershed techniques) of varying equivalent radii (measured in voxels) (Fig. 11b). A dominant contribution of pores at or close to the voxel size would suggest a significant volumetric component of the porosity below the resolution of the image. This analysis showed that for the organic-hosted pore geometry (imaged at the nanometre scale using FIB-SEM) pores around 10 voxels across contributed the most to the total volume of the network (corresponding to an equivalent radius of 12.5 nm). In the sandstone, pores with an equivalent radius of around 13.5 voxels across contributed the most to the total volume of the network (corresponding to an equivalent radius of 121 µm). This indicates that in both systems the pores are well resolved and the lack of connectivity in the shale cannot simply be explained by the finite resolution of the image.
The similarity in connectivity trends between networks generated using different approaches (objectbased approaches and erosion/dilation routines), as well as the statistical similarity (in terms of coordination-number distribution), give confidence that the object-based approaches are effectively capturing the geological processes that are important to the generation of pore-network connectivity. This causal relationship between pore-network generation and connectivity is critical, as it implies that connectivity in authigenic porosity will not occur until much higher porosity values than we would expect in conventional settings. In traditional formulations of multiphase flow, network connectivity is critical to both hydrocarbon recovery (Wardlaw & Cassan 1978;Herring et al. 2013) and single-phase flow (with Navier-Stokes-based fluid flow requiring a well-connected pore structure). These reservoirs, however, do produce in abundance and such a behaviour may be explained by other physical processes such as diffusion through and/or geomechanical deformation of the organic matrix, contributing significantly to overall hydrocarbon transport. The relative impact of these processes would be an interesting target for future work.
The size and extent of organic-hosted pore networks is strongly related to the thermal maturity of the Fig. 12. Connectivity evolution, as measured using the Euler characteristic (a) and the average network coordination number (b), for organic-hosted (red) and intergranular (blue) pore systems. Pore networks generated using erosion/dilation techniques have square markers, and networks generated using object-based techniques have triangular markers. The imaged networks (from which the erosion/dilation series are created) are indicated using circles. host shale, with immature shales (0.5% < R o < 1.5%) displaying few organic-hosted pores, shales of intermediate maturity (1.5% < R o < 3%) displaying a large number of organic-hosted pores and highly mature shales (R o > 3%) displaying a small number of organic-hosted pores (Tissot et al. 1974(Tissot et al. , 1987Wang et al. 2016). As a shale matures, organics are converted into hydrocarbons in situ, with hydrocarbons being principally retained within the organic matrix (Milliken et al. 2013;Tian et al. 2013Tian et al. , 2015. As hydrocarbon generation proceeds, however, a greater proportion of organic matter is converted into hydrocarbon, increasing organic-matter porosity. This gradually increases the connectivity of the hydrocarbon 'bubbles', ultimately causing connected pathways to form and allowing the hydrocarbons to escape from the system. The resulting reduction in pore pressure then causes the collapse and depletion of the organic matter (Zhao et al. 2016;Wang et al. 2018). The relatively poor connectivity of organic-hosted pore networks could explain how relatively high hydrocarbon concentrations can be reached without the hydrocarbons escaping from the systema percolating pathway is not available until high organic-hosted porosities. Once these have been reached (at high maturities) the hydrocarbon is able to migrate away from the system, causing the collapse of the pore network.

Conclusions
The relationship is described between topology, shape and aggregate porosity of two qualitatively different pore networks: intergranular pores, common within simple conventional siliciclastic reservoirs; and organic-hosted pores, common in unconventional shale reservoirs. These pore networks, existing on very different scales, were compared through the use of dimensionless measurements (coordination number, sphericity and normalized curvature). The two networks showed a large difference in connectivity, with a volume-weighted coordination number of 8.14 for the intergranular network and 1.16 for the organic-hosted network. The variation in connectivities is explained by differences in pore shape that are caused by fundamentally differing geological formation processes. This was investigated by examining the connectivity of multiple series of synthetic networks, generated both by the erosion and dilation of the original networks, and using object-based approaches (which implicitly model these processes). The synthetic networks generated using object-based techniques were statistically similar to the original network images. Aggregate connectivity, coordination number and Euler characteristic analyses all showed that the intergranular pore networks were significantly better connected at a given porosity than the organic-hosted network at the same porosity, with intergranular pore structures typically connecting at aggregate porosities of 5-10% and organic-hosted pore structures connecting at aggregate porosities of 35-50%.
Common reservoir rocks exhibit a much wider range of pore systems than the intergranular and relatively simple organic-hosted networks examined in this study, with shales possessing significant additional pore populations that may add significantly to overall connectivity. The techniques for linking real (imaged) pore networks and artificial (simulated) networks could be applied across this broad range of pore geometries to generate a more general understanding of the relationship between a rock's depositional environment, geological origin and diagenetic history.