Document Type : Research Paper
Authors
 Mojtaba Ghaedi ^{1}
 Zoltán E. Heinemann ^{2}
 Mohsen Masihi ^{1}
 Mohammad Hossein Ghazanfari ^{1}
^{1} Department of Chemical and Petroleum Engineering, Sharif University of Technology, Tehran, Iran
^{2} Montanuniversitaet Leoben, FranzJosefStrasse 18, 8700 Leoben, Austria
Abstract
In this paper, a very efficient method, called single matrix block analyzer (SMBA), has been developed to determine relative permeability and capillary pressure curves from spontaneous imbibition (SI) data. SMBA mimics realistically the SI tests by appropriate boundary conditions modeling. In the proposed method, a cuboid with an identical core plug height is considered. The equal dimensions of the cuboid in x and y directions are set such that the cylindrical core plug and the cuboid have the same shape factor. Thus, by avoiding the difficulties of the cylindrical coordinates, a representative model for the core plug is established. Appropriate grid numbers in xy and z directions are specified to the model. Furthermore, the rock and fluid properties of SI test are set in the SMBA. By supposing forms of the oilwater capillary pressure and relative permeability and comparing the oil recovery curves of SMBA and SI data, capillary pressure and relative permeability can be determined. The SMBA is demonstrated using three experimental data with different aging times. Suitable equations are employed to represent the capillary pressure and relative permeability curves. The genetic algorithm is used as the optimization tool. The obtained results, especially for capillary pressure, are in good agreement with the experimental data. Moreover, the Bayesian credible interval (P10 and P90) evaluated by the Neighborhood Bayes Algorithm (NAB) is quite satisfactory.
Keywords
Oilwater capillary pressure and relative permeability are significant data in reservoir engineering (Li et al., 2005; Hou et al., 2012). Centrifuge, steadystate, and unsteady displacement techniques are frequently employed to determine the relative permeability data (Honarpour et al., 1986). Capillary pressure is verified regularly using porousplate, centrifuge, and mercuryinjection methods (Li et al., 2005). Each of these experimental methods has its own advantages and disadvantages. Unsteadystate techniques for low permeability rocks encounter problems and steadystate methods are time consuming (Schembre and Kovscek, 2006; Hou et al., 2012). The porousplate and centrifuge methods sometimes cannot be used for core samples with a very low permeability and mercury injection technique cannot represent the actual reservoir fluids (Li et al., 2005).
For studying matrixfracture communication, SI tests, especially countercurrent SI tests, are of great importance (Haugen et al., 2014; PooladiDarvish, Firoozabadi 2000). It would be quite advantageous to determine the relative permeability and capillary pressure curves from relatively simple, fast, and economical countercurrent SI tests (Li et al., 2005; Haugen et al., 2014; Schembre and Kovscek, 2006). Additionally SI tests are much more close to the reality of the naturally fractured reservoirs; thus the acquired relative permeability and capillary pressure curves based on SI results are more reliable to study oil recovery in water invaded zones.
There are few studies on deriving oilwater capillary pressure and relative permeability from SI data. Bourbiaux and Kalaydjian (1990) used 1D finitedifference numerical model to simulate the flow behaviors during cocurrent and countercurrent SI tests. They reported that the relative permeabilities resulted in a good prediction of countercurrent oil recovery rate are about 30% less than the cocurrent relative permeabilities at a given water saturation. Behbahani et al. (2006) simulated countercurrent imbibition by using a commercial simulator and generating 1D or 2D fine grids. They compared the results with experimental measurements in the literature. Li et al. (2005) proposed a method to infer the capillary pressure and the global mobility data from SI experiments in oilsaturated rock. The presented approach lies on the plot of calculated imbibition. In their method, the relative permeability of one phase has to be known to be able to determine the relative permeability of another phase. Schembre and Kovscek (2006) used recorded water phase saturation profiles during SI by Xray computerized tomography scanning to estimate simultaneously capillary pressure and relative permeability. In their method, the simulated annealing optimization was used and they quantified nonequilibrium effects during SI. Special boundary conditions were described by Haugen et al. (2014), in which cocurrent and countercurrent flow can occur simultaneously. In the suggested method, for the period after the cessation of countercurrent production, the relative permeability behind the front and the effective capillary driving pressure at the front can be estimated. The presented method needs enough production data points to be certain and therefore long samples are required for experiments. Other modeling efforts in this area are based on some simplifying assumptions (Putra et al., 1999; Chen et al., 1995) or restricted to 1D flow (Blair 1964; Kashchiev and Firoozabadi 2003).
A method has been developed herein which determines oilwater relative permeability and capillary pressure functions based on oil recovery curve of the SI test. The proposed method was called SMBA. This approach simulates the SI process in a plug size model. The simulation method is a fully implicit control volume finite difference method. By correct handling of plug boundary conditions during SI experiment, SMBA models the process more accurately. The applicability of the SMBA to three available experimental data with different aging times is also presented. Consistent equations for relative permeability and capillary pressure were employed. In addition, the combination of GA and NAB algorithm were utilized for optimization and uncertainty quantification.
2. Experimental data
Cores used for SI test were cut with the same orientation from a sandstone outcrop. The prepared core plugs have almost similar properties (Table 1). It should be emphasized that there were no fractures in the cores utilized for the imbibition tests. Table 2 illustrates the properties of the crude oil and synthesis brine used as the fluids during the SI experiments. The interfacial tension between the brine and the crude oil was 30 dyne/cm measured by ring tensiometer. The brine contains NaCl, KCL, and CaCl_{2}.H_{2}O mixed with distilled water (Heinemann 2013).
The core plugs were washed and first dried at ambient temperature, and then oven dried at 100 °C. The dried core plugs were vacuumed and saturated with the brine. By measuring the change in weight, the porosity was calculated. Permeability to brine was determined by means of a biaxial core holder. Then, by the injection of the crude oil, the initial water saturation (S_{wi})of about 25% was established in the cores. Afterwards, the core was aged in the crude oil at 90 °C for 72 and 240 hours (Table 1). For each aging time (0, 72, and 240 hours), 3 core plugs were considered as the backup. Standard steadystate (Figure 1) and centrifuge methods (Figure 2) were used to determine relative permeability and capillary pressure respectively. Finally, the SI tests were performed using imbibition apparatus (Figure 3) at room temperature. It should be emphasized that the aging time necessary to reestablish the reservoir differs depending on the crude, brine, and reservoir rock. Generally, it should be aged for 1,000 hours (40 days) at the reservoir temperature (Anderson 1986). Since the purpose of this work was not to restore the reservoir wettability, the aging time was shorter herein.
Table 1
The parameters of cores used in the SI tests.
Sample 
Aging Time (hrs) 
L (cm) 
d (cm) 
K (md) 
φ (%) 
S_{wi} (%) 
1 
0 
8.51 
3.5 
135 
0.21 
0.25 
2 
72 
8.50 
3.5 
134 
0.208 
0.25 
3 
240 
8.50 
3.5 
135 
0.212 
0.26 
Table 2
Properties of fluids used in the SI tests at 20 °C.
Fluid 
Viscosity (cP) 
Density (kg/m^{3}) 
Brine 
1 
1010 
Crude Oil 
2.4 
730.8 
Figure 1
Experimental setup for measuring steadystate relative permeability.
Figure 2
Experimental apparatus setup for determining capillary pressure.
a) 
b) 
Figure 3
(a) Counter current imbibition test setup; (b) A quarter of the simulation model, illustrating matrix blocks (red blocks) and surrounding fractures (green blocks).
3. Description of SMBA
Simulation of a matrix block surrounded with fractures and generation of recovery curves may encounter with several challenges, including numerical instability and unrealistic modeling. Numerical instability problems originate from very small pore volumes of the grid cells, especially for the fracture cells, high permeability contrast between matrix and fracture cells, and the necessity to use bottom hole pressure controlled wells. The difficulty to handle boundary conditions can cause the unrealistic modeling of imbibition test and the determination of recovery curve. A proper simulation of imbibition test requires assuring constant pressure and saturations in the core boundaries.
The SMBA was developed as an auxiliary option in the available research code. SMBA run is similar to a typical reservoir model run. The rock and fluid properties are determined for the model. Then, a simplified model, a single matrix block surrounded by fractures, will automatically be generated and gridded. Selecting the right grid in the numerical modeling is always a challenge. Because the correct balance between the grid resolution and acceptable simulation run time has to be found, Pirker et al. (2007) made a sensitivity analysis and found that the optimal grid block number for the matrix was 8×8×36 (Figure 4). Fracture planes are added to the matrix block to finalize the simulation grid. The fracture planes surrounding the matrix block can be set in all directions or just in the vertical direction. Therefore, different boundary conditions for the matrix block can be simulated. The results of the SMBA for water displacement were verified by comparing the results with analytical BuckleyLeverett theory solution (Pirker et al., 2007). The simplifying assumptions of the BuckleyLeverett theory were applied to the SMBA model (Buckley and Leverett 1942). Figure 5 compares the SMBA and analytical runs.
Figure 4
Sensitivity analysis of grid numbers of the matrix block (Pirker et al., 2007).
Figure 5
Comparison of BuckleyLeverett and SMBA solutions in a simplified water displacement model.
Because of symmetry, only a quarter of the model is calculated in the SMBA; the side length of the cube is reduced in x and ydirection to half. To maintain the gravity force, the dimension in zdirection remains the same. Oil displacement in the matrix block during the SMBA run is shown in Figure 6.
Figure 6
Oil saturation distribution in a quarter of the matrix block with surrounding fractures during the run with SMBA.
Regarding the single matrix block with its boundaries, three types of intercell transmissibilities can be considered:
 Matrixmatrix: this is the singleporosity transmissibility between the matrixes. It should be noted that, since all the cells are of the same size (and with the same parameters), the value of transmissibility in each direction is the same for all the matrix cell pairs.
 Boundaryboundary: since the flow in the boundary blocks is unimportant and the contents of the boundary blocks will always be restored to the initial state, the boundaryboundary transmissibility is by default set to zero to prevent the uninteresting flow between the boundary blocks. The reason is to make the simulation of the imbibition process very close to the reality. Also, in the SI tests, the saturations of the boundaries of the cores change very fast and restore to the unit saturation.
 Matrixboundary: the transmissibility is calculated as in the single porosity model but by only considering the matrix cell length and permeability (since the surrounding cells act just as a boundary to the matrix cells).
The boundary medium has no capillary pressure and straightline relative permeabilities are given to it. The surrounding cells act just as a boundary to the matrix cells and they are not contributed to the flow. The volume of the boundary cells is set large enough so that they can occupy the oil flowing out from the matrix cells.
4. Application of SMBA to simulating SI tests
SMBA is a very efficient tool to simulate the imbibition tests, i.e. to model the process of oil recovery of a core sample saturated by oil and displaced by water. It credibly models the boundary conditions of the core. In SMBA, at the end of each timestep, the oil content of all of the boundary cells is removed and reported as the amount of the extracted oil. Then, the fracture dynamic properties are set back to their initial value. In this way, the problem of constant bottom hole pressure wells is avoided as well. An effective timestep control and NewtonRaphson iteration method has been implemented, which helps achieve a much better stability.
Core plugs are often in cylindrical shape. The simulation of cylindrical models has its own difficulties and challenges. Therefore, to avoid these problems, cuboid matrix block model is generated with the same height as the core plug. The side directions are selected so that the matrix block model and the core plug have the same shape factor value (2.67 cm^{2}). Then, the resulting recovery curves of these two models will be the same. Equation 1 describes the generalized shape factor relation (Kazemi et al., 1992):
(1) 
where, s is the shape factor, V_{b} is the bulk volume of the matrix block, and A_{i} stands for the surface area open to imbibition in the i^{th} direction; x_{i }is the distance from A_{i} to the no flow boundary and n represents the number of the surfaces open to flow. Based on Equation 1, to obtain the same shape factor for the two models, if the matrix block is given the core plug height, h, the size of its sides has to be set to the same core plug diameter (d) (Figure 3).
If the small cell sizes of the single matrix block make numerical problems, higher pore volumes can be given to the cells of the matrix block. To obtain relative permeability and capillary pressure relations, instead of matching the recovery curve versus time, scaled recovery curves can be matched. In the most of the previously developed scaling relations, dimensionless time, t_{D}, is inversely proportional to the square root of the porosity (Ma et al., 1997; Mattax and Kyte, 1962; Zhang et al., 1996; MirzaeiPaiaman and Masihi, 2013; Schmid and Geiger, 2013; Pordel Shahri et al., 2012); for example, the following equation may be used (Ma et al., 1997):
(2) 
In this equation, t, L_{C}, K, and φ are imbibition time, a characteristic length, permeability, and porosity respectively. δ stands for interfacial tension between the phases and μ_{w} and μ_{nw} are the viscosities of the wetting and nonwetting phases respectively. Thus if higher pore volume is assigned to the matrix block to avoid stability problems, except for porosity, the other parameters of the t_{D} for the single matrix block model and core plug are the same. Therefore, if the right capillary pressure and relative permeability functions are set in the SMBA, the plots of oil recovery versus tφ^{0.5} for the simulated and experimental SI are matched.
5. Formulations
5.1. Capillary pressure
Many capillary pressure correlations have been recommended in the literature; Skjaeveland et al. (2000) developed a general capillary pressure correlation:
(3) 
where, S_{w} is water saturation and S_{o} is oil saturation; C_{w} and C_{o} are the entry pressures for oil and water respectively; 1/a_{w} and 1/a_{o} represent the pore size distribution indexes; S_{wR} stands for irreducible water saturation and S_{oR} denotes residual oil saturation. There is one set of C_{w}, a_{w, }C_{o},and a_{o} for imbibition and another for drainage. C_{o }is a negative number and the other constants are positive. In a completely waterwet system, C_{o }is zero(Skjaeveland et al., 2000).
5.2. Relative permeabilities
Two general groups to represent relative permeability curves are the parametric and the nonparametric models. The Coreytype and powerlaw relations are common parametric models (Hou et al., 2012; Skjaeveland et al., 2000). Expressions have also been developed by Skjaeveland et al. (2000) to show relative permeability of water and oil:
(4) 

(5) 
In Equations 4 and 5, k_{rw,ww}is the water relative permeability in a completely waterwet system and k_{ro,ww}represents the corresponding oil relative permeability. Additionally, k_{rw,ow }and k_{ro,ow} are water and oil relative permeability in a completely oilwet medium respectively. Also, is oil relative permeability at irreducible water saturation and stands for water relative permeability at residual oil saturation.
In a completely waterwet system, the following relations can be written for oil and water relative permeabilities (Skjaeveland et al., 2000):
(6) 
where,
(7) 
Similarly in a completely oilwet system, the following equations can describe oil and water relative permeabilities:
(8) 
where,
(9) 
The selected equations for capillary pressure and relative permeabilities have constant parameters in common. This makes the optimization process of finding suitable capillary pressure and relative permeabilities functions more efficient, because less unknown variables have to be determined.
5.3. Objective function
The considered shape of the objective function is illustrated in Equation 10:
(10) 
where, M is a misfit valueand N_{t} is the number of the measured oil recovery curve to be matched; stands for the experimental value of oil recovery at the corresponding t and represents the simulated amount of oil recovery at the equivalent t.
6. Workflow
The workflow of this paper is made of two steps (Figure 7). In the first step, the oil recovery curve output of the SMBA and SI test were matched. In this step, GA was coupled with the SMBA to determine the best and ensemble of capillary pressure and relative permeability functions. The second step involved the statistical inference of the ensemble of the obtained answers.
MATLAB OptimizationToolbox and the NAB method of Sambridge (1999) have been employed as the optimization tool and statistical analysis respectively. NAB is a subset of the Markov chain Monte Carlo method. This technique builds an approximation for the real posterior probability distribution using a Gibbs sampler. Then, by building a multidimensional interpolator in the model space, information is gathered from an ensemble of the forward solutions. In this way, any Bayesian integral can be evaluated. There is no need for more forward simulations, and thus this method is quite fast. For more information about this method see (Sambridge 1999).
Figure 7
Workflow of the paper, oil recovery matching and statistical analysis.
7. Results and discussion
The matching parameters of the capillary pressure and relative permeability functions were C_{w}, a_{w, }C_{o }, a_{o},,, S_{wR}, and S_{oR}. At first, different combinations of GA parameters were tested. Table 3 shows the selected configuration of the GA parameters. The parameters of the best matched models are tabulated in Table 4. During the optimization process, an ensemble of 964 solutions was obtained. They were submitted to the NAB method for uncertainty quantification.
Table 3
GA parameters used for the matching oil recovery curve.
Population 
Maximum no of generation 
Cross over type 
Cross over Rate 
Mutation Rate 
40 
40 
Scattered 
0.75 
0.04 
Table 4
Parameters of the best matched models.
Parameter 
0 hr aging time 
72 hrs aging time 
240 hrs aging time

C_{w} 
1.31 
0.78 
1.25 
a_{w} 
0.21 
0.22 
0.2 
C_{o} 
0.001 
1.01 
1.7 
a_{o} 
0 
0.01 
0.001 
S_{wR} 
0.25 
0.24 
0.27 
S_{oR} 
0.21 
0.31 
0.32 
0.31 
0.35 
0.41 

0.95 
0.92 
0.86 
Figures 814 compare the experimental data with the simulation results. An excellent match between the experimental and the best calculated oil recovery factors for the three core samples can be seen in Figure 8. According to this figure, as expected, the imbibition rate decreases by increasing the aging time (t_{a}). The reason of the strong effect of the aging time on the data can be interpreted as the fast adsorption of the wettabilityaltering compounds (polar compound and/or the deposition of the organic matters originally in the crude oil) during the aging process.
Figures 911 present the calculated capillary pressure for the three core samples. The best calculated capillary pressures are in good agreement with the experimental data. Furthermore, the Bayesian credible interval (P10 and P90) is acceptable. The credible interval was obtained by submitting the 964 solutions to the NAB method. A Bayesian probability interval, also known as credible interval, includes information from the prior distribution into the assessment. In Bayesian inference, a probability interval is a probabilistic area around a posterior moment and using a probability interval is similar to a frequentist confidence interval. For example, if a statistical method approximates the mean of a_{o }of a population to be 0.1, with a 90% credible interval of 0.0080.11, then the estimated posterior probability is that 90% of that population is between 0.008 and 0.11, with a mean of 0.1. These figures also confirm that Equation 3 is general and flexible to represent the capillary pressure curve.
The predicted oilwater relative permeabilities are illustrated in Figures 1214. The best solutions and the predicted intervals are to some extent acceptable. The Nonparametric models such as cubic Bspline model are more general and flexible (Hou et al., 2012). Using these models to represent relative permeability functions may improve the predicted relative permeability. The reason of using Equations 4 and 5 to depict relative permeability relations is the common unknown parameters of Equation 3 utilized for the capillary pressure curve, because the number of the unknowns, which have to be determined during the optimization step, reduces. This may help avoid nonunique solutions for the relative permeabilities and capillary pressure functions, which is the widespread problem of the inverse problems. Another advantage of using setup Equations 35 is the agreement in the shapes of the resulted relative permeabilities and capillary pressure; for example, both relative permeabilities and capillary pressure curves show that the system is waterwet.
The relative permeabilities of cocurrent and countercurrent processes are not the same (Bourbiaux and Kalaydjian, 1990). This could be another source for the differences of the predicted and experimental relative permeabilities. The experimental procedure to calculate the relative permeability has been the standard steadystate method, while the SI experiments are countercurrent imbibition.
Figure 8
Comparison of the experimental and the best calculated oil recovery factors for the three core samples.
Figure 9
Comparison of the experimental capillary pressure data with the calculated capillary pressure functions (the best and uncertainty interval (P10 and P90)) for the sample with 0 hr aging time.
Figure 10
Comparison of the experimental capillary pressure data with the calculated capillary pressure functions (the best and uncertainty interval (P10 and P90)) for the sample with 72 hrs aging time.
Figure 11
Comparison of the experimental capillary pressure data with the calculated capillary pressure functions (the best and uncertainty interval (P10 and P90)) for the sample with 240 hrs aging time.
Figure 12
Comparison of the experimental oilwater relative permeabilities data with the calculated relative permeabilities functions (the best and uncertainty interval (P10 and P90)) for the sample with 0 hr aging time.
Figure 13
Comparison of the experimental oilwater relative permeabilities data with the calculated relative permeabilities functions (the best and uncertainty interval (P10 and P90)) for the sample with 72 hrs aging time.
Figure 14
Comparison of the experimental oilwater relative permeabilities data with the calculated relative permeabilities functions (the best and uncertainty interval (P10 and P90)) for the sample with 240 hrs aging time.
8. Conclusions
This paper presented a very efficient method, called SMBA, to infer relative permeability and capillary pressure functions from SI data. The proposed method handles the boundary conditions more reliably, and therefore it simulates the SI tests realistically. In SMBA, a cuboid model with the same core plug height is generated. By assuming the functions of the oilwater capillary pressure and relative permeability and comparing the oil recovery curves of SMBA and SI data, capillary pressure and relative permeability can be determined. The applicability of the method was illustrated on three core samples with different aging times. Appropriate relations were utilized to represent the capillary pressure and relative permeability curves. Moreover, GA was used as the optimization technique. The predicted forms, especially for the capillary pressure curves, are in very good agreement with the experimental data. The Bayesian credible intervals (P10 and P90) predicted by the NAB method are acceptable as well.
Acknowledgements
We wish to thank Prof. M. Sambridge for providing us with his neighborhood algorithm.
Nomenclatures
1/a_{o} 
: Pore size distribution index for oil 
1/a_{w} 
: Pore size distribution index for water 
A_{i} 
: Surface area open to imbibition in the ith direction 
C_{o} 
: Entry pressures for water 
C_{w} 
: Entry pressures for oil 
d 
: Core plug diameter 
h 
: Core plug height 
K 
: Permeability 
: Oil relative permeability at irreducible water saturation 

: Water relative permeability at residual oil saturation 

k_{ro,ow} 
: Oil relative permeability in a completely oilwet system 
k_{ro,ww} 
: Oil relative permeability in a completely waterwet system 
k_{rw,ow} 
: Water relative permeability in a completely oilwet system 
k_{rw,ww} 
: Water relative permeability in a completely waterwet system 
L_{C} 
: Characteristic length 
M 
: Misfit value 
n 
: Number of the surfaces open to flow 
N_{t} 
: Number of the measured oil recovery curve 
S_{o} 
: Oil saturation 
S_{oR} 
: Residual oil saturation 
S_{w} 
: Water saturation 
S_{wR} 
: Irreducible water saturation 
t 
: Imbibition time 
V_{b} 
: Bulk volume of the matrix block 
x_{i} 
: Distance from A_{i} to the no flow boundary 
: Experimental value of oil recovery at the corresponding t 

: Simulated amount of oil recovery at the equivalent t 

δ 
: Interfacial tension 
s 
: Shape factor 
μ_{nw} 
: Viscosity of the nonwetting phase 
μ_{w} 
: Viscosity of the wetting phase 
φ 
: Porosity 