Regular Article
Consistent prediction of absolute permeability in carbonates without upscaling
1
Deanship of Research, King Fahd University of Petroleum and Minerals, 31261 Dhahran, Saudi Arabia
2
Center for Integrative Petroleum Research, College of Petroleum Engineering and Geosciences, King Fahd University of Petroleum and Minerals, 31261 Dhahran, Saudi Arabia
* Corresponding author: mrkhodja@kfupm.edu.sa
Received:
27
October
2019
Accepted:
14
April
2020
We describe a study focused on the absolute permeability of reservoir carbonate rocks from the Middle East and involving comparison of experimental data and numerical estimates obtained by combining digital-rock and Lattice-Boltzmann Methods (LBM). The question of the “representativeness” of the site at which the simulation is performed is addressed as follows. First, a low-resolution, CT X-ray scan of the core plug is performed to identify regions of large porosity (millimeter-sized vugs, etc.). These regions are then avoided to postselect smaller sites (site volume ~ 1 mm3) which are to be scanned at higher resolutions (voxel size < dominant pore-throat size of the core plug). A “representativeness” criterion based on an empirically-inspired “representativeness” measure (R-measure) is used to eliminate those sites for which R > b, where b is an upper bound (typically, b = 1). Essentially, the measure estimates how well the postselected sites capture the experimental porosity and the dominant pore-throat size of the core plug. This leads to a small set of sites for which the simulations are both computationally manageable and yield a reasonable estimate of the permeability: the experimental and predicted values differ by a factor of about 3 on average, which is a particularly significant result given the challenging heterogeneous pore space of carbonate samples. We believe the suggested methodology to be an adequate and practical way to circumvent upscaling.
© M.R. Khodja et al., published by IFP Energies nouvelles, 2020
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
It is estimated that the Middle East hosts more than half of the world’s proven conventional oil reserves and more than one third of the world’s proven gas reserves. Approximately, 60% of the world’s oil reserves and 90% of its gas reserves reside in carbonate reservoirs. However, production from these reservoirs is rather challenging. Carbonates are characterized by a wide variety of petrophysical and flow properties, oftentimes over very small scales. Understanding how these multimodal small-scale heterogeneities affect the performance of the reservoir is a fundamental and challenging problem that requires the integration of a multitude of multidisciplinary studies. In particular, carbonates are known to exhibit complex pore networks whose structure has a considerable effect on fluid flow and other properties [1–3].
Classically, the petrophysical properties of the reservoir rock are obtained by performing a variety of experimental tests on actual core samples. However, the experimental tests are: (1) costly, (2) time-consuming, and (3) destructive. Digital Rock Physics (DRP) is quickly emerging as a powerful and attractive methodology that complements conventional core analysis. Digital rock analysis is carried out on digitized models of the core [4–8]. Essentially, the digitized model is obtained, for instance, by performing high-resolution (i.e., micrometer to nanometer), three-dimensional (3D), Computed X-ray Tomography (CT) scans of the core. These models are analyzed and turned into realistic representations of the internal structure of the rock. Afterwards, numerical experiments are performed on the 3D models to determine the needed properties. Notably, these experiments can be performed repeatedly under different scenarios since they only involve virtual representations of the rock samples. DRP has been used to study reservoir rock mineralogy [9, 10], elastic properties [11, 12], flow and transport phenomena [13–17], etc.
Recently, several groups have looked into many aspects of the permeability prediction problem. For instance, an attempt has been made at quantifying the permeability anisotropy in carbonates [18], an inverse approach has been adopted in which the measured porosity is used as a starting point to estimate the permeability [19], multiscale microtomography has been used to try and obtain a reasonable representation of the pore-network structure [20–22], further attempts have been made at quantifying the effect of image resolution on simulated permeability [23, 24]. Because of the opaque nature of most porous media, permeability is usually obtained by measuring flow speed under a certain pressure gradient in experiments. Thanks to the development in imaging technologies over the past few decades, it is now possible to visualize the interior pore structures within the rock samples. The void fraction is related to the porosity of the rock, which can be characterized experimentally by mercury porosimetry. Using the measured rock porosity, it is possible to determine a threshold in the grayscale, to discriminate between pores and rock, and use this threshold to filter the attenuation values of the X-ray microtomography image to generate a binary (0 for rock site, 1 for pore site) representation of the rock and the pore space [25]. The binarization process defines all voxels with monochromatic value less than a threshold to be pore space or void, and the rest to be solid and consequently the segmented porosity monotonously decreases with the threshold value [26]. In a combination with the breakthrough in computational methods, we can make direct flow simulations through the pore space of real rock samples to calculate their permeabilities. One of such computational methods is the Lattice-Boltzmann Methods (LBM) [27, 28], which has the advantages of simplicity and flexibility in dealing with complicated geometry over the traditional methods of Computational Fluid Dynamics (CFD), e.g. the Finite-Volume Method (FVM) in which the flow is simulated by solving directly the discretized Navier–Stokes equations. Additionally, the LBM utilizes the nearest neighbor information in an explicit updating algorithm and so it is ideally suited for massively parallel computation.
Here, we describe a study focused on the absolute permeability of reservoir carbonate rocks from the Middle East and involving comparison of experimental data and numerical estimates obtained by combining digital-rock and LBM. This study is comprised of four parts. In the first part, the permeability of three core samples with varied porosity, pore-network structure, and modality is determined experimentally using conventional core analysis methods. Nuclear Magnetic Resonance (NMR) and Mercury-Injection Capillary-Pressure (MICP) measurements are performed to determine the pore-body and pore-throat size distributions, respectively. Also, the mineral composition for the core plugs is determined. In the second part, the core plugs are digitized using low-resolution X-ray scans (voxel volume ~ [30 μm]3) and their pore structure is extracted numerically and compared with NMR-measured pore-size distributions for validation. The purpose of these scans is to identify the regions with large porosity (millimeter-sized vugs, etc.). In the third part, higher-resolution scans (voxel size < dominant pore-throat size of the core sample) are performed on small sites (site volume ~ 1 mm3) semi-randomly selected and located outside the regions with large porosity. In the fourth part, numerical modeling of pore-scale flow is carried out on the sites that are best “representative” of the core sample. “Representativeness” is based on a simple, empirically-inspired measure. Essentially, the measure estimates how well the postselected sites capture the experimental porosity and the dominant pore-throat size of the core sample. The hydrodynamic simulations in this study are performed with our implementation of the LBM.
2 Experimental procedure
Three carbonate core plugs (labeled 1, 2, and 3), obtained from a reservoir in the Middle East, were used in this study. Cylindrical end trims (25 mm in diameter and 6 mm in length) were cut off from the parent core plugs (Figs. 1 and 2). These end trims were used in NMR analysis and in micro-CT X-ray imaging. The remaining 25-mm-by-25-mm cylindrical plugs were used to acquire Mercury-Injection Capillary-Pressure (MICP) measurements. Helium porosity, ϕexp, and helium permeability, κexp, were measured on the core plugs sent for MICP, i.e. without the trim ends, using a Coretest AP-608 Automated Permeameter-Porosimeter. For the end trims, the NMR porosity was determined. The measurements are tabulated in Table 1. We shall assume that the permeabilities of the trim ends are consistent with those of their parent core plugs and the plugs sent for MICP analysis.
Fig. 1 Schematic illustration of the workflow for sample preparation and performed studies. |
Fig. 2 Image showing an example of the end trims used for micro-CT, NMR, and QEMSCAN analyses. The unit on the scale is centimeter. |
Experimental data.
NMR T2-relaxation experiments were performed using the GeoSpec 2-75 (Oxford Instruments) rock-core analyzer with the Green Imaging Technologies (GIT) software. All the experiments were conducted at fixed field strength with a resonating frequency of 2.2 MHz at ambient temperature. The samples were saturated with deionized water for 24 h at a pressure of 1500 psi. The standard Carr-Purcell-Meiboom-Gill (CPMG) pulse sequence was applied to obtain T2-relaxation decay curves using the Echo Time (TE) value constant (TE = 0.053 ms). The anticipated T2max was set to 2000 ms and the Signal-to-Noise Ratio (SNR) for each experiment was retained at 50.
MICP analysis was conducted to determine the pore-throat radius distribution of the core samples (after cutting the trim ends). The instrument used in our analysis is the AutoPore IV 9520 Mercury Porosimeter (Micromeritics). Prior to MICP analysis, the core plugs are dried, weighed, evacuated, and then immersed in mercury at a preset filling pressure of 1.6 psi. About 100 pre-set pressure steps on a log-based scale are programmed and computer-driven over the intrusion pressure range 1.6–60 000 psi, which yields equally spaced, incremental mercury intrusion data points when pore aperture size distribution data is plotted on a semi-log graph. The dominant pore-throat radius corresponds to the peak of the curve.
Scans of the dry (i.e., drained) end trims were conducted using the Versa XRM-500 X-ray micro-CT scanner. The scan recipes involved source voltages in the range 140–160 kV, a power of 10 W, exposure times in the range 1–4 s to obtain the desired flux count, and 1601 projections. First, a low-resolution scan (pixel size ~ 20 μm is performed). The segmentation of the reconstructed 3D models consisted of a combination of two methods: the gradient marker-based watershed methods and the top-hat method. The isolated pores were removed using the axis connectivity tool in PerGeos Software (ThermoFisher Scientific). After segmentation, the regions with well-defined macro-porosity are identified (Fig. 3). To assess the quality of the obtained segmentation, the NMR pore-body radius distribution is compared with the pore-body radius distribution obtained from the micro-CT images. The volume fraction of the segmented pores was then computed and compared with the porosity obtained using other lab measurements. Once the dominant pore-throat radius rmicp is determined from MICP data, one proceeds to the second stage in which higher-resolution scans are performed on small (~1 mm) distinct sites located in the unresolved part of the end trims, while avoiding the regions with millimeter-sized vugs (Fig. 4). These are the postselected sites at which the permeability simulations were carried out. Care was exercised to ensure that fluid flow in the simulations was in the direction that permeability was measured. The connected pore space was used as an input to the pore-network extraction module and a pore-network model of the sample showing the pores and throats was obtained (Fig. 5). The pore-throat distribution was later analyzed to identify the dominant pore-throat size which is one of the factors that affects permeability.
Detailed information on the mineral composition for the plugs was obtained using Quantitative Evaluation of Minerals by Scanning Electron Microscopy (QEMSCAN) on the QEMSCAN 650F (ThermoFisher Scientific) with the iDiscover software. The information is shown in Figure 6.
Fig. 3 Micro-CT image of Sample 2 at 29 μm. Regions with millimeter-sized vugs (inside the yellow circle) were avoided when selecting sites for high-resolution scans. |
Fig. 4 (a) Filtered high-resolution micro-CT image of Sample 2, Site B, at voxel resolution of 2.15 μm. The filter used is a non-local means filter with a search window of 21 and a pixel neighborhood of 5. (b) Binary version of (a) obtained using a combination of gradient marker-based watershed and top-hat segmentations. The pores are shown in black. |
Fig. 5 (a) 3D rendering of Site B, Sample 2, after segmentation. (b) Pore network model of Site B, Sample 2. |
Fig. 6 QEMSCAN analyses of the three carbonate samples: (a) Sample 1, (b) Sample 2, and (c) Sample 3. |
3 Lattice-Boltzmann Method
The ordinary LBM with a single relaxation time is adopted in our pore-scale simulations of single-phase single-component flows to compute the absolute permeabilities of real rock samples. The fundamental idea of the LBM is to construct a kinetic model simplified from the Boltzmann equation at mesoscale that incorporate the essential physical processes to ensure the macroscopic quantities satisfy the desired macroscopic governing equation (i.e., Navier–Stokes equation in the current study). The nonslip boundary condition is imposed by using the simplest bounce-back boundary condition, which is robust and efficient in problems with complicated geometry, although it causes small numerical errors that can be corrected by shifting the solid boundary by half a voxel, as implemented in the mid-plane bounce-back scheme [29]. Thus, the errors are at the same order of magnitude as those inevitably associated with the segmentation of the digital rock.
In the LBM simulations, the grids or lattices are uniformly distributed inside the computational domain and all computational quantities are defined at those discrete points. We use Δx and Δt to denote the lattice distance and time step, respectively. At each point, we have several fixed lattice velocities eα, which are either static with α = 0 or connect the current lattice at x to its neighboring lattices at x + Δteα, where α varies from 1 to Q − 1 and Q is the total number of velocities. The magnitude of eα depends on c = Δx/Δt. The three-dimensional, 19-velocity (D3Q19) model [27] is used in our three-dimensional simulations to determine eα and the weighting factor ωα. The explicit evolution algorithm of the density distribution function fα(x, t) is,
(1)wherein is determined by the kinematic viscosity ν and the equilibrium distribution function is defined as follows,
where and ueq is equal to the flow velocity for flow problems without external force. The transient pressure at each lattice is calculated as p = c2 ρ/3.
After convergence, the volumetric average velocity 〉u〈Ω is used to compute three components of the permeability tensor as follows when the flow is driven in the y direction for example:
where L y is the size of computational domain along the y axis, μ is the dynamic viscosity, pin − pout is the pressure difference between inlet and outlet, uj is the velocity at the void grid j, Ntotal is the total lattice number, and ei is the unit vector along the i-th axis. In our current simulations, the flows are always driven by a pressure difference (i.e., one percent difference between the inlet and outlet). Consequently, only the permeability component along the driving direction is used in the following discussions and will be simply denoted by κ (Figs. 3–5).
4 Numerical results and validations
4.1 Convergence study
In the single-phase single-component simulations using the ordinary LBM, the number of steps required to reach steady state mostly depends on the lattice (or voxel) number particularly along the driven direction and usually also depends on the tortuosity of the rock geometry (namely rock-dependent) and dimensionless relaxation time τ, but is certainly independent of Δx (e.g., voxel size as in the current study), Δt and the pressures used initially and at the boundaries. As usually adopted in the literature, the simulation is deemed as converged when the relative change over a certain number of time steps (e.g., 1000 time steps) is smaller than a preset tolerance (e.g., 10−5). But, the rock samples as well as the specific convergence criteria implemented by different groups are usually different, which leads to different convergence speeds and computational time. Boek and Venturoli used only about 5000–8000 time steps for the simulations of both 1283 lattices and 5123 lattices with τ = 0.6 [25]. Sukop et al. performed the simulations with 1283 and 2563 lattices and more than 25 000 time steps are used to achieve equilibration [30]. Hosa et al. ran the simulation with 1003 lattices for 100 000 time steps to achieve convergence [26]. Li et al. [31] performed the convergence study for both sandstone and tight-sandstone samples using τ = 1 and found that the simulation of sandstone sample with 160 × 160 × 30 voxels takes 35 000 time steps to converge while 162 000 time steps are used for the simulation of tight sandstone with 200 × 200 × 50 voxels.
To show the variation of convergence speed in different simulations, we first study the influence of τ on the convergence speed in simulating a geometry with 5003 voxels from site A of Sample 3 (see Tab. 2). Figures 7a and 7b, show that the simulation with τ = 0.7 takes about 4000 time steps to converge while about 10 000 time steps are needed by the simulation with τ = 1.0, which clearly shows the dependence of convergence speed on the parameter τ. We note that the simulation with smaller τ converges to higher volumetric average velocity because the kinematic viscosity decreases with τ, which does not matter as permeability is independent of viscosity. Then, we repeat the simulation with τ = 0.7 for a different geometry with 5003 voxels as well but from site A of Sample 1 (see Tab. 2) to study the dependence of convergence speed on the specific sample. Figure 7c shows that the number of time steps required for convergence increases to about 400 000 from about 4000 in Figure 7a. Thus, the convergence speed is problem-dependent and it is hard to have a generally valid criterion to automatically stop the simulation in time. Instead, we monitor the evolution curve as in Figure 7 and end the simulation once a plateau is reached.
Fig. 7 Convergence processes of volumetric velocities in the LBM simulations. (a) Convergence for the first geometry (Sample 3, Site A, in Table 2) with τ = 0.7. (b) Convergence for the first geometry (Sample 3, Site A, in Tab. 2) with τ = 1.0. (c) Convergence for the second geometry (Sample 1, Site A, in Tab. 2) with τ = 0.7. |
Numerical results.
4.2 Permeability study with experimental validation
We note that the computed permeability depends on the selection of τ in the ordinary LBM simulations using the BGK model [25, 26, 32–34] which is actually due to the Klinkenberg slippage effect when the Kn number is not close to zero as clearly elucidated in [35]. The higher the Kn number, the larger the ratio of the computed permeability and the intrinsic permeability will be. In ordinary LBM simulations with the Two-Dimensional, 9-Velocity (D2Q9) model [27] or the D3Q19 model, we have,
where N pore is the lattice number used to discretize the dominant pore size [36]. In [35], it is shown that the undesired slippage effect becomes negligible and the computed permeability becomes almost unchanged as desired at τ < 0.7 when Npore is more than 4. In our scans, 〈Npore〉 = 4.2. Additionally, the computed permeability will start dropping due to undesired inertial effect with the increase of the Reynolds number (Re) if we keep on decreasing τ to the limit of 0.5 [35]. Thus, we select τ = 0.7 to make the slippage effect negligible and remove the nonlinear terms of ueq that appears in equation (2) to completely eliminate the inertial effect [35].
As discussed in [35], the multi-relaxation model may provide a parameter range that is wider than that of the single-relaxation model but the basic mechanism remains the same, i.e. in simulations of intrinsic permeability, one should choose appropriate τ and spatial grid resolution values to ensure that Kn is small and thus reduces the rarefaction effect. However, like the single-relaxation model used here, the multi-relaxation model cannot make Kn precisely zero, as the viscosity is nonzero. Hence, we adopt the single-relaxation LBM model for simplicity. Additionally, the single-relaxation model is much more efficient than the multi-relaxation model and is also widely used. The undesired Kn effect in computing the intrinsic permeability can be removed theoretically by using τ = 0.5, however, this choice leads to a vanishing viscosity and thus makes the permeability definition of equation (3) invalid. In fact, the permeability variation with τ due to the Kn effect becomes negligible as (τ − 0.5) ∝ Kn approaches zero. In Figure 8, we provide the permeability variation with (τ − 0.5) for a representative case (i.e., 5003 voxels at Site A of Sample 3, as used in the previous convergence study). The figure shows that the intrinsic permeability extrapolated at τ − 0.5 = 0 is about 2.30 × 10−12, which differs by about 6% from the value of 2.44 × 10−12 obtained by our choice of τ − 0.5 = 0.2. Consequently, we shall use τ = 0.7 in the following applications.
Fig. 8 LBM permeability variation with (τ − 0.5) for Site A of Sample 3. |
Fig. 9 Pore-throat radius distributions. The solid colored curves correspond to the original PNM data. The dashed colored curves correspond to the smoothed PNM data. The smoothing was achieved by using a Gaussian kernel with a radius of 5 pixels. |
How to postselect the small site where the high-resolution scan and the permeability simulation are to be performed? The small size of the site enables us to reduce the computational cost of the simulation to a manageable level. Moreover, the assumption is that if the site is “representative” of the core sample (see Sect. 4), its permeability computed in this manner ought to be of the same order of magnitude as the measured value. The question is how to assess “representativeness”. We define the following heuristic “representativeness measure” (R-measure, for short),
where,
and where ϕ, and r are the numerical porosity and pore-throat radius, respectively, and ϕref is a reference porosity. Typically, ϕref = ϕexp. However, because the core trim ends on which the X-ray scans and NMR analyses were performed in this study were cutoff before ϕexp and κexp were measured, we also need to consider the values of δϕ that correspond to ϕref = ϕnmr. For the sake of comparison, we also define,
where κlbm is the numerical permeability. This particular choice of the R-measure was inspired by the well-established empirical correlation ([37] and references therein),
where m and n are constants. We claim that the best sites to perform the high-resolution scan and the permeability simulation in a core sample are those for which,
where b is an upper bound: for instance, b = 1. (An error-propagation study that shows the stability of the R-measure is presented in the Appendix.) This “representativeness criterion” is not as ideal as requiring that R be the minimum of the R-measure values of all the postselected sites because permeability is not a function solely of porosity and pore-throat radius and variability is inherent to equation (9). In Figure 10, we plot |δκ| versus R for the six postselected sites.
Fig. 10 Plot of the “representativeness measures”, defined by equation (5), using the permeability results of our LBM code. The disks correspond to δϕ values in which ϕref = ϕexp, while the squares correspond to δϕ values in which ϕref = ϕnmr. |
The numerical results are summarized in Table 2. In Table 3, we present a comparison of experimental and computed permeabilities (the experimental values of the permeability are those of Tab. 1 while the computed values are those of Tab. 2). The positive values correspond to cases when the numerical code overestimates the permeability, while the negative values correspond to cases when it underestimates the permeability. The |δκ| values vary from 9.8% to 1314% with an average of 476.7% that corresponds to the mean value 〈κlbm/κexp〉 = 3.1. In Figure 9, we plot the pore-throat radius distributions obtained with MICP, and Pore-Network Modeling (PNM) tool in PerGeos software. The solid black curves represent MICP data, the colored spiky curves represent the PNM data (bin = 0.1 μm). To increase the SNR and simplify the task of determining the dominant pore-throat radius, we applied to the PNM data a Gaussian filter with a radius of 5 pixels (i.e., with standard deviation σ = 5/2 and Full-Width-at-Half-Maximum ). The dashed colored curves represent the results of this filtering. In Table 4, we present a rudimentary comparison with two recently published methods that employ upscaling.
Comparison of experimental and computed permeabilities.
Rudimentary comparison with recent literature.
5 Conclusion
We have described a study focused on the absolute permeability of reservoir carbonate rocks from the Middle East and involving comparison of experimental data and numerical estimates obtained by combining digital-rock and LBM. The permeability simulations are performed with our implementation of the LBM. To assess the “representativeness” of the site at which the simulation is performed, we have defined an empirically inspired measure that estimates how well the site captures the experimental porosity and the dominant pore-throat size of the core sample. This leads to sites for which the simulations are both computationally manageable and yield a reasonable estimate of the permeability, given the simplicity of the method and that upscaling is avoided. The typical workflow is as follows:
-
Determine experimentally the porosity and the dominant pore-throat radius of the core plug.
-
X-ray scan the core plug at low-resolution (voxel volume ~ [30 μm]3) to identify regions with millimeter-sized vugs.
-
Preselect sites outside the vuggy regions and X-ray scan them at high resolution (voxel volume < [pore-throat radius]3).
-
Generate pore-network models of scanned sites and determine the porosities and dominant pore-throat radii of the preselected sites.
-
Use the representativeness criterion to select the “best” site(s) (R ≤ 1) for the permeability simulation(s).
-
Run the permeability simulation(s).
Acknowledgments
The authors gratefully acknowledge partial funding from King Fahd University of Petroleum and Minerals.
Appendix
Error propagation and stability of the R-measure
Assuming that the experimental and digitization errors are random and uncorrelated, one obtains,
wherein,
and,
In the above expressions, Δr and Δϕ are the errors in the numerical pore-throat size and porosity, respectively, due to digitization. In particular, Δr≈ pixel size, while .
Since [40], we also have,
For our MICP system, ΔPmicp/Pmicp ≤ 0.5%. Using arithmetic means for the remaining quantities, yields,
This shows that the representativeness criterion is reasonably robust against measurement and digitization errors.
References
- Moctezuma-Berthier A., Vizika O., Thovert J.F., Adler P.M. (2004) One- and two-phase permeabilities of vugular porous media, Transp. Porous Media 56, 225–244. [Google Scholar]
- Norouzi Apourvari S., Arns C.H. (2015) Image-based relative permeability upscaling from the pore scale. Advanced , Water Resour. 95, 16–175. [Google Scholar]
- Gjetvaj F., Russian A., Gouze P., Dentz M. (2015) Dual control of flow field heterogeneity and immobile porosity on non-fickian transport in berea sandstone, Water Resour. Res. 51, 8273–8293. [Google Scholar]
- Dvorkin J., Armbruster M., Baldwin C., Fang Q. (2008) The future of digital rock physics: Computational methods vs. lab testing, First Break 26, 63–68. [CrossRef] [Google Scholar]
- Kalam M.Z. (2012) New technologies in the oil and gas industry, in: Gomes J.S. (ed.), Digital rock physics for fast and accurate special core analysis in carbonates, Chapter 9, INTECH, Croatia, pp. 201–226. [Google Scholar]
- Andrä H., Combaret N., Dvorkin J., Glatt E., Han J., Kabel M., Keehm Y., Krzikalla F., Lee M., Madonna C., Marsh M., Mukerji T., Saenger E., Sain R., Saxena N., Ricker S., Wiegmann A., Zhan X. (2013) Digital rock physics benchmarks – part I: Imaging and segmentation, Comput. Geosci. 50, 25–32. [Google Scholar]
- Andrä H., Combaret N., Dvorkin J., Glatt E., Han J., Kabel M., Keehm Y., Krzikalla F., Lee M., Madonna C., Marsh M., Mukerji T., Saenger E., Sain R., Saxena N., Ricker S., Wiegmann A., Zhan X. (2013) Digital rock physics benchmarks – part II: Computing effective properties, Comput. Geosci. 50, 33–43. [Google Scholar]
- Cnudde V., Boone M.N. (2013) High-resolution x-ray computed tomography in geosciences: A review of the current technology and applications, Earth-Sci. Rev. 123, 1–17. [Google Scholar]
- Golab A., Knackstedt M.A., Averdunk H., Senden T., Butcher A.R., Jaime P. (2010) 3D porosity and mineralogy characterization in tight gas sandstones, Lead. Edge 29, 1476–1483. [Google Scholar]
- Dernaika M., Uddin Y.N., Koronfol S., Al Jallad O., Sinclair G. (2015) Multi-scale rock analysis for improved rock typing in complex carbonates, in: International Symposium of the Society of Core Analysts, St. John’s Newfoundland and Labrador, Canada, 16–21 August, 2015 [Google Scholar]
- Arns C., Knackstedt M., Pinczewski W., Garboczi E. (2002) Computation of linear elastic properties from microtomographic images: Methodology and agreement between theory and experiment, Geophysics 67, 1396–1405. [CrossRef] [Google Scholar]
- Saenger E., Enzmann F., Keehm Y., Steeb H. (2011) Digital rock physics: Effect of fluid viscosity on effective elastic properties, J. Appl. Geophys. 24, 236–241. [CrossRef] [Google Scholar]
- Blunt M. (2001) Flow in porous media – pore-network models and multiphase flow, Curr. Opin. Colloid Interface Sci. 6, 197–207. [Google Scholar]
- Arns C., Knackstedt M., Pinczewski M., Lindquist W. (2001) Accurate estimation of transport properties from microtomographic images, Geophys. Res. Lett. 28, 17, 3361–3644. [Google Scholar]
- Keehm Y., Mukerji T., Nur A. (2004) Permeability prediction from thin sections: 3D reconstruction and lattice-Boltzmann flow simulation, Geophys. Res. Lett. 31, 4, L04606. [Google Scholar]
- Ovaysi S., Piri M. (2010) Direct pore-level modeling of incompressible fluid flow in porous media, J. Comput. Phys. 229, 19, 7456–7476. [Google Scholar]
- Aryana S., Kovscek A.R. (2013) Nonequilibrium effects and multiphase flow in porous media, Transp. Porous Media 97, 373–394. [Google Scholar]
- Sun H., Vega S., Tao G. (2017) Analysis of heterogeneity and permeability anisotropy in carbonate rock samples using digital rock physics, J. Pet. Sci. Eng. 156, 419–429. [Google Scholar]
- Saenger E.H., Vialle S., Lebedev M., Uribe D., Osorno M., Duda M., Steeb H. (2016) Digital carbonate rock physics, Solid Earth 7, 1185–1197. [CrossRef] [Google Scholar]
- Jouini M.S., Vega S., Al-Ratrout A. (2014) Numerical estimation of carbonate rock properties using multiscale images, Geophys. Prospect. 63, 2, 405–421. [Google Scholar]
- Rahimov K., AlSumaiti A.M., Jouini M.S. (2016) Quantitative analysis of absolute permeability and porosity in carbonate rocks using digital rock physics, in: 22nd Formation Evaluation Symposium of Japan, 29–30 September, Chiba, Japan. Society of Petrophysicists and Well-Log Analysts. SPWLAJFES-2016-J. [Google Scholar]
- Sun H., Vega S., Tao G. (2016) Determination of transport properties in carbonate rock sample using multi-scale ct images, in: 78th EAGE Conference and Exhibition 2016, 30 May–2 June, Vienna, Austria. [Google Scholar]
- Latief F.D.E., Fauzi Z., Irayani U., Dougherty G. (2017) The effect of X-ray micro computed tomography image resolution on flow properties of porous rocks, J. Microsc. 266, 1, 69–88. [CrossRef] [PubMed] [Google Scholar]
- Bazaikin Y., Gurevich B., Iglauer S., Khachkova T., Kolyukhin D., Lebedev M., Lisitsa V., Reshetova G. (2017) Effect of ct image size and resolution on the accuracy of rock property estimates, J. Geophys. Res.: Solid Earth 122, 5, 3635–3647. [CrossRef] [Google Scholar]
- Boek E.S., Venturoli M. (2010) Lattice-Boltzmann studies of fluid flow in porous media with realistic rock geometries, Comput. Math. Appl. 59, 2305–2314. [Google Scholar]
- Hosa A., Curtis A., Wood R. (2016) Calibrating lattice Boltzmann flow simulations and estimating uncertainty in the permeability of complex porous media, Adv. Water Res. 94, 60–74. [CrossRef] [Google Scholar]
- Qian Y.H., d’Humieres D., Lallemand P. (1992) Lattice bgk models for Navier–Stokes equation, Europhys. Lett. 17, 479–484. [Google Scholar]
- Chen H., Chen S.Y., Matthaeus W.H. (1992) Recovery of the Navier–Stokes equations using a lattice-gas Boltzmann method, Phys. Rev. A 45, 5339–5342. [Google Scholar]
- Sukop M.C., Thorne D.T. (2006) Lattice Boltzmann modeling: An introduction for geoscientists and engineers, Springer, Berlin. [CrossRef] [Google Scholar]
- Sukop M.C., Huang H.B., Lin C.L., Deo M.D., Oh K., Miller J.D. (2008) Distribution of multiphase fluids in porous media: Comparison between lattice Boltzmann modeling and micro-x-ray tomography, Phys. Rev. E 77, 026710. [Google Scholar]
- Li R.R., Yang Y.S., Pan J.X., Pereira G.G., Taylor J.A., Clennell B., Zou C.N. (2014) Lattice Boltzmann modeling of permeability in porous materials with partially percolating voxels, Phys. Rev. E 90, 033301. [Google Scholar]
- Ferreol B., Rothman D.H. (1995) Lattice-Boltzmann simulations of flow through Fontainebleau sandstone, Transp. Porous Media 20, 3–20. [Google Scholar]
- Pan C., Luo L.S., Miller C.T. (2006) An evaluation of lattice boltzmann schemes for porous medium flow simulation, Comput. Fluids 35, 898–909. [Google Scholar]
- Prestininzi P., Montessori A., Rocca M.L., Succi S. (2016) Re- assessing the single relaxation time lattice boltzmann method for the simulation of darcy’s flows, Int. J. Mod. Phys. C 27, 1650037. [Google Scholar]
- Li J., Ho M.T., Wu L., Zhang Y.-H. (2018) On the unintentional rarefaction effect in LBM modeling of intrinsic permeability, Adv. Geo-Energy Res. 2, 404–409. doi: 10.26804/ager.2018.04.05. [CrossRef] [Google Scholar]
- Zhang Y.H., Qin R.S., Emerson D.R. (2005) Lattice Boltzmann simulation of rarefied gas flows in microchannels, Phys. Rev. E 71, 047702. [Google Scholar]
- Rezaee M.R., Jafari A., Kazemzadeh E. (2006) Relationships between permeability, porosity and pore throat size in carbonate rocks using regression analysis and neural networks, J. Geophys. Eng. 3, 370–376. [CrossRef] [Google Scholar]
- Sun H., Belhaj H., Tao G., Vega S., Liu L. (2019) Rock properties evaluations for carbonate reservoir characterization with multi- scale digital rock images, J. Pet. Sci. Eng. 173, 175, 654–664. [Google Scholar]
- Islam A., Faisal T.F., Chevalier S., Jouini M.S., Sassi M. (2019) Multi-scale experimental and numerical simulation workflow of absolute permeability in heterogeneous carbonates, J. Pet. Sci. Eng. 173, 173, 326–338. [Google Scholar]
- Wardlaw N.C., Taylor R.P. (1976) Mercury capillary pressure curves and the intepretation of pore structure and capillary behaviour in reservoir rocks, Bull. Can. Pet. Geol. 24, 2, 225–262. [Google Scholar]
All Tables
All Figures
Fig. 1 Schematic illustration of the workflow for sample preparation and performed studies. |
|
In the text |
Fig. 2 Image showing an example of the end trims used for micro-CT, NMR, and QEMSCAN analyses. The unit on the scale is centimeter. |
|
In the text |
Fig. 3 Micro-CT image of Sample 2 at 29 μm. Regions with millimeter-sized vugs (inside the yellow circle) were avoided when selecting sites for high-resolution scans. |
|
In the text |
Fig. 4 (a) Filtered high-resolution micro-CT image of Sample 2, Site B, at voxel resolution of 2.15 μm. The filter used is a non-local means filter with a search window of 21 and a pixel neighborhood of 5. (b) Binary version of (a) obtained using a combination of gradient marker-based watershed and top-hat segmentations. The pores are shown in black. |
|
In the text |
Fig. 5 (a) 3D rendering of Site B, Sample 2, after segmentation. (b) Pore network model of Site B, Sample 2. |
|
In the text |
Fig. 6 QEMSCAN analyses of the three carbonate samples: (a) Sample 1, (b) Sample 2, and (c) Sample 3. |
|
In the text |
Fig. 7 Convergence processes of volumetric velocities in the LBM simulations. (a) Convergence for the first geometry (Sample 3, Site A, in Table 2) with τ = 0.7. (b) Convergence for the first geometry (Sample 3, Site A, in Tab. 2) with τ = 1.0. (c) Convergence for the second geometry (Sample 1, Site A, in Tab. 2) with τ = 0.7. |
|
In the text |
Fig. 8 LBM permeability variation with (τ − 0.5) for Site A of Sample 3. |
|
In the text |
Fig. 9 Pore-throat radius distributions. The solid colored curves correspond to the original PNM data. The dashed colored curves correspond to the smoothed PNM data. The smoothing was achieved by using a Gaussian kernel with a radius of 5 pixels. |
|
In the text |
Fig. 10 Plot of the “representativeness measures”, defined by equation (5), using the permeability results of our LBM code. The disks correspond to δϕ values in which ϕref = ϕexp, while the squares correspond to δϕ values in which ϕref = ϕnmr. |
|
In the text |