Energy & Fuels 1996,9, 894-900
894
An Artificial Neural Network for the Prediction of Immiscible Flood Performance Ridha Gharbi," Mansour Karkoub, and Ali ElKamel College of Engineering and Petroleum, Kuwait University, P.O. Box 5969, Safat 13060, Kuwait Received May 4, 1995. Revised Manuscript Received July 5, 1995@
Despite several decades of artificial neural network reserach in other engineering disciplines, only recently work has been reported on its use as a prediction tool in petroleum engineering applications. Existing methods for the prediction of fluid flow in porous medium include numerical simulation techniques and laboratory core flood experiments. Both of these methods are generally expensive and time consuming. However, neural networks, once successfully trained, can be used t o predict reservoir performance in a short time with a personal computer. An artificial neural network was developed using data obtained from fine-mesh numerical simulation to predict the breakthrough oil recovery of immiscible displacement of oil by water in a two-dimensional vertical cross section. The network is able t o predict the results of the fine-mesh numerical simulations without actually performing these simulation runs. Various neural network connections were investigated using the back-propagation with momentum algorithm for error minimization. This paper describes the design, development, and testing of the neural network.
1. Introduction
the displacement behavior. Instead, adjustable parameters are used in a simple one-dimensional model to Most enhanced oil recovery (EOR) processes involve match the recovery performance of the displacement. the displacement of one fluid by another fluid in a One of the problems with the use of such models is that heterogeneous porous medium. The importance of EOR they are not based on the actual physics of flow and, as displacements in the petroleum industry makes its a result, the adjustable parameters often have no understanding and prediction critical in the decisions physical significance. on the applicability of certain recovery techniques. A second approach for reservoir performance involves Unfortunately, heterogeneity and the interaction of the use of numerical simulations. This approach has several forces, namely viscous, capillary, gravity, and enjoyed a resurgence of interest in recent years due to dispersive forces, can conspire t o make many EOR the availability of larger and faster vector computers displacement processes difficult to predict. than were available previously. The advantage of this Empirical models have been developed to predict the prediction approach is that it is based on sound physical recovery performance of fluid d i s p l a ~ e m e n t . ~ , ~ , ~ Jprinciples ~J~ of multiphase flow in porous media. HowThese models make no attempt to predict the details of ever, in order to obtain accurate predictions, the number of grid blocks used in the simulation must be large, thus * Author to whom correspondence should be addressed. increasing the computational requirements. @Abstractpublished in Advance ACS Abstracts, August 15, 1995. (1)Baldwin, J. L. Using a Simulated Bi-Directional Associative In addition to these predictive tools, laboratory core Neural Network Memory with Incomplete Prototype Memories to floods have been used in the petroleum industry as a Identify Facies from Intermittent Logging Data Acquired in a Silicismall-scale replica of petroleum reservoirs. Expericlastic Depositional Sequence. A case Study. Proc. SPE Annu. Tech. Conf. Exhibition 1991,273-286. ments are then performed on these replicas to predict (2)Briones, M. F.;Rojas, G . A.; Moreno, J. A,; Hidaigo, 0. Thermothe effectiveness of enhanced oil recovery processes dynamic Characterization of Volatile Hydrocarbon Reservoirs by before they are tried in the field. The disadvantage of Neural Networks. Proc. SPE Latin A m . Caribbean Pet. Eng. Conf. Buenos Aires, Argentina 1994,235-243. core flood experiments is that they are expensive and (3) Fayers, F. J. An Approximate Model with Physically Interprettime consuming. able Parameters for Representing Miscible Viscous Fingering SPE Res. Eng. 1988,551-558. The current prediction methods for reservoir perfor(4jFayers, F. J.; Newley, T. M. J. Detailed Validation of an mance are either accurate but time consuming, or fast Empirical Model for Fingering With Gravity Effects. SPE Res. Eng. but not necessarily accurate. Artificial neural networks 1988,3,542-550. (5) Garcia, G . ; Whitman, W. W. Inversion of a Lateral Log Using offer a cost effective and a reliable prediction tool. The Neural Networks. Proc. SPE Pet. Comput. Conf. 1992,295-304. objective of this paper is to develop an artificial neural (6) Habiballah, W. A.; Startzman, R. A,; Barrufet, M. A. Use of Neural Networks for the Prediction of Vapor-Liquid Equilibrium network for the prediction of the breakthrough oil K-Values. In-Situ 1993,17, 227-242. recovery in the displacement of oil by water in two(7)Juniardi, I. R.; Ershaghi, I. Complexities of Using Neural dimensional porous medium. Neural networks can Network in Well Test Analysis of Faulted Reservoirs. P m . SPEAnnu. Western Regional Meeting, Anchorage, AK 1993,711-721. effectively predict outputs from given inputs due to their
( 8 ) Koval, E.J . A Method for Predicting the Performance of Unstable Miscible Displacement in Heterogeneous Media. J.Pet. Technol. 1963, 145-154. (9) Kumoluyi, A. 0.; Daltaban, T. S. Higher-Order Neural Networks in Petroleum Engineering. P m . SPE Annu. Western Reghnal Meeting, Long Beach, CA 1994,555-570. (10)Lake, L. W. Enhanced Oil Recouery; Prentice Hall: Englewood Cliffs, NJ, 1989.
(11)Mohaghegh, S.; Arefi, R.; Ameri, S. Design and Development of an Artificial Neural Network for the Prediction of Formation Permeability. Proc. SPE Pet. Comput. Conf., Dallas, TX 1994. (12)Odeh, A. S.A Proposed Technique for Simulation of Viscous Fingering in One-Dimensional Immiscible Flow. SPE Res. Eng. 1989, 4 , 304-38.
0887-0624/95/2509-0894$09.00/00 1995 American Chemical Society
Energy & Fuels, Vol.9,No. 5, 1995 896
Prediction of Immiscible Flood Performance
Sigmoidal F d o n
7
1
0.9
_
I
......................... i........................... i..........
.j
.................). .................
.+ .........................
.:
.......................i........................
j ......................
:/
0.6 . . . . . . . . . .+.....................
0.5 ....................
Figure 1. An individual processing unit.
0 2 ......................
i
,..... ................................ i ~
L
.........................
j
......................
......................... i........................ 4 ........................... i ......................... inherent capability of capturing nonlinear functional behaviors. Neural networks are now being applied in -6 -4 -2 0 2 4 6 several areas of petroleum operations, such as the prediction of fluid properties,2,6well logging,1~5~11i16~19~20 Figure 2. Sigmoid function: f(A) = 141 + e-A). well t e ~ t i n g ,and ~ , ~drilling.18 Artificial neural networks are likely to become the preferred tool of the future for The processing elements within an artificial neural the design and performance prediction of enhanced oil network are arranged in a certain topology such as recovery processes. shown in Figure 3. The massive interconnections in the network give rise t o a large number of degrees of 2. Neural Networks freedom. This allows it to capture the nonlinearities in the systems’behavior, consequently, making it a good 2.1. Background. Neural networks are computing predictive tool. This is one of the reasons that artificial tools that consist of large number of simple, highly neural networks are superior to more typical regression interconnected processors called neurons or processing models. Neural networks do not depend on a single elements (PE). A neuron processes a pattern (input node as empirical models do. The neurons in a neural to give an vector) I with components 11,1 2 , 13, .... network operate in parallel and each node affects the output a (Figure 11, which can serve as input to other final output only slightly. A second reason is that neurons. Several factors, such as the weight vector w artificial neural networks have the ability to learn and and the bias B, determine the output a. The output of if back-propagation is used, the weights in the neural the PE is obtained by applying a transfer function on networks are periodically adjusted whenever a situation B). Several the total activation, TA, (TA = C W J i occurs for which the neural networks’ performance is functional forms can be used for the transfer function, inadequate. Finally, artificial neural networks can be which include linear, sigmoidal, and exponential. Howused to process multiple-input multiple-output data ever, it has been shown that sigmoidal (S-shaped) unlike other empirical modeling tools which can map functions are advantageous, especially in supervised one dependent variable only. neural nets; i.e., the input and the corresponding As shown in Figure 3, the artificial neural network outputs are known a ~ r i 0 r i . lA~typical sigmoid function consists of an input layer, an output layer, and if is given by necessarily one or several hidden layers. The input layer receives information from an external source and passes it to the hidden layers. These in turn quietly process this information and pass the results to the output layer. The hidden layer never communicates and is graphically shown in Figure 2. with the outside world. Its purpose is to serve as a (13) Pollard, J. F.; Broussard, M. R.; Garrison, D. B.; San, K. Y. repository for the information contained in the network. Process Identification Using Neural Networks. Comput. Chem. Eng. The procedure for picking the number of layers in a 1992,16,253-270. (14) Rumelhart, D. E.; Hinton, G. E.; Williams, R. J. Learning network is trial and error. First, only one layer is used Internal Representation by Error Propagation. Parallel Data Processand the number of neurons as well as the initial guess ing; The M.I.T. Press: Cambridge, MA, 1986 Vol. 1, Chapter 8, pp are varied. The network is repeatedly exposed to 318-362. (15) Shook, M.; Li, D.; Lake, W. L. Scaling Immiscible Flow Through training data until it Yearns” the correct input-output Permeable Media by Inspectional Analysis. In-Situ 1992,16,311-349. behavior. If the error between the desired and predicted (16)Smith, M.; Carmichael, N.; Reid, I.; Bruce, C. Lithofacies output is still unacceptable for large number of neurons, Determination from Wire-Line Loe Data Usine a Distributed Neural Network. Proc. IEEE Workshop xeural NetLorks Signal Process., the number of layers is increased by one and the same Princeton, NJ 1991, 482-292. procedure is repeated. The errors from the networks (17) Todd, M. R.; Longstaff, W. J. The Development, Testing, and Application of a Numerical Simulator for Predicting Miscible Flood using 1, 2 , .... and n layers are compared and the Performance. Trans. AIME 1972,253,874-882. connection that leads t o the lowest error is retained. It (18) Waller, M.D.; Rowsell, P. J. Intelligent Well Control. Trans. is worth noting that increasing the number of layers Inst. Mining Metall. 1994, 103, a47-a51. (19)Zhou, C. D.; Wu, X. L. Neural Network-Based Formation does not necessarily mean better results. Parameters Estimation from Well Logs in Quantitative Log Analysis: Training a network is an important factor for the A Comparative Study. Proc. Asia Pacific Oil Gas Conf., Singapore - . 1993,357-364. success of the neural network. There are many common (20) Zhou, C. D.; Wu, X. L.; Cheng, J. A. Determining Reservoir algorithms for training neural networks, and the method Properties in Reservoir Studies Using a Fuzzy Neural Network. Proc. we used is the back-propagation with momentum. SPE Annu. Tech. Conf. Exhibition, Houston, TX 1993, 141-150.
+
Gharbi et al.
896 Energy & Fuels, Vol. 9, No. 5, 1995
Back-propagation14minimizes a cost function such as the sum-squared error between the predicted and the desired values. In the minimization process, the weights of all the neurons' connections are adjusted until the error goal is achieved. The momentum procedure is used to help the search process not to get stuck in a local minima. A general methodology for the backpropagation method with momentum is now described. 2.2. Back-Propagation with Momentum. The first step in performing back-propagation is to calculate the output of the last layer of the network or output layer, n. In general, the output from neuron j in layer k is calculated by the following:
j = 1,2,
...,N k and k = 1, 2, ...,n
(2)
where for k = 0
N o = m and a i o = l i
(3)
The output of the output layer can be obtained by setting k = n in eq 2. The sum-squared error, sse, of the network is given by
cc P Nn
sse =
(dj -
(4)
p=lj=l
=(mn-k)p
where q is the learning rate and y is the momentum constant. The weights can be updated as follows:
The above procedure is repeated until the sum squared error is within the desired limits. In order to improve the search process, an adaptive learning procedure is used in the following fashion: if the new sse is greater than the previous one by more than E , then the new weights are rejected and recalculated using a new learning rate (see Figure 4). The new learning rate is given by
E and qo are two values that can be fxed at the beginning of the computations. Typically, E is chosen to be 5% and qo is 1/&.
3. Process Modeling
where d j and aj, are the desired and actual outputs from neuron j in the output layer n, respectively. The subscriptp represents a specific input pattern. The next step is to calculate the gradient descent of the output layer. The gradient descent of layer k is given by (dn-k)p
one by one going backwards until the input layer is reached. When all the gradient descents for all the layers are calculated, the adjustments for the weights are then calculated:
( w T - k ) p (dn-k)p
lz = 1,2,
...,n - 1 ( 5 )
The superscript T represents the transpose of the matrix. The neuron gradient matrix for an input pattern p is given by
(6) and the corresponding weight matrix is given by
In order for the neural network to be successfully trained to predict the breakthrough oil recovery, enough input data is required. Due to the lack of quantitative experimental data, the input and the output of the finemesh numerical simulations are used to train the network. The ultimate goal is to show that neural networks can learn the functional relationship between the inputs and the outputs of the simulator. This is demonstrated by considering an immiscible displacement of oil by water in a two-dimensional,vertical cross section. The porous medium is assumed homogeneous and anisotropic (diagonal permeability tensor with kx * k,) with constant porosity and dip angle (a)measured from the horizontal. The fluids and the pore space are assumed incompressible. The x-axis is the principal flow direction, while the z-axisis the direction transverse to the flow. Figure 5 shows the geometrical setting for the displacement. Under these assumptions, the displacement can be modeled by the following set of equations:
as,
au,
*+=+==o
au,
(13)
Note that the third subscript of all the connection weights, k , in eq 7, is removed to avoid cumbersome subscripts. The values of the gradient descent for the output layer is given by
where
krorpo
u, = -kz-
Po
+
]
- ,ogsin(a) ax
(16)
Prediction of Immiscible Flood Performance Pwz
+
= -k,h[%
cos(a)]
Energy &Fuels, Vol. 9, No. 5, 1995 897
(17)
PW
U,
= -k,-
krorpo +
Po
]
- eg cos(a)
(18)
[y
(19)
+0-
Po= Pw
j(S,)
Constant injection rate, constant production pressure, and no-flow boundary conditions a t z = 0 and z = H are considered. In addition, the porous medium initially has an initial water saturation, S w i and a residual oil saturation, Sor. Shook et al.15presented a generalized form of inspectional analysis which derives the minimum number of linearly independent dimensionless scaling groups for this displacement. The authors found that there are five dimensionless scaling groups that describe the immiscible displacement of oil by water in two dimensions:
(20)
strongly at low displacement rates, when the viscous and pressure forces are low. The magnitude of the gravity number determines the degree of gravity instability, whereas the sign determines whether the instability is gravity override or underride. If the gravity number is negative, the problem is override whereas if the gravity number is positive, the problem is underride. The gravity number is an important scaling group which should always be considered because immiscible displacement in the laboratory or in the reservoir occurs in the earth’s gravitational field. The scaling group, RL, is the effective length-thickness aspect ratio, which controls the amount of vertical mixing. The effective aspect ratio has long been recognized as the sole controlling factor governing the approach to vertical equilibrium.1° A reservoir is said to be in vertical equilibrium when the transverse driving forces are zero. In other words, vertical equilibrium is a state of maximum transverse communications. The vertical equilibrium condition is obtained when the effective aspect ratio is large. The scaling group, Na, is the dip angle number, which is purely geometric. The number has been shown to play a sigdicant role in the effect of the gravity forces.15 4. Discussion of Results
(21)
(22)
--[-]L
k,
RL - H k , L
112
Na = - tan(a) H
(23) (24)
where A@ = ew - eo. Equations 20-24 show the dimensionless scaling groups that must be the same in the laboratory model and in the field in order to scale immiscible displacements in porous media. For any geometrically similar porous medium, matching these dimensionless scaling groups and the petrophysical properties will yield identical dimensionless results. The scaling group, W , the end-point mobility ratio, has long been recognized as controlling the areal sweep efficiency of an immiscible displacement. When this group is less than unity, the displacement is said to be favorable and the areal sweep efficiency will be high. When this group is greater than unity, the displacement is said to be unfavorable and the areal sweep efficiency will be low. The scaling group, Npc,the ratio of capillary to viscous forces, is a measure of the influence of surface or capillary forces on the displacement. This influence is strongest at low displacement rates in which the surface forces promote the displacement efficiency of the nonwetting phase by the wetting phase. The gravity number, Ng,is the ratio of gravitational to viscous forces in the displacement. Like the capillary force, the gravity force is independent of the displacement rate and will therefore manifest itself most
In order to generate date t o train and test the neural network, the five dimensionless scaling groups are systematically varied and the corresponding breakthrough oil recovery is obtained from fine-mesh simulations of the displacement. The purpose is not to establish when a particular displacement is gravity, capillary pressure, or viscous force dominated; rather it is to generate sufficient data points to successfully train the neural network. The simulations were done using a two-dimensional simulator which solves eqs 1319 using an explicit finite-difference approximation. A 300 x 300 grid was used for all simulation runs. This grid was found to be accurate enough and a further increase in the number of grid blocks yielded the same simulation results. A limited number of simulation runs were performed using dimensionless scaling groups with ranges listed in Table 1. Some of these data points were used to train the networks, while the others were used to test the effectiveness of the training process. Several neural network architectures are used to predict the breakthrough oil recovery. A one-hidden-layernetwork (i.e., a sigmoid layer and a linear output layer) was initially utilized. The number of neurons in the hidden layer was systematically varied to minimize the sse, defined in eq 4, using 50 000 epochs. Table 2 shows that the sse decreases as the number of neurons increases. However, even for a large number of neurons (greater than 501, the sse is till not acceptable (sse > Therefore, an additional sigmoid hidden layer is added to the network and the number of neurons in each layer is varied until an acceptable sse is reached (see Table 2). Several architectures with two hidden layers yielded acceptable sse. However, most of these designs were rejected after the cross-validation process was carried out. The cross-validationprocess involves two steps. First, a new sse, based on data points not used in the training process, is calculated and only the networks that led to a new sse less than low3are retained. Consequently,
898 Energy & Fuels, Vol. 9, No. 5, 1995
.
Gharbi et al.
n 0
U
n 0
U
I -u
\
I
1st Hidden Layer
~~
2nd Hidden Layer
Output Layer
h p u i Layer Figure 3. Multilayer neural network diagram.
Table 1. Ranges of Training Data scaling group Ng NCl
range 0.01-10.0 0.0 -2 0.0 0.10-12.0 0.50-3.0 0.0-10.0
RL Mo NP,
Table 2. Results of the Training Process for Various
Neural Network Architectures
no. of layers 1 1 1 1 1 2 2 2 2 2 2 2
4 8 12 16 20 50 5 7 7 7 7 8 20
-
2
2
1
YES
Figure 4. Flow chart for adaptive learning rate.
Figure 5. Displacement of oil by water in a two-dimensional porous medium. the two-hidden-layer network with 20 neurons in the first layer and 6 in the second was rejected. Second, it is made sure that none of the retained architectures lead to overfitting. The overfitting problem in these networks is checked using the results obtained from in-
no. of neurons first second layer layer
sse
new sse
5 4 5 6
-
7 8 2
0.56328 0.41979 0.22086 0.25090 0.20581 0.18345 0.04308 0.01187 0.00704 0.00521 0.00214 0.00098 0.00041
20
4
0.00026
0.0005
20
6
0.00095
0.004
-
-
0.0002 0.0004
termediate data points. If oscillation occurs, then the architecture is rejected. The two-hidden-layer network with 20 neurons in the first layer and 4 in the second, for example, is rejected because of this problem as clearly illustrated in Figure 6. All other retained architectures were also rejected except the two-hiddenlayer network with 8 neurons in the first layer and 8 in the second layer. This network architecture led to an excellent fit as shown in Figure 7. Therefore, this network was chosen for the prediction of the breakthrough oil recovery.
Energy & Fuels, Vol. 9, No. 5, 1995 899
Prediction of Immiscible Flood Performance 0.45,
0
2
1
.
.
. .
. . .,
I
/
*,
0.35.
t,
0.3.
+.,
1
0,251
1
0 05'
21 0015
.
. . __....
.
lo-'
10,
1 OD
__
, 1oo
.
i
;;.A.j 10'
Ng
Ng
Figure 6. Neural network estimated recovery using two hidden layers with 20 neurons in the first layer and 4 in the second layer: actual (solid), estimated (dotted).
Figure 9. Breakthrough oil recovery as a function of Ng a t various RL for MO = 0.5, N , = 0.0, Np, = 0.0; actual (pluses), estimated for RL = 0.1 (solid),estimated for RL= 1.0 (dashed), estimated for RL = 5.0 (dotted).
a 05
103
10:
w Figure 7. Neural network estimated recovery compared to fine-mesh simulation using two hidden layers with 8 neurons in the first layer and 8 in the second: actual (solid),estimated (dotted). 1
1
0.08' 10"
10.'
1OD
Figure 10. Breakthrough oil recovery as a function of Ng a t various RL for MO = 2.0, N , = 0.0,Np, = 0.0; actual (pluses), estimated for RL = 0.1 (solid),estimated for RL= 1.0 (dashed), estimated for RL = 5.0 (dotted). 0.5
0.051 10''
Figure 8. Breakthrough oil recovery as a function of Ng a t various MO for RL = 8.0, N , = 0.0, Np, = 0.0; actual (pluses), estimated for MO = 0.5 (solid),estimated for MO = 1.0 (dashed), estimated for MO = 2.0 (dotted).
Figures 8-12 show a comparison of the results obtained from the fine-mesh simulation with those obtained from the neural network model. The agreement between the two methods is good over the entire range of the simulated data. The fit is further tested on data not used in the training process as illustrated
10'
Ng
1on
10" Ng
Figure 11. Breakthrough oil recovery as a function of Ng at various N , for RL = 10.0, MO = 2.0, Np, = 0.0;actual (pluses), estimated for N , = 20 (solid),estimated for N , = 1.2 (dashed), estimated for N , = 0.3 (dotted), estimated for N , = 0.0 (dashdotted).
in Figures 13 and 14. The agreement is still good even with these data. This shows that, for this type of displacement, the artificial neural network has been trained successfully and that the model has good
900 Energy & Fuels, Vol. 9, No. 5, 1995
Gharbi et al.
of dimensionless scaling groups, thus eliminating the need for the expensive fine-mesh simulations. 5. Concluding Remarks
0
4
An artificial neural network has been developed to predict the breakthrough oil recovery of immiscible displacement from data obtained using two-dimensional fine-mesh numerical simulations. Various network architectures were investigated using a back-propagation with momentum algorithm for error minimization. For the case studied, the most effective network architecture consist of two-sigmoid hidden layers and a linear output layer with 8 neurons in the first hidden layer and 8 in the second. This latter network was then used to predict the breakthrough oil recovery for new situations which were not included in the training phase. These predictions compared well with the results of the fine-mesh simulations. Although the reservoir model we used was a relatively simple one (Le., does not consider heterogeneity), the main purpose of this work was to demonstrate the capability of artificial neural networks to learn the functional relationship between the input and the output of fluid displacements in porous media. This work establishes a foundation for applying neural networks in the prediction of oil recovery in heterogeneous reservoirs, which will be the subject of future study.
1
Nomenclature aJk = output of neuron j in layer k Fk = transfer function of the neurons in layer k
01 104
io-'
1 o0
10'
Ng
Figure 13. Predicted breakthrough oil recovery as a function of Ng at various RL for MO = 1.0, N , = 0.0, Np, = 0.0; actual (pluses), estimated for RL = 0.10 (solid), estimated for RL = 0.3 (dashed) estimated for RL = 1.0 (dotted), estimated for RL = 10.0 (dashed-dotted). 0 45,
a
8
I
-
0.3
0.25 -
I
"'1
0.15l lo-*
1on
lo-'
*. 10'
Ng
Figure 14. Predicted breakthrough oil recovery as a function of Ng at various RL for MO = 0.5, N , = 0.0, Np, = 0.0; actual (pluses), estimated for RL = 1.0 (solid), estimated for RL = 3.5 (dashed).
interpolative properties. Therefore, it can be used to predict the breakthrough oil recovery from other sets
g = gravitational acceleration H = system height Z = input pattern J(S,) = Leverett's J-function k,, k, = permeability in the x and z direction, respectively k, k,, = relative permeability to water and to oil, respectively L = system length MO = end-point mobility ratio n = total number of layers Nk = number of neurons in layer k N , = dip angle number Np, = capillary number Ng = gravity number P = number of input patterns Po = pressure in the oil phase P, = pressure in the water phase RL = effective aspect ratio So, = residual oil saturation S, = water saturation S , = connate water saturation t = time uWx,uw2= water velocity in x and z direction, respectively uar, u m = oil velocity in x and z direction, respectively U T = total system velocity Wyk = weight for input i, neuronj, and layer k x , z = coordinates a = dip angle measured fron the horizontal /?Jk = bias for neuronj in layer k = learning rate y = momentum constant = end-point water and oil relative mobility, respec,:A tively p,, po = water and oil viscosity, respectively C$ = porosity ew,eo = water and oil density, respectively u = interfacial tension EF950087N
