An introduction for non-experts on using X-ray micro computed tomography as a tool for pore scale digital subsurface characterisation of siliciclastic materials

Abstract

This paper presents an overview of using X-ray micro computed tomography (μCT) as a valuable tool for micro scale investigation of siliciclastic materials. When processed using digital image analysis (DIA), valuable quantitative data can be extracted from μCT 3D images. Subsurface reservoirs are of great importance to society as fluid-bearing formations, but also as storage reservoirs for carbon dioxide. μCT imaging has the capability to perform preliminary, highly detailed investigations of potential reservoirs. This approach has a range of benefits when compared to traditional 2D techniques, such as optical and scanning electron microscopy (SEM). Key advantages include the technique being non-destructive and capable of 3D and 4D visualisation. This facilitates rapid repeated digital measurements and experiments on microstructures. Digital samples can also be readily shared within the scientific community to replicate results and quickly launch new investigations. However, limitations still exist, posing challenges to the wider application of such a methodology. Such limitations include the identification of a representative elementary volume (REV), computational cost, and suitable processing of the output image data. Here, we highlight the value of using μCT and DIA, from our first experiences, to facilitate pore scale siliciclastic reservoir characterisation, but also highlight our perceived limitations and barriers to its much wider application. This paper introduces the key processing stages, opportunities and limitations of these techniques.

Lay summary

Subsurface characterization of potential reservoirs is a key aspect of the current energy transition phase. A digital approach based upon μCT imaging, digital image analysis (DIA) and digital rock physics (DRP) are techniques to make the initial assessment of reservoir deposits at the microscale to help inform subsequent decisions on further investigation. Application of the proposed digital methodology allows for investigation of the relationship between porosity and permeability, the upper percolation threshold, and a range of other pore body and throat characteristics. The aim of this contribution is to provide non-experts with a basic guide on the potentiality and limitations of these techniques.

1. Overview

X-ray micro computed tomography (μCT) is an imaging technique that uses the differing X-ray attenuations of materials to visualise them. The technique involves firing X-rays at a sample through a range of cross sections to map the X-ray attenuation throughout a sample (Cnudde & Boone, 2013 and references therein). These attenuation maps are reconstructed as a stack of equally spaced 2D greyscale images, effectively representing a 3D sample volume, with individual 2D pixels becoming effective 3D voxels. The differences in X-ray attenuation of a phase are directly related to its density, which is used to distinguish between materials (Ketcham & Carlson 2001; Blunt et al., 2013). Darker pixels in the image represent phases of lower density, attenuating fewer X-rays, whilst brighter pixels represent phases of greater density, attenuating more X-rays. To extract valuable data from these images, a series of steps must be performed through the technique of Digital Image Analysis (DIA). Over the past decade, lab-based micro computed tomography has gained popularity for tackling pore scale problems in a geological context (Blunt et al., 2013; Cnudde & Boone, 2013, Wei et al., 2014; Bultreys et al., 2015, 2016; Menke et al., 2016, 2019; Singh et al., 2017; Thomson et al., 2018, 2020a, b; Hertel et al., 2018; Yu et al., 2019; Fei et al., 2019; Garfi et al., 2020; Payton et al., 2022b, a, 2021, 2022c). To begin with, images are often cropped to remove voxels not belonging to the study sample to reduce the effects of processing artefacts and limit computation times. Cropping can also be employed to remove extremities of the 3D volume, which are especially prone to beam hardening artefacts. Beam hardening causes the outer pixels of an image to appear erroneously brighter than the centre. This is due to the average strengthening, or hardening, of the polychromatic X-ray beam as it passes through the sample (Ketcham & Carlson, 2001). Images are then typically processed, using a variety of filters (Ketcham & Carlson 2001; Buades et al., 2008; Thomson et al., 2018; Garfi et al., 2020; Payton et al., 2021) in order to remove image artefacts and noise (Ketcham and Carlson 2001; Cnudde & Boone 2013), which improves the sharpness of phases and features (Figure 1). Images in a stack may also be normalised so that segmentation is reliable throughout the sample.

Figure 1. Example of a raw μCT image alongside the filtered equivalents. A non-local means filter is effective in removing most of the background noise and will suffice in many cases. The use of a median filter may be beneficial in other instances, particularly where grain boundaries are of significance.

These pre-processing steps are an important stage in the workflow to enable reliable data extraction. From here, the process of allocating image voxels to discrete phases (i.e., segmentation) is performed, facilitating quantitative analysis (Ketcham & Carlson 2001; Iassonov et al., 2009; Campbell et al., 2018; Thomson et al., 2018). Given the direct impact that segmentation has on the output data, it is often considered the most important step in DIA.

Following segmentation, a range of processes and challenges can be investigated, including reactive transport (Menke et al., 2016; Payton et al., 2022a), microstructural analysis (Wei et al., 2014; Menke et al., 2016; Hertel et al., 2018; Thomson et al., 2018, 2020b; Fei et al., 2019; Payton et al., 2021, 2022c), single phase flow (Wei et al., 2014; Bultreys et al., 2015; Thomson et al., 2018, 2020b; Menke et al., 2019; Payton et al., 2021, 2022a, c) and multiphase flow (Blunt et al., 2013; Bultreys et al., 2015; Singh et al., 2017; Menke et al., 2019; Garfi et al., 2020). Microstructures extracted from μCT images of natural samples can be used as the basis for a range of analyses and investigations into properties and features such as porosity, permeability, grain shape, pore geometry and tortuosity. Significantly, this approach allows for digital representations of microstructures to be used as physical domains for numerical modelling of processes such as reactive transport and multiphase flow (Figure 2).

Figure 2. An example of a 3D pore structure extracted from μCT images. The triangular faces of tetrahedral volume elements are shown, which facilitate use as a modelling domain. The red-orange colours in the domain are an example from a numerical model of carbon precipitation from carbonated pore water, where red and white represent the greater and lesser presence of precipitated carbon, respectively.

μCT imaging and subsequent image analysis are especially effective for quantitative porosity and permeability investigation. Traditional imaging approaches include optical microscopy, scanning electron microscopy (SEM) and electron backscatter diffusion (EBSD), whilst the traditional, go-to methods for quantitative analyses include helium pycnometry and mercury intrusion. Traditional imaging techniques can provide valuable measures of porosity but have limitations in that they are 2D in nature and are destructive to a certain degree. Meanwhile, pycnometry and mercury intrusion have the limitations of requiring laboratory facilities, provide no visual representation of the pore structure, and are time-consuming when running many and repeated experiments. With that said, they do sample larger representative volumes than μCT imaging. Further to this, the use of mercury is dangerous and of environmental concern.

Each of the limitations highlighted in non-digital or traditional approaches are overcome by μCT coupled with DIA. μCT imaging is non-destructive, which enables easy repeat experimentation whilst preserving the material for use with other study techniques (Payton et al., 2022d). Furthermore, μCT imaging enables 3D analysis is not only quantitative but also qualitative, as study volumes are easily and effectively visualised. In order to conduct a preliminary pore scale characterisation, no further laboratory time or equipment is required after imaging. Instead, a typical laboratory workstation with a suitable amount of memory, usually around 2–3 times more than the data volume to be processed, can effectively perform basic analysis quicker and cheaper in the long term than repeated laboratory experiments.

μCT alongside DIA is a valuable technique for micro scale siliciclastic reservoir characterisation but comes with some limitations. There are a number of barriers stopping this method from addressing more challenges. These barriers include identifying a representative elementary volume (REV) (Al-Raoush and Papadopoulos 2010; Mostaghimi et al., 2013; Jackson et al., 2020; Huang et al., 2021) for upscaling of observations. A REV is the smallest possible volume that allows for suitable representation of larger volumes of the same material. However, due to the small scale of μCT study volumes, determining a suitable REV is challenging; a challenge that exists across all characterisation techniques at this scale. Furthermore, the processing of the output data (Iassonov et al., 2009; Shi & Yan 2015; Chauhan et al., 2016; Furat et al., 2019; Leonti et al., 2020) is often contentious due to the potential for user bias. The process of segmentation is directly reliant on filtering and user input in terms of either choice of automatic thresholding algorithm (Iassonov et al., 2009; Bultreys et al., 2016) or implementation of a manual threshold (Bultreys et al., 2016; Thomson et al., 2020a, b; Payton et al., 2021). This manifests from the image acquisition process itself that can introduce artefacts, such as beam hardening, but also the voxelised nature of digital images and the size that these voxels can be. Whilst there is an increasing availability of automated techniques to perform segmentation based on unsupervised machine learning (Chauhan et al., 2016; Andrew 2018; Furat et al., 2019; Purswani et al., 2020; Alqahtani et al., 2022), in our experience of working with quartz-dominated samples, we found that a manual approach was effective due to clear phase separation. In such cases, user bias may be considered constant across samples and therefore, results are considered relative rather than absolute (Cnudde & Boone, 2013).

Whilst μCT imaging has opened up the opportunity to use real microstructures as physical modelling domains (Jiang & Tsuji 2014; Payton et al., 2022a), computational cost needs to be taken into account. A trade-off exists between the level of detail and accuracy in a model, and the amount of computational time and power required.

As these issues are yet to be satisfactorily addressed in their entirety, the application of μCT to some scientific challenges is still limited. Such challenges include characterising microporosity (Thomson et al., 2019), using true pore geometries for continuum scale modelling and upscaling micro scale observations.

The goal of this article is not to provide a comprehensive review of μCT imaging and DIA, but instead to share the experience of a non-expert with these techniques and how they can be used in sedimentology. We aim to highlight their value and limitations without entering into too much detail, which can be found instead from specified references, providing a first-aid-style reference for researchers approaching this technique – not for specialists in the topic. We hope that our experiences will help other new users with their approach to adopting these methodologies.

2. X-ray micro computed tomography images

To allow for a digital representation of a geological sample to be as useful in terms of measurements and observations as a physical sample, all details at all scales need to be accounted for. The quality of the imaging output is limited by the size of the pixels or voxels that can be acquired to make up the image itself. The voxelised nature of the acquired images means that features will never display perfect curved boundaries, e.g., owing to the cubic nature of the voxels (Fei et al., 2019; Payton et al., 2022a, b). The higher the level of detail acquired through smaller voxels, the more information can be acquired about the physical sample, tending toward a completely accurate digital representation (Zhan et al., 2010) (Figure 3). However, greater image resolutions result in much larger datasets and can only be acquired for much smaller sample volumes (Razavifar et al., 2021). Consequently, a trade-off must be found whereby an acceptable resolution is acquired for the purpose of the study without generating an unmanageable amount of image data to be stored and processed.

Figure 3. A schematic representation of how voxel size directly influences the accuracy of feature edge digitisation in μCT images.

μCT images acquired from siliciclastic sandstone samples typically have voxel sizes of between 1 and 5 μm³ (Payton et al., 2021, 2022c, b, a). For such sample types, it is reasonable to use voxel sizes in this range given that sand grains are those classified as having a diameter of >63 μm (Wentworth 1922). Consequently, the work presented on porosity and permeability must be considered to represent sample macro porosity rather than the true total porosity. As highlighted by Thomson et al. (2019, 2020b), in the Brae Formation sandstone, microporosity contributing to connectivity can be present beyond the resolving capability of μCT imaging.

3. Digital image analysis and processing

A common approach to measure permeability is to use a numerical model simulating single-phase flow (Thomson et al., 2018, 2019, 2020b; Payton et al., 2021, 2022c). The computation time for these simulations can be substantial, at minimum taking five days and at maximum three weeks on a workstation, running 6 hyperthreaded Intel Xeon W-2133 CPUs with clock speeds of 3.60 GHz supported by 120 GB of memory, performed by Payton et al. (2021, 2022c, b). A direct modelling approach, such as the application of the finite element method on a detailed domain mesh, provides the greatest degree of accuracy in a result. This is opposed to other indirect techniques, such as pore network modelling, where a less detailed domain is used but at the cost of computation time (Blunt et al., 2013; Bultreys et al., 2016). This means that running multiple test cases per sample at varying scales may not be feasible due to time constraints. Consequently, exploration of indirect methods could be beneficial, such as using pore network models (PNMs) that have run times far smaller (Blunt et al., 2013; Varloteaux et al., 2013a, b). PNMs have been shown to offer a good degree of agreement in their permeability results when compared with direct numerical modelling approaches (Varloteaux et al., 2013b). Therefore, using a PNM-based approach could allow for further permeability investigation in a given time period, requiring a fraction of the time for a small decrease in accuracy.

The process of segmentation is of great importance owing to the direct impact that it has on the final output results acquired from μCT images. Segmentation is the process by which image voxels are assigned to a discrete phase. There are a number of approaches that may be taken to carry out this task, including thresholding, edge-based, region-based, watershed, clustering-based and machine learning segmentation. Here we focus on the simplest method, thresholding segmentation for initial phase identification (Kaur & Kaur 2014). This process is always prone to user bias or fails to work consistently across different materials in the cases of manual and automatic segmentation, respectively (Bultreys et al., 2016; Campbell et al., 2018). Automatic segmentation, performed by an algorithm, is used to provide comparability between samples and reliability (Payton et al., 2021, 2022a), however, it is apparent that this possibly comes at the cost of accuracy. Meanwhile, manual segmentation, performed by a user, provides greater accuracy at the cost of reliability, which changes between samples and users (Thomson et al., 2019, 2020b; Payton et al., 2021, 2022c, b). A wide range of automatic segmentation algorithms exist to tackle images with different types of features and levels of noise (Ridler & Calvard 1978; Otsu 1979; Huang & Wang 1995; Iassonov et al., 2009), However, it is unrealistic to find the best automatic algorithm for each study sample, as the definition of ’best’ is difficult. The accuracy of segmentation will always be questioned owing to the direct impact it has on any output measurements (Bultreys et al., 2016; Garfi et al., 2020). Steps can be taken to estimate and quantify segmentation error as described by Iassonov et al. (2019), however, often results within a study are considered to have constant bias and are therefore, relative rather than absolute (Cnudde & Boon, 2013).

With specific regard to the digital analysis of individual grains in sedimentary materials, the ability of watershed segmentation to both accurately and reliably distinguish between individual features is of great significance (Figure 4). Watershed segmentation conceptually converts a 2D greyscale image into a 3D topographic surface based on greyscale values. This allows for topographic lows to be ‘filled’ with water and the point at which two lows spill into one another, the watershed, is where a feature boundary is defined. Payton et al. (2022b) explored processes that can aid the accurate and reliable segmentation of individual grains, highlighting the benefit of a combination of Non-local Means (NLM) (Buades et al., 2008) and median filter (Gupta 2011). The NLM filter operates by searching the wider image for pixels that are not just similar to the target pixel, but also have similar windows of pixels around them. The target pixel is then averaged according to the other similar pixels and windows found. Meanwhile, the median filter replaces the greyscale value of each pixel with the median value of the pixels in a given local neighbourhood area of that pixel. Whilst different filtering algorithms are available (Mirmozaffari, 2020), we found that the use of these image filters resulted in a significant improvement in grain segmentation quality compared to prior work (Thomson et al., 2020a). However, when comparing the segmented grain volume with the raw μCT images, the generated grain boundaries can still be considered inaccurate in some instances. For example, some grain boundaries (see Figures 4b and 4d) appear to split up what could be seen as a single grain in Figure 4b. Whilst the combination of a NLM and median filter provides acceptable results, it is apparent that accuracy could be further improved. This stems from a combination of image voxel size, image filtering and available segmentation algorithms.

Figure 4. The sequence of processes required to perform effective segmentation of individual grains within a single phase. (A) and (B) show the effect of the non-local means and median filters, respectively. The watershed algorithm is employed on the granular phase in (C) and those grains that are not intersecting the volume boundaries are highlighted in (D). (E) Shows the 3D output following individual grain segmentation.

4. Numerical modelling

To achieve an acceptable degree of accuracy without excessive computational cost, simplification is required. μCT images can be used to produce highly detailed physical numerical 3D modelling domains (Madonna et al., 2012; Bultreys et al., 2016; Mostaghimi et al., 2016; Shams et al., 2021; Payton et al., 2022a). Using such a detailed modelling domain results in a large number of mesh elements on which to perform numerical computations at a high cost. Consequently, investigations may be simplified in terms of the chemical system by considering fewer reaction pathways (Payton et al., 2022a) or fewer physical flow processes. In doing so, a model becomes less accurate, highlighting the barrier that computational resources pose to extracting the most value from this approach.

The process of mesh generation from μCT images to produce a modelling domain often involves a degree of modification to the triangulated grid representing the microstructure, especially in 3D (Schöberl 1997, as shown in Figure 2). This can result in the intersection of mesh elements, which causes issues with the convergence of a numerical model. These usually arise from sharp and irregular clusters of elements on the mesh surface, which are not effectively dealt with during mesh building (Lo 2002). Measures that can be taken to address issues with the mesh include manual intervention, mesh decimation or simplification and adaptive mesh refinement. Payton et al. (2022a) highlighted that a degree of simplification through smoothing of the surface mesh was effective in allowing their model to successfully converge on a solution. The smoothing process fundamentally changes the volume that has been segmented to a small degree, introducing a source of uncertainty regarding the accuracy of the pore structure. As with image filtering, which also changes the data, some is necessary to improve the quality of the output, but too much can result in degradation of the output quality.

5. Representative elementary volumes

Determination of a representative elementary volume (REV) is a limitation with regards to upscaling, which manifests in many pieces of micro scale work. Identification of a REV allows for observations and measurements to be made, which are then applicable throughout a given area, volume or geological formation (Al-Raoush & Papadopoulos 2010; Mostaghimi et al., 2013; Jackson et al., 2020; Huang et al., 2021). Owing to the heterogeneous nature of the subsurface environment, a REV is usually quite large, in order to factor in these different heterogeneities, and cannot always be found within a typical μCT volume (Al-Raoush & Papadopoulos 2010). Due to the scale on which μCT imaging takes place, finding a REV of a substantial area is challenging.

6. Summary and future outlook

X-ray micro computed tomography (μCT) imaging paired with digital image analysis (DIA) is an effective tool in addressing a range of pore scale challenges in Earth Science. Despite the proven success in applications, including microstructural analysis, flow modelling and geophysical property investigation, there remain limitations.

At the core of using this technique to address far-reaching geoscience problems is the issue of upscaling from a representative elementary volume (REV). This limits the wider applicability of results acquired using this technique at present. Further work is necessary to allow for micro scale studies to produce value far beyond pore scale understanding, up to the continuum scale. This is an active area of current research, with a range of methods to achieve upscaling currently proposed and developed. Such approaches include use of fractal-scaling (Munawar et al., 2021), implementation of the Brinkman equation (Wei et al., 2021), a combined Darcy-Brinkman-Stokes approach (Menke et al., 2021), pore network model stitching (Kohanpur & Valocchi 2020), correlation with whole core imaging (Hertel et al., 2018), and percolation theory (Liu and Regenauer-Lieb 2011). The focus of these studies on tackling relatively small upscaling from the pore scale highlights the size of this challenge. If scaling relationships based on any of these approaches can be derived, which are widely applicable and allow the continuum scale to be accessed, this would add huge value to μCT imaging.

Another barrier to expanding the use of this technique is the amount of data acquired and the associated computational limitations. If larger volumes of material are to be imaged in such high levels of detail, then there must be resources available to store and process this data. With the current rapid growth in cloud computing use in research, it stands to reason that the widening accessibility to high performance computing facilities and large volumes of cloud storage could allow for this issue to be addressed. Increasing availability of computational resources means that image pre-processing and analysis can be carried out more rapidly, expediting the process of answering research questions. This encompasses the ability to grow the complexity of numerical models based upon μCT microstructures. At present, the physical grid itself requires substantial computing power to deal with, even before the complexity of the model itself is considered. Widening use of cloud computing facilities and larger supercomputing facilities may overcome this barrier, allowing for more detailed and larger grids to be combined with less simplified mathematical approaches to simulating reactive transport or multiphase flow.

Finally, the importance of accurate and repeatable image segmentation is very clear. At present, manual segmentation approaches appear to dominate in the geosciences, owing to the user’s ability to tackle difficult areas resulting from natural heterogeneity. Automatic segmentation algorithms are often not effective and may only be applicable to a certain type of material. To address this topic, a range of machine learning (ML) approaches are emerging and continue to be improved (Andrew 2018; Purswani et al., 2020). The ability to produce effective ML models rapidly and easily would allow for a combination of the benefits of automatic and manual segmentation. This would mean models that are accurate, directly comparable and repeatable between samples could be applied for segmentation.

Acknowledgement

The authors thank Andrew Tye (BGS) for his comments on an early version of the manuscript, and two anonymous reviewers and the AE for constructive comments improving our contribution. Authors would like to thank Frank Lehane (RHUL) for his computational support throughout the project and Brett Clark (NHM) for carrying out μCT analysis. This manuscript was published with the permission of the Executive Director of the British Geological Survey. RLP and DC acknowledge support from a NERC DTP studentship (Grant Number NE/L002485/1) as well as further financial support through a CASE partnership with the British Geological Survey as part of their British University Funding Initiative. AK acknowledges support from the British Geological Survey via NERC national capability.

Conflict of interest

The authors declare no conflict of interest.

Data Availability Statement

Data sharing is not applicable to this article as no datasets were generated or analysed.

Authors Contributions

Conceptualisation: R.L.P., D.C., A.K.; Methodology: R.L.P.; Formal analysis and investigation: R.P.; Writing—original draft preparation: R.P.; Writing—review and editing: R.L.P., D.C., A.K.; Supervision: D.C., A.K.