Document Type: Research Paper
Authors
^{1} 1Graduate Student, Institute of Petroleum Engineering, Faculty of Chemical Engineering, College of Engineering, University of Tehran, Tehran, Iran
^{2} Associate Professor, Institute of Petroleum Engineering, Faculty of Chemical Engineering, College of Engineering, University of Tehran, Tehran, Iran
^{3} Researcher, IHS Global Canada Limited, Alberta, Calgary, Canada
Abstract
Keywords
Main Subjects
In describing the fluid flow in oil and gas reservoirs as porous media, relative permeability is a key petrophysical property with a significant effect on the evaluation and forecasting of the reservoir performance (Xu et al., 2013). The relative permeability curves are essential input data in oil and gas reservoir simulators and the accurate determination of their values is of great importance for various types of reservoir studies (MohamadiBaghmolaei et al., 2016). Corescale displacement experiments, including steadystate, unsteadystate, and centrifuge experiments, are performed primarily in order to generate relative permeability curves. The steadystate method is a direct method for obtaining relative permeability from experimental data. This method has the disadvantage of being timeconsuming especially for low permeability material (Kamath et al., 1995; Schembre and Kovscek, 2003).
The essence of unsteadystate methods is to estimate relative permeability curves based on measured production (oil and water) histories. Unsteadystate experiments are rapid and relatively easy to perform. The technique of Johnson et al.(1959) and that of Jones and Roszelle (Jones and Roszelle, 1978) are some known unsteadystate methods that neglect capillary pressure based on the BuckleyLeverett theory (Buckley and Leverett, 1942). Simplification assumptions make these methods generally unreliable for low permeability rocks (Akin and Kovscek, 1999). History matching methods are other types of unsteadystate methods that are based on matching experimental data. However, these methods have the disadvantage of being complicated and relatively time consuming. Furthermore, these methods need accurate experimental data, which are not easily obtained, such as a saturation profile and pressure gradient to calculate relative permeability. They need at least an Xray computerized tomography (CT) scanning method (Schembre and Kovscek, 2006). Table 1 summarizes various available major methods currently being used in the literature to estimate relative permeability.
Table 1
Summary of various available major methods to estimate relative permeability.
Steadystate methods (Christiansen and Howarth, 1995) 


Unsteadystate methods (Hussain et al., 2010) 
Analytical methods 

Semianalytical methods 


History matching 

The mentioned unsteadystate methods have been recently challenged because they assume localequilibrium in porous media during fluid flow modes. Localequilibrium leads to a selfsimilar behavior in spontaneous imbibition experiments (Li et al., 2003). In other words, the previous unsteadystate models assume capillary pressure and relative permeability to be solely functions of water saturation that lead to a selfsimilarity condition. In this regard, based on the theory that the redistribution of two phases in porous media is not instantaneous when the media is filled by water, some recent studies consider nonequilibrium effects on reservoir recovery and its properties (Schembre and Kovscek, 2006; Barenblatt et al., 2003; Silin and Patzek, 2004; Le Guen and Kovscek, 2006). Mainly focused on countercurrent spontaneous imbibition, these studies attempted to show that capillary pressure and relative permeability are not solely functions of instantaneous water saturation. Schembre and Kovscek (Schembre and Kovscek, 2006) provided a methodology to involve nonequilibrium effects in the interpretation of unsteadystate laboratory data in order to obtain steadystate relative permeability curves. Saturation profiles in their method are measured as a function of time by Xray CT during the experiments.
Empirical correlations and analytical mathematical models (Purcell, 1949; Corey, 1954; Honarpour et al., 1982; Ibrahim and Koederitz, 2000; Lomeland et al., 2005; Mosavat et al., 2013; Xu et al., 2015) are used for forecasting relative permeability curves in two different ways: first, predicting relative permeability as a function of some media properties such as capillary pressure, absolute permeability, porosity, fluid saturation etc.; second, defining relative permeability as a parametric model in which required parameters are obtained by implementing measured experimental data (MohamadiBaghmolaei et al., 2016).
In this work, we introduce an analytical approach which uses laboratory measured data as input data in an effective way to predict water and oil relative permeability curves in imbibition experiments. As mentioned earlier, taking into account the nonequilibrium effects leads to the need of recording saturation profile history as a function of time. However, our new approach uses input data to obtain a saturationdistance graph which does not change with time; therefore, there is no need to monitor the saturation distribution in the laboratory that may be time consuming and costly. This approach makes the experiment a fast way to estimate imbibition relative permeability curves, at least for nontight cores. Considering the relative permeability of a media property that may be affected by nonequilibrium effects, the saturation parameter itself is set to be timedependent in our approach, which leads to the inclusion of time effects in the interpretation. We first present the mathematical derivation of the new approach and evaluate it by related experimental data extracted from open literature. Here, any approximations and assumptions of the method are discussed. Then, we investigate the method’s applicability range, and ten synthetic cases are modeled by a numerical simulator to define the effects of various involved variables on the accuracy of the proposed method.
2. Mathematical derivation
We describe a new method that can be used to calculate water and oil relative permeabilities for waterflooding test. The four using fundamental equations are: (i) a continuity equation for a 1D fluid flow procedure in porous media, (ii) an extended form of Darcy’s law for twophase fluid flow, (iii) the cumulative oil production equation, and (iv) a water fractional flow equation for a system producing water and oil simultaneously.
In a coreflooding process, a wetting phase displaces a nonwetting phase, and this leads to the same movement direction of the two phases. We relate the inlet and outlet rates through the media by the continuity equation then combine it with the extended form of Darcy’s law for twophase fluid flow:
(1) 
To calculate water saturation at the end of the core, S_{wL}, we call the cumulative oil production equation:
(2) 
Considering , assuming that is constant, and then taking integration from the above equation leads to:
(3) 
Recalling the continuity equation while neglecting fluid compressibility with the assumption that x is not a function of and taking integration from that, we have:
(4) 
We rearrange the above equation for the end of the core, where x is replaced by core length, L, and then combine it with Equation 3:
(5) 
By taking integration from the second term of Equation 5 and then rearranging the resulting equation, S_{wL} can be calculated as:
(6) 
where, water fractional flow, f_{w}, can be obtained from production data:
(7) 
Thus, water saturation at the end of the core can be measured at any time. This timedependent parameter plays a critical role in the calculation of relative permeabilities in our approach. Water does not appear before the breakthrough time; therefore, t in Equation 6 starts from the breakthrough time and is followed by the afterbreakthrough times. Using and assuming that the pressure gradient in the oil phase is negligible lead to rearranging Equation 1 as:
(8) 
Considering that P_{c} is a function of S_{w} only, the water relative permeability is then obtained from:
(9) 
This equation suggests that water relative permeability can be calculated by finding the three derivative terms in this equation, including , , and . Because the afterbreakthrough times are used to calculate S_{wL}, we use the breakthrough time to calculate each x term related to each water saturation term at the end of the core:
(10) 
We assume that is approximated by:
(11) 
in which, is timedependent, but the graph of versus does not change with time. When we replace it in Equation 9:
(12) 
Applying the general expression for water fractional flow in a horizontal system of wateroil:
(13) 
and then rearranging the above equation under the same assumptions for , can be obtained as:
(14) 
3. Application of the method
We use recovery data of an experimental work by Qadeer et al. (2002) to demonstrate the application of the new method. These data were obtained from waterflooding into Berea Quarry sandstone, a material known as relatively high permeability (720 md) media which is also homogeneous and strongly water wet with relatively uniform properties. The average porosity of the samples was 0.2. The total designed experiments by Qadeer et al. are presented in Table 2.
Table 2
Sequence of displacement experiments by Qadeer et al., 2002.
Displacement 
Flow Rate (cc/min) 
Primary drainage Imbibition1 Drainage1 Imbibition2 Drainage2 Imbibition3 Drainage3 
8.59 8.43 8.60 3.77 4.28 1.04 1.07 
We only use the Imbibition3 experimental data in this study due to the fact that there is little recovery after the breakthrough time in the other imbibition tests. The experiments were conducted using 8% sodium bromide as the wetting phase (with initial saturation of 0.46) and cyclohexane as the nonwetting phase. The measured imbibition capillary pressure curve by the centrifuge method for similar Berea sandstone cores (Qadeer, 1988) was used as shown in Figure 1(a).
Qadeer et al. used CT imaging to measure the saturations during their displacement experiments. They obtained relative permeability curves based on the automated history matching procedure of LevenbergMarquardt as described by Press et al. (1990).
The measured recovery data by Qadeer et al. are presented in Figure 1(b) for the Imbibition3 test.
All parameters needed to estimate water and oil relative permeabilities using our approach are provided by these experimental data. Recovery data, Figure 1(b), is used to calculate f_{w}, t, and N_{p}. An area of 20.258 cm^{2} and a length of 25.4 cm are used for A and L respectively. Water saturation at the end of the core is, therefore, calculated by Equation 6 for all recorded time data. Breakthrough time is obtained from Figure 1(b), where oil recovery is changed immediately; then, can be calculated by using Equation 10. The capillary pressure curve, Figure 1(a), is used to calculate . The resulting and , respectively defined by Equations 12 and 14, are compared with the history matching results by Qadeer et al. as displayed in Figure 1(c). As shown in this figure, good agreement is obtained between the two methods.
a) 

b) 

c) 

Figure 1
a) Capillary pressure curve by the centrifuge method for Berea sandstone cores (Qadeer, 1988); b) Recovery match for the Imbibition3 test by Qadeer et al. (2002); c) Comparison of relative permeability curves for Imbibition3 test obtained by HM (Qadeer et al., 2002) and the new approach.
4. Approximations
The new estimation method in this study makes several approximations. First, the method is appropriate for experiments in which recovery is considerable after the breakthrough time because t starts from the breakthrough time in Equation 6. Second, we suppose that capillary pressure is previously estimated by another procedure and use it directly in our new method; this makes the method dependent on the input capillary pressure. After fixing the capillary pressure, however, there is no need to use and fix parametric equations to define a unique solution for relative permeability curves. This means that we do not need prefiguration of the relative permeability curves. Third, we neglect the pressure gradient in the oil phase in our derivation of relative permeabilities. Generally, this makes our method especially suitable for water/gas flows or for water into light hydrocarbons. This assumption was previously used by other methods in the literature (Hughes and Blunt, 2001; Schembre and Kovscek, 2003). Although the oil phase pressure gradient is neglected when finding water relative permeability, the method also yields nonwetting phase relative permeability according to Equation 14. Fourth, instead of any recorded water saturation distribution in the core, we assume that the slope of the graph of water saturation at the end of the core versus x, approximated by Equation 11, is used in our calculations. This assumption can be true because the trend (slope) of water saturation distribution versus the position of the fluid front, S_{w} versus x, at the breakthrough time is regenerated by versus . Figure 2 shows the similar trends of these two curves for the experimental data of Qadeer et al. (2002) at an injection rate of 1.04 cc/min. Fifth, due to the fact that we assumed unique porosity and permeability in the derivations of the method’s equations, the recovery data are obtained for a homogeneous media. Finally, because capillary end effects emerge due to the discontinuity of capillarity in the wetting phase at the outlet of the core and this effect is prominent in drainage process where the wetting phase is being displaced by nonwetting phase (Huang and Honarpour, 1998), we neglect this effect for the derivation of relative permeability in our imbibition tests.
Figure 2
Water saturation profile for the Imbibition3 test.
Besides the assumptions made above, it should be mentioned that we consider the relative permeability a property that may be affected by nonequilibrium effects in our approach where we used a timedependent saturation parameter. In (MirzaeiPaiaman et al., 2011) various contradictory claims about the existence and contribution of nonequilibrium effects to multiphase flow models was investigated, and the results show that nonequilibrium effects exist in all imbibition modes (cocurrent and countercurrent imbibition). They stated that nonequilibrium effects exist during the spontaneous imbibition process but their effect may not be significant enough to affect oil recovery versus time, while water saturation versus the similarity variable for different times will be affected.
All these approximations in combination together allow the model to successfully predict water and oil relative permeability curves in twophase fluid flow.
5. Verification of the method
The applicability of the new method in determining relative permeability curves is examined by analyzing the effects of changing some included parameters on the accuracy of the developed approach. We use a 1D numerical simulator in order to design appropriate cases to study these effects.
5.1. Simulated example
A 1D numerical simulator is used to show the efficiency of our derivations. The descriptions and details of the simulator are presented in Appendix. In order to show the application of the data generated by the simulator in our proposed method, we apply an imbibition experiment as given by Schembre and Kovscek (Schembre and Kovscek, 2003). This experiment is simulated in the numerical simulator; then, its resulted recovery data are applied to the new method. Schembre and Kovscek used an Xray CT method to obtain saturation profiles along the cores. In their method, a previously measured capillary pressure curve was used as the input data. ndecane was the nonwetting fluid, and the initial water saturation was zero. The properties of the porous media and the fluid data are presented in Table 3. The input capillary pressure and the relative permeability curves of their method are also used as the input data to the numerical simulator. These curves are shown in Figure 3(a) and Figure 3(b) respectively.
Table 3
Characteristics of the porous media and the fluid data of the case study (Schembre and Kovscek, 2003).
Item 
Value 
Permeability, k (md) Porosity,(fraction) Length, L (cm) Compressibility, c (kPa^{1}) Brine density, ρ_{w} (g/cc) Oil density, ρ_{o }(g/cc) Brine viscosity, µ_{w }(cP) Oil viscosity, µ_{o}(cP) 
6.6 0.65 8.8 7.25×10^{6} 1.0 0.8 1.0 0.86 
a) 

b) 
Figure 3
a) Input capillary pressure curves for Schembre and Kovscek’s work (Schembre and Kovscek, 2003); b) The new approach results in comparison with the input curves and Schembre and Kovscek’s (Schembre and Kovscek, 2003) relative permeability curves.
The resulting recovery data from the simulator is used in our new proposed method to calculate water and oil relative permeabilities by Equations 12 and 14 respectively. Figure 3(b) presents input relative permeability curves, the curves obtained by our proposed method and the ones obtained by Schembre and Kovscek in comparison with each other. Like the results of Schembre and Kovscek’s method, the results of our approach are well matched with the input curves in an acceptable range of water saturation.
5.2. Sensitivity analysis
We vary some involved variables in the proposed method, including the water injection rate, threshold pressure, oil viscosity, input capillary pressure, and input relative permeability curves to study the applicability of the new approach in various displacement experimental designs. As the base case, a 1D porous media with an absolute permeability of 100 md and a length of 60 cm is considered. The input relative permeability and imbibition capillary pressure functions are assumed to be expressed by BrooksCorey correlations (Brooks and Corey, 1966):
(15) 

(16) 
where, is the BrooksCorey pore size distribution index, and is the normalized water saturation expressed by:
(17) 
Table 4 tabulates the properties used for the porous media and the fluids in the base case example. The input data are used to build the base case that is modeled by the numerical simulator; then, its resulted recovery data are applied directly to the new method to regenerate the input (true) curves of relative permeabilities, independent of them. Good matched results are obtained by the comparison of the input and resulting relative permeability curves as shown in Figure 4.
Table 4
Base case media and fluid properties.
Item 
Value 
Absolute Permeability, k (md) Area, A (cm^{2}) Length, L (cm) Porosity, (fraction) Oil viscosity, µ_{o}(cP) Water viscosity, µ_{w}(cP) Entry Pressure, P_{th}(kPa) k_{rw}^{0} k_{ro}^{0} λ Water injection rate, q_{inj} (cc/min) 
100 9 60 0.4 0.5 1 20 1 1 4 0.006 
Figure 4
Comparison between the input curve (solid line) and the proposed method’s result for the base case.
We designed nine different cases to cover different ranges of involved parameters. As with the base case, by applying the numerical simulator and using its resulting recovery data as the inputs to the new proposed method, water and oil relative permeabilities are calculated. Then, in comparison with the input relative permeabilities, the accuracy of the new method is indicated. A list of the designed cases is tabulated in Table 5. The parameter studied in this table presents the change of the case against the base case.
Table 5
List of cases studied in the estimation of water and oil relative permeabilities.
Case ID 
Parameter studied 
Value 
1 2 3 4 5 6 7 8 9 
Water injection rate Water injection rate Water injection rate Threshold pressure Threshold pressure Oil viscosity Oil viscosity λ λ 
0.003 cc/min 0.008 cc/min 0.01 cc/min 10 kPa 30 kPa 1 cP 3 cP 2 ∞ 
a. Injection rate
Cases 1–3 in Table 5 are designed to study the effect of water injection rate on the accuracy of the new proposed method as shown in Figure 5(a).
It is revealed from this figure that by increasing the injection rate the accuracy of the approach increases. This can be described by considering the changes in the Capillary Number (N_{c}) as a result of changes in the flow rate. N_{c} is usually large for highspeed flows and low for lowspeed flows; thus, it is expected to increase as a result of increasing the injection rate:
(18) 
Oil is connected from the inlet to the outlet at any time during the displacement processes in the case of low N_{c} and during the commencement of the displacement processes in the case of intermediate to high N_{c}. Progressing displacement in the latter leads to decreasing oil saturation, and the oil connected region startpoint, known as the oil subnetwork, progressively moves towards the outlet (Idowu and Blunt, 2009). Since no pressure drop occurs outside the oil subnetwork, and we assume negligible oil pressure gradient in the proposed method, it is expected that the results show improved accuracy when the oil subnetwork moves faster in the region. The assumption for a much longer time during the displacement process seems more desirable. We confirm the effect of N_{c} by changing the threshold pressure, P_{th}, in cases 4–5 as presented in Table 5. For a cylindrical throat of radius r with no contact angle hysteresis, we have:
(19) 
where, P_{th} is in direct relation with interfacial tension, . Assuming no changes in the pore’s radius, r, and contact angle, , an increase in leads to an increase in P_{th} and a decrease in N_{c}. Considering the results in Figure 5(b) and the base case plot shows improved accuracy for the cases having smaller P_{th}’s.
b. Oil viscosity
The effect of changing oil viscosity is studied in cases 6–7 as shown in Figure 5(c). It is clear that by increasing oil viscosity and having it approach water viscosity the method makes no changes in the accuracy of the results. Then, by increasing oil viscosity to a value higher than that of water viscosity, the results show some errors. This refers to the approximation that makes the method accurately describe displacement processes for gases and relatively light hydrocarbons. It can be inferred from this analysis that cases with a lower wateroil mobility ratio can be better described by the approach.
c. Properties of the porous media
Since we consider Coreytype correlations to describe input relative permeability and capillary pressure curves for the numerical simulator, we can study the effects of these properties on the applicability of the new proposed method by changing λ in Equations 15 and 16. The related designed cases are cases 8–9 in Table 5. The results are shown in Figures 5(d) through 5(e) and show that by increasing λ from 2 to ∞ the accuracy of the method increases. Brooks and Corey (Brooks and Corey, 1966) stated that the values of 2, 4, and ∞ express the situations of a wide range of pore size, a medium range of pore size, and a single uniform pore size respectively. In reality, we are changing the porous media heterogeneity by changing λ. By increasing λ from 2 to ∞ to decrease the heterogeneity of the porous media, the accuracy of the method increases. However, due to the fact that we made some approximations regarding media heterogeneity in our numerical simulator, like assuming unique porosity and permeability, the recovery data are obtained for a homogeneous media. Applying these data to the proposed method may also affect the errors related to heterogeneity.
a) 

b) 

c) 

d) 

e) 
Figure 5
a) Comparison between the input curve (solid line) and the proposed method result for cases 1–3; b) Comparison between the input curve (solid line) and the proposed method result for cases 4–5; c) Comparison between the input curve (solid line) and the proposed method result for cases 6–7; d) Comparison between the input curve (solid line) and the proposed method result for case 8; e) Comparison between the input curve (solid line) and the proposed method result for case 9.
6. Conclusions
A saturation profile should be recorded as a function of time through unsteadystate experiments in order to be used in multiphase flow properties assessment, especially in the new methods, which have been introduced in the last decade and take into account nonequilibrium effects. A new interpretation method for the recovery data from laboratory experiments of waterflooding processes was introduced providing a saturationdistance graph that does not change with time and can be used instead of the common saturation distribution profile to obtain relative permeability curves. The inclusion of time in the saturation calculation in this method may be considered a way to capture the time effects in the interpretation. The experimental data from the literature was used to validate the method. The comparison of the developed model with other accepted methods revealed a good match between the results. The effects of capillary number, oil viscosity, and the heterogeneity of the porous media on the applicability of the method were investigated by a numerical simulator, and the following results were obtained.
1 A relatively high capillary number yields very good matched results.
2 Due to the method approximations, better results were achieved for cases with a relatively low wateroil mobility ratio.
3 The results of the method for the porous media with a uniform pore size show higher accuracy than the results of the method for the porous media with a wider range of pore sizes.
Acknowledgments
The authors are deeply grateful to Prof. Mehran PooladiDarvish for his valuable suggestions and recommendations that were essential in the derivation of the analytical formulations at every stage of this research.
Nomenclature
A 
: Cross sectional area 
K 
: Absolute permeability 
k_{r} 
: Relative permeability 
: End point relative permeability at residual saturation 

L 
: Length of core 
P 
: Pressure 
P_{cow} 
: Capillary pressure 
P_{th} 
: Threshold pressure 
Q 
: Rate 
R 
: Production time 
S 
: Saturation 
S_{r} 
: Residual saturation 
Greek Letters 

: Porosity 

: Viscosity 