Numerical investigation of oil-water drainage and imbibition in digitized sandstones

– We simulated the primary drainage and ﬁrst imbibition processes in a numerical porous-plate experiment using a lattice-Boltzmann based pore scale simulator. The pore space representations of the Berea, Fontainebleau and Bentheimer sandstones were obtained by high resolution micro-CT and subsequent segmentation. The pore body and throat distributions were identiﬁed using an open-map based method. The rock surface was assumed completely water wet. The capillary pressure as function of water saturation and the oil and water distributions were investigated for diﬀerent geometries represented by the diﬀerent sandstone samples. The calculated residual oil saturation after imbibition compares well with recently published high quality laboratory data. The dependence of the residual oil saturation on the initial oil saturation was calculated and the results are consistent with published experimental data and correlations. Finally the drainage and imbibition processes were calculated for several Fontainebleau samples with diﬀerent porosity. It was found that the residual oil saturation decreases with porosity and the variation agrees well with published experimental data.


Introduction
Oil-water drainage and imbibition processes in rocks are of fundamental interest to oil industry.For long time, the applied capillary pressure at different water saturations was measured experimentally [1][2][3][4][5].Due to the emergence of advanced computerized tomography (CT) technology, high-performance computing and advanced computational methods, the numerical study of drainage and imbibition has become feasible.The rock geometry can be digitized into a 3D voxel matrix representing solids and pores.In order to calculate the transport of multiphase flow, either pore network model (PNM) or direct numerical simulation (DNS) can be used.PNM replaces real pore bodies and throat shapes with geometric shapes such as ducts with circular, rectangular and triangular or more complex cross sections [6,7].The fluid transport and distribution are determined by employing analytical or semi-analytical solutions for the individual elements that are coupled together.DNS directly simulates two phase flow including two phase interaction characterized by interface tension and flow-wall interaction a Corresponding author: xiaobo@ingrainrocks.com represented by contact angle in complicated rock surfaces.We use a DNS method, the lattice-Boltzmann method (LBM) which directly simulates the fundamental equations of multiphase flow in the pore space using simplified kinetic equations [8][9][10][11][12].The LBM has emerged as an alternative numerical technique for simulating fluid flow for three decades.It can simulate the fluid equations accurately and efficiently, and still preserve some of the advantages of kinetic equations, such as being very suitable for parallel processing and accurately simulating flows with complicated wall geometries and interfaces.This study employs a LBM method to investigate the oil-water drainage and imbibition process in three kinds of sandstone rocks samples and effects of the geometry, initial oil situation and porosity on the results.We digitized a total of 10 rock samples including one Berea, one Bentheimer and eight Fontainebleau using micro-CT.The pore body and throat size distributions were identified using an open-map based approach [13,14].We discuss the features of the PC curves, the water and oil distributions and specially the residual oil saturation (S or ) after imbibition.
Article published by EDP Sciences

Material and simulation methods
As listed in Table 1, 10 sand stone samples are considered.We used the first 3 samples, Berea, Bentheimer and Fontainebleau1 to investigate the dependence of S or on rock geometry.In order to compare with recently published experimental data [1], we selected samples which have similar porosity, permeability and formation factor as those used in the experiment.The samples Fontainebleau 1-8 (FB 1-8) have different porosity and were used to investigate the variation of S or with porosity.
Figure 1 shows 3D pore structures and 2D slices of the Berea, Bentheimer and FB1.Their pore body shapes look different.We used a morphological approach [13] to compute pore body size distributions and a digital mercury injection experiment to compute pore throat distributions [14].As shown in Figure 2, the pore body radius distributions have small differences while the three rock samples have obviously different pore throat size distributions.The Bentheimer sample has the largest throat size and the FB1 has the smallest.The average throat sizes of the Berea, Bentheimer and Fontainebleau samples are 9.6, 12 and 7.7 µm.
A LBM solver has been developed to simulate multiphase flow in the imbibition and drainage processes.The LBM solver uses a D3Q13 multiphase model that has multiple relaxation times [15].The solver was parallelized for multiple CPUs and GPUs using MPI and CUDA.To cope with the large computational demand the simulations were run on high performance computing hardware.In the LBM the fluid wall interaction is modeled by wall wettability and fluid density fraction.The contact angle can be set by assigning certain wettability to the wall surface.A water wet surface may change to less water wet or even oil wet when oil stays in contact with that surface long enough, such as the surfaces of pores occupied by oil in a reservoir or the oil-aged rock surfaces in experiments.In our LBM simulations only the intrinsic contact angle is assigned if the wall roughness is resolved.The receding and advancing contact angles are determined by the wall wettability and roughness through the fluid-wall interaction.If the wall roughness is not resolved, receding and advancing contact angles are input parameters.In pore network modeling, the receding and advancing contact angles are always input parameters.In this study the wettability was assumed to be strongly water wet to compare to the water-wet experiments.We developed special boundary conditions at the solid wall to model thin water-or oil-film flows that allow transport of fluids through those thin films.The water film appears on the water wet solid surface and the oil film on the oil wet surface.We simulated a digital porous plate experiment as described in reference [12].In the digital porous plate experiment, the one surface of the rock is attached to an oil reservoir and the opposite surface is attached to a water reservoir.If a water film is connected to the water reservoir through bulk water or a water film, the pressure boundary condition is applied to the water that contacts with the film.Otherwise the bounce back boundary condition is applied to the water solid surface.A similar boundary condition is applied to the oil.These boundary conditions prevent e.g.: unrealistic high values of trapped water in the drainage and allow snap-off effects through the films in imbibition.In the digital porous plate experiment, the water pressure is kept constant throughout the simulation.In the drainage process the oil pressure increases step by step till the desired residual water saturation is reached.In the imbibition process the oil pressure decreases till it reaches the water pressure (only controlled spontaneous imbibition).The pressure difference between oil and water is called the applied capillary pressure.

Results and discussions
The drainage and imbibition processes were calculated for the samples listed in Table 2.A zero contact angle was assigned to the water.The interfacial tension was set to 0.03 N.m −1 .In Figure 3 the applied capillary pressure for three different rock samples is shown as function of water saturation.The water saturation at the end of drainage in the simulation was matched with the experimental data in reference [1].The petrophysical properties of rock sample and the calculated residual oil saturation for the experiments and the LBM simulations are listed in Table 2.  Reference [1] reports two sets of results for plugs that were about 8 cm in length and 3.3 cm in diameter and minplugs that were 6 mm in length and 5.9 mm in diameter.Our results for Berea and Bentheimer are close to the experiment results.The simulation result of Fontainebleau is much larger than the experimental one.However it was found that the simulation result is consistent with more systematic experimental results reported in reference [5] as shown in Figure 11.The experiment results of S or at the porosity of Fontainebleau 1 (0.114) spread from 0.52 to 0.62 and the simulation result is 0.525.It was pointed out in reference [5] that the observed smaller value of S or for Fontainebleau rocks in experiments might result from gas diffusion.Since the contact angle was the same in the simulations, the difference of S or resulted from different geometries characterized by the pore structure of the rocks.It is noticed that the residual oil satura-tion correlates proportionally with the ratio of pore body to pore throat size ratio of the rock samples shown in Figure 2.
Figure 4 shows the 3D structures of residual oil and 2D slices of oil and water distributions after imbibition for the three rock samples.The oil is trapped in the center of the pore because the rock surface is water wet.The water stays in pore throats and on solid surfaces.The shape of small oil cluster is sphere-like and the large clusters occupy several pore bodies resembling the pore structure of the rock which is consistent with the reported experimental results in reference [16].
In the drainage and imbibition processes the mean pressure gradient is applied along z direction.It is found that the trapped oil is uniform along the other two transverse directions and the trapped oil fraction increases with z coordinate.The less oil presents near the end at z = 0 202-page 4   where the water enters the rock in the imbibition process.However, as further identified by the slice averaged oil saturation in Figure 5, the trapped oil is approximately uniformly distributed in the middle region that spans about 400 voxels.The uniformly distributed oil saturation indicates that the rock sample is large enough to be a representative volume for the imbibition process.

Oil saturation Porosity
The pore size distribution of trapped oil cluster is compared to the pore body size distribution in the rock sample in Figure 6.The oil cluster size distribution is very close to the pore body size distribution and very slightly shifts to larger size.It indicates that oil is trapped in pores with different size of the rock sample.It is consistent with the oil clusters shown in Figure 4.
We investigated the dependence of residual oil saturation S or on the initial oil saturation S oi using the Berea sample.Figure 7 shows the applied capillary pressure as function of water saturation S w for the Berea sample for different initial oil saturations S oi = 1 − S w .The residual oil saturation increases with initial oil saturation.The calculated oil saturation agrees well with experimental results as shown in Figure 8 [3].In order to reduce the end effects (see Fig. 5) in the simulation the oil saturation was calculated by removing part of the sample that is located within 100 voxels from the both ends.The calculated data also agree well with the Land correlation, a widely used empirical model in literature [17]: Here S wc is the connate wetting phase saturation.It has been assumed S wc = 0 in the simulation.In Figure 9 the pore body and trapped oil cluster size distributions are plotted.The oil cluster size distribution gradually shifts to the larger size with decreasing initial oil saturation S oi .For smaller initial oil saturation, the oil occupies the larger pore at beginning of the imbibition.Therefore the trapped oil occupied the larger pores after imbibition.
We simulated the drainage and imbibition process for eight Fontainebleau samples to investigate the effect of porosity.Figure 10 shows the applied capillary pressure as function of water saturation for different Fontainebleau samples.The dependence of residual oil saturation is plotted in Figure 11.The results of residual oil saturation are consistent with two sets of published experimental results in references [4,5].Fontainebleau is a mostly clay-free sandstone.The larger residual oil saturation at low porosity may result from larger pore body to throat size ratio.Aissaoui noticed that the pore throats sizes diminish with   porosity [18].This is also true for our Fontainebleau rock samples evident in Figures 12 and 13.The pore body sizes for samples with different porosity are close while the pore throat size changes obviously with the porosity.Furthermore the pore body to throat size ratio was calculated.As shown in Figure 14, the pore body to throat size ratio and residual oil saturation S or change with porosity in a similar manner.The correlation coefficients between S or and R1 and between S or and R2 are 0.94 and 0.93 respectively.Here two kinds of ratios are calculated by R1 = r pore / r throat and R2 = r pore /r throat .The brackets indicate an average.
Figure 15 shows pore body and trapped oil size distributions for the eight Fontainebleau samples.It is similar to what is observed in Figure 6   oil is trapped in different sizes of pores and also is trapped in clusters of pores.

Conclusion
A lattice Boltzmann solver has been developed to investigate drainage and imbibition processes for different rock geometry, porosities and initial oil saturations in water wet conditions.The trapped oil after imbibition uniformly distributes in the middle of the sample.The oil is trapped in the center of pore bodies of different sizes while water stays in pore throats and on solid surfaces.The small oil clusters are sphere-like.The larger clusters occupies several pore bodies resembling the pore structure of the rock.The residual oil saturation and trapped oil cluster distributions after imbibition are consistent with recently reported experiments.The dependence of residual oil saturation on the initial oil saturation agrees well with the Land correlation.The residual oil saturation in Fontainebleau decreases with porosity and is highly correlated with the pore body to throat size ratio.

Fig. 3 .
Fig. 3. Applied capillary pressure as function of water saturation for Berea, Bentheimer and Fontainebleau sandstones.

XFig. 4 .
Fig. 4. 3D structures of residual oil and 2D slices of oil and water profiles after imbibition for Berea, Bentheimer and FB1 samples.In the 2D slices water is indicated by black, oil as dark gray and solid as light gray.In the 2D slices the z-axis is vertical and the y-axis points to the right.

Fig. 6 .
Fig. 6.Pore body and trapped oil drop radius distributions in the Berea, Bentheimer and Fontainebleau samples.

Fig. 7 .Fig. 8 .
Fig. 7. Applied capillary pressure as function of water saturation Sw for the Berea sample for different initial oil saturations Soi = 1 − Sw.

Fig. 9 .
Fig. 9. Pore body and oil drop radius distributions after imbibition for different initial oil saturations Soi.

Fig. 10 .
Fig. 10.Applied capillary pressure as function of water saturation for Fontainebleau samples.

Fig. 11 .
Fig. 11.Residual oil saturation as function of porosity for Fontainebleau samples.

Table 1 .
X. Nie et al.: Mechanics & Industry 17, 202 (2016) Rock samples used in the simulation and digital properties.Here Δx is the voxel length, L is the dimension of rock sample, φ is the porosity, Ka is the absolute permeability and FF is the formation factor.