STRUCTURE OF SATURATION JUMPS IN NONEQUILIBRIUM DISPLACEMENT FROM POROUS MEDIA O. B. Bocharov, O. V. Vitovskii, and V. V. Kuznetsov

UDC 532.546

The nonequilibrium displacement characteristics are investigated experimentally and by numerical calculation of the nonequilibrium percolation equations under various conditions of wetting of the porous medium by the displacing and displaced fluids. A scheme for calculating the disequilibrium parameter is proposed. The classical theory of the equilibrium displacement of oil by water is based on the relative phase permeability functions ki(s) determined experimentally from steady displacement tests and the Leverett function J(s) characterizing the capillary pressure jump at the interphase menisci. The unsteady displacement process is characterized by the presence of zones of sharp change in the water saturation s in the porous medium (water saturation jumps), in which the Muskat--Leverett equilibrium model of the two-phase flow of immiscible fluids may break down [i, 2]. In the simplest scheme for taking disequilibrium into account [I, 2] it is assumed that the functional parameters J(s) and ki(s) are the same in the nonequilibrium as in the equilibrium flow, but depend not on the true water saturation s but on the fictitious saturation ~(s). In [i, 2] it is proposed that ~ be determined from the kinetic equation 65

where t ~ is time, and ~ is the characteristic time of fluid substitution in the porous medium. The solutions of the nonequilibrium percolation equations were investigated in [2] for small r/tl, where t I is the displacement time, using the method of matched asymptotic expansions and in [3] numerically for intermediate ~/t I < i. It was shown that in this case the length of the stabilized zone increases with increase in the displacement rate, as confirmed by the experimental data of [4]. At the same time, we still lack experimental data on the structure of the saturation jumps for intermediate values of ~/t I < 1 and, moreover, we have no means of determining the characteristic substitution time 9 for nonequilibrium displacement. i. Scheme for Calculatin ~ the Disequilibrium Parameter The Muskat--Leverett equilibrium model of two-phase flow through porous media is based on the assumption that the relative phase permeabilities ki(s I) and the Leverett function J(sl) are determined only by the saturation of the pore channel by the displacing fluid s I. In the case of steady two-phase flow experiments to determine the ki(s l) for various porous media have shown that for capillary numbers Ncl less than 10 -4 they are in fact mainly determined by the structure of the pore channel and the fluid saturations [4, 5]. The capillary number Ncl = v0~i/o is defined in terms of the overall filter velocity v0, the viscosity of the displacing fluid Ul, and the interfacial tension o. In this case the capillary forces are large as compared with the hydrodynamic pressure difference on scales of the order of the pore length and determine the part of the pore channel which can be occupied by a fluid wetting the porous medium, water, for example, with saturation s I. In unsteady displacement for a given value of s I in a volume element comprising many pores an equilibrium two-phase flow structure is formed only at a low rate of change of saturation 8sl/St ~ In this case the filling of the pores takes place under equilibrium conditions. At first, the small channels are occupied, together with the regions Novosibirsk. Translated from izvestiya Akademii Nauk SSSR, Mekhanika Zhidkosti i Gaza, No. 6, pp. 97-104; November-December, 1990. Original article submitted July 3, 1989.

0015-4628/90/2506-0891512.50 9 1991 Plenum Publishing Corporation

891

on the walls of the pores with small average radii of curvature of the interphase menisci r !. The resulting capillary pressure P2 -- Pl = ~ prevents the larger channels from being filled until the hydrodynamic pressure differentials on scales of the order of the length of the displacement zone lead to an increase in the phase pressure difference to P2 -- Pl > o/rl, which makes possible the occupation of the larger channels. Since any pore space is formed of a system of channels, the latter are occupied in steps, and the step interval t I is determined by the capillary pressure PG (sl) = o/ s I) and by the hydraulic resistance of the channel for the fluids k~(sl). Using Darcy's law to write the equations of motion of the fluids along the capiilary under the influence of capillary pressure only, we can determine the time, equal to the step interval, taken by the displacing fluid to fill the channel:

ti

=

"

2 o(ffmj(s,)k,*(s,)

( 1+~o

"

where ~0 = ~i/~2, ~i is the viscosity of the displacing fluid (water), ~2 is that of the displaced fluid, k~(s I) are the relative permeabilities of the occupied channel, is the length of the pore channel, and k is the absolute permeability of the porous medium. For a given value of s~ the permeabilities of the pore channels may differ substantially from the relative phase permeabilities ki(sl), since to a considerable extent the latter are determined by the number of pore channels through which the fluids pass. In particular, this is important near the limiting values of the saturations. For a high rate of change of the water saturation in the volume element even at low capillary numbers Ncl < i0 -4 the stepwise equilibrium occupation of the pore channels cannot ensure the necessary value of 8sl/St ~ determined by the external conditions. This leads to a local increase in the pressure differences between the fluids and to the simultaneous opening in one step of several pore channels with similar breakdown capillary pressures. In this case for a given value of s I the relative phase permeabilities for the displacing fluid increase as compared with equilibrium displacement, while those for the displaced fluid decrease and will correspond to a value of the water saturation al greater than the real value s I. The effective water saturation ~i is determined by s I and the number of pores n s simultaneously opened in one step in unit volume of the porous system: a I = a1(sz, ns). For small n s near the equilibrium state

(1.2 s

The number of pores n s simultaneously opened in a single step is related to the rate of change of water saturation per unit volume by the equation

~ Ot ~

~'g(s'____!)

(1.3)

tt

where t I is the step interval determined from (I.i), and $(s l) is the fraction of the pore channel occupied by water in one step (determined by the equilibrium shape of the interphase menisci in the pore channel for a given capillary pressure). Determining n s from (1.3) using (i.i), we obtain the expression for the effective saturation a~=s~+

O~t[ I~2 Ost On, ,=o2oYkmlF*k~*~(sO Ot~

(1.4)

where F* = (1 + u 0 k ~ / k ~ ) i s t h e B u c k l e y - - L e v e r e t t f u n c t i o n f o r t h e r e l a t i v e p e r m e a b i l i t i e s of the pore channels. I n t h e c a s e o f g r a n u l a r p o r o u s m e d i a s u c h as q u a r t z sand o r s a n d s t o n e w i t h a low c e m e n t c o n t e n t t h e l e n g t h o f t h e i n d i v i d u a l p o r e ~ i s c l o s e t o t h e a v e r a g e g r a i n s i z e d. U s i n g t h e w e l l - k n o w n Kozeny--K~rmSn e x p r e s s i o n f o r t h e p e r m e a b i l i t y o f t h e s e p o r o u s m e d i a k = d2m2/180(1 -- m) 2 and n o n d i m e n s i o n a l i z i n g ( 1 . 4 ) w i t h r e s p e c t to the characteristic d i s p l a c e m e n t t i m e t e = mLSS/v0, we o b t a i n 038

a=s+Dtl ( s ) - -

892

( 1.5 )

D=

90(t--m) 2~. Yk

-~--max[,

~

m~m

R(s)= F*~*~On: 31--81,i Si,l--Si,i '

L

~=o

[

Og I ] J F * k 2 * ~On~ .~=o

max~ ~t--81 i 31j--S1,i

1

-

~=o

On:

Yo~2 G

where D is the disequilibrium parameter, R(s) is the disequilibrium function normalized on unity, s is the normalized water saturation, a is the normalized effective saturation, and s~, i and sl, f are the initial and final limiting water saturations. In view of the complexity of the structure of the pore space relations (1.6) can only be estimates and, in reality, the disequilibrium function is an empirical functional parameter like the relative phase permeabilities and the Leverett function. At the same time, it follows from (1.6) that the disequilibrium function is determined by the form of the functions J(s), F*(s), k~(s), ~(s), and 8a/Sn s, which depend on the conditions of wetting of the porous medium by the fluids. We will assume that the function 8a/On s depends only weakly on the water saturation over the entire range of s and is close to zero near the limiting water saturations. This makes it possible to estimate the form of the function R(s) for various conditions of wetting of the porous medium by the displacing and displaced fluids. Figures la and ib show the equilibrium phase permeabilities in quartz sand with k = 6"10 -12 m 2 for wetting (i, 2) and nonwetting (3, 4) fluids (Fig. la) and the Leverett functions (Fig. ib, i, 2) for imbibition with advancing contact angles 8 equal to -20 ~ and -75 ~ , respectively. The contact angles were measured on a rough plate with roughness similar to that of quartz sand. As the fluids we used water and a liquid hydrocarbon -transformer oil. The contact angle was varied by adding to the water violet ink containing ethylene glycol. The Leverett function was measured by the semi-impermeable membrane method. The phase permeabilities were measured at the limit points and then calculated from the equilibrium displacement data using the water saturation distribution structure measured b$ the radio isotope method. The relative permeabilities of the water-filled channels k~(s I) were estimated from the phase permeabilities ki(s l) taking into account the change in the number of fluid-conducting channels with increase in the water saturation s I. The theory of percolation in regular lattices [5] shows that the dependences of the phase permeabilities of a lattice with equal pore channel permeabilities on the pore channel concentration are almost linear (curves 5 and 6 in Fig. la), except for the intervals near the limiting saturations. Then the ratios of the phase permeabilities for curves 5, 1 and 6, 2 give an estimate of the relative permeabilities of the channels ki(s z) for wetting and nonwetting fluids in a strictly hydrophilic medium. Figure ib also shows the dependence for the fraction of the volume of a pore channel occupied by water in a single step in a strictly hydrophilic medium (curve 3) and for almost neutral wetting (curve 4). These curves are only estimates and take into account the properties of the measured phase permeabilities (Fig. la). At small values of s I the change in the volume of the water in the pores is due mainly to the saturation of t h e n e c k s of the constrictions of the pore channels, whose volume is small, and the value of ~(s I) does not exceed 0.i. At saturation values s I > 0.3 the pore expansions also begin to be saturated and ~(s z) increases, which corresponds to an increase in km(s 1) (Fig. la). Near the limiting water saturation much of the nonwetting fluid has already been displaced from the pore channels in the previous steps and estimates show that ~(s I) cannot exceed 0.5--0.6. For almost neutral wetting the phase permeabilities for water lie much lower than when 8 = 20 ~ , and the fraction of the volume of the pore channel occupied by water in a single step is less than for a strictly hydrophilic medium (Fig. Ib, curve 4). The values of the disequilibrium function for a strictly hydrophilic porous medium (curve I) and for almost neutral wetting (curve 2), calculated from (1.6) on the basis of the data of Figs. la and Ib, are presented in Fig. 2. The form of the disequilibrium

893

/

ki

a

7

/

oA

o,,

!\

,/,

0.5

~Z

o3

0.6

s7 1.0

0.4

0,8

s7

0

Y

0.5

Fig. 2

Fig. 1

functions for different wetting conditions is essentially different, which is chiefly associated with the very low values of the capillary pressure when s I > 0.5 for almost neutral wetting (Fig. ib). Relations (1.5) and (1.6) close the system of equations of nonequilibrium flow through the porous medium, which in the variables saturation s and overall filter velocity v has the form [2]:

m - as, - ~ + div [vf(a,)

+

o~km F(~,) kz (~l)grad J ( a i ) ] : 0 div v=O

(1.7) (1.8)

In the one-dimensional case for a constant water flow velocity at the inlet v = v 0 at x = 0 Eqs. (1.7), (1.8) reduce to a single evolution equation describing the saturation distribution in space: as o o ao~ - ~ + ~ x F ( Cz) + e -~-x a ( a ) .~x

ec=s+DR(s)-~, t = - -m LrotA~s '

:

x~ x = --~ ,

(1.9)

0 ol/krn e = Lvol~---~

( 1.10 )

where the disequilibrium parameter D and the disequilibrium function R(s) are determined from (1.6). In the problem there are two small parameters associated with higher derivatives s and D and only one of them, the disequilibrium parameter D, is contained in the strongly nonlinear function F(a), which complicates the analysis of the solutions of (1.9). At the same time, it is clear that disequilibrium will mainly be manifested in regions with high rates of change of s, i.e., near the displacement front, as is confirmed by numerical calculation [3]. Therefore the nonequilibrium displacement may be esentially different in nature for different forms of the disequilibrium functions, in particular on the saturation interval 0 < s < s c, where s c is the frontal saturation in the solution of the Buckley--Leverett problem (Eq. (1.9) with D = 0, s = 0). This should be most clearly expressed when nonequilibrium displacement in a strictly hydrophilic porous medium is compared with the case of almost neutral wetting, since in these two cases the disequilibrium functions are essentially different (Fig. 2). In order to investigate the nature of nonequilibrium displacement at different values of the disequilibrium parameter and the disequilibrium function we carried out displacement experiments and compared the experimental data with the results of a numerical analysis of the system of equations (1.5), (1.6), (1.9). 2. Experimental Apparatus and Measurement Technique The experiments were carried out on horizontal cylindrical and rectangular models of the porous medium. The diameter of the working section of the cylindrical model was 5.2.10 -2 m, and its length 1.7 m; the dimensions of the rectangular working section were 0.01 x 0.2 • 0.6 m. The working section was filled with vibrocompacted quartz sand

894

with a permeability ~i0-I0 -12 m = and a porosity m x 0.4. In order to obtain a uniform porous medium the working section was filled with water before introducing and vibrocompacting the sand. The working section was then dried, vibrocompacted, evacuated and filled with an aqueous solution of NaCI at a concentration of 2 g/liter or with a 50% solution of violet ink subsequently displaced by a liquid hydrocarbon of the necessary viscosity. This preparation of the porous medium made it possible to obtain a uniform and easily reproducible wettability of the porous medium by the displacing fluid and to simulate the initial water saturation characteristic of oil reservoirs. The liquid hydrocarbon was displaced by the aqueous solutions at a constant flow rate varied over a broad interval, which made it possible to obtain both equilibrium and nonequilibrium displacement regimes. All the experiments were carried out at a pressure at the working section outlet P0 = 0.i MPa and a temperature =20~ During displacement the structure of the saturation jumps was recorded by means of the gamma indicator method. Promeran containing a radioactive mercury isotope was introduced into the displacing fluid. During the process of displacement from the rectangular model the spatial distribution of the gamma radiation was measured with a SEGAMS (Hungary) gamma camera, whose measuring element is a Na/l scintillation crystal 0.3 m in diameter, which was displaced along the length of the working section. The use of a standard high-resolution collimator made it possible to obtain a resolution of 5 • 5 mm with respect to area and 14% with respect to energy. In processing the image on the computer we took into account the nonuniform sensitivity of the detector photomultiplier and the statistical nature of the gamma radiation. The pressure difference along the length of the working section was measured with a slide-wire pressure gauge, and during the experiments we also measured the volume flow rates of the fluids at the outlet from the working section. The phase permeabilities corresponding to the experimental conditions and the capillary pressure are shown in Figs. la and lb. 3. Calculation Scheme In order to determine the evolution of the saturation distribution along the length under nonequilibrium displacement conditions we solved the system of equations (1.9), (i.i0) numerically for the initial condition s(x, 0) = 0, the boundary condition at the inlet x = 0 a5 8a (5) ~

- F (5) = - i

corresponding to the absence of a displacing fluid flow rate and for fluid take off at the outlet proportional to the mobilities (no end effect) a5 ea (5)--

=

0

Ox

The equation for s(x, t) is of the third order and has not been solved for the higher derivative 8s/St; therefore, differentiating (i.i0) with respect to t and substituting for 8s/St from (1.9), we first obtain the equation for

at

Ox

I

at

at

By virtue of the fact that R(0) = 0 the initial value for ~ is taken equal to zero, For (3.1) by the integrointerpolationmethod we constructed an implicit conservative difference scheme of the second order in the space coordinate. In the n-th layer the term R(s)Ss/St was approximated using Eq. (I.i0) in terms of a n and s n, which made it possible to preserve a two-layer scheme. By the method of quasilinearization the nonlinear equation for an+1 was reduced to a linear equation, which was solved by the threepoint sweep method, After determining ~n+l we solved Eq. (i. I0) for s n+l in accordance with an implicit scheme with iterations in the nonlinearity using Newton's method. Then employing the values of an+i and s n+1 obtained, we improve the accuracy of the solution until the specified degree of accuracy was achieved. The initial approximation was chosen using linear extrapolation with respect to time.

895

9

az

9

9

e, e o

I

Oo o

0,5

s

t~Lxz,t,t,~r,?

~ ~t, Lx~t,-

0.~

Fig. 3 1,gk~zi O

,Z a3 o

o~ e O

0.~

I

7i

.

:/

OO

/

0

~176176 oz tl

Fig. 4

m

0,8

0,7

Fig.

0,8

s~

5

4. Discussion of the Results The experiments showed that in a strictly hydrophilic porous medium, OH = 20 ~ , when the length of the working section L = 0.6 m and P0 = 0.058 for capillary numbers Nc2 less than 1.7"i0 -s the experimental data on the structure of the saturation jumps are well generalized by the data of the numerical calculations for D = 0 and the length of stabilized zone decreases as the displacement rate increases. For capillary numbers Nc2 greater than 1.7"i0 -s we observed an increase in the length of the stabilized zone with increase in the displacement rate, which indicates nonequilibrium displacement [2]. In Fig. 3 we have plotted the experimental profiles of the saturation jumps for U0 = 0.058, o = 30"10 -3 N/m, k = 6.10 -12 m 2, Nc2 = 4.7-10-5, and ~ = 0.061 (points i) for the four moments of time t = 0.13, 0.33, 0.53, and 2.13. The points 2 in Fig. 3 represent the initial distribution of the water saturation along the length of the porous medium. Curve 3 in Fig. 3 shows the structure of the saturation jumps for equilibrium displacement. Under the experimental conditions of Fig. 3 the calculations for the equilibrium model show that near the displacement front the structure of the solution is essentially different from the experimental data. Curve 4 in Fig. 3 represents the results of the nonequilibrium displacement calculations made for a disequilibrium function corresponding to a strictly hydrophilic porous medium (curve i in Fig. 2). As in the experiments, the increase in the disequilibrium parameter D with increase in the capillary number leads to the flattening out of the leading front. The closest agreement between calculation and experiment was observed when D = 8"10 -3, which corresponds to a value of ~a/~n s in (1.6) equal to 300. Almost the same value of ~ / ~ n s for a strictly hydrophilic porous medium was obtained by comparing the experimental water saturation profiles with the calculations for U0 = 0.i and Nc2 = 1.8-10 -5 The displacement experiments on a working section with L = 1.7 m showed that for the same capillary number the disequilibrium is much more weakly expressed, as follows directly from (1.6). As the length of the working section increases, the value of ~ / L decreases, the value of the disequilibrium parameter D falls and, consequently, for nonequilibrium displacement the length of the stabilized zone decreases. The experiments on displacement in a strictly hydrophilic porous medium showed that the nonequilbirium effect is strongest in the stabilized zone at the displacement front, whereas behind the displacement front the effect is weak (Fig. 3). In the case of almost neutral wetting the displacement experiments showed that when L = 0~ m the nonequilibrium effect is also observed at Nc2 = 8.8"10 -7 , which is lower than for a strictly hydrophilic porous medium. In Fig. 4 we have plotted the distribution 896

of the water saturation along the length of the working section for U0 = 0.058 and a capillary number Nc2 = 5.4-10 -7 , s = 5.95, t = 0.56 (points i) and Nc2 = i.i'i0 -~, s = 0.248 (points 2) for the moments of time t = 0.2 and 0.42. The points 3 in Fig. 4 represent the initial water saturation distribution. In the experiments with Nc2 = 5.410 -7 the experimental profile is well generalized by the calculations made in accordance with the equilibrium model (curve 4), D = 0. As distinct from the strictly hydrophilic porous medium as the displacement rate increases (Nc2 = I.i.i0 -s) in the experiments a considerable decrease in the frontal saturation was observed and the displacement coefficient at the moment of breakthrough of the water fell from 0.58 to 0.42. The curves 5 in Fig. 4 represent the nonequilibrium displacement calculations for a disequilibrium parameter D = 0.2 determined in accordance with (1.6) for 8a/3n s = 300, and the disequilibrium function corresponding to curve 2 in Fig. 2. As in the experi7 ments, the calculation results also give a considerable fall in the frontal saturation. This is associated with the high degree of nonuniformity of the disequilibrium function on the saturation range from zero to the frontal value in equilibrium displacement: s c = 0.6. On the saturation interval 0 < s < 0.4 the disequilibrium function is small and only slightly flattens out the leading front (Fig. 4). When s > 0.4 the disequilibrium function increases sharply, which leads when D = 0.2 to the strong flattening of the saturation profile and the formation of a new frontal saturation s~. The value s~ is close to the saturation value s ~ 0.5, at which the disequilibrium function reaches a maximum. The experiments and calculations on the fluid viscosity ratio interval 0.25 < ~0 < 0.023 for almost neutral wetting of the porous medium by the water showed that in all cases during nonequilibrium displacement there is a considerable fall in the frontal water saturation and the displacement coefficient at the moment of breakthrough of the water. A fundamental feature of nonequilibrium displacement is the fact that the functional dependences of the phase permeabilities remain the same as in equilibrium displacement. Figure 5 shows the dependence of the phase permeability ratio kz2 = kz/k 2 for the displacing and displaced fluids on the water saturation s I for equilibrium (curve I) and nonequilibrium (points 2) displacement for ~0 = 0.058 and the experimental conditions of Fig. 4. For nonequilibrium displacement the phase permeability ratio was determined from the experimental data on the ratio of the fluid flow rates at the outlet from the working section after breakthrough of the water, when the saturation gradients are small over the entire displacement zone (Fig. 4). The data for equilibrium and nonequilibrium displacement in Fig. 5 are similar, which confirms the basic hypothesis of the nonequilibrium displacement model [i, 2] concerning the constancy of the functional dependences of the relative phase permeabilities under the experimental conditions. LITERATURE CITED i. G. I. Barenblatt, "Flow of two immiscible fluids in a homogeneous porous medium," Izv. Akad. Nauk SSSR, Mekh. Zhidk. Gaza, No. 5, 144 (1971). 2. Go I. Barenblatt and A. P. Vinnichenko, "Nonequilibrium flow of immiscible fluids through porous media," Usp. Mekh., ~, 35 (1980). 3. O. B. Bocharov, V. V. Kuznetsov, and Yu. V. Chekhovich, "Numerical investigation of the nonequilibrium flow of immiscible fluids through porous media," Inzh-Fiz. Zh.,

5J_7, 91 (1989). 4. G. I. Barenvlatt, V. M. Entov, and V. M. Ryzhik, Theory of Unsteady Flow of Liquids and Gases Through Porous Media [in Russian], Nedra, Moscow (1972). 5. S. Kirkpatrick, "Percolation and conduction," Rev. Mod. Phys., 45, 574 (1973).

897