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 digitalrock and LatticeBoltzmann Methods (LBM). The question of the “representativeness” of the site at which the simulation is performed is addressed as follows. First, a lowresolution, CT Xray scan of the core plug is performed to identify regions of large porosity (millimetersized vugs, etc.). These regions are then avoided to postselect smaller sites (site volume ~ 1 mm^{3}) which are to be scanned at higher resolutions (voxel size < dominant porethroat size of the core plug). A “representativeness” criterion based on an empiricallyinspired “representativeness” measure (Rmeasure) 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 porethroat 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 smallscale 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) timeconsuming, 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 highresolution (i.e., micrometer to nanometer), threedimensional (3D), Computed Xray 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 porenetwork 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 Xray 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 LatticeBoltzmann 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 FiniteVolume 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 digitalrock and LBM. This study is comprised of four parts. In the first part, the permeability of three core samples with varied porosity, porenetwork structure, and modality is determined experimentally using conventional core analysis methods. Nuclear Magnetic Resonance (NMR) and MercuryInjection CapillaryPressure (MICP) measurements are performed to determine the porebody and porethroat size distributions, respectively. Also, the mineral composition for the core plugs is determined. In the second part, the core plugs are digitized using lowresolution Xray scans (voxel volume ~ [30 μm]^{3}) and their pore structure is extracted numerically and compared with NMRmeasured poresize distributions for validation. The purpose of these scans is to identify the regions with large porosity (millimetersized vugs, etc.). In the third part, higherresolution scans (voxel size < dominant porethroat size of the core sample) are performed on small sites (site volume ~ 1 mm^{3}) semirandomly selected and located outside the regions with large porosity. In the fourth part, numerical modeling of porescale flow is carried out on the sites that are best “representative” of the core sample. “Representativeness” is based on a simple, empiricallyinspired measure. Essentially, the measure estimates how well the postselected sites capture the experimental porosity and the dominant porethroat 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 microCT Xray imaging. The remaining 25mmby25mm cylindrical plugs were used to acquire MercuryInjection CapillaryPressure (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 AP608 Automated PermeameterPorosimeter. 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 microCT, NMR, and QEMSCAN analyses. The unit on the scale is centimeter. 
Experimental data.
NMR T2relaxation experiments were performed using the GeoSpec 275 (Oxford Instruments) rockcore 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 CarrPurcellMeiboomGill (CPMG) pulse sequence was applied to obtain T2relaxation decay curves using the Echo Time (TE) value constant (TE = 0.053 ms). The anticipated T2_{max} was set to 2000 ms and the SignaltoNoise Ratio (SNR) for each experiment was retained at 50.
MICP analysis was conducted to determine the porethroat 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 preset pressure steps on a logbased scale are programmed and computerdriven 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 semilog graph. The dominant porethroat radius corresponds to the peak of the curve.
Scans of the dry (i.e., drained) end trims were conducted using the Versa XRM500 Xray microCT 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 lowresolution scan (pixel size ~ 20 μm is performed). The segmentation of the reconstructed 3D models consisted of a combination of two methods: the gradient markerbased watershed methods and the tophat method. The isolated pores were removed using the axis connectivity tool in PerGeos Software (ThermoFisher Scientific). After segmentation, the regions with welldefined macroporosity are identified (Fig. 3). To assess the quality of the obtained segmentation, the NMR porebody radius distribution is compared with the porebody radius distribution obtained from the microCT images. The volume fraction of the segmented pores was then computed and compared with the porosity obtained using other lab measurements. Once the dominant porethroat radius r_{micp} is determined from MICP data, one proceeds to the second stage in which higherresolution scans are performed on small (~1 mm) distinct sites located in the unresolved part of the end trims, while avoiding the regions with millimetersized 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 porenetwork extraction module and a porenetwork model of the sample showing the pores and throats was obtained (Fig. 5). The porethroat distribution was later analyzed to identify the dominant porethroat 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 MicroCT image of Sample 2 at 29 μm. Regions with millimetersized vugs (inside the yellow circle) were avoided when selecting sites for highresolution scans. 
Fig. 4 (a) Filtered highresolution microCT image of Sample 2, Site B, at voxel resolution of 2.15 μm. The filter used is a nonlocal 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 markerbased watershed and tophat 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 LatticeBoltzmann Method
The ordinary LBM with a single relaxation time is adopted in our porescale simulations of singlephase singlecomponent 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 bounceback 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 midplane bounceback 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 threedimensional, 19velocity (D3Q19) model [27] is used in our threedimensional 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 u^{eq} is equal to the flow velocity for flow problems without external force. The transient pressure at each lattice is calculated as p = c^{2} ρ/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, p_{in} − p_{out} is the pressure difference between inlet and outlet, u_{j} is the velocity at the void grid j, N_{total} is the total lattice number, and e_{i} is the unit vector along the ith 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 singlephase singlecomponent 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 rockdependent) 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 128^{3} lattices and 512^{3} lattices with τ = 0.6 [25]. Sukop et al. performed the simulations with 128^{3} and 256^{3} lattices and more than 25 000 time steps are used to achieve equilibration [30]. Hosa et al. ran the simulation with 100^{3} lattices for 100 000 time steps to achieve convergence [26]. Li et al. [31] performed the convergence study for both sandstone and tightsandstone 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 500^{3} 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 500^{3} 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 problemdependent 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 TwoDimensional, 9Velocity (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 N_{pore} is more than 4. In our scans, 〈N_{pore}〉 = 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 u^{eq} that appears in equation (2) to completely eliminate the inertial effect [35].
As discussed in [35], the multirelaxation model may provide a parameter range that is wider than that of the singlerelaxation 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 singlerelaxation model used here, the multirelaxation model cannot make Kn precisely zero, as the viscosity is nonzero. Hence, we adopt the singlerelaxation LBM model for simplicity. Additionally, the singlerelaxation model is much more efficient than the multirelaxation 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., 500^{3} 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 Porethroat 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 highresolution 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” (Rmeasure, for short),
where,
and where ϕ, and r are the numerical porosity and porethroat radius, respectively, and ϕ_{ref} is a reference porosity. Typically, ϕ_{ref} = ϕ_{exp}. However, because the core trim ends on which the Xray 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 Rmeasure was inspired by the wellestablished empirical correlation ([37] and references therein),
where m and n are constants. We claim that the best sites to perform the highresolution scan and the permeability simulation in a core sample are those for which,
where b is an upper bound: for instance, b = 1. (An errorpropagation study that shows the stability of the Rmeasure is presented in the Appendix.) This “representativeness criterion” is not as ideal as requiring that R be the minimum of the Rmeasure values of all the postselected sites because permeability is not a function solely of porosity and porethroat 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 porethroat radius distributions obtained with MICP, and PoreNetwork 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 porethroat radius, we applied to the PNM data a Gaussian filter with a radius of 5 pixels (i.e., with standard deviation σ = 5/2 and FullWidthatHalfMaximum ). 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 digitalrock 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 porethroat 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 porethroat radius of the core plug.

Xray scan the core plug at lowresolution (voxel volume ~ [30 μm]^{3}) to identify regions with millimetersized vugs.

Preselect sites outside the vuggy regions and Xray scan them at high resolution (voxel volume < [porethroat radius]^{3}).

Generate porenetwork models of scanned sites and determine the porosities and dominant porethroat 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 Rmeasure
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 porethroat size and porosity, respectively, due to digitization. In particular, Δr≈ pixel size, while .
Since [40], we also have,
For our MICP system, ΔP_{micp}/P_{micp} ≤ 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
 MoctezumaBerthier A., Vizika O., Thovert J.F., Adler P.M. (2004) One and twophase permeabilities of vugular porous media, Transp. Porous Media 56, 225–244. [Google Scholar]
 Norouzi Apourvari S., Arns C.H. (2015) Imagebased 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 nonfickian 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) Highresolution xray computed tomography in geosciences: A review of the current technology and applications, EarthSci. Rev. 123, 1–17. [CrossRef] [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) Multiscale 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 – porenetwork 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 latticeBoltzmann flow simulation, Geophys. Res. Lett. 31, 4, L04606. [Google Scholar]
 Ovaysi S., Piri M. (2010) Direct porelevel 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., AlRatrout 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 WellLog Analysts. SPWLAJFES2016J. [Google Scholar]
 Sun H., Vega S., Tao G. (2016) Determination of transport properties in carbonate rock sample using multiscale 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 Xray 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) LatticeBoltzmann 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 latticegas 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 microxray 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) LatticeBoltzmann 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. GeoEnergy 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) Multiscale 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 microCT, NMR, and QEMSCAN analyses. The unit on the scale is centimeter. 

In the text 
Fig. 3 MicroCT image of Sample 2 at 29 μm. Regions with millimetersized vugs (inside the yellow circle) were avoided when selecting sites for highresolution scans. 

In the text 
Fig. 4 (a) Filtered highresolution microCT image of Sample 2, Site B, at voxel resolution of 2.15 μm. The filter used is a nonlocal 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 markerbased watershed and tophat 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 Porethroat 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 