Consistent prediction of absolute permeability in carbonates without upscaling

. 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 mm 3 ) 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 signi ﬁ cant result given the challenging heterogeneous pore space of carbonate samples. We believe the sug-gested methodology to be an adequate and practical way to circumvent upscaling.


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][2][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][5][6][7][8]. Essentially, the digitized model is obtained, for instance, by performing high-resolution (i.e., micrometer to nanometer), threedimensional (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][14][15][16][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][21][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 porethroat 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 lm] 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 mm 3 ) 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.

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-25mm cylindrical plugs were used to acquire Mercury-Injection Capillary-Pressure (MICP) measurements. Helium porosity, / exp , and helium permeability, j 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.
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 T2 max 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 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 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 lm 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 r micp 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.

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 Dx and Dt to denote the lattice distance and time step, respectively. At each point, we have several fixed lattice velocities e a , which are either static with a = 0 or connect the current lattice at x to its neighboring lattices at x þ Áte a , where a varies from 1 to Q À 1 and Q is the / exp and j exp are respectively, the porosity and permeability of the core samples after cutting off the end trims, / nmr is the NMR porosity of the end trims, r micp is the MICP dominant pore-throat radius, and q g is the grain density. Sample 1 is bimodal (see Fig. 9). total number of velocities. The magnitude of e a depends on c = Dx/Dt. The three-dimensional, 19-velocity (D3Q19) model [27] is used in our three-dimensional simulations to determine e a and the weighting factor x a . The explicit evolution algorithm of the density distribution function f a ðx; tÞ is, f a ðx þ Át e a ; t þ ÁtÞ ¼ f a ðx; tÞ þ f eq a ðx; tÞ À f a ðx; tÞ s ; wherein s ¼ 3m Át Áx À2 þ 0:5 is determined by the kinematic viscosity m and the equilibrium distribution function is defined as follows, f eq a ðx; tÞ ¼ qx a 1 þ 3 c 2 e a Á u eq þ 9 2c 4 e a Á u eq ð Þ 2 À 3 2c 2 u eq Á u eq ! ; where q ¼ P QÀ1 a¼0 f a and u eq is equal to the flow velocity a¼0 e a f a for flow problems without external force. The transient pressure at each lattice is calculated as After convergence, the volumetric average velocity u h i X is used to compute three components of the permeability tensor j 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, l 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 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 j (Figs. 3-5).

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 s, but is certainly independent of Dx (e.g., voxel size as in the current study), Dt 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  To show the variation of convergence speed in different simulations, we first study the influence of s 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 s = 0.7 takes about 4000 time steps to converge while about 10 000 time steps are needed by the simulation with s = 1.0, which clearly shows the dependence of convergence speed on the parameter s. We note that the simulation with smaller s converges to higher volumetric average velocity because the kinematic viscosity decreases with s, which does not matter as permeability is independent of viscosity. Then, we repeat the simulation with s = 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 / and r are the numerical porosity and dominant porethroat radius of the postselected site, respectively. j lbm is the permeability prediction of our LBM code.  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.

Permeability study with experimental validation
We note that the computed permeability depends on the selection of s in the ordinary LBM simulations using the BGK model [25,26,[32][33][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 s < 0.7 when N pore is more than 4. In our scans, hN pore i ¼ 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 s to the limit of 0.5 [35]. Thus, we select s = 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 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 s 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 s = 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 s due to the Kn effect becomes negligible as ðs À 0:5Þ / Kn approaches zero. In Figure 8, we provide the permeability variation with (s À 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 s À 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 s À 0.5 = 0.2. Consequently, we shall use s = 0.7 in the following applications.
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 porethroat 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 j exp were measured, we also need to consider the values of d/ that correspond to / ref ¼ / nmr . For the sake of comparison, we also define, where j lbm is the numerical permeability. This particular choice of the R-measure was inspired by the well-established empirical correlation ( [37] and references therein), j / / m r n ; ð9Þ 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 errorpropagation 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 dj j j versus R for the six postselected sites.
The numerical results are summarized in Table 2. In Table 3, we present a comparison of experimental and  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 lm). 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 r = 5/2 and Full-Width-at-Half-Maximum ðFWHMÞ ¼ 2 ffiffiffiffiffiffiffiffiffiffi ffi 2 ln 2 p r % 5:9). 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.

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: 1. Determine experimentally the porosity and the dominant pore-throat radius of the core plug. 2. X-ray scan the core plug at low-resolution (voxel volume~ [30 lm] 3 ) to identify regions with millimetersized vugs. 3. Preselect sites outside the vuggy regions and X-ray scan them at high resolution (voxel volume < [porethroat radius] 3 ).
4. Generate pore-network models of scanned sites and determine the porosities and dominant pore-throat radii of the preselected sites. 5. Use the representativeness criterion to select the "best" site(s) (R 1) for the permeability simulation(s). 6. Run the permeability simulation(s).  [38] Yes 6À16~120~400 3 À900 3 [39] Yes 5À1400~300~700 3 À1000 3 T refers to the estimated time it would take to arrive at the results using the resources currently available in our laboratory. V is the approximate volume of the sites on which the permeability simulation was run.