Abstract
Microcomputed tomography (μCT) and Lattice Boltzmann Method (LBM) simulations were applied to continental carbonates to quantify fluid flow. Fluid flow characteristics in these complex carbonates with multiscale pore networks are unique and the applied method allows studying their heterogeneity and anisotropy. 3D pore network models were introduced to single-phase flow simulations in Palabos, a software tool for particle-based modelling of classic computational fluid dynamics. In addition, permeability simulations were also performed on rock models generated with multiple-point geostatistics (MPS). This allowed assessing the applicability of MPS in upscaling high-resolution porosity patterns into large rock models that exceed the volume limitations of the μCT. Porosity and tortuosity control fluid flow in these porous media. Micro- and mesopores influence flow properties at larger scales in continental carbonates. Upscaling with MPS is therefore necessary to overcome volume-resolution problems of CT scanning equipment. The presented LBM-MPS workflow is applicable to other lithologies, comprising different pore types, shapes, and pore networks altogether. The lack of straightforward porosity-permeability relationships in complex carbonates highlights the necessity for a 3D approach. 3D fluid flow studies provide the best understanding of flow through porous media, which is of crucial importance in reservoir modelling.
1. Introduction
Porosity and permeability control the storage and fluid flow in reservoir rocks. An example of potential reservoir rocks is continental carbonates, such as travertines (a term here used sensu lato [1]), which are highly heterogeneous as a result of their geological evolution, influenced by sedimentary origin, diagenetic processes, and burial history. The latter processes influence the size and shape of pores, producing some of the most complex pore networks recorded in sedimentary rocks. Recent hydrocarbon discoveries highlighted the continental carbonate reservoir potential in the presalt exploration, offshore Brazil [2, 3] and in the Namibe basin, Angola [4–6]. Quantitative data about lithofacies’ occurrence, distributions, and their related porosity and permeability are key to understanding the reservoir behavior. The sedimentology of continental carbonates has been widely studied [7–13], but recently the focus of continental carbonate studies shifted towards the rocks’ petrophysical properties like porosity, permeability, and acoustic velocities [14–22]. Noteworthy is the fact that the permeability inside rocks is strongly dependent on the geometric and topological properties of the porous medium at microscopic scales [23].
Many efforts have been made to reconstruct digital pore networks based on computer tomography (CT) to study the fluid flow in rock samples [24–29]. In this study, the 3D pore networks were acquired with different CT systems, which allowed scanning a range of sample volumes at different resolutions. After segmentation, the obtained pore networks were loaded into Palabos, a software tool for classic computational fluid dynamics (CFD). The source code and scripts are available at the Palabos website [30].
In order to evaluate the reservoir potential of continental carbonates, a critical assessment has to be made regarding the scale at which petrophysical measurements should be performed. This requires covering several spatial scales. Multiscale flow modelling is defined as any method used to explicitly represent the flow properties at more than one scale within a reservoir [31]. Reservoir models typically cover at least twelve orders of magnitude, ranging from pore to core to interwell to full field simulations. Four general scale orders can be recognized, which are pore to lithofacies, lithofacies to geomodel, and geomodel to reservoir model. Multiple-point geostatistics (MPS) gained importance in the recent years for modelling geological structures at different spatial scales, from μm to m [32, 33]. For example, in the field of hydrology, MPS allows building models on decameter scale [34], while Okabe and Blunt [35] used the technique to model pore space on a micrometer scale. MPS enables simulating complex interconnected pore structures by directly inferring the patterns from training images and furthermore allows modelling heterogeneity. This allows simulation of complex interconnected structures [36, 37] and the upscaling of complex pore systems to scales relevant to flow property considerations.
In this study, the CFD were for the first time applied to representative continental carbonate samples, taken from different reservoir analogues, in order to study micro- to mesoscale fluid flow in lithofacies-specific pore networks. A major benefit to these models is that the samples are differentiated based on sedimentological interpretation. Furthermore, it was verified that scan resolution and also porosity and tortuosity influence the simulated permeability. In heterogeneous carbonates, pores sizes range over several orders of magnitude and pores at different scale influence the interconnectivity. It is thus necessary to study the pore network at different scales. Not only voxel and pore network based stochastic reconstruction methods [38, 39] but also simulated annealing and Markov Chain Monte Carlo methods [40] have been proposed to merge data obtained from multiscale imaging. Here, a MPS scale-independent workflow was proposed which integrates small- and large-scale pore patterns, thus generating large volume porosity models at high resolution. MPS contributes to the comprehension of reservoir behavior in rocks by solving the common reservoir upscaling problem. To verify whether generated pore patterns were accurate, a comparison between simulated and measured permeabilities on matching artificial and original samples was made. The latter assured that MPS can in future studies be used to optimize digital pore networks and the accuracy of simulated permeabilities.
2. Research Material
Continental carbonate samples from outcrops in the Ballık area (Turkey) and Süttő and Budakalász (Hungary) were selected for analysis. The continental carbonates are of Quaternary age. The dataset covered the four dominant facies types present in the quarries, that is, the subhorizontal, reed, cascade, and waterfall facies (Figure 1), as defined by Claes et al. [1]. All of these facies types have pores at least varying from nanometer to centimeter scale and could have relevant contributions to fluid flow in this material. For a sedimentological background of the samples (out of the scope of this study), the reader is referred to earlier published literature [1, 14, 41, 42].
(a)
(b)
(c)
(d)
The facies types were usually characterized by different meso- and macropore networks, here defined as pores that have sizes of 1–100 μm and >100 μm, respectively. Pseudofenestral and interpeloidal pores, aligned along the layering, dominate the subhorizontal (Figure 1(a)) and cascade (Figure 1(b)) facies. In the latter, interlayer and shelter pores were also present between shrub crusts [1]. The reed facies (Figure 1(c)) was named after its characteristic phyto-moldic porosity. In the waterfall facies (Figure 1(d)), hanging plants on steep slopes became encrusted and generated highly porous zones in which well-connected phyto-moldic pores and high interstitial pore spaces formed framework porosity. It has to be noted that pore sizes of vugs, caverns, and shelters in these continental carbonates can exceed the sample size of classical plugs or cores.
3. Methodology
3.1. Core Analysis
Petrophysical measurements in this study were conducted on plugs with 2.5 cm and 3.4 cm diameter. These cylindrical plugs were taken in the core laboratory with a Hilti water-cooled diamond coring drill. After plugging, the samples were cut and polished to make the ends flat and parallel. In the first step of the research, samples were sent to Panterra Geoconsultants (Leiderdorp, Netherlands) for porosity-permeability analyses. The effective porosity in the plugs is measured by means of helium expansion porosimetry. Gas permeability, in this case with nitrogen gas (N2), is measured in a steady-state permeameter.
3.2. Computer Tomography (CT)
CT imagery is frequently applied in recent geomaterials research and industrial applications [43–46]. In this study, CT was used to simulate petrophysical properties in porous continental carbonate media. An inherent characteristic of CT is the relationship between resolution and sample size. The resolution of the scan () is a combination of the pixel size () of the detector, the magnification of the object (), and the size of the X-ray focal spot (). Equation (1) gives the relationship between these parameters:where SSD is the source-detector distance and SOD is the source-object distance.
Different and new CT scanners generate scans at different and ever improving scan resolutions, which can influence the estimated porosity. For example, the HECTOR scanner developed at UGCT [47] has a 240 kV setup which allows the scanning of large core samples with a diameter of 10 cm at a 28 μm resolution. This resolution is normally only reachable on smaller plug samples using standard micro-CT scanners. HECTOR CT scans were processed in Octopus [48]. In addition to the HECTOR CT, samples in this study were also scanned with a Phoenix Nanotom S instrument (GE Measurement and Control Solutions, Wunstorf, Germany), equipped with a 180 kV/15 W high-performance nanofocus X-ray tube and a 2304 × 2304 pixel Hamamatsu detector [21]. Radiographs were reconstructed with the volume processing software Phoenix datosx (GE Measurement and Control Solutions, Wunstorf, Germany) and images with an isotropic voxel size of 2, 4, 12, or 16 μm were exported. Slices were segmented in Matlab. Volume of Interest (VOI) selections and further analysis were done in Avizo Fire.
3.3. Representative Elementary Volume
The Representative Elementary Volume (REV) is a crucial concept when evaluating petrophysical properties of reservoir rocks. The REV or unit cell is the smallest volume over which a porosity measurement can be made which yields a value representative of the whole. In this study, porosity and the pore network are used as starting point to perform permeability simulations. Such simulations could only be correct if the REV for porosity was reached. When the REV for the porosity and pore types present is not reached, the pore network will yield variable porosity and thus unrepresentative simulated permeabilities. The REV of 10 samples, with a scan resolution of 16 μm, was determined according to the statistical REV calculation method described in Claes [49], which follows the approach of Bear [23]. The latter author used the chi-square criterion () as a measure of porosity fluctuation inside a selected sample volume. The smaller this value was, the closer the correspondence with the volume of the REV was. For each sample, the size of the REV was calculated 100 times, according to the procedure of Claes [49]. QQ plots of these 100 REV simulations and value calculations of the chi-square goodness of fit test are used to verify whether the REV simulations are log-normally distributed. 95% confidence intervals for the REV were subsequently calculated. This method allowed determining objectively whether the characteristic meso- and macropores in the selected samples were representative.
3.4. Multiple-Point Geostatistics (MPS)
The Single Normal Equation Simulation (SNESim) Algorithm. Traditional geostatistical simulation algorithms can be subdivided into two groups based on their principal method: pixel based group and object based group. The former simulates one pixel at a time, while the latter fits an object or pattern onto the simulation grid in one step. One of the important differences between both methods is the ease with which they can handle conditioning data. Pixel based methods are easily adapted to incorporate hard conditioning data. In contrast, object based methods faithfully fit the large-scale patterns but are difficult to condition to local data if that data is abundant. The multiple-point concept proposed by Journel [50] and Guardiano and Srivastava [51] combines the strengths of both methods. The main difference between MPS and traditional geostatistical methods is the use of a training image (TI). This TI is a representation of the expected geometry and spatial distribution of the objects present in the actual field and does not need to be a real image of the field. The simulation is performed pixelwise with the conditional probabilities respecting the conditional proportion of the TI. The TI must at least be the size of the REV to capture pore geometries and connections.
Figure 2 provides a schematic overview of the different elements used in a MPS approach, as well as its strengths for simulating complex geological structures. Figure 2(a) represents the original image, in which three facies types can be recognized. In Figure 2(b), the principal components of the original image are retained and a training image is created. Details concerning facies 2 (grey), with the lowest occurrence probability in Figure 2(a), are not retained in the training image and hence are not present in the simulated images (Figures 2(d) and 2(e)). Figure 2(c) shows the conditioning data (e.g., well data) used in the simulations. Figures 2(d) and 2(e) show the result of MPS simulations using the SNESim algorithm and the classical Sequential Indicator Simulation (SISim) algorithm, respectively. In Figure 2(d), the connectivity of facies 3 (black) is better preserved compared to the results in Figure 2(e). This observation indicates the advantage of a multiple-point approach in contrast to two-point correlations used in SISim.
(a)
(b)
(c)
(d)
(e)
In this study, the SNESim algorithm implemented by Strebelle [52] was used. Using this algorithm, the TI was scanned once and all conditional proportions were stored in a search tree data structure. In the next step, these proportions were used to create the simulated values. At each simulation node , the search template was used to retrieve the conditional data event , which is defined aswhere is a filled-in nodal value.
SNESIM Algorithm(1)Define search template and construct search tree specific to template .(2)Relocate hard data to the nearest simulation grid node and freeze them during simulation.(3)Define a random path visiting all locations to be simulated and for each location do the following:(a)Find conditioning data event defined by template τJ.(b)Retrieve conditional probability distribution function (CPDF) from .(c)Draw a simulated value from the conditional distribution and add it to the dataset.
Several parameters play an important role in the obtained simulation results. They were varied in order to improve the results of the simulations. Liu [53] and Meerschman et al. [54] provide an extensive analysis of the influence of the different parameters on the simulations. The most important parameters for upscaling using MPS are described below.
In order to capture the large-scale structural information in MPS simulations, the original algorithm was adapted to allow the use of multigrids. Tran [55] introduced this technique to create a large-scale template with a reasonably small number of nodes. It was important to keep the number of nodes limited because otherwise the search tree would become too large, resulting in an exponential increase of calculation time. In a first step, nodes on the coarsest grid were simulated using the rescaled template. Subsequently, the nodes on the second coarsest grid were simulated and so on. The relationship between “” different grids was expressed in (3). Figure 3 shows an example in which three multiple grids are used ():In MPS, the TI is used to provide detailed information on which the simulation is based. This information is derived from the patterns present in the TI. However, often the marginal distributions of the different categories in the TI are not equal to those of the real sample, which needs to be simulated; that is, the percentages of occurrence of the different categories diverge between the TI and the desired simulation. Strebelle and Journel [32] implemented a method to adapt this discrepancy of the marginal category distribution. However, if the marginal probability of the TI and the desired marginal probability are not close enough to each other, the algorithm would not reproduce the desired proportions of each category present in the simulated images. The servo corrector, introduced by Strebelle [52], bends the running simulated marginal probability further towards the target proportions by introducing a parameter λ. The larger λ is, the stronger the impact of the applied correction is (see the following equation):where is the proportion which was calculated based on the original sample and all previously simulated nodes. It is, however, important to note that and the MPS cannot completely be decoupled; hence, the target proportions should not be too different from the TI [53].
MPS in Upscaling. Figure 4 provides a schematic outline of the workflow. Figure 4(a) is a detailed image (TI) obtained using μCT. Hence, this dataset, with an isotropic voxel resolution of 4 μm, holds detailed information about micrometer scale porosity, present in a specific facies type. Figure 4(b) originated from a medical CT dataset (200 by 200 by 500 μm3) of the same sample. Medical CT usually incorporates a much larger sample volume, as inferable from the difference in scale in Figures 4(a) and 4(b). The medical CT datasets provided information about the pore network on a larger spatial scale than the μCT and was used as conditioning data in the simulations. Thus, the TI (2D slice in Figure 4(c)), derived from the μCT scans, is representative of the nonresolved matrix of the medical CT (2D slice in Figure 4(d)). The TI porosity patterns were used to simulate micrometer scale matrix porosity, which is below the resolution of the medical CT. This approach allowed retaining the connectivity of the pore network over a larger spatial scale, resulting in more reliable pore network reconstructions (Figure 4(e)) as well as more accurate permeability simulations, at least if samples obeyed the REV criterion. In order to determine the conditioning data, the assumption that the centers of the larger pores stayed the same and were recognizable in all used datasets was made, that is, in TI, condition data, and simulated results, regardless of the resolution. Centers of the largest pores were obtained by applying the workflow described below.
The dataset with the largest voxel size was segmented (Figures 4(b) and 5(a)) and the distance map of the pore facies was calculated (Figure 5(b)). This resulted in a dataset in which the center of the pore had the highest value. Subsequently, the position of the pore center was calculated by using a regional maximum algorithm. The approach not only resulted in the position of the pore centers but also preserved information about the pore sizes. In order to retain more information about the larger pores in the dataset, a threshold distance was introduced, which was calculated as σ times the maximum distance.
(a)
(b)
(c)
(d)
Only distance map values higher than σ times the maximum distance were replaced by an index representative for the property under investigation. The effect of the newly introduced parameter σ in the simulation workflow can be interpreted as follows: a small σ value will assign the property index to a larger amount of pore pixels around the pore center (Figure 5(c)), while a large σ value preserves more details about the center of the pores itself (Figure 5(d)).
By combining the patterns of the TI and the conditioning data, computer-generated rock samples were retrieved. Hence, this allowed investigating several other petrophysical parameters such as total and effective porosity, tortuosity, and permeability. Training images used in this upscaling workflow do not meet the ideal REV criterion that was previously mentioned. This option was willingly chosen because of computational cost, which is linearly related to the cube of the TI edge length. In the SNESim algorithm, all patterns are stored in the RAM of the computer and additional different patterns increase the CPU requirements. Therefore, an arbitrary size of 2503 voxels has been chosen for training images. This corresponds to a size of 1 mm3 at a resolution of 4 μm. These TI volumes have been subjectively chosen in zones that are relevant for the simulated properties. The difference with the concept of an REV is that the volume of the REV should theoretically be able to select a random zone in the sample which should always be representative. This has now been changed to an operator controlled step, which should not influence the ultimate results.
3.5. Lattice Boltzmann Method (LBM) Simulations
Stationary, pressure-driven fluid flow through porous media was simulated in Palabos by imposing a constant pressure gradient between the inlet and outlet of the pore network. Instead of solving Navier-Stokes equations, LBM solves the discrete Boltzmann equation by simulating flow for Newtonian fluids with collision models. In the applied single-phase fluid simulations under laminar flow conditions, the permeability was described by Darcy’s law (see the following equation):where is the permeability; is the fluid viscosity; is the flow velocity; and is the pressure difference over the length of the pore network.
The network geometry was stored in an input or geometry file and was based on a stack of binary images (Figure 6(a)). In the input file, value zero (blue) was assigned to every fluid voxel. Value two (red) described internal rock voxels, neighboring only other rock voxels, and value one (green) was assigned to boundary voxels that touched both fluid and rock voxels (Figure 6(b)). In the geometry file, fluid pathways had to be sufficiently resolved to avoid nonhydrodynamic effects [25, 29]. This was ensured by stepwise refining of the binary image stacks. First, pore objects with an equivalent diameter shorter than two times the length of a voxel side were deleted. In a second step, only the effective porosity, that is, all the interconnected pore space at the scanned resolution, was selected.
(a)
(b)
In this study, the standard Bhatnagar-Gross-Krook (BGK) collision operator was applied, together with the D3Q19 lattice. Flow direction during the simulations was limited to a finite vector set, representing the particle travel directions. In case of the D3Q19 lattice, there are 18 discrete lattice velocities for a fluid particle at rest (Figure 7). Palabos comes with several lattice models, including the D3Q15, D3Q19, and D3Q27 lattices. The D3Q15 lattice was not used here because of the less accurate results and numerical instabilities at high Reynold’s numbers [56, 57]. The D3Q27 lattice is more complex, when compared to the D3Q19 lattice. It provides eight additional particle transport directions, namely, to the corners of the cubic lattice nodes. Because of this, the accuracy of the model will increase, but so will the computational needs per iteration step. Both models were tested on the continental carbonates and yielded similar results, in agreement with earlier publications [56, 57]. Hence, continuing working with the D3Q19 lattice to limit the computational cost of the permeability simulations was decided.
The initial fluid velocity in the simulation was set to zero. The simulated flow accelerated during the iteration steps under influence of the fixed pressure gradient between the inlet and outlet of the pore network. Fluid velocity vectors with terminal points that coincided with the fluid/rock interface underwent a bounce back with no slip boundary conditions. This means that particle displacement from a fluid node to a solid surface (Figure 8(a)) resulted in a bounce back vector with the same initial point and direction but the opposite sense (Figure 8(b)). The initial vectors (A, B, and C) and bounce back vectors (D, E, and F) neutralized each other and resulted in a zero velocity for the displacement towards the solid surface. The standard deviation of the average energy was calculated while running simulations in Palabos. As also mentioned by Degruyter et al. [25], steady state was reached when the standard deviation of the average energy fell below a given threshold value. In this study, the threshold value was set to 10−3, with a maximum of 20000 iteration steps. Similar to Degruyter et al. [25], it was checked whether the permeability stayed constant when the applied pressure gradient was varied over several orders of magnitude. This ensured that laminar flow was established during the simulations [58, 59].
(a)
(b)
During the LBM simulations in Palabos, dimensionless lattice units are used to describe fluid flow properties, which are easily converted to any kind of physical system [57, 60]. In this study, the dimensionless lattice units were converted to Darcy units by multiplying with the square of the effective length of a voxel side. Knowing that the permeability in Darcy’s law is proportional to the ratio between fluid flow rate and the applied pressure gradient () between the inlet and outlet () of the pore network, the permeability () was described as follows:where is the effective length of a voxel side; is the lattice viscosity; and is the fluid flow rate.
Permeability was simulated along the length axis (-axis) of the plugs and spatial transformation of the binary images of nine samples allowed simulating permeability in the horizontal and directions, orthogonal to -axis. Permeability measurements were carried out on CT-derived pore networks and rock models generated in the MPS approach with cubic (500 px3) or cuboid volumes of 500 by 500 by 1000 pixels or 800 by 800 by 450 pixels, respectively.
3.6. Tortuosity
The porosity in rocks is defined by the amount of void spaces present in a sample. These pores could be connected in permeable pathways. Then, complexity (sinuosity) of the pathways for fluids through these pores is described by a property known as the tortuosity. This property provides information on the interconnectedness of the pore objects as part of a pore network and is used in transport models for porous media [61–63]. The tortuosity () is defined by the length of the flow path () and the shortest trajectory () between a defined inlet and outlet plane in the direction of the applied pressure gradient (see the following equation, [63]):The tortuosity equals 1 for a straight path through the sample () and will be infinite for a cyclic path ().
4. Results
The results below are presented in order of the subsequent steps in the simulation approach. First, the LBM method simulations were conducted on pore networks derived from natural samples and the effect of the scan resolution on CT based porosity and simulated permeability was checked. Subsequently, MPS was used to integrate high-resolution porosity patterns into lower-resolution, large-volume CT scans. LBM simulations were used to calibrate several MPS parameters and to verify whether the generated pore networks are a good representation of the real pore system.
4.1. Spatial Resolution Effect on the Recorded Pore Network
A comparison was made between datasets obtained by medical CT with 230 μm resolution (Figures 9(a) and 9(e)) and the HECTOR scans with 28 μm resolution (Figures 9(b) and 9(f)) of two different continental carbonate facies types: a reed facies and a subhorizontal sample. It is clear that much more detail is visible in the HECTOR scans. Figures 9(c) and 9(d) show the comparison of the visible porosity in both scans of both facies types. The mean porosity difference between both scans of the reed facies type was 1.5% (Figure 9(c)). The subhorizontal facies type was characterized by a larger mean porosity difference of 2.4% (Figure 9(d)). This was explained by pore volumes below the resolution of medical CT, which were more abundant in the latter facies type.
(a)
(b)
(c)
(d)
(e)
(f)
In order to assess the influence of resolution on smaller samples, a plug with 7 mm diameter was scanned at 3 different resolutions: 4 μm (Figure 10(a)), 12 μm (Figure 10(b)), and 16 μm (Figure 10(c)). Figure 10(d) shows a comparison of the calculated porosity in function of plug height. The porosity difference between the 4 μm and 16 μm scan varied between 1.4 and 2.35% along the height of the sample. These observations illustrated the importance of resolution in correctly characterizing the pore networks. Below, a workflow is proposed to combine information of different datasets which have a different resolution in order to improve the accuracy of larger samples.
(a)
(b)
(c)
(d)
4.2. Representative Elementary Volume
Pore REVs of ten μCT-scanned plugs (3.4 cm diameter) at 16 μm resolution were calculated (Table 1). For all the samples, the REV was reached. The QQ plots of 100 REV simulations per sample (not shown) had determination coefficients of at least 0.94. The value of the chi-square goodness of fit test was always well above the 0.05 hypothesis threshold. The latter tests implied that the null hypothesis was never rejected, that the calculated REVs were lognormally distributed, and that upper and lower bounds of the 95% confidence intervals for the mean size of the REV could be calculated.
The confidence intervals for the side length and volume of the cubic REV are listed in Table 1. The scanned sample volumes were all representative for the pore types and sizes that they contained. The largest REVs were found in samples originating from the reed lithofacies, which is in agreement with their heterogeneous pore networks that were observed by Claes [49]. The smallest REVs were observed in the subhorizontal and cascade facies. The pore network had a more homogeneous distribution throughout these lithofacies. The volumes over which the permeability was simulated were chosen with the size of the REVs kept in mind. Simulation cell volumes for samples from the reed facies, for example, which was the facies with the largest REV size, were always 1180 mm3, while for the subhorizontal facies, which in general had a smaller REV size, simulation cells of 512 mm3 could be used. The simulation cells were always at least over two times the REV volume.
4.3. Permeability Simulations
In total, 18 VOIs were prepared from continental carbonate plugs. The simulated permeabilities along the vertical and orthogonal horizontal axes are given together with other sample characteristics, such as lithofacies types and μCT calculated porosities (Table 2). The permeabilities obtained from LBM simulations are in good agreement with physical laboratory measurements conducted on continental carbonate samples (Figure 11). The laboratory core analyses were conducted on 5 cm3 or 10 cm3 cylindrical plugs, while volumes analyzed with LBM varied between 0.5 cm3 and 1.2 cm3.
Despite the volume difference in order of magnitude, a similar spread for porosity and permeability was observed for the core measurements and LBM simulations, with the latter plotting along the best fit curve to the core measurement data. The laboratory plug measurements include the data published in Soete [22], which indicated that the subhorizontal facies data plot mainly below 15% porosity and 100 mD. This is in agreement with the findings for both vertical and horizontal permeability in this study, where only one simulation for the subhorizontal facies yielded permeability > 100 mD and a median of 12 mD was found, while never exceeding 12.3% porosity.
For the cascade facies, phyto-rich samples were avoided and the simulations focused on samples containing shrub crust lithotypes. Dominant cascade pore types were therefore similar to those observed in the subhorizontal facies and yielded simulation results that were in the same range. The reed and waterfall facies were both characterized by a strong increase in calculated μCT porosity. In the reed facies, this increase was due to the presence of reed moldic porosity [1, 21], which formed elongated tubes that locally ran through the simulation volumes. These tubes provided high permeability pathways. The simulated permeability in the reed facies exceeded core lab results from continental carbonates from the Ballık area (Turkey), published in Soete [22], but that was expected, since highly permeable reed samples from Süttő and Budakalász in Hungary [42] were added to the dataset of this study. Waterfall pore types captured in the VOIs possessed grass and bryophyte moldic pores that together formed plug scale framework porosity.
Framework porosity, formed by plants larger than grass and bryophytes, was not included in this study. This explained the somewhat lowered permeability results for the LBM approach when compared to laboratory-measured permeability for the waterfall facies [22]. The similarities between measured and simulated permeabilities lead to the verification that the analyzed tomographic volumes qualify as permeability REVs. The μCT scans were thus at high enough spatial resolution to capture the primary permeability contributors of the pore network.
Miniplugs (7 mm diameter) were taken from some of the samples to further investigate the effect of sample size and spatial resolution on simulation results. Miniplugs were μCT-scanned at 4 μm resolution and permeability was simulated over a volume of 18.4 mm3, which was in most cases below the calculated pore REV of the samples. It was shown in Figure 10 that porosity increased for a resolution step from 16 to 4 μm. Based on the increasing porosity, higher permeabilities were expected for scans at 4 μm resolution. Despite the higher porosity and sharper μCT images, the volumes of 4 μm scans were too limited to capture the variability of the pore network. The 4 μm miniplug simulations (no REV reached), although within the order of magnitude, yielded permeabilities that overestimated and underestimated the 16 μm plug simulations, in which the REV was reached (Table 3(a)).
In addition, two miniplug VOI scans were conducted at 2 μm resolution, the highest achievable resolution for the Phoenix Nanotom S instrument in these rocks. The 2 μm VOI scans were projected back in their respective 4 μm miniplug scans. Permeability simulations for the latter VOIs were not expected to be representative, since both 2 μm and 4 μm scans were based on volumes below the REV size for porosity. The latter simulations were purely conducted to demonstrate the effect of spatial resolution. The simulations that were conducted over the exact same volumes but at different spatial scales yielded similar results (Table 3(b)). The latter demonstrated that an increase in resolution from 4 to 2 μm did not significantly improve the accuracy of the digital pore network. Based on the above observations, it was concluded that the best digital representations of the pore network were in this case achieved in 4 μm resolution CT scans but with VOIs that equal at least . Achieving 4 μm resolution while scanning large VOIs is at this point impossible and would need significant improvement of CT scanning equipment. An MPS workflow (see Section 3.6) that integrates high-resolution pore network details into large volume datasets and that captures the variability of the pore network can provide a solution and is thus proposed in this paper.
The total porosity () and connected porosity in direction () were calculated from μCT-scanned samples. For some samples, was slightly smaller, which meant that the majority of recorded pore objects in the μCT scans contributed to the connectivity. Other samples displayed large differences between total and connected porosity, for example, in simulation 14 (12.4 versus 1.4%). This indicated pore network heterogeneity and isolated porosity in the simulated flow direction at the resolution of the acquired CT images. The simulated permeability along -axis () of the samples was plotted against and (Figure 12). Porosity was in this case plotted on the vertical axis of the graph to clearly show the lowering from to . A power-law relationship was observed between and , with a determination coefficient of 0.75. This implied that good permeability estimates could be obtained from μCT porosity calculations at 16 μm resolution. However, it has to be kept in mind that pores below the 16 μm resolution could influence flow properties in continental carbonates but were not considered in these permeability estimations. Two high-permeability outliers with relatively low porosities can be observed and will be further discussed in the next section (Section 4.4).
4.4. Tortuosity
The tortuosity (τ) was calculated and plotted against permeability (Figure 13). Tortuosity controlled permeability in the studied complex pore networks. The observed power-law relationship between and was consistent with findings of studies treating volcanic rocks, where increasing tortuosity results in lower permeabilities [25, 26, 65]. Highest tortuosities were reported for the subhorizontal and cascade facies, while the waterfall and reed facies had the lowest tortuosities (Table 2). Tortuosity of the pore network helped to understand discrepancies between porosity and permeability. Two outliers from the porosity-simulated permeability regression line (sample 14 and 15 in Figure 12) were highlighted in the tortuosity-simulated permeability cross-plot (Figure 13). Both samples are part of the reed facies and consists of only a few reed molds. The limited, patchy reed mold presence did not result in high overall porosities for the samples. Despite the limited porosity, both samples are characterized by low tortuosities, indicative for the straight, high velocity flow paths that the reed molds provide through the sample.
4.5. Simulating Permeability in Orthogonal and Directions
For nine samples, including four from the subhorizontal and cascade facies and five from the reed and waterfall facies, the permeability was also simulated in the horizontal and directions, orthogonal to the vertical simulations along -axis, totaling 18 horizontal simulations. These operations were limited to nine samples because of the high computational cost of the rotation and simulation procedures. The simulations , , and for the nine samples were plotted and the influence of different facies types and related pore types on the results was checked (Figure 14). For the subhorizontal and cascade facies, the highest permeabilities were observed in the horizontal direction, often being over an order of magnitude higher than the vertical permeability of the same sample. Even in sample Al18, where an open network in the vertical direction was absent, high horizontal permeabilities were encountered. Pseudofenestral, interpeloidal, interlayer, and shelter porosities in the subhorizontal and cascade facies were described as pore types that appear aligned with the horizontal or inclined laminations. In the vertical direction, these pore types were usually discontinuous and overgrown by younger rock fabrics, which formed barriers and resulted in low vertical permeabilities (<2 mD). A pore network volume rendering of a subhorizontal facies sample, which included gastropod moldic pores (Figure 15(a)), demonstrated this low connectivity in direction. The alignment of the pores along the laminations provided better horizontal connectivity.
(a)
(b)
(c)
For the reed facies, permeability was strongly dependent on reed moldic pores. The 3D visualization of flow paths inside a reed sample (Figure 15(b)) demonstrated the heterogeneous nature of fluid flow through these porous media. High velocities were achieved within only a couple of reed tubes, which dominated the fluid flow. Depending on the orientation of the reed molds, that is, eroded (horizontal) or in-growth position (vertical), the direction of the highest permeability differed; for example, in sample Al09, which contained eroded reed, the highest permeability was achieved in the horizontal direction (Figure 14).
Grass and bryophyte framework pores in the waterfall facies formed vertically oriented pore networks. The small framework pipes yielded high permeabilities in the vertical direction, while the horizontal connectivity was dependent on the connections between individual pipes. The horizontal permeability was therefore lowered in samples from the waterfall facies, for example, samples CK31 and CA31 (Figure 14). The streamline rendering (Figure 15(c)) shows the low tortuous path in direction, with relatively straight connections between the bottom and top plane of the sample.
4.6. Permeability Simulations on the MPS Generated Rock Models
In Section 3.3, the spatial resolution effect on simulated permeability on the scanned sample set was shown. The results demonstrated that (1) high resolution scans were needed to capture micrometer scale pores which impacted flow properties and that (2) lower-resolution CT scans should be performed on sufficiently large volumes to fully capture the heterogeneity of the pore network. In order to overcome the volume versus resolution problem of CT scanners, MPS-generated rock models are used as input for permability simulations. The latter allowed assessing the applicability of the MPS workflow in generating larger-volume rock models at higher resolutions.
In order to assess the influence of the above-introduced conditioning data MPS parameter σ on the rock models, a series of simulations were performed using different σ values (Figures 16(a)–16(d)), while keeping all other parameters equal. The simulations were performed on continental carbonate samples from the reed facies. This facies type was chosen because of its typical abundant open porosity network and the characteristic elongated shape of the pore bodies. The latter resulted in permeability values of 50 mD and higher [41]. The TI had a resolution of 4 μm and the conditioning dataset had a resolution of 16 μm. Hence, the resulting simulated dataset contained the volume of the larger 16 μm CT scans, with a higher 4 μm resolution. The TI had a size of 2503 pixels or 10003 μm. The search template distance is 25 pixels. The servo corrector was kept zero and six multigrids were used. The TI volume was chosen with extra care in order to represent the typical patterns or in this case pore structures. The resulting simulated samples are shown in Figure 16. The shape of the pores corresponded well with pore shapes that were present in the TI. Rod- and blade-shaped pores made up 50% of the pores. Because permeability was measured in the direction of the aligned pore shapes, a permeability value of around 1000 mD was expected based on petrophysical measurements indicated in Figure 11.
(a)
(b)
(c)
(d)
Analyzing the shape of the pores, according to the classification of Claes et al. [66], provided an excellent tool for the quality control of the simulation. Similar pore shape occurrence percentages between the physical sample to be simulated and the MPS simulated sample were indicative of good quality simulations. The results also indicated high dependency on conditioning data in order to retain the desired larger-scale connectivity in the simulated sample. The connectivity of the pore network had a direct influence on the simulated permeability values of the simulated rock samples. The permeability values increased when more conditioning data, that is, lower σ values, were used in the simulation. Simulated permeability values were in the expected range for reed continental carbonate samples, but the most representative permeabilities were found for σ = 0.25, that is, the largest amount of conditioning data. Large amounts of conditioning data were needed to simulate complex carbonate rocks.
In order to further test the applicability of the workflow, three adjacent volumes of a sample from the subhorizontal facies were simulated using a TI (2503 pixels) with a resolution of 4 μm. Figure 17(a) depicts the conditioning data (1600 × 1600 × 800 pixels) used in a simulation of this facies type, which had a resolution of 16 μm and σ value of 0.25. Three different zones were generated using the same TI and the results (8003 pixels) were depicted in Figures 17(b), 17(c), and 17(d), respectively.
(a)
(b)
(c)
(d)
The TI was selected inside a porous layer characteristic for this facies type. The resulting subvolumes had a side length of 3.2 mm and hence depicted only one sedimentary lamina. This explained why the typical horizontal layering of these samples, with porous and less porous layers, was not clearly visible in the simulated volumes. They would only have become visible when the spatial scale of both the TI and simulation grid was increased. The resulting simulated permeability values ranged from 0 mD, in case no connected network was present, to 203 mD. Hence, the permeability values were highly dependent on the provided conditioning data, as was already demonstrated in the reed facies samples in Figure 16. Consistent with the laboratory permeability measurements, the overall values were lower compared to the reed facies (1367 mD for σ = 0.25, Figure 16(d)). The lower permeabilities in the vertical direction for the subhorizontal facies were in line with the findings in Section 3.5, where permeability was simulated in horizontal and vertical directions. The simulated rock models preserved the characteristic laminated pore structure of the subhorizontal facies, with limited connectivity in the vertical direction.
The unique datasets generated by the HECTOR scans provided an excellent case study to evaluate the potential of the proposed workflow. Because of their detailed resolution (28 μm), HECTOR scans delivered excellent TI. The medical CT scan data provided the larger-scale conditioning data and had a resolution of 230 μm in and directions and 500 μm in direction. This allowed performing petrophysical simulations on core samples. In this case, the difference in TI resolution and conditioning data was more than doubled compared to the first set of simulations. This allowed simulating artificial samples at a larger spatial scale with a resolution of 28 μm. Because of the porosity difference between the TI and the generated rock sample, a servo factor of 0.75 was used. If a lower value was used, the porosity value of the simulated rock samples became too high, resulting in visually unrealistic models and unrealistic simulated permeability values. For example, when the reed sample (Figure 18(c)) was simulated using a servo factor of 0, the calculated permeability became 22 D. This observation illustrated the importance of the servo factor, which should be taken higher when the resolution difference increases. A clear difference between the simulated porosity networks of both the reed and subhorizontal facies types was observed. The moldic reed pores were retained in the simulation. In contrast, typical sequences of porous and less porous layers in subhorizontal samples were observed in simulated samples.
The simulation-obtained permeabilities were in the expected order of magnitude compared to measured values. Moreover, the anticipated difference between the subhorizontal and reed facies type pore networks was clearly apparent. In the subhorizontal sample, no open network was present in direction ( shown in Figure 18(a)). Figures 18(b) and 18(d) show the simulation results in which the same TI was used as in Figures 18(a) and 18(c), respectively, but the conditioning data were derived from a different sample of a similar facies type. Simulations with the same TI but different CD demonstrated that the quality of the HECTOR based TI was sufficient to be used on multiple samples.
5. Discussion
CT scanning and computational fluid dynamics (CFD) proved to be excellent tools for the investigation of flow properties in complex continental carbonate pore networks. Different flow properties and paths were observed for different facies types and their characteristic pore types. Consistent with De Boever et al. [16] and Claes et al. [66], the permeability was higher for samples from the reed and waterfall facies. The subhorizontal and cascade facies were characterized by lower permeabilities and often yielded higher permeabilities horizontally in contrast to the vertical direction.
The main factors influencing the permeability in continental carbonates were connected porosity () and tortuosity (). For both variables, a power-law relationship with the permeability was established. However, the determination coefficients (0.75 and 0.66, resp.) revealed data scatter; hence, permeability predictions based on either of these variables alone were not very accurate. and counteracted each other at several occasions, like when porosity would overestimate simulated permeability of a sample and when tortuosity would underestimate it and vice versa. Therefore, more accurate permeability predictions for a given sample were obtained by using (8), which includes both and :A cross-plot of predicted versus LBM-simulated permeability (determination coefficient of 0.85) illustrates how and interact in the determination of the permeability of the samples under investigation (Figure 19). LBM-simulated permeability was slightly overestimated by the predicted permeability at low values. The predicted values tend to eliminate the most extreme values, limiting the 0.1–10000 mD simulated permeability range to 0.35–4500 mD in the estimations.
The REV measurements and CFD simulations demonstrated the effect of sample scale and spatial resolution on the permeabilities. Claes (2015) confirmed that different REV sizes can be found at different scales, corresponding to different geological length scales. Here, only pores below plug size were considered in REV analyses. For this pore scale, the porosity REV was reached in all 16 μm VOIs, over which permeability simulations were performed. Large-scale framework pores, caverns, and decimeter-sized vugs were, however, also encountered in the field but do not form a part of this research. Increasing the sample scale to include decimeter-sized vugs would drastically increase the REV size. In general, with increasing sample size, new porosity REVs would be encountered whenever new large-scale pore types would be introduced. Macroscopic observations also showed that large-scale framework, cavern, and vug porosity were mainly interconnected through the pores encountered within the studied plugs. The main pore scales governing connectivity and fluid flow at the scale of the carbonate body were therefore expected to be included in this study.
In general, the following can be concluded: the larger the CT analyzed volume and the finer the scan resolution, the more realistic the representation of the pore network. The MPS workflow proposed in this manuscript allows generating artificial rock samples which can be used to investigate the influence of rock heterogeneity on fluid flow dynamics. Moreover, this technique permits bridging the gap between different datasets with different resolutions. The workflow was used on datasets where the resolution of the conditioning data was four and eight times the resolution of the TI. In both cases, the simulations yielded realistic results.
The amount of conditioning data, determined by the σ value, and the use of the servo factor proved to be important parameters to obtain realistic simulations. Comparing measured and simulated permeabilities showed that σ = 0.25 yielded the most realistic pore networks models for the investigated continental carbonate samples. Large amounts of conditioning data were necessary to model pore networks of continental carbonate samples in a representative way. The MPS workflow was applied to upscale from miniplug (diameter 7 mm) to plug (diameter of 3.81 cm) and to large core volumes (diameter of 9.7 cm), but Zhang [67] also proved the applicability of MPS datasets on larger spatial scales such as borehole imaging and facies distribution models in reservoir modelling.
The LBM simulations for generated rock models are plotted on top of the measured relationship between porosity and permeability for continental carbonate plug samples (Figure 20). The simulated volumes were indicated in red. The simulated values follow the general trend and are therefore considered to be realistic.
6. Drawbacks and Future Perspectives
As with all modelling techniques, there are some drawbacks to this research even though efforts are made to reduce them. In this paper, the choice for fluid flow modelling techniques and upscaling approaches is based on the need for correctly representing the porous network. This is achieved by avoiding simplifications as much as possible. In first instance, this is evidenced by the chosen CFD method. The LBM approach is chosen as it offers a good trade-off between computational power and representativeness. In LBM, the fluids are represented as small packages in individual voxels rather than as individual molecules. The actual pore network is used in the simulations rather than simplifications that pore topology models typically use. In order to obtain realistic permeability simulations, representative pore networks are preferred, even though simplifications would decrease the computational cost drastically.
Secondly, the choice was made to improve resolution of spatially larger models rather than simplifying the problem through classical upscaling techniques that smoothen out small-scale variations. This choice is made based on the rationale that small-scale features, variations, and heterogeneities hold the information that can determine large-scale reservoir properties. This rationale is not only relevant in, for example, unconventional reservoirs, but also in carbon capture and storage (CCS), in which the smallest pores are used for capillary capture of fluids.
The use of multiple-point statistics in upscaling is still a relatively new approach and as such a methodology which is not yet perfected. Known issues include the determination of input parameters in the SNESim algorithm and its computational requirements. A lot of parameters, like the number of multigrids, search template size, and λ parameter, have to be adjusted in the SNESim algorithm, which introduces subjectivity into the models. To avoid subjectivity as much as possible, the workflow is here calibrated with physical laboratory experiments. Visual inspection is always applied on the results to ensure that (1) a good fit between simulated properties and physical measurements is obtained and that (2) models are geologically relevant.
Apart from visual inspection, it is observed that the σ-parameter for conditioning data is best kept constant for specific upscaling problems. A σ-parameter of 0.25 is an ideal value when the resolution for the TI and conditioning data differs with a factor of 4. A σ-parameter of 0.75 proved to be more appropriate for a factor of 8 resolution difference. This difference in the σ-parameter is logical given the difference in the amount of overlap in the samples. Another drawback is that even when applying the σ-factor, it sometimes occurs that new larger pores are unnecessarily generated. The larger pores are, however, expected to be sufficiently represented in the lower-resolution scans. The training images are chosen at 2503 voxels, which is smaller than the REVs, to limit the computational cost. The selection of these smaller training images was user-controlled in order to select highly representative zones. The quality of the TI therefore will not be compromised by the smaller size.
The final drawback of any upscaling technique that aims to increase resolution, while also increasing scale, is computational cost. For example, increasing the resolution by a factor of 4 increases the required storage space by a factor of 43. Such volumes might easily exceed the abilities of standard workstations to model fluid flow. Because of the computer power limitation, this research focused on μCT and HECTOR scans, which yield sample volumes smaller than the outcrop scale macropores that are regularly present in carbonates. Such oversized pores potentially influence large-scale fluid flow and can thus be important components in reservoir models. Even though these pores have not been treated here, they could be incorporated into models with the proposed reservoir scale upscaling purpose, as it is intrinsically scale-invariant. There is, ad hoc, no method that allows imaging of oversized pores three-dimensionally. MPS input data such as training images or conditioning data could nevertheless be created through photogrammetry or LiDAR [68, 69] or through methods that use 2D input data like field pictures to generate 3D volumes, similar to Okabe and Blunt [35]. Generating models that include macropores observed at outcrop scale, that is, vugs, framework porosity, caverns, and fractures, can significantly increase the sensitivity towards reservoir scale fluid flow. Outcrop scale and plug to core scale models, however, have to be considered separately, since incorporating all elements of a multiscale pore system would make the models extremely large and computationally expensive.
The future of MPS-driven upscaling techniques should aim to reduce subjective parameters as much as possible. Both this research and the research by Zhang [67] indicate that MPS is a potent method in upscaling of reservoir properties, but neither of these approaches provides the most solid solution. For fluid flow experiments on increasingly larger models, it might be required to use intensive parallel computing, apply modified versions of the LBM, or use pore network models with microlinks [24] in order not to nullify the increased resolution.
7. Conclusion
Transport properties are fundamental in the characterization of reservoir bodies. In this study, computed tomography of continental carbonate samples and Lattice Boltzmann Method (LBM) permeability simulations were applied to obtain 3D quantifications of the pore network and flow paths. Working with 3D datasets was an absolute necessity to understand flow in complex continental carbonates, because accurate estimates of properties which govern fluid flow in a rock, that is, connectivity, tortuosity, and so forth, could not be obtained from 2D observations. The simulated permeabilities were in good agreement with physical plug and core permeabilities. Power-law relationships were observed between permeability and connected porosity on one hand and permeability and tortuosity on the other hand. Estimates of rock permeability for a sample were strongly improved by including both tortuosity and connected porosity. Both parameters balance one another, and so their integration allows better prediction of the fluid flow through the complex porous network of continental carbonate samples.
The subhorizontal and cascade facies were associated with pseudofenestral and interpeloidal pores and yielded the lowest porosities and permeabilities. The reed and waterfall facies, with, respectively, reed moldic and framework porosity, were characterized by much better reservoir properties. The resulting porosities easily surpassed 30% and exceeded 50 D, in both measured and simulated permeabilities. Facies and pore type-dependent heterogeneity was observed. Better horizontal connectivity and permeability for the subhorizontal and cascade facies were in contrast to the generally higher vertical permeabilities that were observed for samples from the reed and waterfall facies.
The proposed multiple-point geostatistics (MPS) workflow is an additional tool to solve common problems in the upscaling of reservoir properties. An additional advantage of this scale-independent approach is its applicability on 2D as well as 3D datasets. MPS does not introduce general averaging that is commonly used in upscaling.
The workflow allowed upscaling to realistic models at different scales by starting from detailed, high-resolution information and incorporating it into larger-scale models. Moreover, using recent advances in fluid flow models, this workflow has the potential to replace expensive permeability measurements on large core samples and can serve as input for simulations of multiphase fluid flow. The results have also indicated that the TI image can be used for different samples of the same facies type. This would allow scanning samples more quickly using conventional medical CT scanners. When a more detailed description of a part of the pore network is desired, a pore model can be generated using this workflow.
The here-presented innovative and integrated LBM-MPS workflow was applied in continental carbonate rocks with complex pore networks. The methodology is directly applicable to other lithologies, comprising different pore types, pore shapes, and pore networks altogether. The lack of straightforward porosity-permeability relationships in complex carbonates highlights the necessity for a 3D approach. Studying fluid pathways in 3D provides the best possible understanding of flow through porous media and will be of crucial importance in reservoir modelling.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The authors thank the owners and managers of the respective quarries for allowing them to work in actively excavated environments, which provided insight into the 3D pore networks of the continental carbonate rocks. They would like to thank Herman Nijs for aiding in practical matters during the sample preparation processes. The authors would like to acknowledge the Hercules Foundation (Flanders) for founding the micro- and nano-CT project for the hierarchical analysis of materials. J. Soete was funded by a Ph.D. grant from “Agentschap voor Innovatie door Wetenschap en Technologie” (IWT), Flanders, Belgium. Last but not least, they would like to thank Professor Mehmet Özkul and Sándor Kele for their help and the interesting discussions during the field trips in Turkey and Hungary.