Computer Simulation and Comparison of the Efficiency of Conventional, Polymer and Hydrogel Waterflooding of Inhomogeneous Oil Reservoirs

The oil displacement in a layered inhomogeneous reservoir using two types of physical-chemical technologies (polymer flooding and hydrogel flooding) is the subject of this research. In the first case the aqueous polymer solution of the desired concentration is injected into the porous reservoir creating the high-viscous moving fields. Unlike this technology, the hydrogel flooding is characterized by creation and evolution of the moving hydrogel field directly in porous medium in result of chemical reaction between the water solutions of two gel-forming components which one after another are injected into the oil reservoir with given time interruption. The first component is sorbed more intensively and moves slower than the second one, so when it gradually overtakes the first solution, they begin chemically react with creation of hydrogel. Special numerical methods, algorithms and computer software are developed to solve these systems of nonlinear equations, study and compare an efficiency of the oil field development at the different type of waterflooding. It is shown that creations of the moving polymer or hydrogel fields significantly increases the uniformity of oil displacement in all layers of reservoir and improve their basic exploitation parameters due to the cross-flows between layers and creation of the moving structures in the velocity field of two-phase flow. In doing so, hydrogel technology may be much more effectiveness in comparison with polymer flooding.


Mathematical theory
Development of reservoirs with complex structure by usual waterflooding is characterized by a low oil recovery due to uneven water-oil displacement of low-permeable areas of the porous medium. To improve this technology highly viscous agents are used to stimulate of the displacing ability of water and for the enlargement of the flooding zone within the oil reservoirs.
One of the ways to redirect the filtrational flow is the creation of motionless (stationary) barriers in the oil reservoir. However, field experiments [1][2][3][4] show that these barriers have low efficiency in the layered inhomogeneous reservoirs, consisting of hydraulically interrelated layers with different physical properties. This is because fluid flows round the stationary barriers, so even at a short distance from them these barriers have weak effect on distribution of fluxes within layers. In its turn, the creation of controlled moving areas of the thickening agent within the inhomogeneous reservoirs can significantly intensify their development [1][2][3][4][5][6]. Therefore these technologies are preferential for increasing oil extraction efficiency.
It is possible to distinguish two types of physical-chemical technologies for stimulation of oil field development. In the first of them the thickening agent is injected into the reservoir as a finished mixture that is typical for the polymer waterflooding. The filtration theory of finished mixtures is not required to describe their formation. The mixture composition is completely determined by the value of the thickening agent concentration in water [7][8][9][10]. By second technology both the filtration and the creation of the highly viscous mixture are due to chemical reactions are interrelated and occur in the porous media at the same time. Mathematical model of this phenomenon includes two groups of differential equations. The first, hydrodynamic group includes the system equations of two-phase filtration [10]: , Here τ is the time, S is the water saturation, P is the pressure, C is mass concentration of the thickener; m is the porosity, K is the absolute permeability of porous media, v w , v o are the velocities of water and oil, respectively; v is the total filtration velocity of two-phase mixture; μ w , μ 0 are viscosities of water and oil phases; is viscosity of the thickener water solution of at C = 0, α g is an empirical coefficient; , are the relative phase permeabilities: where S * is the irreducible water saturation; S * is the limiting water saturation.
The second group consists of the mass transfer equations of the mixture components ( ) ( ) These equations are written on basis of general principles of hydrogeochemistry [11] and should be used together with Eqs. (1) -(3). Concretization of (4), (5) was done in [10] at the following assumptions. Firstly, the thickening agent R is created as a result of the chemical reaction v 1 R 1 + v 2 R 2 = vR in the mixture between its two components R 1 , R 2 , where C, C 1 , C 2 are their mass concentrations. Secondly, the composition of the mixture (water phase) has an effect on two-phase filtration only through dependence (2) of viscosity μ w on the concentration C of the thickener: μ w = μ w (C). In Eqs. (4), (5) ω is difference between the numbers of acts of direct ("+") and reverse ("-") chemical reactions per unit time; M i , M are molecular masses; [C] = C/M is the molar concentration of the thickening agent per unit volume of fluid; α i , α are mass concentrations of the i-th gel-forming component (i = 1,2) and thickener in the absorbed (immovable) state. Thirdly, the chemical reactions are the local equilibrium ones, so hereinafter the value ω is equal zero in the last relation (5), hence, or , where is k R the constant of the chemical reaction. This means that Eqs. (1) -(5) describe the filtrational and chemical processes in the porous media when the mixture composition is very close to the equilibrium. In this case it is generally agreed [11] that in the differential Eqs. (4), (5) ω is unknown function. The dependencies α = α(C), α i = α i (C i ), i = 1, 2 on the concentrations in the water flow are determined by the sorption isotherms. The term «sorption» combines the physical and chemical adsorption of admixture on the surface of the porous medium skeleton, dissolution of admixture in the material of skeleton grains and the mechanical retention of admixture in the narrowing of porous channels. Note also that in the present paper we consider the situation when the thickener (hydrogel) is formed by chemical reaction between two gelling components which are successively injected into the reservoir. The first component R 1 is sorbed more intensively, and therefore, moves slower, so the second component R 2 comes close to the first one and begins to chemically react with it at the rear edge of the moving area that R 1 occupies in the oil reservoir at this time point.
It is necessary to note that this paper focuses on the operating regimes of the oil reservoir development when convection dominates diffusion so in modeling we do not take into account the last factor. Diffusion becomes prevail process only in situation when filtration is practically stopped.
The model (1) -(5) of two-phase five-component filtration allows us to study the influence of the controllable hydrogel field on the flow in the layered inhomogeneous reservoir. Hereinafter we consider a plane-parallel filtration in the vertical cross section of reservoir, consisted of N layers, see Fig. 1. Bottom holes of injection and production wells can include all layers on thickness of reservoir at its left-and right-hand boundary, respectively, or some parts γ L , γ R of these boundaries. Injection well can operate at given pressure P L Note that the second relation (7) is valid only at the permeable part γ L of the left-hand boundary of reservoir and this is a prerequisite for uniqueness of the system solution (1), (2).
At the bottom-hole γ L of injection well we set also the boundary conditions for water saturation S and concentrations C 1 , C 2 , C: Here and -the initial time and the end moment of injection the gel-forming component R m , -its initial concentration in water phase, m = 1,2; .
Obviously, if the water solution of the first component R 1 is injected from the beginning of the reservoir operation. Otherwise, the reagent R 1 can be entered into the reservoir at some point when the water content in the flow rate of producing well reaches a predetermined value , or after injection a specified amount of water that determined by the ratio , where V m is the total volume of the pore space of reservoir; L is its length. Moment comes after injection into the reservoir the water solution of component R 1 in quantity , measured in share of V m . Parameters and of R 2 -injection are defined by analogy.
The producing well, located at the boundary x = L, is operated at the given bot-tom-hole pres- Absolute permeability K has discontinuities of the first kind at the boundaries γ l of layers, so that, the matching conditions are satisfied at these boundaries: The solution of the problem (1) -(11) with the initial conditions for the functions S, C 1 , is to be computed in the area Ω = {0<x<L; 0<z<H} for 0<τ≤T, where T is the time of the reservoir exploitation, S 0 (x,y) is the water saturation distribution in the reservoir at τ = 0. The functions S, C 1 , C 2 , C are the one-valued piecewise continuous bounded functions and pressure P is continuous function with piecewise continuous derivatives.
Note that equations of two-phase three-component filtration (water -oil -polymer) describing the water-oil displacement using a polymer as a thickener, are a particular case of system (1) -(3) when C 1 = 0, C 2 = 0, ω = 0. The viscosity of the [ ] 0 where is the jump of a function y, ; is the fraction of water in the two-phase flow rate (the Buckley-Leverett function); . Condition (10) implies that the function S is continuous if the phase permeabilities and are determined by the identical formulas like (3) in every layer.
The roof, the bottom of the reservoir at z = 0, z = H and the non-perforated pieces Г L , Г R of the side boundaries γ L , γ R are impermeable, so: polymer solution depends on concentration C, absolute permeability K of the reservoir and mass amount α of thickener that adsorbed in pores. This dependence is determined from experimental data and can be written as following [10]: where C max is the maximum value of the thickener concentration; is the viscosity of water without the thickener at C = 0; A, B, and α p are empirical coefficients that must be determined throughout the range of values of absolute permeability of the oil reservoir. The expression in Eq. (8), enclosed in square brackets, is characterized by the resistance factor and the residual resistance factor. It takes into account the decrease in mobility of the polymer solution in the presence of the oil phase in a porous medium. In this paper we assume the adsorption α = α (S,C) is in equilibrium and irreversible, so that this dependence is determined by the Henry sorption isotherm α(S,C) = ГSC at , where Г is the Henry coefficient which depends on the absolute permeability and can be obtained experimentally.
Initial time τ b and end moment τ e of injection the water polymer solution can be set in much the same way as for the gel-forming component R 1 and The Eqs. (1) - (12) are solved by the finite-difference methods. We note here only some principal features in discretization of the solution region Ω and development of the completely conservative scheme [10,12,13]. The grid Ω h is non-uniform: it is built in such a manner that its step h x along the horizontal x-axis is constant and in the vertical direction z steps h z in the layers are different. Furthermore, the points of the additional mesh are also appended in the every layer by means of displacement of initial grid per half a step in the xand z -directions. Both grids Ω h and are used to improve approximation of the filtration flows of phases in the mass-transfer Eqs. (1), (4), (5) which are used for calculation of the functions S, C 1 , C 2 , C. These equations are approximated in the points of these grids by the up-wind finite-difference scheme of the second order. To provide for stability of solution the Eqs. (1), (4), (5) the time step of the scheme is determined in according to the Courant-Friedrichs-Lewy criterion. The pressure P is calculated from the quasi-linear elliptic equa- , that would not be difficult to obtain from Eqs. (1), (2). This equation is approximated by the five-point difference scheme of the second order.
The appropriate system of algebraic equations in respect to P is solved by the iterative method [10,12,13] of a high rate convergence (3-5 iterations) at every time step. Its convergence rate does not depend on the total number of the grid points and matrix coefficients whose values are deter-mined by function k(S,C) and can differ by several orders. Convergence of the numerical solution was checked using larger number of grid points and is also confirmed comparing the numerical and analytical solutions of one-and two-dimensional test problems. The numerical model is implemented in the software using the C# programming language, Task Parallel Library (TPL), a set of libraries NVIDIA CUDA and OpenCL. This package allows us to simulate two-phase multicomponent filtration with simultaneous visualization of results of computations which are presented and analyzed below. ; S (x,z,0) = S * = 0.2; S * = 0.8. At the formation of hydrogel the viscosity of its water solution even for relatively low values of concentration is significantly higher than the viscosity of the polymer solution, so that in the formulas (2) and (13) α g =1500 and α p =10; A = 4; B = 1.1, respectively.

Results of the computations
Water solution of the first component R 1 is injected into the reservoir from the beginning of its exploitation ( = 0), and the second reagent R 2almost at once upon after finishing the injection of R 1 . The quantities , of these water solutions are equal to 0.15 and 0.2 in shares of the reservoir pore volume V m . The flow rate Q L is specified at the injection well. Calculations are carried out till 96 per cent of the water cut of the flow rate Q R of the production well.
The results of the computations are presented below in tables and figures. Variants 1, 2 and 3 correspond to conventional, hydrogel and polymer waterflooding. Each color of palette at the bottom of figures corresponds to the certain value range of the displayed two-dimensional function. The dotted lines of the double thickness are shown the boundaries between layers and side boundaries of the oil reservoir. Herein-after, in all figures C 1 , C 2 and C are the dimensionless concentrations of the gel-forming components R 1 , R 2 and hydrogel (or polymer) that are defined as the ratio of the actual dimensional concentration of the corresponding reagent to its maximum concentration.

Waterless period of reservoir development
Oil displacement by conventional waterflooding has a relatively even character in the high-permeable middle layers, but the low-permeable upper layers are practically not involved in the development, see to 0.216 shares of the pore volume V m , and its oil recovery factor reaches 27%, see Table 1. Layered heterogeneity produces a very intricate character of filtration in the reservoir due to cross-flows between its layers and creation of the moving structures [10] (the leading and rear waves) in the velocity field of two-phase flow. These structures can form at the boundaries of layers in result of interaction the moving jumps of the piecewise-continuous function S and the discontinuity lines of absolute permeability K.
The profile of such a structure in the vicinity of the leading edge of water saturation S at time τ = 16 days is shown in Fig. 2b. For clarity, the directional field of vector V of two-phase flow is plotted only in the left part of the oil reservoir where the moving structure is located. The vectors V are shown by colored arrows of the same length. Dark colored arrows correspond to a greater flow velocity.

Table 1
Development parameters of the oil reservoir and its layers at the moments and T of water breakthrough in the producing well and completion oil production Variant → ῆ at τ = at τ = T at τ = T As easy to see in Fig. 2b, the specific feature of the S -structure is a change in direction of the vertical component v z to the opposite, so the flow pattern in the S -structure is similar to the conventional wave. The S -structure moves along the 1st, 3rd and 5th layers very slowly because of their lower permeability in comparison with the 2nd and 4th layers. However, its prolonged exposure at the slow movement of phases along the reservoir has a significant effect on redistribution of water saturation. As a result, the cone-shape profile of function S is formed in the middle layer, see Fig. 2b. Namely owing to action of such moving S -structure, a round shape zone of high oil content is formed in the 3-rd layer at 50 < x < 100 m at the time point = 107 days when waterless period operation of the producing well is being finished, see Fig. 2a.
Water saturation S is about 0.34 ÷ 0.5 in this zone.
In case of hydrogel flooding the injection of the water solution R 1 is completed at the time point ≈ 50 days. Almost immediately, at ≈ 55 days, the injection of the reagent R 2 begins, which continues more than 3 months until ≈ 155 days. After that at τ > only water is received into the oil reservoir.
At the stage 0 ≤ τ ≤ of the R 1 -injection the process of the water-oil displacement is the same as at conventional waterflooding because the concentration C 1 does not affect the viscosity of the water solution R 1 and the filtration characteristics (velocity field V, distribution of water saturation S, etc.) of the oil reservoir development.
Similarly, the reagent R 2 at the stage of its injection does not influence filtration as long, until at ≈ 65 days its leading edge does not reach the rear edge of the component R 1 . At τ < the distributions of concentrations C 1 , C 2 are formed only under action of the same S -structures that are created at conventional waterflooding of the layered inhomogeneous reservoir.
Thereafter the moving field of the thickener R appears in the reservoir as a result of the chemical reaction between gel-forming reagents R 1 and R 2 in the common intersection region of the fields C 1 and C 2 . Maps of concentrations C 1 , C 2 and C are shown in Fig. 3 at the time of water breakthrough in the producing well ( = 107 days, Ṽ gel = 0.216). To aid the visualization, the figure shows those parts of the oil reservoir, in which concentrations C 1 , C 2 and C differ from zero. Comparison of the water saturation maps in Figs. 2a and 3 shows that hydrogel already begins influence significantly on the phase redistribution in the reservoir. However, for the period 0 < τ < thickener has no time to influence significantly the integral characteristics of reservoir development.
As it is easy to see in the Table 1, the oil recovery factors ῆ from the reservoir are practically identical both for conventional and hydrogel waterflooding, though the values ῆ 1 already begin to differ and their small distinction is caused by redistribution of water saturation S between the layers, compare Figs. 2a and 3. As result the oil recovery factors ῆ 1 , ῆ 3 , ῆ 5 of the less-permeable layers are increased but inverse effect takes place in the high-permeable layers. Small impact of hydrogel on the reservoir development is caused by short duration (about two months) of waterless stage 0 < τ < , so the moving field of the thickener has a time to significantly equalize water oil displacement in the layers due to redirection of the phase flows between them only in a small part of the reservoir. In consequence, the operational efficiency of the producing well remains the same for both type of flooding. However, as illustrated in Fig. 3 below the leading edge of the hydrogel field already starts forming a high-saturated petroleum zone -«the oil billow». Water saturation S in this zone is close to constant, the value of which is equal to 0.5 in the case being considered.
It is interesting to note that in the non-uniform porous medium thickener R additionally generates its own moving structures [10] caused by the motion of jumps of the piecewise-continuous function C. Their initiation is associated with dependence of hydrogel viscosity and, as a result, the relative conductivity of two-phase flow from the magnitude of C. Such C -structures are similar to S -structures, and they also are formed at the intersection points of the leading and rear edges of concentration C and the discontinuity lines of absolute permeability K. The most intense C -structure, which has almost vertical vectors of the total flow V, is formed in the vicinity of the intersection point of the leading edge of the gel concentration C and boundaries γ 1 of layers. It is important to note that structures of both types may interact with each other. In multi-layered oil reservoirs this effect has an additional impact on the filtration flow. Typical forms of C -structures which arise in the directional field of the vector V in the vicinities of the leading edge of the function C at τ = are demonstrated in Fig. 3. This front is shown here as solid line.
In polymer waterflooding, the influence of viscous thickener on the movement of phases in the reservoir starts at τ = τ b = 0 and lasts continuously to the water breakthrough into the producing well at the time point = 127 days. By this time a large part of the polymer solution (Ṽ = 0.254) has been already injected into the reservoir, so that unlike hydrogel technology, the thickener does affect the filtration during all period of waterless operation of the producing well. This results in the growth of the oil recovery factor ῆ of the reservoir, see Table 1. Fields of water saturation S in Figs. 3 and 4 illustrate distinction between two filtration processes with the use of the polymer waterflooding and hydrogel technology. It is easy to see that the lower intensity of Cstructures in the motion of the polymer solution in porous media, which is caused by its smaller viscosity in comparison with the viscosity of the hydrogel, leads to the less uniform water-oil displacement in layers and to the smaller saturated oil billow.

Discontinuity in thickener field
In case of hydrogel waterflooding, the reagent R 2 is injected into the reservoir even more 1.5 months after water breakthrough to the production well ( ≈ 155 days).Throughout this process, creation and motion of the hydrogel barrier in the high permeable second, fourth and middle layers leads to flow around a barrier by liquids through the less permeable first and fifth layers. These liquids flow into 1-th and 5-th layers at the rear edge of hydrogel field and flow out them at the leading edge. Eventually, the first component R 1 is pressed to the boundaries γ 1 and γ 4 between these low and high permeable layers that, in turn, affects the amount of the created thickener in these areas. Within this process, the hydrogel field also migrates into these layers and accumulates at the vicinities of boundaries γ 1 and γ 4 , see Fig. 5. The narrow barrier, located in the second, third and fourth layers, pro-motes improvement the oil displacement from the less permeable upper and lower layers.
Up to a certain point in time the quantities of reagents R 1 and R 2 in the common intersection region of their fields are sufficient for creation the gel τ gel V pol τ pol V singly-connected region of hydrogel, and its moving field is able to withstand and redirect the filtration flows from the injection well to the top and bottom of the oil reservoir.
However, by the time point = 245 days ( = 0.489) values of concentration, C 1 , C 2 of the reagents R 1 , R 2 in the intersection region become too small to provide the singly-connectedness of the hydrogel barrier. Its breakdown begins in the high-permeable layers with formation of the «corridor» for water-oil flow. This corridor is clearly demonstrated in the enlarged fragment of the oil reservoir in right part of Fig. 5 which shows a place of discontinuity of the hydrogel field, filled in blue. With widening of corridor an intensity of the oil displacement from the less permeable layers significantly decreases.
In the polymer technology, injection of the thickener solution of quantity V pol = 0.3 is being completed at = 151 days. Then only water enters into the reservoir. As illustrated in Fig. 6, the breakdown of the polymer barrier under the pressure of water flow also begins in the second most permeable layer at the time point = 227 days ( = 0.453).

The end of reservoir development
Under conventional waterflooding, the water content in the flow rate of the producing well reaches a limit value = 96% when T wat = 945 days, and its operation stops thereafter. The total amount of fluid extracted from the reservoir during its development, is equal to the volume = 1.886 of the injected water. Characteristics of the filtration process are shown in Fig. 7. Additionally, Table  2 demonstrates the values of the total cross-flows along the whole boundary γ k between the k-th and (k+1)-th layers in forward and reverse direction for the three considered types of flooding during the period of growth the water cut θ of well production. The dimensionless values q k→k+1 are given in shares of the k-th layer pore volume (e.g., q k→k+1 = 1 is the amount of liquid that is equal to one such volume).
As it is easy to see in Table 2, for this period of conventional waterflooding the cross-flows between all layers are significantly increased that results to their more uniform development. However, in spite this, these cross-flows θ wat V η η gel V q k→k+1 , q k+1→k are not so intensive, to align the saturation of phases throughout the thickness of the reservoir: at time point τ = T wat heterogeneity of the layered structure is detected in distribution of S as easily as it does at , see Figs. 2 and 7. By the end of the conventional waterflooding there are no S -structures in the oil reservoir, so the vectors of velocity field become practically horizontal, see Fig. 7a. Such simple filtration process leads to very considerable nonuniformity of water-oil displacement. Owing to this, at τ = T wat the oil recovery factors ῆ i of layers with high-and low permeability differ approximately in 1.5 to 2 times, and the total recovery factor of the reservoir is enough small: = 48.3%, see Table 1. With the use of the hydrogel flooding the values T gel = 390 days and = 0.78 are 2.42 times less than in previous case, but cross-flows q k→k+1 , q k+1→k , between the layers are almost 3 ÷ 10 times more, see Table 2. As a result, the development uniformity of reservoir becomes so high, that lamination of its structure is not practically reflected in the map of water saturation S, see Fig. 8.

Table 2
Total cross-flows between the layers of the reservoir at τ = and τ = T: 1, 2, 3 -conventional, hydrogel and polymer waterflooding, respectively At the time point τ = T gel a large part of the hydrogel field is characterized by very low concentration C, see Fig. 8, however, it continues to affect the structure of the velocity field V. By example, the fragment of the reservoir in the vicinity of the injection well is shown in Fig. 7b. The areas of residual low concentration continue to create Cstructures with the transverse flows of the phases in the reservoir. Unlike the conventional waterflooding, C -structures stimulate mass transfer between layers until the completion of the reservoir development at τ = T gel .
In case of polymer waterflooding the values T pol = 510 days and = 1.021 are 1.3 times less, and the cross-flows q k→k+1 , q k+1→k are approximately 1.5 to 5 times more intense compared to similar values of conventional waterflooding, but about two times lower than with the use of the hydrogel technology, see Table 2. Because of this, the polymer waterflooding usually does not allow to reach the same high uni-formity of oil displacement in the layers, as in case of hydrogel flooding, so that maps of water saturation S quite accurately reproduce the layered heterogeneity of oil reservoir, see As we said earlier, oil billow is being created ahead of the moving hydrogel (or polymer) barrier. When this wave comes to the right-hand boundary γ R , the water cut θ of well production fast decreases and stabilizes at the lower level. As a result, the rate of the oil recovery begins to grow, see Fig. 10. With further evolution of the filtration process, the field of thickener also reaches the boundary γ R and the water cut of the well production quickly rises to the maximal value . With this, zones of maximal concentrations C are mainly located in the less permeable layers, see Figs. 8 and 9. Such character of filtration process causes to the comparatively short exploitation periods T gel = 390 and T pol = 510 days of the producing well at a high level of the water cut of well production in comparison with the duration T wat = 945 days of conventional waterflooding. Moreover, in case of the hydrogel flooding, the production quantity of the associated water and duration T gel of the oil reservoir exploitation significantly decrease (see Table 1 and Fig. 10). It is also easy to see in this table that the creation of the hydrogel field enhances the ultimate oil recovery of the reservoir (of 14.1%) and every layer (from 8 to 24%). It is interesting, that the incremental oil production is much more (24%) in the first less-permeable layer in comparison with the second high-permeable seam (8%). Furthermore, the oil recovery factor of the third layer is only by 1.4% smaller than , although absolute permeability K 2 of the second layer is four times as large ( is about 9.8% and 4.9% as less as in case of the conventional and polymer waterflooding, respectively.

Conclusions
The principal results of the simulations can be briefly presented as follows. The oil displacement in the layered inhomogeneous reservoir using the moving hydrogel fields can greatly differ from conventional and polymer waterflooding. A controllable hydrogel field allows us to increase significantly the uniformity of oil displacement in all layers of reservoir and improve their basic exploitation parameters. However, the use of hydrogel technology requires a detailed study of the structure and properties of the reservoir. Optimal technological parameters (properties of the injected gelling components, their amounts and etc.) can be determined in each concrete case by computer simulation.