Simulating permeability reduction by clay mineral nanopores in a tight sandstone by combining computer X-ray microtomography and focussed ion beam scanning electron microscopy imaging

Computer X-ray microtomography (µXCT) represents a powerful tool for investigating the physical properties of porous rocks. While calculated porosities determined by this method typically match experimental measurements, computed permeabilities are often overestimated by more than 1 order of magnitude. This effect increases towards smaller pore sizes, as shown in this study, in which nanostructural features related to clay minerals reduce the permeability of tight reservoir sandstone samples. Focussed ion beam scanning electron microscopy (FIB-SEM) tomography was applied to determine the permeability effects of illites at the nanometre scale, and Navier–Stokes equations were applied to calculate the permeability of these domains. With these data, microporous domains (porous voxels) were defined using microtomography images of a tight reservoir sample. The distribution of these domains could be extrapolated by calibration against size distributions measured in FIB-SEM images. For this, we assumed a mean permeability for the dominant clay mineral (illite) in the rock and assigned it to the microporous domains within the structure. The results prove the applicability of our novel approach by combining FIB-SEM with X-ray tomographic rock core scans to achieve a good correspondence between measured and simulated permeabilities. This methodology results in a more accurate representation of reservoir rock permeability in comparison to that estimated purely based on µXCT images.

1 Introduction

Depositional environment and subsequent diagenetic alterations are two key factors that influence the bulk mineralogical composition and the authigenic clay mineral inventory of a reservoir (Wilson and Pittman, 1977; Worden and Morad, 1999), and therefore the fluid-flow properties of the porous rock. A well-established technique to image and analyse rapidly the 3D physical properties of porous rocks is computer X-ray microtomography (µXCT) combined with the concept of digital rock physics (Andrä et al., 2013a, b; Okabe and Blunt, 2004). By applying monochromatic synchrotron radiation, it is possible to overcome conventional µXCT artefacts like beam hardening and problems that arise due to limited phase contrast, lack of resolution and edge preservation, as well as low signal-to-noise ratios (Brunke et al., 2008; Kling et al., 2018; Lindquist et al., 2000; Mayo et al., 2015; Spanne et al., 1994). Synchrotron-based µXCT scans with image resolutions in the order of 1 µm can provide a sound basis for flow and transport modelling of tight sandstones as suggested by Peng et al. (2014). They found that synchrotron µXCT imaging is necessary for tight sandstones when the connectivity of the pore space is low and pore throats cannot be resolved using a conventional µXCT scanner. They further concluded that a high abundance of the smallest resolvable pores falsifies modelled permeabilities due to an overestimation of actual pore sizes. Several studies have shown sub-micrometre pore structures to be a frequent feature of tight reservoir rocks (Jiang, 2012; Shah et al., 2016; Soulaine et al., 2016). Most of these nanostructures are related to different types of clay minerals; most commonly illite, kaolinite, chlorite and smectite (e.g. Wilson and Pittman, 1977; Worden and Morad, 1999; Desbois et al., 2016).

Although known for years, considering such structural features below µXCT resolution in pore-scale models remains challenging (Alyafei et al., 2015; Guan et al., 2019; Liu and Mostaghimi, 2018; Liu et al., 2020; Menke et al., 2019; Peng et al., 2012). Soulaine et al. (2016) systematically analysed the effect of sub-resolution domains with varying permeabilities on the simulated permeabilities of Berea sandstone (20 vol % porosity, 2 vol % sub-resolution domains) and found that calculated permeabilities can be reduced by up to 50 %, if microporous domain permeabilities converge towards zero. Thus, it is evident that neglecting sub-resolution information can lead to a significant overestimation of rock permeability in such simulations (e.g. Saxena et al., 2018, 2017). Menke et al. (2019) utilized a multi-scale Brinkman area approach, applying different permeabilities for each microporous domain to simulate flow in mono-mineralic carbonate rock. They showed that Stokes–Brinkman models are in good agreement with experimental data, whereas Stokes and/or Navier–Stokes models alone were not able to predict permeability in a conventional flow scenario. They also demonstrated that for pure carbonates, a direct correlation can be established between observed density contrasts and specific physical properties, such as porosity and permeability. However, this approach is not applicable in a system with more than one rock-forming mineral, such as a tight sandstone, where density contrasts relate to different mineral phases as well as sub-resolution porosities.

The lack of distinct material information for a voxel is often ascribed to the partial volume effect (e.g. Kessler et al., 1984; Ketcham and Carlson, 2001). To overcome the impact of this effect on physical properties of a rock, imaging techniques that can resolve the pore structure at different length scales have to be applied. For estimating the permeability of reservoir rocks, the resolution achieved by synchrotron radiation imaging lies within an acceptable range (Saxena et al., 2018). Several studies have demonstrated that, by combining X-ray and scanning electron imaging, the pore space of tight clay-bearing rocks can be spatially resolved from the millimetre down to the nanometre scale (e.g. Desbois et al., 2016; Hemes et al., 2015; Markussen et al., 2019).

In this study, we aim to demonstrate a new approach by combining synchrotron-based µXCT imaging with focussed ion beam scanning electron microscopy (FIB-SEM) to improve flow simulations in a tight sandstone formation with high clay mineral content. First, we used machine learning-based image segmentation to enhance pore space segmentations of artefact-rich FIB-SEM topologies. Then, we conducted Navier–Stokes simulations on FIB-SEM topologies. Finally, we subsequently used these simulations as input data for sub-resolution domains in µXCT-based Stokes–Brinkman models (Brinkman, 1949; Neale and Nader, 1974). This novel morphology-based approach for sub-resolution-rich materials results in simulated permeabilities that fit experimental results significantly better than using Navier–Stokes simulations alone.

2 Sample preparation and characterization

In this study, a well-characterized illite-bearing Upper Rotliegend sandstone from Germany was used, which was sampled from the quarry Schwentesius near Bebertal village (Heidsiek et al., 2020). This well-known location exposes an analogue of the Permian gas reservoir sandstone of the Flechtingen High, which formed part of the North German Basin. Thin section analyses of this reservoir rock have shown that illite is the main clay mineral, which primarily occurs as a coating along detrital grains, illite meshworks grown on coatings, diagenetically altered K-feldspar and illitized lithoclasts (Fischer et al., 2012).

Samples were taken from a large sandstone block displaying a variety of different aeolian and fluvial depositional facies (Fig. 1a). We identified several facies denoted A to J from four different types of sedimentary depositional layers. The samples of the different facies were numbered from left to right. Samples were drilled out and extracted from marked locations in the form of plugs with a diameter of 2.5 cm. The plugs were cut into three segments which were used for X-ray diffraction analyses, helium porosimetry measurements and FIB-SEM imaging. Mini plugs were drilled directly beside the plugs with a diameter of 3 mm and a length between 10 and 20 mm. The larger plugs were used for mineralogical and geochemical analysis, and to measure permeabilities experimentally. The mini plugs were scanned at the synchrotron radiation-based µXCT imaging beamline P05 at the PETRA of DESY Hamburg (Germany). After the synchrotron measuring campaign, the mini plugs were additionally examined by FIB-SEM and energy dispersive X-ray spectroscopy (EDX) imaging to obtain qualitative and quantitative information about the clay mineral particles found within the rock pore space. Microporous structures in the Rotliegend sandstone sample could be resolved by comparing µXCT and FIB-SEM images. The term “microporous” refers to the definition of sub-micrometre porosity by Soulaine et al. (2016), who differentiated between void, solid and microporous voxels in µXCT images. According to for example Tinet et al. (2020), the term “nanoporous” is used to describe structures with predominant pore sizes in the nanometre range (0.2–1000 nm). In this range, the clay minerals often appear as a greyish pore-filling phase in µXCT images.

Figure 1 (a) Upper Rotliegend sandstone block showing four main deposition facies. (b) Sampling locations and sizes of the extracted plugs (green) and the mini plugs (red) extracted.

3 Analytical procedures

3.1 Mineralogical characterization

X-ray diffraction analyses using the Rietveld analysis program Profex 4.0 (Doebelin and Kleeberg, 2015; Ufer et al., 2012) have shown a homogenous mineralogical composition along the layers of the sample block with only slight variations in content (Fig. 2).

Figure 2 XRD analysis of the main minerals in the plugs based on Rietveld calculations, determined for several samples of facies A to J of the sandstone block.

The main components are quartz (58 wt %–69 wt %), authigenic and diagenetic feldspars (12 wt %–20 wt %), calcite (1 wt %–18 wt %), and illite (10 wt %–17 wt %). In the <0.2 µm size fraction, we observed traces of swelling smectite interpreted as contaminants from surface weathering. Accessory hematite and barite particles also occur with abundances of <1 wt %. In the Bebertal tight sandstone, illite is by far the most abundant clay mineral and makes up ∼ 95 wt % of the <2 µm fraction. Commonly, the nanometre-scale microstructural features of reservoir rocks containing significant amounts of clay minerals are usually not detectable at µXCT resolution due to a low absorption contrast (Ahmad et al., 2018).

3.2 Permeability measurements and pore size distribution

A helium gas-driven permeameter under steady-state conditions was used to experimentally obtain permeability values for the individual plugs (e.g. Filomena et al., 2014). Compressed air was used to apply a pressure of 10 bar to the core samples that were coated by a latex membrane. A Bronkhorst F-111C-RBB-33-V was used as a flowmeter. It was calibrated by polynomial calibration and has a measuring accuracy of ±0.5 % plus ±0.1 % of the full scale. To measure the differential pressure, an Emerson Rosemount 3051S Ultra was used with an accuracy of ±0.045 %. To control the inflow and outflow, two pressure sensors (Brooks 5866) with an accuracy of ±0.5 % were installed. For permeability measurements, the inflow and outflow pressures of the helium flux were sequentially increased in up to six pressure steps. The differential pressure was kept constant at 200–500 mbar depending on the sample properties. The intrinsic sample permeability was derived from the apparent gas permeability, K_g, determined for each pressure step using Darcy's law (Liu et al., 2017):

$$\begin{array}{}\text{(1)}& K_{\mathrm{g}}={\displaystyle \frac{\mathrm{2}Q{p}_{\mathrm{2}}\mathit{\eta }L}{A(p_{\mathrm{1}}^{\mathrm{2}}-p_{\mathrm{2}}^{\mathrm{2}})}},\end{array}$$

where Q is the measured gas flow rate, η is the dynamic viscosity of the permeant, L is the sample length, A is the sample cross section, and p₁ and p₂ are the inflow and outflow pressures. By plotting K_g against $\frac{\mathrm{1}}{(p_{\mathrm{1}}+p_{\mathrm{2}})/\mathrm{2}}$, the data can be fitted by a straight line. The intercept of the best-fit line at the K_g axis corresponds to the intrinsic sample permeability, K_int (Gao and Li, 2016; Klinkenberg, 1941). Also, MICP measurements were conducted with an Autopore IV Series (Micromeritic Instrument Corp.) to determine the pore size distribution of a dried sub-sample with a weight of ∼ 2.5 g, which was taken from a cross-bedded aeolian layer of the sandstone block. Based on the capillary law, MICP enables the analysis of a wide spectrum of pore sizes (3 nm to >900 µm), corresponding to a pressure range of 0–414 MPa. As a non-wetting liquid with a high contact angle (130–140^∘), mercury only penetrates a pore when pressure is applied. Under the assumption of cylindrical pores, the applied pressure is directly proportional to the pore throat diameter as described by the Washburn equation (Washburn, 1921):

$$\begin{array}{}\text{(2)}& D={\displaystyle \frac{-\mathrm{4}\mathit{\gamma }\cdot \cos \mathit{\theta }}{P}},\end{array}$$

where D is pore throat diameter, γ is the surface tension, θ is the contact angle and P is the applied pressure.

3.3 Synchrotron-based µXCT

For synchrotron tomography, a beam energy of 29.87 KeV was used while the sample-detector distance was 1.2 cm. The effective image resolution of the detector equipped with a CCD camera was 1.22 µm per pixel, while the image size was 3056×3056 pixels. We used an advanced reconstruction script described in Moosmann et al. (2014) with the MATLAB^® software and binned the images by a factor of 2 before reconstruction to increase the signal-to-noise ratio. This decreased the voxel resolution to 2.43 µm and changed the image size to 1528 × 1528 pixels. The number of projections was 1200, with the information of five subsequent images used to calculate an average for every projection image. After reconstruction of the 3D image stacks, the scans were denoised using the non-local means filter of the GeoDict 2020 software package (Buades et al., 2011). Image segmentation of the mini plugs was realized by conventional greyscale thresholding. The greyscale values for each phase on a 16 bit intensity range were 0–8550 for pores, 8551–12 880 for grains, 12 881–25 000 for cements and 25 001–65 535 for high-density cements and oxides. Since the best threshold values for each phase varied slightly from sample to sample, these values were adjusted for each segmentation within a small range of ±50. A qualitative comparison with machine learning segmentation revealed a better pore-to-solid segmentation and resolving of small pore throats by thresholding (Fig. 3). Since the main goal was to achieve the best possible permeability estimation, the differentiation between pore and solid is more important for permeability estimation than the accurate segmentation into different phases (Khan et al., 2016). Leu et al. (2014) point out that even a small variation in pore throat morphology can have a large impact on the estimation of permeability.

Figure 3 (a) Greyscale 2D µXCT image of a mini plug. (b) Segmented phases using multi thresholding. (c) Segmented phases using the machine learning image classification module of the software ilastik (Version 1.3.3) by Berg et al. (2019).

For permeability simulations, we either used the fast Fourier transformation (SIMPLE-FFT) or the left-identity-right (LIR) solvers, both implemented in the FlowDict module of the GeoDict software package (Linden et al., 2015; Moulinec and Suquet, 1995). While the SIMPLE-FFT solver is fast for calculating low-porosity domains, the LIR solver is better suited for high-porosity domains and requires less memory. Both iterative finite volume solvers can apply Navier–Stokes and Navier–Stokes–Brinkman equations. Derived from these equations, Darcy's law (Eq. 1) is used to calculate the permeability of a material (Darcy, 1856):

$$\begin{array}{}\text{(3)}& \mathit{u}=-{\displaystyle \frac{K}{\mathit{\eta }}}(\mathrm{\nabla }p-\mathit{f}).\end{array}$$

In Eq. (3), u is the three-dimensional average fluid-flow velocity, K is the permeability, η is the fluid viscosity, p is the intrinsic average pressure tensor and f is the force density field, which was defined using the Navier–Stokes conservation of momentum equation for all three dimensions (Eq. 4):

$$\begin{array}{}\text{(4)}& -\mathit{\eta }\mathrm{\Delta }\mathit{u}+\left(\mathit{u}\cdot \mathrm{\nabla }\right)\mathit{u}+\mathrm{\nabla }p=\mathit{f}.\end{array}$$

The Brinkman term can be added to the Navier–Stokes equation where porous voxels are required. These voxels include the nanoporous flow resistivity:

$$\begin{array}{}\text{(5)}& -\mathit{\eta }\mathrm{\Delta }\mathit{u}+\left(\mathit{u}\cdot \mathrm{\nabla }\right)\mathit{u}+\mathit{\eta }{K}^{-\mathrm{1}}\mathit{u}+\mathrm{\nabla }p=\mathit{f},\end{array}$$

where K⁻¹ is the inverse of the permeability tensor and ηK⁻¹ the flow resistivity. The applicability and robustness of combining Navier–Stokes equations with the Brinkman term has been validated by Iliev and Laptev (2004).

We calculated the permeability with symmetric boundary conditions in tangential and flow direction with a pressure drop of 20 370 Pa. The symmetric boundary conditions are valid for low-porosity structures with non-periodic pore throat geometries. The differential pressure value was set to be able to compare the results with helium permeation flux measurements where similar values have been used. As a convergence stopping criterion, a low error bound of 0.05 was chosen. The error bound criterion uses the result of a previous iteration to predict the final solution by linear and quadratic extrapolation. This stops the iteration when the relative difference regarding the prediction is smaller than 5 %. During this process, the solver recognizes oscillations and local minima or maxima which prevents a premature stopping of the solver during the iteration process.

3.4 FIB-SEM measurements

In this study, a Zeiss Auriga crossbeam electron microscope equipped with a Gemini electron column and an Orsay Physics ion beam was used. SEM images were taken at 1 kV with an in-lens secondary electron (SE2) detector, and FIB slicing was executed with a beam current between 0.5 and 2 nA and a voltage of 30 kV. This resulted in a slice thickness of 25 nm. A large FOV of ∼ 20 µm could be reached. To derive structural information from the FIB-SEM images, extensive post-processing of the data was required. Following image alignment and cropping, stripes and shadow artefacts were filtered out before image segmentation. The Slice Alignment operation of the module ImportGeo of the GeoDict 2020 software package was used to align the images, while the Curtaining Filter was used for stripes correction (Fig. 4). In general, the segmentation of pores in FIB-SEM images was not straightforward since scans of porous polished sections are pseudo-2D and contain information from behind the actual imaging plane (De Boever et al., 2015).

Figure 4 (a) Raw SEM image of secondary illite growth in a porous feldspar with streaks and shadow artefacts. (b) Images with aligned slices of a cropped region of interest. (c) Filtered greyscale SEM image of the illite meshworks. (d) The binarized result after image segmentation with pixel classification algorithms of the software ilastik.

As multi-thresholding and watershed segmentation algorithms have problems with shine-through artefacts (Prill et al., 2013), capturing the correct 3D pore space geometry is of crucial importance for the determination of a realistic permeability. Recent advances have shown that machine learning image segmentation software can successfully be utilized to segment pore space in CT scans (Berg et al., 2018).

The software ilastik, an interactive learning and segmentation toolkit by Berg et al. (2019), was used for the interactively supervised machine learning segmentation of the phases in our FIB-SEM images (Fig. 4). The built-in pixel classification module groups the probability of each pixel to be assigned to a phase according to their different imaged features. The software provides a set of different images based on features derived from the original image (filtering and smoothing, edge detection, etc.). Examples of the correct class (phase) provided by the user builds a decision surface in feature space. Based on these features, a random forest classifier assigns each pixel or voxel to a certain phase. This method works for both 2D and 3D data. The accuracy of this method increases with the number of user-provided training data. In a manually controlled workflow, it was possible to reach high segmentation accuracies with only minor over- or underestimations of the pore space (Fig. 4c, d).

3.5 Defining microporous domains

The need to define microporous domains results from the mismatch of permeability between µXCT simulations and gas-driven permeameter tests. While the simulation of permeability in structures with high permeability and porosity obtained by µXCT scans is precise, the effect of nanoporosity below resolution on permeability increases with decreasing permeability (Pittman and Thomas, 1979; Saxena et al., 2018). When comparing backscattered electron (BSE) images with µXCT images of the same slice, it becomes apparent that the smallest pores in µXCT images simplify the real pore structure (Fig. 5). Furthermore, SEM and EDX images revealed that most of the pores are filled with clay minerals. Since both void and microporous regions share similar greyscale values, it is impossible to correctly differentiate upon segmentation. In this approach, the segmentation of the pores includes clay minerals with a low absorption contrast. To determine the distribution of illite in µXCT scans, we used correlative µXCT, SEM and EDX measurements. For this, a 3 mm plug was embedded into epoxy resin and then ground and polished until the region of interest was reached. Two sites were chosen for EDX mappings (Fig. 5). Comparing EDX mappings with CT images shows that the distribution of illite agrees with the textural findings of Fischer et al. (2012).

Figure 5 (a) µXCT and corresponding (b) BSE image of the polished section including the close-up regions S1 and S2 that were used for EDX mapping to analyse the distribution of illites. (c) Ternary plot of Ca, K and Al content of each pixel in the mappings. Calcite, illite and K-feldspar plot in three distinct regions. By selecting these areas, a phase map overlay was created that visualizes the location of these phases in SE images. (d) Region S1 and S2 with phase map overlays. Illite is predominantly found in the vicinity of altered feldspars, as coating along grain boundaries and as authigenic pore filling.

Furthermore, it becomes apparent that illite enrichments coincide with regions that are usually referred to as void pore space in µXCT images. The mismatch between real pore structure and segmented pore space is highest in small pores and throats. An average size of ∼ 5 µm was measured for illites which cover the surfaces of the quartz and feldspar grains. This corresponds to approximately two voxels in the µXCT images. To refine flow and reduce the influence of overestimated pore sizes in these specific regions, we define all pores with a diameter ≤2 voxels as microporous domains. The Brinkman term accounts for small pores, where grains are porous themselves (Brinkman, 1949). To extract these regions from the initial pore space segmentation F_p∗ we calculated Euclidian distance maps as used by Maurer et al. (2003):

$$\begin{array}{}\text{(6)}& d=\sqrt{\left(x_{\mathrm{2}}-x_{\mathrm{1}}\right)+\left(y_{\mathrm{2}}-y_{\mathrm{1}}\right)+(z_{\mathrm{2}}-z_{\mathrm{1}})},\end{array}$$

where d is the distance between two points; x₁, y₁ and z₁ are the coordinates of the first point; and x₂, y₂ and z₂ are the coordinates of the second point. By selecting d(≤2) all pore voxels with the nearest distance of ≤2 voxels to the next solid surface are masked, including the outer rim of larger pores (Fig. 6b). To obtain a mask with pores ≤2 voxels, this rim layer surplus must be removed. This was realized by performing a morphological dilation operation on d(>2) using a structure element SE2 of 2³ voxels. By subtracting the dilated image from d(>2), the microporous domain data F_μ2 (Fig. 6c) can be obtained.

$$\begin{array}{}\text{(7)}& F_{\mathit{\mu }\mathrm{2}}=d\left(\le \mathrm{2}\right)-{\mathit{\delta }}_{\mathrm{SE}\mathrm{2}}\left(d(>\mathrm{2})\right)\end{array}$$

Figure 6 (a) Segmented µXCT image of the pore space of a mini plug. (b) Image with applied microporous domain using the Euclidian distance map approach. (c) Image with removed outer rim layer of larger pores.

In GeoDict the permeabilities of the microporous domains are calculated with the Brinkman term. Based on flow simulations on 3D FIB-SEM images of illite meshworks, an isotropic permeability was assigned to the microporous domains (Fig. 6). Yoon and Dewers (2013) confirmed the validity of the approach of extrapolating structural features of clays measured by FIB-SEM to the pore scale. The complete workflow for a precise permeability simulation is illustrated in Fig. 7.

Figure 7 Workflow for microporous domain modelling combining FIB-SEM nanoporosity with µXCT scans. Sample E9 of the aeolian layer with well-distributed flow paths was selected for visualization. Upon the modelling of the microporous domains, the total fluid-flow velocity decreases. Displayed simulated permeabilities include the mean values of five samples with the best distribution of percolation paths. Fluid-flow velocities close to zero are transparent.

4 Results and discussion

4.1 Permeability measurements and pore size distribution

Four to six pressure stages were applied to the plug samples to determine their gas permeability. In Table 1, parameters which were used for the permeameter test are shown for sample F8. These values yield a Klinkenberg factor of 0.56 and a quadratic regression coefficient of 0.94.

Table 1 Parameters of the gas permeameter flux experiment (sample F8), including measured helium gas flow rate, recorded inflow and outflow pressures as well as calculated apparent gas permeability for each pressure stage. For evaluating the intrinsic sample permeability, the sample cross section (0.025 m) and length (0.05 m) as well as the dynamic viscosity of helium (18.4 × 10⁻⁶ Pas) is also needed.

Overall, the measured intrinsic permeabilities range between 1.1 and 5.4 mD with an average value of 2.9 mD for both aeolian facies (Table 2; Fig. 1, sample series E and F). The quadratic regression coefficient ranges from 0.732 to 0.998, which indicates a high accuracy of the permeameter test. The small variation of the physical rock properties between the main layers is induced by different angles of the grain layering in each plug sample and small-scale variations in grain and pore size distributions. The observed variance is in the typical range of observed permeability fluctuations for tight reservoir sandstones (Lis-Śledziona, 2019). Considering the observed permeability, the studied rock samples are at the lower end of the permeability range known for sedimentary reservoir rocks (Gluyas and Swarbrick, 2004).

Table 2 Intrinsic permeabilities and porosities of the measured sandstone plug samples.

The pore size distribution of the sub-sample analysed by mercury intrusion porosimetry (Sect. 3.2) was obtained by a semi-logarithmic representation of the normalized intrusion volume achieved per pressure interval (Fig. 8). The consumed capillary stem intrusion volume of 27 % was within the acceptable range needed for precise measurements. We discovered a major peak at approximately 1 µm, representing the most common pore throat diameter.

Figure 8 Distribution of the pore throat diameters of a sample of the aeolian facies in layer J obtained from mercury intrusion porosimetry (MICP).

In addition to the mercury porosimetry results, calculations of the pore size distribution of FIB-SEM and µXCT mini plug images were conducted using the GeoDict module PoroDict. This module provides an algorithm that virtually pushes spheres of different sizes into a medium to determine the 3D pore size distribution (Münch and Holzer, 2008). FIB-SEM measurements of the pore size distribution cover the small pores related to the illite nanoporosity, while µXCT porosimetry illustrates the pore size distribution of the larger-scale intergranular pore space skeleton.

The µXCT pore volume analysis shows a constant rise towards the smallest observed pore diameters (2.44 µm; Fig. 9). The steepest slope is observed between 7.3 and 12.2 µm, which indicates that the largest volume of pore bodies is observed in that range. The calculated pore size distribution from nanoporous 3D structures of illite meshworks obtained by FIB-SEM imaging shows three distinct peaks at 75, 125 and 250 nm. This illustrates that the most common inscribed pore diameters are relatively small compared to the actual extent of the pores, which is considered typical for slit-shaped pore systems where the pore axial ratio is high (Desbois et al., 2016). Ultimately, this results from the different thicknesses of clay platelets and heterogeneity in the alignment of the illite nanostructure (Aplin et al., 2006).

Figure 9 Pore throat diameters from µXCT (a) and FIB-SEM images (b).

An apparent gap exists between diameters observed by MICP and 3D imaging with FIB-SEM and µXCT. While MICP peaks at around 1 µm, diameters observed by µXCT start at 2.4 µm and the largest inscribed diameters observed by FIB-SEM are below 700 nm. However, it must be noted that pressure-controlled MICP generally gives information about the dimension of pore throats, whereas pore size distributions provided by imaging techniques give information about pore body volumes (Zhao et al., 2015). Furthermore, pore shielding may cause an underestimation of larger pores for MICP (Abell et al., 1999; Gane et al., 2004). The occurrence of authigenic illites is the likely cause of this effect as they are commonly found in pore throat areas. Since the Washburn equation assumes ideal pore throats of cylindrical shape, the underestimation of larger pores becomes more evident with the increasing complexity of the pore throat system at both the millimetre and nanometre scale (Washburn, 1921). Therefore, a direct comparison between the used methods is unlikely to result in compatible results.

4.2 Permeability simulation

A calculated isotropic permeability of $\mathrm{2.53}\times {\mathrm{10}}^{-\mathrm{3}}$ mD was used for the microporous domains based on the Navier–Stokes fluid-flow simulations of permeability of FIB-SEM scans of the illite meshworks. The number of porous voxels resulting from clay mineral modelling ranges between 3.3 vol % and 7.1 vol % of the total structure volume for all µXCT scans considered. A comparison of the modelled clay minerals in µXCT scans with XRD mass balancing highlights a large difference between the mineral abundances measured (Table 3, Fig. 2).

Table 3 Modelled clay mineral content within the microporous domains in µXCT scans.

While the mean amount of clay minerals based on XRD measurements was 12.7 wt % (about 11.3 vol % within a structure with 8 % porosity), an average amount of 5.4 vol % was modelled by the distance map algorithm. This is expected since the illite content inside grains was not modelled since it has no effect on permeability. We simulated permeability of 12 mini plug samples that were scanned by µXCT and compared them with measurements from gas permeameter experiments. As a first step, we extracted and illustrated the 10 largest and 100 largest open flow paths through the pore space of all mini plug cores before the modelling of the microporous domains (Fig. 10). This yields information concerning the heterogeneity of the flow fields and allows the validity of the Navier–Stokes simulations to be checked. Structures with flow impingement often cause numerical problems, which results in an artificial underestimation of the permeability simulations. Significant underestimations of permeability after clay mineral modelling were also found in areas where the percolation paths in the samples were limited to a few voxels in the structure before modelling.

Figure 10 Comparison of percolation paths in two µXCT reconstructions of the 3D imaging data without microporous domain modelling. (a) Well distributed percolation paths in a µXCT reconstruction of mini plug sample E9 without microporous domain modelling. (b) Constricted percolation paths with flow impingement limited to a small region of the structure in mini plug sample I3 without microporous domain modelling. The 10 largest percolation paths are coloured in purple; the 100 largest percolation paths are coloured in cyan.

This effect leads to an artificial permeability drop, which renders calculations less precise. Since evenly distributed flow paths are necessary to determine the true permeability of a volume of a rock, we considered only permeability calculations of samples which show no flow impingement for modelling (Bear, 1972; Leu et al., 2014; Zhang et al., 2000). The validity of the simulations was ensured by applying the variance algorithm on the flow fields discussed in Jacob et al. (2019).

The mean value of the experimentally obtained intrinsic gas permeabilities was 2.9 mD (Table 2, Fig. 11a). Samples with no flow impingement had an initial mean permeability of 26.5 mD, whereas the mean permeability simulated on µXCT images was 26.5 mD (Fig. 11b). With applied microporous domains, a mean permeability of 1.9 mD was calculated (Fig. 11c).

The mismatch between measured and simulated permeabilities could be decreased to 1 mD (−34.5 %) compared to 23.6 mD (+813.8 %) prior to the modelling. Hence, our approach significantly improved the match between measured and simulated permeabilities and lies within the standard deviation of the permeabilities measured by the gas permeameter (Table 2). Histograms of the simulated fluid-flow velocities indicate a strong decrease in velocities due to the modelling of the microporous domains in sample E9 (Fig. 12). Correlation histograms of the velocities depict differences induced by the modelling. While identical fluid-flow field velocities would plot as a single straight line, the velocities, in this case, show straight lines with varying slopes. This indicates a general decrease in the fluid-flow velocities with a splitting of different flow path velocities. Points which plot as a line represent a main fluid-flow path with a direct correlation of the velocities between the original and the modelled structure's flow field. Furthermore, a wider spread of the distribution of high velocities results from microporous domain modelling. This indicates narrow pore throats where fluid-flow velocities are locally enhanced compared to the structure before the modelling.

Figure 11 Comparison of measured and simulated permeabilities. The boxes indicate the upper and lower quartile of derived permeability values. Median values are coloured in red; mean values are coloured in black. (a) Measured permeabilities using the gas permeameter test. (b) Simulated permeabilities using µXCT data. (c) Simulated permeabilities using µXCT data with modelled microporous domains.

Figure 12 Comparison of unmodelled µXCT images and modelled microporous domain images (a). Correlation histogram prior to and after the modelling that shows a decrease in the main fluid-flow velocities, while few velocities after the modelling were increased with a wide spread of distribution (b).

Based on our combined analytical and numerical study, further research may help to increase the accuracy of simulated permeabilities even further. Since isotropic permeabilities of the microporous domains were applied to µXCT images, the accuracy can be improved by taking the anisotropy of clay mineral fabrics and surface topology into account. This can be done by applying anisotropic permeabilities in the calculation of the microporous domains. While this study showed a good match between the experimental and simulated permeability, the need to include heterogeneities of clay mineral layering to improve the simulations was depicted in Villiéras et al. (1997).

The presented approach might also be applicable to other rock types such as mudrocks, carbonates, etc., if the preconditions are met. The examined rock properties (e.g. permeability) should be based on the same scale dependencies as the rock analysed in this work (e.g. Grathoff et al., 2016).

5 Conclusions

Overall, the outcome of this study shows that combining µXCT and FIB-SEM imaging with numerical models constitutes a valuable and novel approach for determining physical properties of clay-bearing tight reservoir rocks. Considering the high number of accessible pores in the scans, the phenomenon of flow impingement was mainly attributed to the unresolvable nanoporosity. While permeability, which is one of the most important reservoir properties, is often determined by simulations based on µXCT scans of small samples taken from a field-scale reservoir, we were able to demonstrate that an accurate estimation for clay-rich and low-permeability rocks is only possible if nanoscale porosity is also included. Thus, our simulations using the Euclidian distance map approach resulted in an improved match with stationary gas permeameter measurements in contrast to permeability simulations merely based on unmodified µXCT images. Adopting this multi-method approach, we increased the accuracy of simulated permeabilities of samples measured by µXCT. These results have important implications for improved modelling of reservoirs relevant to gas and water applications. A realistic simulated permeability of a tight reservoir sandstone could only be achieved by appropriate modelling of the nanoporosity related to matrix clay minerals (illite) that occur below the µXCT resolution. The simulated permeability based on combined µXCT and FIB-SEM images and modelled microporous domains showed good agreement with the experimental results. Obtaining an even distribution of the simulated fluid-flow paths through the sample without flow impingement was necessary to obtain an accurate permeability estimation from 3D imaging. Resolving the nanopore structure and distribution of clay mineral-related features by the combined analytical and numerical modelling approach represents a tool for achieving a more accurate understanding of the fluid-flow behaviour within tight sandstones, with direct relevance to predicting the injection, storage or extraction of gas or water in a reservoir rock. Our multi-method approach can be applied to determine more accurate permeability values and flow paths for reservoir rocks with high clay mineral content if direct experimental measurements are not successful. Hence, future studies should focus on distinguishing the different morphologies of clay minerals and their related anisotropic effect on rock permeability. The permeability of the nanoporous structures depends highly on the layering of the clays and their spatial orientation on the grain surfaces and within feldspars. This approach should include a variety of various sedimentological facies, also with high porosity and permeability, to investigate whether clay mineral modelling is also a valid tool for such sedimentary rocks.