Study on influence of two-phase filtration transformation on formation of zones of undeveloped oil reserves

In order to study the process of fluid filtration during flooding of an oil field, article uses Rapoport – Lis model of non-piston oil displacement by water. During plane-radial filtration in a homogeneous formation, radii of disturbance zones are determined with and without taking into account the end effect. Influence of changes in value of capillary pressure gradient on distribution of water saturation coefficient in the non-piston displacement zone for high and low permeability reservoirs is revealed. Application of an element model for a five-point injection and production well placement system showed that, using traditional flooding technology, flat-radial fluid filtration is transformed into rectilinear-parallel. At solving equation of water saturation, Barenblatt method of integral relations was used, which allows determining the transformation time. By solving the saturation equation for rectilinear-parallel filtration, change in the value of water saturation coefficient at bottomhole of production well for an unlimited and closed deposit is determined. It is shown that an increase in water cut coefficient of a production well is possible only for a closed formation. To determine coefficient of water saturation in a closed deposit, a differential equation with variable coefficients is obtained, an iterative solution method is proposed. In the element of the five-point system, oilsaturated zones not covered by development were identified. For channels of low filtration resistance, conditions for their location in horizontal and vertical planes are established. It is shown that, at maintaining formation pressure, there is an isobar line in formation, corresponding to initial formation pressure, location of which determines direction of fluid crossflow rates. Intensity of crossflows affects application efficiency of hydrodynamic, physical and chemical, thermal and other methods of enhanced oil recovery.

Introduction. Model of oil displacement by water has been well studied in classical works on underground hydrodynamics, nevertheless, study of oil recovery process during flooding is one of the most difficult tasks in design and control of oil field development. This is due to fact that when water is injected in bottomhole formation zone (BHFZ) of injection well, an end effect (EE) takes place. Water accumulates in a region of radius r k , in which water saturation coefficient increases to a certain value of s k . When this value is exceeded, water enters the formation, process of two-phase filtration and non-piston oil displacement by water begins. Water flows at internal boundary of EE region, and at external boundary water does not flow out, but oil does. Study of filtration process is complicated by a nonlinear boundary condition on external region of EE [14]. Because of it end effect is neglected, coefficient of water saturation at bottomhole of injection well takes limit (maximum) value s * . Porosity and permeability properties (PPP) of the formation, depending on saturation coefficients, geological structure of deposit, location of wells affect change in nature of the fluid flow in formation and, therefore, development of mobile oil reserves.
Statement of the problem. Existing problem of undeveloped residual mobile reserves (RMR) of oil is due to physical properties of the fluids, implemented development system, macroheterogeneity of the filtration parameters in lateral and thickness of oil reservoir [2,8,9,12,15,16,18]. Surface phenomena occurring at phase boundary, capillary pressure gradient, hydrophilicity or hydrophobicity of reservoir (microheterogeneity) significantly affect deviation of actual development indicators from design technological ones [7,11]. To eliminate discrepancy between the indicators, technologies (methods) of enhanced oil recovery (EOR) are used. When selecting and evaluating the effectiveness of corresponding method, it is necessary to take into account the change (transformation) of the filtration processes and influence of the fluid saturation degree in reservoir capacitive space.
Methodology. Physical models for multiphase filtration of non-piston oil displacement are considered in [4,6,10,13]. From equations of continuity, motion and rheological equations, a system of partial differential equations is obtained, solution of which is reduced to determining the distribution of saturation and phase pressure coefficients in disturbance region -two-phase filtration. To connect water saturation coefficient with pressure, an additional closing equation is introduced. All considered models can be divided into two groups: taking into account and not taking into account the capillary pressure gradient. Solution of differential equations systems, as a rule, is carried out by approximate numerical methods. In the present work, to solve Rapoport -Lis equations, Barenblatt method of integral relations is used [1]. Transformation of filtration zones in process of developing a production object (PO) is considered. Obtained approximate analytical solutions of water saturation coefficient values allow identifying the zones with existence of two-phase filtration and residual oil.
Discussion. Results of the research. To study transformation of filtration processes, a model of flooding a homogeneous reservoir is considered. With commissioning of injection well, a radial non-piston displacement of oil by water occurs. Pressure at the bottomhole of injection well is greater than initial formation and bottomhole pressure of production well. Over time, flat-radial filtration transforms into rectilinear-parallel filtration between injection and production wells [17]. To determine values of water saturation coefficient in the model proposed by Rapoport -Lis, a special case of which is Buckley -Leverett model, formation and fluids are considered incompressible. Task is reduced to solving the differential equation for saturation coefficient taking into account capillary pressure at phase boundary and corresponds to a rigid elastic water drive regime: where m 0 -open porosity coefficient; s -water saturation coefficient in a two-phase filtration zone; W(t) -total phases' filtration rate; x -coordinate, for radial filtration x is replaced by r; f(s) -Buckley -Leverett function; k 0 -absolute permeability coefficient; k 2 * (s) -relative phase permeability (RPP) of oil reservoir; p k -capillary pressure; Δρ = ρ 2 -ρ 1 , ρ 2 and ρ 1 -oil and water densities; g -acceleration of gravity; Δg 0 = g 2 -g 1 , g 2 and g 1 -initial pressure gradients during oil displacement by water in a two-phase filtration zone.
Buckley -Leverett function is determined experimentally: where  1 ,  2 -dynamic viscosity coefficients of water and oil; ) ( * 1 s k , ) ( * 2 s k -reservoir RPP for water and oil. Equation (1) contains time derivative of the function s(x,t), and operator div corresponds to an equation of parabolic type. Therefore, to solve (1), Barenblatt method of integral relations can be used [1].
Origin of coordinates will be compatible with bottomhole of injection well. Radius of disturbance zone of two-phase filtration is denoted by R(r, t). Since the change in water saturation coefficient corresponds to change in pressure of water injected into reservoir, solution is found in the form [6]  where s 0 -coefficient of residual water saturation (bound water); s * -water saturation limit value; r w -well radius. Heaviside unit function (t) equals 0 at t = 0 and 1 at t > 1: Note that f 1 (r, R(t)) < 0 for r w ≤ r ≤ R(t).
If the end effect is neglected, then function s(r, t) satisfies following boundary conditions: Pressure at the boundary of two-phase filtration zone R(t) is equal to initial pressure in deposit p 0 .
With a flat radial displacement of oil by water, equation (1) takes form , Q(t) -flow rate of injection well; h -formation thickness.
To determine the radius of two-phase filtration zone R(t), (6) is integrated over r in the interval from r w to R(t). Thus It follows from boundary conditions (5) that, for r = R(t) water saturation coefficient s = s 0 . RPP of reservoir takes maximum value ) ( 0 If q s = const after separation of variables and integration, a transcendental equation for determining R(t) is obtained: Consequently, under given boundary conditions (5), capillary pressure and RPP do not explicitly affect the radius of two-phase filtration zone.
However, if end effect of radius r k is not neglected, then second boundary condition (5) at bottomhole of injection well will take the form After integration in the interval [r k , R(t)] right part of equation (7) takes following form After substituting (10) in right part of equation (7), a differential equation is obtained that differs from (8) by its right part: From formula (10) it can be seen that parameter A k depends on technological parameter q s , phase permeability of oil, derivative of capillary pressure with respect to s at point s k , radius of a two-phase filtration zone R. According to the data given in [3], capillary pressure derivative with respect to s is equal to tangent of inclination angle to axis s (tg) in interval s 0 < s < s * . The highest values of capillary pressure derivative are observed in vicinity of water saturation coefficient s 0 corresponding to front of oil displacement by water, distanced from bottomhole of injection well ( Fig.1). Derivative values vary in range -10 < tg < -8. At the bottomhole of injection well, a change in capillary pressure for a high permeable reservoir (curve 1) has practically no effect on oil displacement: tg  -0.1. For a low permeable reservoir (curve 2) tg  -1.
Let us consider an element of operational object consisting of injection and production wells (Fig.2). While inequality R(t) + R 2 (t) < L holds, where L -distance between wells; R 2 (t) -radius of disturbance zone 2 for production well, in which only oil is filtered, there is flat radial filtration. Since injection pressure is greater than initial formation pressure and bottomhole pressure of production well, bulk of injected water will move to perforation zone of production well. Filtration becomes rectilinear-parallel. Pressure at interface of zones 1 and 2 (points С i ) is equal to initial formation pressure p 0 . Another part of water will flow at a speed V 1 to zone 1. In disturbance region 2 of production well, oil will flow into linear filtration channel at a speed V 2 . Moreover, flow speeds will decrease over time. Note that, under given boundary conditions, points С i are motionless.
For n = 2, distance from injection well to transformation boundary is R = 80 m, pressure at the boundary is p 0 . In zone 1, two-phase filtration takes place -oil displacement by water according to Rapoport -Lis model. Disturbance zone radius of production well, in which single-phase oil filtration exists, is 320 m. Since commissioning of injection well, time of changing the nature of the filtration process for given PPP and specified technological parameters does not exceed 4 days. Since pressure at bottomhole of production wells is less than pressure at other points in deposit, flat-radial filtration is transformed into rectilinear-parallel. Flooding of production wells will begin through channels A 1 B i .
To describe linear displacement of oil by water, equation (1) is used. Origin of coordinates coincides with bottomhole of injection well. Boundary conditions in this case change, since on interval [x w , R] a change in coefficient of water saturation has already occurred in formation, and at bottomhole of injection well for x = x w = r w coefficient of water saturation s = s * . Over time at border x = R coefficient of water saturation will increase: s > s 0 . Therefore consider s(x, t) on two intervals: For Here R -constant, corresponding to n (Table 2); с, b, a 0 , a 1, a n -unknown parameters determined from boundary conditions and conditions for merging of functions s 1 and s 2 .
Boundary conditions for function s 2 are set at the boundary of displacement front: Conditions for merging of functions s 1 and s 2 for x = R will be From a joint solution of system (16), (17)   After substitution of (18) into (15) we get For practical testing of determining the coefficients of water saturation and watering at bottomhole of production well, let us consider the case R < x ≤ l(t), s = s 2 .
After substitution of s 2 into equation (1) and subsequent integration over x in range from R to l(t), we obtain a differential equation for determining linear size of two-phase filtration zone l(t): Here s is value of water saturation coefficient corresponding to x = R. It should be noted, that s  s 0 is variable value, which will increase over time. Buckley -Leverett function varies in the interval   s f and also increases. Verification calculation of parameter A kl shows that for various values included in formula, its order does not exceed 10 -4 , and from (20) with increasing l(t) the second term in square bracket tends to zero. Therefore, equation for determining l(t) has the form From where cubic equation is obtained Using the Cardano method, l(t) is found: Value s is determined from formula (19), which for x = R will be To obtain an approximate analytical solution, we estimate the function (t), included in right part of (21). Let total phase filtration rate W(t) be constant. Then, using mean value theorem, we can write It follows from (24) that value s depends on time, consequently, s  also increases over time t.
Limit value For s 0 = 0.2; Δs = 0.6 and R = 80 m value of s = 0.2485. Thus, Buckley -Leverett function f (s) in the integrand expression changes insignificantly (a weakly changing function) and its removal beyond the sign of integral is justified. In the future we will take s s   . Note, that at t →  formula (23) will take a simple form: From (23) and (27) Table 3  Note, that if x = L -distance between wells, then from formula (19) it follows that for l(t) → , s(L,t) = s 0 + a 2 , х -constant value. This is contrary to practice, since water cut increases, therefore, at bottomhole of production well, proportion of water in fluid flow should increase. Accordingly, size of filtration zone should be limited to l k = l(t k ), which is determined from geometric structure of the deposit.
Pressure distribution in confined deposits was considered in [1,16]. By analogy, change in coefficient of water saturation in a closed limited deposit is found in the form Unknown function D(t) is defined after substitution (29) in equation (1) and after integration in the interval [R,l k ] we get differential equation: where Here s R corresponds to value of water saturation coefficient at x = R, s k corresponds to x = l k , s * > s R > s k > s 0 . For the Buckley -Leverett functions, relation f (s k ) < f (s R )  1 is true.
It follows from (29) that function D(t) is limited above, its maximum value 0 * max /s s D  . Expanding the exponent in a row, we obtain an approximate solution: For calculation, following iterative scheme is proposed: 1. By known values l k and R from first formula (31) parameter H is found.  Setting a new value of time, we determine increase in coefficient of water saturation s 3 (L,t) at bottomhole of production well.
To determine water cut ratio  on channel A 1 B 1 we use the formula , ; where v -water filtration rate; v 2 -oil filtration rate; W -total fluid flow rate at bottomhole of production well. If W is constant, then terms change, oil filtration rate decreases, and water filtration rate increases.
Water filtration rate is expressed through total rate and depends on capillary pressure gradient [15]:  Filtration zone size l k depends on geological structure of deposit and significantly affects time of watering. In Rapoport -Lis model, porosity coefficient is constant, fluids and reservoir are incompressible, and therefore, linear filtration channels can be neglected. Consequently, in areas 3 (Fig.2, a) there will remain oil reserves not involved in development.
To determine oil reserves not involved in development in zones 3 we define t 22 -contact time of disturbance zones of production wells B 1 B 2 in point D 1 .
Distance between production wells B 1 B 2 = L 2 = 2 L . Radius of oil filtration zone R 22 for homogeneous reservoir Denote by t 22 time corresponding to reaching of point D 1 , -interface between zones of filtration of production wells: where χ 22 -coefficient of oil piezoconductivity; С -numerical coefficient. Radius R 2 and transformation start time t * satisfy the relation From formulas (36) and (37) we get We use data of example 2 and calculation results of Table 2. For L = 400 m, R = 80 m, t * = 3.3 days and R 2 = 320 m from (38) and (39) we get t 22 = 2.6 days, R 22 = 283 m. Consequently, disturbance zones of two production wells are in contact at a point D 1 earlier, than they connect at point C 1 of disturbance zone of injection and production wells.
We make an essential remark. Formation pressure in point D 1 at t 22 = 2.6 days is p 0 , which will further decrease. Therefore, in time interval t 22 ≤ t ≤ t * in zone 3 oil filtration will begin in direction of point D 1. If due to "proximity" of values time t 22 and t * is neglected, it is possible to determine oil reserves in areas not covered by flooding. So, at time t * three zones with different saturation patterns form in the PO element of the fivepoint system (Fig.2).
1. In vicinity of injection well first zone has the shape of a circle with radius R. Saturation pattern -water and oil.
2. Four zones 2, each of which has the shape of a closed region bounded by curved lines С i D i . Saturation pattern -oil.
3. Four zones 3, bounded by curved lines С i С j D i , not covered by development, saturated with oil, pressure in which is equal to initial pressure p 0 .
Denote the area of zone 3 by S 3 . Consider an isosceles right triangle B 1 B 2 A 1 , whose legs are equal L. From Fig.2 determine the area of zone 3: where S Δ -area of a triangle Consider several options for calculating the area of a figure B 1 D 1 C 1 (Fig.3). As can be seen from the figure, arc D 1 C 1 is convex with respect to the center B 1 . At changing the angle in the interval 0 ≤  ≤ /4 radius R  increases from R 22 = 283 m to R 2 = 320 m. This condition is satisfied by function For interval 0 ≤  ≤ /4 distance from injection well А 1 to point D  is determined by ratio It follows from formulas (42) and (43), that curvature radius R  and distance А 1 D  to injection well A 1 depends on angle . Calculation results are presented in Table 4.
If, to determine area S  , we take average radius r av , equal to half-sum of radii R 2 = 320 m and R 22 = 283 m, then Thus, approximation of section B 1 D 1 C 1 to sector of constant average radius gives an overestimated value of area, and, as follows from formula (40), a smaller value of area of zone 3, which will lead to an incorrect estimate of reserves.
Element of PO B 1 B 2 B 3 B 4 (see Fig.2) contains four zones 3, each with area S 3 . Volume of oil (geological reserves) To determine volume of recoverable oil reserves in zone 3 we take coefficient of oil displacement by water where  2 -coefficient of formation coverage by flooding. Example 5. Take m 0 = 0.2; h = 10 m; s 0 = 0.2; s * = 0.8;  2 = 1. Remaining data will be taken from example 3.
From formulas (40), (45), (46) we get S 3 = 17613 m 2 , V 3 = 12721 m 3 , V 3iz = 84540 m 3 . If input flow rate of well that penetrated zone 3 is 50 m 3 /day, then time spent on extraction of mobile reserves V 3iz , is 4 years, if boundaries of zone 3 do not change. Volume of geological oil reserves in PO element is 512000 m 3 . Share of unextracted reserves will be 0.22. However, over time, there will be an intrusion of water in zone 3, in which two-phase filtration will occur, the oil flow rate will fall.
Let us consider how filtration processes change over time with simultaneous operation of wells. Following remark is appropriate. Due to given boundary conditions in PO element, following inequality is true: where arbitrary point G   zone 3. Condition (47) is satisfied from the moment the zone formation begins. Consequently, flooding of zones 2 and 3 will occur unevenly. With growth of time t > t * > t 22 in zone 3, two-phase filtration begins, boundary C 1 G  of which is limited by an isobar with initial pressure p 0 . Single-phase filtration zone 2 is deformed, water injected into formation will begin to flow through the boundary along channel C 1 D  .
Thus, in a reservoir, which is permeability homogeneous in Oxy plane, linear filtration will occur in zones 3 and 2. Values of saturation coefficients and pressures in filtration channels will differ from initial values at beginning of operation.
Consider the filtering process in element of PO in form of a channels system connecting injection A 1 and production B 1 well.
Article [5] shows that watering of well fluid happens in presence of low filtration resistance (LFR) channels in formation. Cases of LFR channels' location along the thickness of reservoir and in the planes Oxy are considered. To determine channel length, parameters  j , were introduced, depending on concentration of indicator peaks and fixation time, j is indicator fixation number. For the element of the five-point development system considered above, we have From which Angle  = 45 corresponds to a straight line A 1 B 1 = L -the shortest distance between injection and production wells, parameter  j = 1. Angle  = 0 с corresponds to the greatest distance A 1 D 1 B 1 = 2 L , parameter  j = 2 . Therefore, LFR channels that satisfy the conditions lie in horizontal Oxy planes. If conditions (50) are not satisfied, then LFR channels are located in a vertical plane.
In Tables 3-5 of article [5] cases of horizontal arrangement of LFR channels L j for three given parameters  j are considered. Table 5 shows values  j and L j from article [5, see Table 3] and dependencies  j and A 1 D  B 1 on angle . Obviously, for a set parameter there are two channels in horizontal plane L j : j = 1 and j = 2. Remaining LFR channels are located in a vertical plane. Note that methodology for calculating the volume of pore channels, necessary to apply alignment of injectivity profile is given in the article [5].
Each model has its own limitations, its own field of application, including the Buckley -Leverett model. Confirmation of reliability is compliance of calculated and actual field indicators [7]. It is known that watering of a production well does not occur "abruptly", but gradually increases with increasing water saturation coefficient.
If formation is layered-inhomogeneous, consists of interlayers (seams) of different permeability, then for each ith seam, its own zones 1, 2, 3 are formed. Filtration transformation time for each interlayer will be different, determined by formulas (9), (12), (13). In absence of crossflows between seams (within framework of accepted Rapoport -Lis model), following boundary conditions are possible: flow rate of water injected into formation should be distributed in proportion to conductivity coefficients k i h i of seams [8,15] or it should be assumed that bottomhole pressure is constant in all seams [11]. Note that each case requires a separate study. LFR channels correspond to highly permeable seams (HPS). For stationary filtration, presence of water "tongues" in the most permeable sections of formation, which is carried out through HPS, has been established. Change in water saturation coefficient is determined by formulas (29), (31), (33).
Original method for studying the distribution of residual mobile oil reserves in a laterally heterogeneous formation was proposed by M.I.Maksimov in a monograph [8]. It is shown that, depending on location of zones with permeability of 40, 200, and 1000 mD, a "development selfregulation" takes place in an area heterogeneous formation. In our opinion, this is a consequence of fluid flows (see Fig.2, b), individual for each zone.
In the 50s of the twentieth century. V.N.Shchelkachev [18] formulated conditions (criteria) for the use of forced fluid extraction (FFE) recommended for PO at a late stage of development: fluid water cut of at least 90 %, high productivity coefficients, collector stability, etc. We believe that large value of water cut coefficient is not a decisive criterion for determining the late fourth stage of PO.
As it is known, FFE technology consists in a phased change in operating modes of production and injection wells, in creation of high-pressure gradients in formation, allowing low-permeability interlayers, dead end and stagnant zones to be involved in development. This leads to a change in directions of the filtration flows, and presence of a hydrodynamic connection between differently permeable seams involves developing low-permeability differences, i.e. cyclic flooding. Therefore, there is a combined unstationary flooding. Existence of regions 3, formed as a result of filtration flows transformation and containing mobile unextracted oil reserves, indicates the need for oil recovery enhancement methods (EOR) that were not previously used in this PO. There are techniques that determine density of oil reserves in areas not covered by development, in particular the Voronoi method. In numerous scientific works of recent years [15,16], commissioning of second wellbores with a shallow and horizontal end and multi-stage hydraulic fracturing are recommended. Length of the wellbore in reservoir, as follows from above calculations, should not exceed the radius of drainage zone R 2 of production well (Fig.2, b). Methods for calculating the horizontal drain in formation are quite fully considered in literature [15,16].
Thus, a change in direction of the filtration flows leads to creation of linear filtration channels between production wells and formation of zones not covered by development (Fig.3). Obviously, this must be taken into account at designing and arranging second wellbores for efficient development of oil reserves.
Since bottomhole pressure of injection wells exceeds initial formation pressure p 0 , and bottomhole pressure of production wells, on the contrary, is less than p 0 , then isobar line corresponding to it -G 0 will exist in the formation. Moreover, existence of G 0 line does not depend on chosen physical model. Location of this line in space depends on boundary conditions and affects formation of filtration flows. On the other hand, isobar line G 0 limits efficient effect of all existing EOR: hydrodynamic, physical and chemical, gas, heat, fracturing, zones of efficient operation of horizontal wells endings, etc. Presence of fluid crossflows in horizontal plane (see Fig.2, b), as well as flows in hydrodynamically interconnected seams in layered-heterogeneous formations affect efficiency of the technologies used.

Conclusion
1. Study of traditional flooding technology using Rapoport -Lis model of non-piston oil displacement by water for plane-radial filtration taking and taking not into account the end effect found that changing the capillary pressure gradient practically does not affect distribution of saturation coefficients in two-phase filtration zone.