© M.R. Khodja et al., published by IFP Energies nouvelles, 2020

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]³) 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³) 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.

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 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 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.

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,

$$f_{\alpha }(\mathbf{x}+\Delta t {\mathbf{e}}_{\alpha }, t+\Delta t)=f_{\alpha }(\mathbf{x},t)+\frac{f_{\alpha }^{\mathrm{eq}}(\mathbf{x},t)-f_{\alpha }(\mathbf{x},t)}{\tau },$$(1)wherein $\tau =3\nu \Delta t \Delta x^{-2}+0.5$ is determined by the kinematic viscosity ν and the equilibrium distribution function is defined as follows,

$$f_{\alpha }^{\mathrm{eq}}(\mathbf{x},t)=\rho {\omega }_{\alpha }\left[1+\frac{3}{c^2}{\mathbf{e}}_{\alpha }\cdot {\mathbf{u}}^{\mathrm{eq}}+\frac{9}{2c^4}{\left({\mathbf{e}}_{\alpha }\cdot {\mathbf{u}}^{\mathrm{eq}}\right)}^2-\frac{3}{2c^2}{\mathbf{u}}^{\mathrm{eq}}\cdot {\mathbf{u}}^{\mathrm{eq}}\right],$$(2)

where $\rho ={\sum }_{\alpha =0}^{Q-1}\mathrm{}{f}_{\alpha }$ and u^eq is equal to the flow velocity $\mathbf{u}=\frac{1}{\rho }{\sum }_{\alpha =0}^{Q-1}\mathrm{}{\mathbf{e}}_{\alpha }{f}_{\alpha }$ for flow problems without external force. The transient pressure at each lattice is calculated as p = c² ρ/3.

After convergence, the volumetric average velocity 〉u〈_Ω is used to compute three components of the permeability tensor $\overline{\stackrel{̅}{\kappa }}$ as follows when the flow is driven in the y direction for example:

$${\kappa }_{\mathrm{iy}}=\frac{\mu L_y}{p_{\mathrm{in}}-p_{\mathrm{out}}}{\langle \mathbf{u}\rangle }_{\mathrm{\Omega }}\cdot {\mathbf{e}}_i=\frac{\mu L_y}{p_{\mathrm{in}}-p_{\mathrm{out}}}\frac{1}{N_{\mathrm{total}}}\sum _{j\in \mathrm{void}}\mathrm{}{\mathbf{u}}_j\cdot {\mathbf{e}}_i;\hspace{1em}i=x,y,z,$$(3)

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 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⁻⁵). 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³ lattices and 512³ lattices with τ = 0.6 [25]. Sukop et al. performed the simulations with 128³ and 256³ lattices and more than 25 000 time steps are used to achieve equilibration [30]. Hosa et al. ran the simulation with 100³ 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 500³ 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³ 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.

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,

$$\mathrm{Kn}={\left(\frac{\pi }{6}\right)}^{1/2} \frac{\tau -0.5}{N_{\mathrm{pore}}},$$(4)

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 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., 500³ 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⁻¹², which differs by about 6% from the value of 2.44 × 10⁻¹² obtained by our choice of τ − 0.5 = 0.2. Consequently, we shall use τ = 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),

$$R\equiv {\left(\delta r^2+\delta {\varphi }^2\right)}^{1/2},$$(5)

where,

$$\delta \varphi \equiv \frac{\varphi -{\varphi }_{\mathrm{ref}}}{{\varphi }_{\mathrm{ref}}},$$(6)

$$\delta r\equiv \frac{r-r_{\mathrm{micp}}}{r_{\mathrm{micp}}},$$(7)

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,

$$\delta \kappa \equiv \frac{{\kappa }_{\mathrm{lbm}}-{\kappa }_{\exp }}{{\kappa }_{\exp }},$$(8)

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),

$$\kappa \propto {\varphi }^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,

$$R<b,$$(10)

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.

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 $\left(\mathrm{FWHM}\right)=2\sqrt{2\ln 2} \sigma \approx 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.

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:

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 μm]³) to identify regions with millimeter-sized vugs.

3. Preselect sites outside the vuggy regions and X-ray scan them at high resolution (voxel volume < [pore-throat radius]³).

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).

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,

$$\frac{\Delta R}{R}=\frac{\sqrt{{\left(\frac{\Delta \left(\delta \varphi \right)}{\delta \varphi }\right)}^2+{\left(\frac{\Delta \left(\delta r\right)}{\delta r}\right)}^2{\left(\frac{\delta r}{\delta \varphi }\right)}^4}}{1+{\left(\frac{\delta r}{\delta \varphi }\right)}^2},$$(A.1)

wherein,

$$\frac{\Delta \left(\delta \varphi \right)}{\delta \varphi }=\frac{\Delta {\varphi }_{\mathrm{ref}}}{{\varphi }_{\mathrm{ref}}}\sqrt{1+\frac{1+{\left(\frac{\Delta \varphi }{\Delta {\varphi }_{\mathrm{ref}}}\right)}^2}{{\left(1-\frac{\varphi }{{\varphi }_{\mathrm{ref}}}\right)}^2}}$$(A.2)

and,

$$\frac{\Delta \left(\delta r\right)}{\delta r}=\frac{\Delta r_{\mathrm{micp}}}{r_{\mathrm{micp}}}\sqrt{1+\frac{1+{\left(\frac{\Delta r}{\Delta r_{\mathrm{micp}}}\right)}^2}{{\left(1-\frac{r}{r_{\mathrm{micp}}}\right)}^2}}.$$(A.3)

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 $\Delta \varphi \approx \Delta {\varphi }_{\mathrm{ref}}\approx 10^{-3}$.

Since $r_{\mathrm{micp}}\propto P_{\mathrm{micp}}^{-1}$ [40], we also have,

$$\frac{\Delta r_{\mathrm{micp}}}{r_{\mathrm{micp}}}=\frac{\Delta P_{\mathrm{micp}}}{P_{\mathrm{micp}}}.$$(A.4)

For our MICP system, ΔP_micp/P_micp ≤ 0.5%. Using arithmetic means for the remaining quantities, yields,

$$\frac{\Delta R}{R}\approx 0.15.$$(A.5)

This shows that the representativeness criterion is reasonably robust against measurement and digitization errors.

References

All Tables

All Figures