A Simulation of Core Displacement Experiments for the Determination of the Relative Permeability

ABCTRACT Computations of the relative permeability curves were made through their representation by two functions for wetting and nonwetting phases. Each function contains one parameter that controls the shape of the relative permeability curves. The values of these parameters are chosen to minimize an objective function, that is represented as a weighted sum of the squared differences between experimentally measured data and the corresponding data calculated by a mathematical model simulating the experiment. These data comprise the pressure drop across core samples and the recovery response of the displacing phase. Two mathematical models are constructed in this study to simulate incompressible, one-dimensional, two-phase flow. The first model describes the imbibition process and the other describes the drainage process. The values of the relative permeability parameters are calculated by employing Rosenbrock optimization procedure. The reliability of this procedure has been confirmed by applying it to four displacement cases. The optimum values of the relative permeability parameters, which reflect the final shape of the relative permeability curves, are achieved at the minimum value of the objective function. All the above processes are be embodied in relative permeability package RPP which is constructed in this study using FORTRAN language.


INTRODUCTION
Relative permeability curves are important for many reservoir engineering calculations.Frequently, the relative permeability curves are constructed on the basis of laboratory displacement experiments.They are not measured directly, but instead are inferred from flow data measured during the displacement experiments through the use of some interpretive methods.The dynamic fluiddisplacement experiments data of core samples are often used in the calculation of the relative permeability curves.called implicit interpretive approach, was developed to overstep the limitations that are associated with the explicit approach .
In the implicit methods, the relative permeability curves are computed by representing them by two functions, each of which contains certain coefficient that controls the shape accuracy of the relative permeability curves.Relative permeability curves are adjusted until the values computed with a mathematical simulation of the laboratory experiment.Match in same sense the measured data (Kerig and Watson,1987, Sigmund and McCafery,1997.Archer and Wong (1973) used a reservoir simulator to model laboratory tests.The relative permeability curves are adjusted by trail-and-error procedure until calculated oil recovery and relative injectivity curves match those obtained from the laboratory displacement tests.Sigmund and McCaffery (1979) suggested that the determination of relative permeability curves be executed by representing it with two functions, each of which contains one coefficient.The values of these coefficients are calculated for different rocks by employing nonlinear-least-squares optimization procedure.
In this study, estimation of relative permeability curves from two-phase displacement experiments are considered.This is achieved through preparing a package RPP by which relative permeability curves may be determined.Employing numerical techniques solves the mathematical models, which are used to develop the package.The implicit method is used for this purpose.

MATHEMATICAL FORMULATION
The relative permeability curves are to be computed by assuming a function of two parameters, then use is made of optimization technique to calculate the optimum value that fulfill a preassigned conditions.To estimate relative permeabilities with a parameter estimation approach, an objective function is constructed as a weighted sum of squared differences between the measured and the calculated data from a mathematical model of coreflood experiment.The objective function may be determined from : where, J : is the objective function to be minimized; F → (β →) : is the (NPx 1) vector of the unknown parameters; Y → : is the (N x 1) vector of the observed (measured ) data; W : a (NxN) weighting matrix; N : total number of data point; NP: number of parameters; and T : transpose operation.
The above form of the objective function is adopted by numerous workers (Batycky et al.,1981,Tao andWatson,1984).For a typical displacement experiment, the measured data consists of the pressure drop a cross the core sample, ∆P obs , and the cumulative volume of displacing phase recovery response, E R obs Two mathematical models have been considered; one for imbibition process and the other for drainage process.The quantities, ∆P cal ,and E R cal may be obtained by solving these models.

FUNCTIONAL REPRESENTATION OF THE RELATIVE PERMEABILITY CURVES
The relative permeability curves for both wetting and nonwetting phases may be written as follow: in which; where, a1, b1 are constants .The values of the constants a1, and b1 were suggested to be 0.01 for computational purpose to linearize the relative permeability

[
] [ ] ...( 1) curves as these curves approach zero; but which otherwise do not influence the shapes of curves.The functional relationships between relative permeabilities and saturation defined by Eqs.2 with their two adjustable saturation exponents, and the values used for the constants a1 and b1 appeared to describe adequately the flow characteristics of the porous media that were studied Sigmund and McCaffery (1979).
Capillary pressure effects are expressed in terms of the following equation 3 : in which, and Theoreatically, λ could have any positive value greater than zero (Brooks and Corey,1964); being small for media having a wide range of pore size and large for media with a relatively uniform pore size.
In the current project and based on Brooks and Corey study it is found that the values of the parameter range from 1.0 to 7.5 for different porous media.In this context an alternative approximation for the value of λ is made, that depends upon the degree of homogeneity for the estimation of λ.The range of the parameter λ may be subdivided into four categories: A-Very homogeneous porous media, λ ≥ 7.5 B-Homogeneous porous media, 7.5 >λ ≥ 5.0 C-Non -homogeneous porous media, 5.0 >λ ≥ 3.0 D-Highly non homogeneous porous media, 3.0 >λ ≥1.7 Application of these approximations is found to give good results The end points relative permeabilities k rwe.and k rnwe.are calculated using Darcy's Law as follow: 1-In the case of imbibition.

IMBIBITION AND DRAINAGE MODELS
The flow equations for unsteady state displacement of incompressible, onedimensional, two phases fluid flow (which include capillary pressure and ignore gravity effects) may be expressed as (Azaz and Settari,1979,Nolen and Berry,1979,Peaceman,1977). The derivation of above equations was given by Dake (1978).To solve these equations, the finite difference form of these equations has been used in the imbibition and drainage simulator.The methods of solution were used for treating the nonlinear terms of Eqs. 6, 7, and 8 are: 1-One -point upstream transmissibility weighting; 2-Fully implicit transmissibilities using chord slope method to estimate derivatives; 3-Newton's-Raphson's method to handle nonlinearities resulting from the use of capillary pressure.The Boundary and initial conditions for Eqs.6 through 8 have been discussed by Aziz and Settari, and Sigmuned and McCaffery among others.The initial conditions are a uniform distribution of saturation S winit , and pressure P winit .The inlet flow condition at X=0 is constant flow rate for the wetting or nonwetting phases ( Q wi or Q nwi ) for t>0.The outlet boundary condition at X=L must follow Darcy's law, and at the same time allows the pressure for both fluids in the porous medium to be continuous with the pressure just outside the outlet.In the case of imbibition, Aziz and Settari 1 stated these conditions mathematically at the outlet as: and where, Qwo, Qnwo are the flow rate of the produced wetting and nonwetting phases respectively.Where the outlet breakthrough of the wetting phase occurs when the wetting phase saturation increases from its initial value S winit , to the value S wo corresponding to zero capillary pressure.In the case of drainage, S wo remains constant until the value (1-S winit ) is reached at the outlet face, which represent the maximum nonwetting phase saturation.Thus, the outlet condition is: The above conditions embodied the frontal advance of fluid-flow concept, which was suggested by Buckly and Leverett (1942) The curved sides of the simulated core sample are assigned no flow boundaries, as actually is treated in the simulator models of the laboratory tests.

EVALUATION OF THE RELATIVE PERMEABILITY PACKAGE(RPP)
The dimensionless cumulative injection Qi, the recovery response E R , and the dimensionless pressure response ∆P D for both the imbibition and drainage cases are required for comparison purpose.These may be calculated at the end of each time step from the following equations: Recovery response is define for imbibition as: and for drainage as: The dimensionless pressure response for both imbibition and drainage are calculated from the following equations: Four runs were made on example problems which were presented by Sigmund andMcCafgfery (1979), andBatycky et al. (1981), to confirm the reliability of the Resenbrock optimization procedure for obtaining the parameters Ew and Enw rather than the leastsquares fitting procedure, which was adopted by a number of workers in this field.(Watson et al.,1988,Yong andWatson,1991).
In all the four runs bad initial estimates were used for the parameters Ew and Enw to confirm the efficiency of the RPP in computing the optimum values of these parameters.
In the subsequent sections the fourdisplacement case are discussed.The first case is the imbibition displacement for Swan Hills core, (limestone core type).The second is the drainage displacement for Swan Hills core.The third case is the imbibition displacement for Rainbow core, (dolomite core type); and the fourth case is the imbibition displacement for core sample of Carig, 1971.

CASE1-IMBIBITION DISPLACEMENT FOR SWAN HILLS CORE
Initial estimates of the parameters, Ew and Enw and initial step size were introduced to the RPP, the RPP will give the optimum values of these parameters.The preceding descriptions are summarized in Table 2.
Core flood parameters for this case are given in Table 1.
Fig. 1 shows the observed dimensionless pressure drop ∆P D , and recovery response E R , data obtained from the experimental work, as functions of the dimensionless cumulative injection, Qi compared with simulator -calculated values which are obtained by using the Rosenbrock optimization procedure which is listed in Table 2 .The simulated values of Ew and Enw determined from the values are seen to match the experimental results remarkably excellent.The two parameters relative permeability curves characterized by the values of Ew and Enw which are given in Table 2 is shown in Fig. 2, for both the imbibition and the drainage displacement (case 2).The optimization constrains used in this case for the parameters Ew and Enw in the RPP are [0.1 ≤ Ew ≤8] and [1≤ Enw ≤10.5].

CASE 2-DRAINAGE DISPLACEMENT FOR SWAN HILLS CORE
The coreflood parameters necessary for the calculation are given in Table1.Fig. 3 shows the comparison between the observed dimensionless pressure drop and recovery response data with those calculated by the simulator approach at the optimum values of Ew and Enw were shown in Table2.The simulated values are seen to match the experimental data reasonably well, with exception of the pressure at the end of the flood.These differences in the latter situation may be attributed to calculate the end -point relative permeability to the nonwetting pahase.The relative permeability curves as a function of wetting phase saturation is shown in Fig. 2, where the comparison between drainage and imbibition values (case1), is made.The wide divergence between the relative permeability curves makes it mandatory to choose the correct one for predicting flooding behavior.Because the wetting and nonwetting phases relative permeability is history dependent, wetting phase flood performance should be predicted more accurately by unsteady-state.The optimization constrained used in this case for the parameters Ew and Enw in the RPP are [0.3 ≤ Ew≤8] and [1 ≤Enw ≤10.5], respectively.

CASE 3-IMBIBITION DISPLACEMENT FOR RAINBOW CORE
In this case the analysis of transient response curves for displacement of nonwetting by wetting phase in a heterogeneous, water-wet, Rainbow core is considered Core-flood parameters for this case are given in Table1.The observed values are seen to be less than simulated values before breakthrough, this effect is due to the high degree of heterogeneity.Therefore, the application of the Johnson -Bossler -Naumann method to the data plotted in Fig. 4 is impossible because of the high nonlinearity before breakthrough .
The concave shape presented for the wetting phase relative permeability in Fig. 5 is a consequence of having only a single curvature presenting the relative permeability curves in Eqs. 2. This could be overcome by using a more general relationship between saturation and relative permeability or by employing a lower value of Swr in the simulator such as would be achieved experimentally with a highrate drainage displacement.In This case, the core properties data which were presented in Batycky et al. (1981) are considered.By employing RPP conformable results of relative permeability curves have been obtained.Fig. 6 illustrates the response curves for Ew = 1.525 and Enw = 7.759, the other core -flood parameters are given in Table1.The response curves of the results obtained from RPP are seen to match the experimental data reasonably well with the exception of the pressure response before breakthrough.This problem was treated in RPP by using weighting factor.The choice of this weighting scheme partly reflects the spacing of data.Fig. 7 shows the relative permeability curves of this case.The initial guesses of the parameters Ew and Enw, and initial step sizes were presented in Table 2.The constrained used in this case for the parameters Ew and Enw in the RPP are [0.1≤Ew≤15] and [0.1 ≤Enw ≤15].

CONCLUSIONS
The following conclusions are drawn from this study 1-The RPP can be adopted in the two-phase relative permeability curves computations, for the unsteady-state data, with high reliability.2-The parameter estimation approach, which adopted in this work, overcomes significant limitations of the classic calculation procedure of the Johnson -Bossler -Naumann method and related methods.3-Rosenbrock optimization procedure that has been utilized in this study is more efficient than the Gauss-Newton.Nonlinear least squares technique that was used by Sigmund and Mc Caffery.4-The application of the RPP give an occasion to calculate relative permeability curves for heterogeneous cores rather than the explicit methods, which applied for homogeneous cores only.NOMENCLATURE A = cross sectional area of the core sample, L 2 E R = recovery efficiency, percent of movable fluid recovery, dimensionless E nw = parameter in nonwetting phase relative permeability expression, dimensionless E w = parameter in wetting phase relative permeability Expression, dimensionless k = absolute permeability, L 2 k r = relative permeability, dimensionless k rwe = end-point relative permeability of wetting phase, dimensionless k rnwe =end-point relative permeability of nonwetting phase, dimensionless L = length of the core sample, ft.P cb = capillary pressure scaling coefficient, m/L 2 P c = capillary pressure, m/L 2 ∆ P D = dimensionless pressure drop across core ∆P= transient pressure drop across core, m/L 2 P nw = pressure in nonwetting phase, m/L 2 P w = pressure in wetting phase, m/L 2 Q = flow rate, L 3 /t q = flow rate, 1 / t Qi = dimensionless cumulative injection R D = dimensionless measure of viscous-tocapillary force ratio.Se = saturation, normalized with respect to Sw min and Sw max.S pc = stauration, dimensionless normalized with respect to S wo and S wr S w = wetting -phase saturation, dimensioless.
Fig 4 shows a comparison between the computed dimension pressure drop and recovery response with the observed results using the values of Ew and Enw which are given in Table2; as determined by the Rosenbrock procedure.Initial estimates of Ew= 2.0 and Enw =3.5 which are bad relative to the initial estimates that were given by Sigmund and McCaffery (Ew = 0.5, Enw = 3.0) for least -squares procedure, were deliberately chosen to confirm the efficiency and superiority of the developed optimization procedure than that development by Sigmund and Mc Caffery.The last procedure requires several preliminary simulator runs to make good initial guesses of the parameters Ew and Enw, which are necessary for the Gauss -Newton nonlinear least-squares that was a dopted by Sigmund and Mc Caffery.The relative permeability curves characterized by the Ew and Enw given in Table2 are shown in Fig 5.
S nw = nonwetting -phase saturation, dimensionless S wo = wetting -phase saturation corresponding to zero capillary pressure, dimensionless S winit =initial wetting-phase saturation, dimensionless S wmin = is the minimum wetting phase saturation established by a drainage displacement, dimensionless S wr =irreducible wetting phase saturation, dimensionless S wmax = is the maximum wetting phase saturation established by an imbibition displacement, dimensionless S w,av =average wetting phase saturation, dimensionless t = time, t ∆t =time increment, t x= distance, L ∆x = distance increment, L φ = porosity, dimensionless λ = measure of pore-size distribution, dimensionless µ = viscosity, m/Lt σ = interfacial tension, m/L 2 SUBSCRIPTS nw =nonwetting phase m = wetting or nonwetting phase w = wetting phase SUPERSCRIPTS cal.=calculated data value obs.=observed or measured data value → = vector notation