Co-estimation of saturation functions (kr and Pc) from unsteady-state core-flood experiment in tight carbonate rocks

The production from oil and gas reservoirs is greatly affected by rock and fluid properties of the porous rock. Capillary pressure (Pc) and relative permeability (kr) are two important properties employed in the mathematical simulation reservoirs for predicting oil recovery from underground hydrocarbon resources. In this study, various core-flood experiments were performed using different tight carbonate rock samples for oil–water and oil–gas systems. The objective of this research is to investigate the multi-phase flow functions (kr and Pc) in tight formations. The kr curves of each sample were obtained by two different mathematical methods: the history-matching (ant colony optimization) technique and analytical method (JBN). The comparison between the relative permeability of the history-matching technique with that of the JBN method revealed a significant discrepancy between them. The modeling of an experiment using kr of JBN revealed a significant difference between experimental and simulation oil production, whereas the relative permeability of history matching accurately reproduced the experimental oil recovery. This observation highlights the inadequacy of the JBN technique for determination of relative permeability in particular in the tight rock where capillary forces are important. In addition to the relative permeability, the capillary pressure values as a function of saturation were estimated from core-flood tests using a history-matching technique. The comparison between oil–water capillary pressures obtained from centrifuge tests with those of core-flood experiments depicted good agreement, whereas the capillary pressure of oil–gas system measured from core-flood tests was considerably different from centrifuge experiment results. This outcome demonstrated that the capillary pressure obtained from centrifuge experiments in some cases may not be representative of dynamic capillary pressure governing the multi-phase flow in porous media.


Introduction
Due to the great importance of relative permeability (k r ) and capillary pressure (P c ) curves for the mathematical modeling of the multi-phase flow in porous media (e.g., petroleum reservoir, carbon dioxide geological storage, and underground water resources), many types of research have been directed toward the prediction of these parameters. Over the past years, various numerical and experimental methods have been developed for the direct measurement and estimation of both k r and P c (O'Meara Jr and Lease 1983; Watson et al. 1988;Nordtvedt et al. 1993;Shahverdi et al. 2011;Shahverdi 2012;Pini and Benson 2013). Shahverdi (2012) developed a new method for simultaneous estimation of three-phase k r and P c values based on history matching of unsteady-state experimental data using the genetic algorithm. He verified his method by making use of two sets of three-phase experimental data published in the literatures and then applied his methodology to calculate three-phase k r values of different WAG cycles. Pini and Benson (2013) used three different fluid pairs and measured k r and P c drainage curves simultaneously on a single Berea Sandstone core. They proved reliability of the measured k r curves by comparing their measured data on Berea Sandstone cores and with various gas/liquid pairs and those obtained from the literature obtained.
The most common experimental methods for the measurement of relative permeability are steady-state and unsteady-state experiments. In the steady-state method, two immiscible fluids are injected co-currently at a specific ratio through the core until the same production ratios are achieved from the outlet of the core. In this method, the relative permeability values can be calculated directly using Darcy's law. Although a steady-state test facilitates the direct determination of k r curves, it is a very time-consuming and costly method with some limitations such as neglecting capillary effects, one-dimensional flow in the core, and isothermal and incompressible fluid (Cinar et al. 2007). Moreover, the ratio of the injected fluids is not technically allowed to be less than a certain value, hence the relative permeability at low saturations cannot be measured. An alternative technique for determining the relative permeability values is the unsteady-state method which is much quicker than the steady-state test. In the unsteady-state test, one fluid (which is immiscible with the fluid in the core) is injected through the core to displace the resident fluid inside the core. Then the fluid production and pressure across the core are measured against time. Since the flow in the porous media is occurring under the unsteady-state condition, the relative permeability data cannot be obtained directly from Darcy's equation. Hence, two general methods have been devised for the measurement of relative permeability from unsteady-state tests: the analytical (explicit) method and the optimization (implicit) method. In the analytical method, the Buckley-Leverett theory (Buckley and Leverett 1942) is used to derive the relative permeability from the production and pressure profile (Johnson et al. 1959). However, in the optimization technique, the relative permeability curve is tuned in a simulation-optimization procedure to reach an acceptable match between the laboratory results of fluid production and pressure with the corresponding values calculated by the model (Shahverdi et al. 2011). These methods are described in detail later in this manuscript.
Another experimental method for estimating relative permeability that was first introduced by Purcell (1949) involves using the capillary pressure data. In this method, the porous medium is considered as a bundle of capillary tubes with different diameters that are not connected to each other. Purcell (1949) derived an integral formula representing the relative permeabilities as a function of capillary pressure data. This method was later modified by other researchers (Fatt and Dykstra 1951;Burdine 1953).
In addition to the experimental methods, modeling approaches such as pore-network modeling and empirical correlations have been developed for predicting relative permeability (Rajaram et al. 1997;Laroche and Vizika 2005;Sinha and Wang 2007). In the pore-network modeling approach, the void space of a rock is described as a network of pores connected by the throats. The pores and throats are assigned some idealized geometry and rules are developed to determine the multi-phase fluid configurations and transport in these elements. The rules are combined in the network to compute effective transport properties on a macroscopic scale some tens of pores across (Blunt et al. 2002). The triumph of pore-scale network models depends strongly upon the quality of representing real pore spaces based on their topographical and geometrical characteristics (Xiong et al. 2016).
Various empirical correlations were developed for predicting two-phase and three-phase relative permeability as a function of fluid saturation (Corey 1954;Hirasaki 1973;Sigmund and McCaffery 1979;Honarpour et al. 1982;Chierici 1984;Lomeland et al. 2005;Shahverdi and Sohrabi 2013;Shahverdi and Sohrabi 2014). These correlations contain some parameters depending on the rock and fluid condition that should be adjusted using experimental data.
The geological specification of the rock is one of the key factors that significantly affect multi-phase flow in porous rocks. Many studies have been directed toward investigations of flow functions (k r and P c ) in sandstone rocks (Oak 1991;Dana and Skoczylas 1999;Haugen et al. 2014;Krause and Benson 2015;Moghadasi et al. 2015;Kianinejad et al. 2016), while much fewer works have been carried out to investigate k r in carbonate rocks (Spearing et al.;Schneider and Owens 1970;Sola et al. 2007). In general, the carbonate rocks have complex textures and pore networks with high degrees of heterogeneity, unlike sandstone rocks which are more homogeneous in terms of pore structure. In addition to the above features, the capillary forces play a substantial role in governing the fluid flow in carbonate rocks because they typically have low permeability. Hence capillary pressure should be taken into account in simulation (Roehl and Choquette 2012).
In this study, various core-flood experiments are performed on carbonate rocks (with the order of 1 mD permeability) using the oil-water and oil-gas system in order to investigate the performance of the relative permeability in the very low permeability formations. Different mathematical methods are used to estimate the relative permeability and capillary pressure curves from the core-flood experiments. Then the accuracy of each method is discussed. The capillary pressure curves measured from core-flood experiments are also compared to those obtained from centrifuge tests based on which some conclusions are drawn.

Experiment
In this research, four different carbonate rocks taken from an Iranian oil reservoir have been used to perform the unsteadystate core-flood experiments. The physical properties of the core samples are given in Table 1. As can be seen, the order of the permeability of the rock samples used in this study is around 1 mD which is much less than the permeabilities found in conventional reservoir rocks. The fluid properties used in the experiments are given in Table 2.
Water injection and gas injection tests were carried out into the cores which were initially saturated with the oil and irreducible water (connate water). The experiments were performed at reservoir conditions and under constant pressure drop across the core; hence, the injection rate varies with time. The individual fluid production from the outlet of the core and the injection rate against the time were recorded as the experimental results. The detailed procedure for the experiments is described below: 1. The core plugs were cleaned with toluene and methanol through a Soxhlet extraction apparatus at 60-70 °C to remove residual hydrocarbons, formation brine, salts and other contaminants. Then they were dried in an oven at 70 °C for 8 h. A helium porosimeter and liquid permeameter were used to measure the porosity and permeability of the core. 2. Each core plug was placed in the core holder and then flooded continuously with brine using a constant pressure pump until the core was fully saturated with the brine. As fluid injection proceeds, brine injection rate is measured at different time steps; if it reaches a stable constant value, the only produced fluid is brine which means the core is fully saturated with brine. 3. Crude oil was injected into the core sample to displace the brine and establish the connate water saturation (S wi ) in the core. The oil was injected at the representative fluid velocity in the reservoir and the flooding continued until no further brine was produced, i.e., any remaining water in the core is immobile. 4. The core was aged in the crude oil at the temperature of 65 °C for 15 days to restore reservoir wettability condition.
5. Water or gas was then injected to displace the oil in the core. The cumulative oil production and fluid injection rates were recorded versus time. The fluid injection continued until no further oil was produced from the core, i.e., residual oil saturation (S orw for water injection test and S org for gas injection test) was attained in the core. It should be noted that during the whole process the pressure drop across the core was held constant.
The schematic of the above procedure and the state of saturation in the core at different stages for the unsteady-state water-flooding are depicted in Fig. 1. It should be noted that a similar process was performed for the gas injection test,  where nitrogen was injected into the core with a constant pressure drop across the core, and the oil production and gas injection rates were recorded during the injection process. The initial and boundary conditions employed in the experiment are presented in Table 3.

Mathematical method
The most common method for estimating the k r values from the unsteady-state test is an analytical method that implements the Buckley-Leverett theory (Buckley and Leverett 1942) to derive the relative permeability curve from the results of core-flood test (e.g., production and pressure data). This method was originally devised by Johnson-Bossler-Naumann (Johnson et al. 1959). Later in 1998, (Grader and O'Meara Jr 1988) modified the JBN method to present a more straightforward formulation and procedure for estimating relative permeability curve from the experimental results. In this method, first the saturation of the phase j at the outlet of the core is calculated by the following equation: where S o j is the initial saturation of phase j, PV is the amount of the injected fluid as a fraction of pore volume, R j is the fluid production of phase j as a fraction of pore volume, and S j is the saturation of phase j at the outlet of the core.
Having calculated S j , the corresponding relative permeability at this saturation can be computed using: where k rj is the relative permeability of phase j, q t is the total production flow rate (cc/h), f j is the fractional flow of phase j, which is equal to the ratio of q j to q t , k is absolute permeability (mD), A is the cross section of the core (cm 2 ), j is the viscosity of phase j (cp), L is the core length (cm), ΔP is the pressure drop across the core (atm), PV is pore volume of the core (cm 3 ).
Despite the simplicity of the JBN method, the formative assumptions of this method (e.g., numerical differentiation error, neglecting the effect of capillary forces and flow in one direction) are very limiting and lead to inaccurate results in many cases. The numerical differentiation of the experimental data (in Eqs. 1-3) can lead to erroneous results, in particular, when the experimental data points are affected by fluctuations. Also ignoring the effect of capillary forces can result in a significant inaccuracy of the relative permeability data especially for the tight rocks. A more accurate method for determining of the relative permeability from the unsteady-state test overcoming the aforementioned limitations is the implementation of an optimization technique. In this approach, a relative permeability and capillary pressure model are tuned in an optimization   procedure to reach an acceptable match between the laboratory results of the fluid production, pressure drop and the corresponding values calculated by the model. The algorithm of the optimization technique is given in Fig. 2. Firstly, the rock and fluid properties (i.e., the size of the core, porosity, permeability, fluid density, and viscosity) measured in the laboratory, together with an initial guess for relative permeability and capillary pressure, are used in a core-flood simulator to obtain the fluid production and pressure against time. Having simulated the core, the difference between experimental production, pressure and those of the simulation data are calculated as the objective function. Then an optimization algorithm is used to minimize the objective function by adjusting the unknown parameters (k r and P c ). The core-flood simulator in the optimization algorithm (Yaralidarani and Shahverdi 2016) is responsible for simulation of the multi-phase flow in the rock accounting all the governing forces (i.e., capillary, viscous and gravity). The outputs of the simulation are saturation and pressure distribution along the core as well as the fluid production. The mass conservation law (continuity equation) and Darcy's law are two fundamental equations used in the numerical simulation. For the oil-water system, the set of equations are as follows: where ws is water density (g/cc) at standard condition, w is water viscosity (cp), B w is water formation volume factor (cc/Scc), k is absolute permeability (mD), k rwo is two-phase water relative permeability, k row is two-phase oil relative permeability, ∇P w is pressure gradient of water (atm/cm), q Sw is source/sink term related to water at the standard condition (cc/h), ∅ is porosity, S w is water saturation, t is time (h), wo is oil density at standard condition (g/cc), o is oil viscosity (cp), B o is oil formation volume factor (cc/Scc), S o is oil saturation and q So is source/sink term related to water at standard condition. Capillary pressure is defined as follows: where P o and P w are water and oil phase pressure (atm), respectively. In addition to the above equation, considering the definition of saturation leads to the following equation: By solving Eqs. (4-7) simultaneously, the unknown variables (e.g., P o , P w , S w , and S o ) and fluid production are obtained. It should be noted that similar equations are employed for the oil-gas system. In this study, a modified version of the ant colony optimization (ACO) technique is used as an optimization tool in conjunction with a core-flood simulator (as described above) to obtain the relative permeability and capillary pressure simultaneously from the unsteady-state tests (Yaralidarani and Shahverdi 2016).
The aim of the optimization technique presented in Fig. 2 is to obtain the relative permeability and capillary (6) P cow = P o − P w , pressure curves that reproduce a reasonable match of fluid production and pressure with the corresponding experimental data. The objective function (shown in Fig. 2) is defined as follows: In the above equation, o and I are standard deviations of experimental values of the oil production and fluid injection rates, respectively. These quantities represent the errors associated with the core-flood experiments, which can be obtained by repetition of the experiment. The terms W o and W I are weighted coefficients compensating for the magnitude difference between the two quantities (oil production and fluid injection rate). Finally, Q o represents the oil production and q f is the fluid injection rate. The proper adjustment of these factors allows us to match fluid production and pressure to the same extent. (Yaralidarani and Shahverdi 2016).
The relative permeability and capillary pressure curves are represented by a flexible model in the optimization f,i 2 . Fig. 4 Experimental results of cumulative oil production compared to the simulation results using k r of HM and k r of JBN method loop. In this research, the relative permeability and capillary pressure data versus water saturation are produced from the LET model (Eqs. 9-11) (Lomeland et al. 2005) and Skjaeveland's model (Eq. 12) (Skjaeveland et al. 1998), respectively.
where S w is water saturation, S o is oil saturation, S * w is normalized water saturation, S wi is irreducible water saturation (connate water), S orw is residual oil saturation, k rwo and k row are water and oil relative permeability in oil-water system, respectively and P cow is oil-water capillary pressure. In Eqs. (9) and (10), parameter T describes the upper part of the curve, L parameter describes the lower part, and parameter E controls the curvature or elevation of the k r curves. The parameters, k o rwo and k o row have physical meaning and represent the maximum two-phase relative permeability of water and oil, respectively. In Eq. (12), c w and a w account for the positive part of capillary pressure curve while c o and a o control the negative part.
Similarly, LET equations may be used for producing gas-oil relative permeability values. Contrary to the water injection tests which are among imbibition processes, gas injection is a drainage process, and hence gas-oil capillary pressure values are produced using the following equation which is a shortened form of the Skjaeveland's capillary pressure model:

Fig. 5
Experimental results of water injection rate compared to the simulation results using k r of HM and k r of JBN method where S * g is normalized gas saturation, S g is gas saturation, and S org represents residual gas saturation.
The initial and boundary conditions of the model are given in Table 3. The parameters of the relative permeability and capillary pressure model (Eqs. 9-12) are accounted as unknown parameters and the modified ant colony optimization (ACO) algorithm is used to adjust these functions (Yaralidarani and Shahverdi 2016). In other words, the ACO algorithm changes these parameters to obtain an acceptable agreement between experimental data (i.e., cumulative oil production and water injection rate) and those of simulation. Hence, the unknown variables can be represented as the following vector for oil-water tests: and as follows for gas-oil tests:

Fig. 6
Oil-water capillary pressure versus normalized water saturation obtained from the centrifuge test and that resulting from history matching of the water injection tests Range of initial guesses within which unknown variables are changed automatically by the program to match experimental data is presented in Table 4. It is worth noting that these values have been obtained using sensitivity analysis. It should also be added that 'c w ' and 'a w ' are set equal to zero manually, because they are representative of positive part (spontaneous imbibition) of P c curve which does not apply to our oil-water experiments. Also note that ending criterion has been defined as 5% discrepancy between the experimental and history-match data.

(a) Water injection test
The experimental results of the water injection tests of different cores (Table 1) were history-matched by the (16) ACO method as described in the previous section. In this method, the relative permeability and capillary pressure were tuned simultaneously. Also, another set of relative permeability curves were obtained from the water injection test by the JBN technique. The oil-water relative permeabilities versus water saturation obtained from these two methods are presented in Fig. 3. As it was highlighted earlier, the k r obtained from the history matching is much more realistic than that of the analytical method. The history-matching technique takes into account different kinds of forces (i.e., capillary, viscous and gravity) in the flow calculation, whereas the JBN method ignores the capillary and gravitational forces. The results shown in Fig. 3 depict that the k r obtained from the JBN method cannot adequately predict the k r of history matching. This significant disagreement is mainly attributed to the impact of capillary pressure which is pronounced for the tight rocks. This observation is very important because it suggests that the routine approach (JBN method) to deriving the relative permeability curve from the unsteady-state experiments for use in reservoir engineering calculations can produce erroneous results in particular for the tight rocks (carbonate). In addition, JBN method assumes 1D flow and The k r obtained from three-dimensional HM and JBN techniques was employed in the core simulation to reproduce the oil recovery and injection rate of the experiment.
The cumulative oil production, as well as water injection rate versus time obtained from experiment and those resulting from simulation by JBN and HM relative permeability curves, is depicted in Figs. 4 and 5, respectively. As can be seen, using k r from the JBN method substantially overestimates or underestimates the oil production and injection rates of all experiments. These results confirm the previous conclusion that the JBN method cannot accurately estimate relative permeability from core-flood tests.
The oil-water capillary pressure versus normalized water saturation obtained from centrifuge experiment and history matching of the water injection test is presented in Fig. 6. As can be seen, there is good agreement between experimental and history-matched P c curves.
Estimated parameters of the LET relative permeability model and Skjaeveland capillary pressure model are provided in Table 5. Figure 7 shows the relative permeability of the oil-gas system obtained by the history matching (ACO method) and implementation of the JBN technique on the gas injection test. Relative permeabilities in Fig. 7 show large discrepancies between the JBN and HM methods, especially for the gas phase. In order to realize which method (JBN or HM) estimates the relative permeability more accurately, the relative permeability of both techniques was employed in the simulation of gas injection tests to reproduce the oil recovery and injection rate. Figure 8 presents the cumulative oil Fig. 8 Experimental results of cumulative oil production compared to the simulation results using k r of HM and k r of JBN method production versus time obtained from the gas injection test and simulations using k r from the HM and the JBN methods. Figure 9 shows the cumulative gas production versus time obtained from experiment and simulations using k r from the HM and the JBN methods. As shown in these figures (Figs. 8 and 9), similar to the oil-water system, the simulation results using k r from the JBN technique are largely different from the experimental results, demonstrating the inaccuracy of the relative permeability data estimated using the JBN technique.

(b) Gas injection test
In order to ensure accurate performance of the ACO algorithm, evolution of objective function (Eq. 8) versus iteration number is depicted in Fig. 10 for water injection tests. The declining trend obviously shows that objective function is being minimized meaning a good match between experiment and simulation data is achieved.
These results indicate the significant impact of oil-gas capillary on the relative permeability. In other words, neglecting capillary pressure in the JBN technique leads to the erroneous oil-gas relative permeability curve. Figure 11 presents the oil-gas capillary pressure versus normalized oil saturation for various cores estimated from gas injection tests using the history-matching technique and centrifuge experiment. As can be seen, the oil-gas capillary pressure predicted by the history-matching method is considerably different from those of the centrifuge test. This figure illustrates that centrifuge test may not measure the oil-gas capillary pressure precisely. This discrepancy may be attributed to two reasons. First, the centrifuge capillary pressures are measured for an oil-air system and then scaled for the oil-gas system using an interfacial tension (IFT) correction factor. Hence, the scaling of capillary pressure data from air-oil to the oil-gas may be associated with some errors. The second reason is that the capillary pressures of the core-flood test are tuned simultaneously with the relative permeability. In other words, the impact of relative permeability is accounted for in the estimation of capillary pressure from core-flood test, whereas the capillary pressure of the centrifuge test is calculated without consideration for the effect of relative permeability. The authors believe that the capillary pressure of the core-flood test is more representative of the dynamic capillary forces governing flow in the reservoir. Hence, it is advisable for reservoir Fig. 9 Experimental results of cumulative gas production compared to the simulation results using k r of HM and k r of JBN method engineers to use core-flood capillary pressure for oil reservoir simulations.
The parameters of the LET model, as well as the Skjaeveland capillary model, calculated using the history-matching method are provided in Table 6.

Conclusion
The saturation functions (k r and P c ) were successfully coestimated from unsteady-state core-flood experiments performed in tight carbonate rocks, using ant colony optimization (ACO) technique. In order to accurately incorporate the physic of the flow, the simulation model was built in three dimensions (X, Y, and Z) with fine grids considering the porosity and permeability distribution (to account for rock heterogeneity). Moreover, the good agreement between simulation and experimental results was obtained at reasonably low iteration number in optimization technique.
The results depicted that neglecting capillary pressure in core-flood simulation can lead to erroneous relative permeability data. Although adding the capillary pressure function in history-matching process can increase the uncertainty of estimated functions, having an approximate value for P c (from other tests such centrifuge) as initial guess can modify the estimated k r and P c .
The comparison between the results obtained from the analytical method (JBN) and those calculated from the optimization technique revealed significant difference between the two methods. The difference is primarily because the JBN technique ignores the effect of capillary and gravitational forces. In addition to this, the analytical method assumes one-dimensional flow with the The oil-water capillary pressure derived from the coreflood experiment reasonably agreed with that obtained from the centrifuge test, while the capillary pressure of the oil-gas system for the core-flood test was significantly different from the centrifuge test. This has serious implications for reservoir simulation studies that use centrifuge data as simulation inputs.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons .org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Fig. 11
Oil-gas capillary pressure versus normalized oil saturation obtained from the centrifuge test and those resulting from history matching (HM) of the gas injection tests