The Numerical Simulation Study of the Oil–Water Seepage Behavior Dependent on the Polymer Concentration in Polymer Flooding

It is well acknowledged that due to the polymer component, the oil–water relative permeability curve in polymer flooding is different from the curve in waterflooding. As the viscoelastic properties and the trapping number are presented for modifying the oil–water relative permeability curve, the integration of these two factors for the convenience of simulation processes has become a key issue. In this paper, an interpolation factor Ω that depends on the normalized polymer concentration is firstly proposed for simplification. Then, the numerical calculations in the self-developed simulator are performed to discuss the effects of the interpolation factor on the well performances and the applications in field history matching. The results indicate that compared with the results of the commercial simulator, the simulation with the interpolation factor Ω could more accurately describe the effect of the injected polymer solution in controlling water production, and more efficiently simplify the combination of factors on relative permeability curves in polymer flooding. Additionally, for polymer flooding history matching, the interpolation factor Ω is set as an adjustment parameter based on core flooding results to dynamically consider the change of the relative permeability curves, and has been successfully applied in the water cut matching of the two wells in Y oilfield. This investigation provides an efficient method to evaluate the seepage behavior variation of polymer flooding.


Introduction
After the primary depletion and the waterflooding developmental processes, there is still much oil remaining in reservoirs [1][2][3]. Since the water cut of oil blocks gains over 90%, resulting in no longer remarkable waterflooding economic benefits, polymer flooding has been widely applied in oilfields as one of the main enhanced oil recovery (EOR) methods [4][5][6]. In order to clearly characterize the oil-water flow capacity in the polymer flooding, the oil and water relative permeabilities are the focus of this study.
In past decades, numerous researchers have reported that for polymer flooding, the oil-water relative permeability curve is regarded as an important role in describing seepage behavior, well performance matching and developmental plans establishment. To reflect the effect of polymer solutions on relative permeability curves, various numerical methods have been employed based on experiments or field data. In 1989, pore-network modeling (PNM) for describing the polymer shear in capillaries was presented [7]. Based on this research, dynamic PNM and adaptive dynamic PNM were respectively applied to consider the complicated influence factors of the relative permeabilities 2 of 19 in polymer flooding, such as polymer properties [8], interfacial tension and capillary number [9], rock wettability, and the hysteresis phenomenon [10,11]. Although the PNM has the ability to intensively describe the interaction between solid surface and polymer solution, the scale of the PNM model could hardly achieve the representative elementary volume (REV) of the actual field. For the field-scale application, automatic-history-matching methods have been used to inversely calculate the oil-water relative permeability of polymer flooding [12]. With more robust inversion methods based on the iterative ensemble Kalman filter (EnKF) [13] or the Levenberg-Marquardt (LM) algorithm [14], relative permeabilities and rock properties were modified by large actual data and massive simulation computation [15,16]. To further improve the predictive capability of numerical evaluation, the simulator UTCHEM combined the Corey-type relative permeability model with the dimensionless trapping number [17,18], which depends on the interface tension (IFT), rock permeability, buoyancy force and viscous force [19][20][21]. The residual saturation in the relative permeability models also could be dynamically recalculated by the trapping number in the process of polymer injection [22]. Besides the trapping number, the viscoelastic properties of polymers were paid attention to improve the microscopic oil displacement efficiency [23,24]. In the presented polymer flooding simulator, the viscoelastic properties, which depend on the polymer concentration and the relative molecular mass of polymers [25], were combined with the trapping number to describe the oil saturation variation and the increase of the oil's relative permeability [26]. To avoid the fluid modeling failure caused by the missing date of polymer properties and to improve the computational efficiency, the oil-water relative permeability curves for polymer flooding in the commercial simulator were modified by polymer absorption without considering the factors, viscoelastic properties or trapping number [27].
In summary, many calculated methods have been emphasized to describe the oil-water relative permeabilities in polymer solutions. Caused by the viscoelastic properties and trapping number, the variation of the fluid flow capability has been fully and carefully considered in non-commercial scientific simulators; however, few researchers have explored the synthesis and simplification of the two main factors for their convenient application in simulation processes, such as history matching adjustment.
To fill this gap, a polymer concentration interpolation factor Ω is presented to simplify the two main factors and conveniently reflect the flow behavior in polymer solutions. Then, compared with the results of the commercial simulator, the case study is launched to discuss the effect of the factor Ω on well performances. Finally, with the experimental results, the factor Ω is set as a flexible parameter to successfully match the polymer flooding water cut of two producers in Y oilfield.

Assumption
To construct the numerical calculation of the polymer flooding, some assumptions are made in this work: (1) The oil phase, water phase and the polymers are considered as three components; (2) The percolation of each component in the formation is isothermal and obeys the Darcy flow model; (3) The solubility of the polymer in the oil phase is neglected; (4) The oil-water flow capability could be affected by the polymer concentration.

Polymer Concentration Interpolation Model of Oil-Water Relative Permeabilities
To integrate the effect of the viscoelastic properties and the trapping number into oil-water relative permeabilities, the polymer concentration interpolation model is presented in this work. In the interpolation model, polymer concentration is selected to concisely reflect the viscoelastic properties and the interface tension (IFT) based on the relative permeability curves of polymer flooding experiments [28][29][30], since the relative permeability curves are for a given polymer type. Further, with large intervals of the experimental relative permeabilities between different polymer Energies 2020, 13, 5125 3 of 19 concentrations, the interpolation method of the relative permeability curves is consulted to dynamically evaluate oil and water mobility based on the change in the polymer concentration.
Based on the above analysis, a table interpolation factor Ω related to the normalized polymer concentration R is proposed to describe the changing degree of the relative permeability curves (see Table 1). The normalized polymer concentration is defined as the ratio of the polymer concentration to the maximum polymer concentration c pmax . In Table 1, the c 1 p and c 2 p are the polymer concentrations in the experiments. In the first row of Table 1, R is set as 0 and Ω = 0, since this row corresponds to the relative permeability curve at the zero polymer concentration (waterflooding case). In the last row of Table 1, R is 1 constantly and Ω = 1, since this row corresponds to the relative permeability curve at the maximum polymer concentration. With the interpolation factor Ω, a new relative permeability curve at the objective normalized polymer concentration can be obtained through the following steps: Table 1. The interpolation factors.

Normalized Polymer Concentration (R)
Ω Step 1: According to the table function, the specific interpolation factor Ω * of the objective normalized polymer concentration (R = c * p /c pmax ) can be calculated. Then, based on the Ω * and taking the equal proportional points from l 1 (water-phase relative permeability curve at zero polymer concentration) and l 2 (water-phase relative permeability curve at the maximum polymer concentration), the corresponding scaled permeability curves l 3 and l 4 can be obtained (Figure 1a).
Step 2: Selecting a certain water saturation S w , one can get points A and B on the scaled permeability curves l 3 and l 4 ( Figure 1a).
Step 3: Based on the points A and B, point C can be calculated by a similar method to Step 1.
Step 4: By selecting multiple water saturation values and repeating Step 2 and Step 3, an interpolated water relative permeability at the objective polymer concentration can be calculated.
Further, oil-relative permeability curves can be achieved via the above interpolation method. Figure 1b is a schematic diagram showing the use of the interpolation method proposed in this paper to calculate the relative permeability curves at any normalized polymer concentration. The blue line and the red line are regarded as the basis relative permeability curves reflecting, respectively, the cases of waterflooding and the maximum polymer concentration. Based on the basis curves, the interpolated relative permeability curves for different polymer concentrations are shown as the dot lines in Figure 1b. Step 3: Based on the points A and B, point C can be calculated by a similar method to Step 1.
Step 4: By selecting multiple water saturation values and repeating Step 2 and Step 3, an interpolated water relative permeability at the objective polymer concentration can be calculated.   Further, oil-relative permeability curves can be achieved via the above interpolation method. Figure 1b is a schematic diagram showing the use of the interpolation method proposed in this paper to calculate the relative permeability curves at any normalized polymer concentration. The blue line and the red line are regarded as the basis relative permeability curves reflecting, respectively, the cases of waterflooding and the maximum polymer concentration. Based on the basis curves, the interpolated relative permeability curves for different polymer concentrations are shown as the dot lines in Figure 1b.

Mathematical Model
The mathematical model contains the polymer component, the oil phase and the water phase.

Mathematical Model
The mathematical model contains the polymer component, the oil phase and the water phase.
The key difference from the black oil model is the polymer's properties. In this paper, the polymer flow mechanisms are considered as follows: (1) Polymer adsorption effect The interactions between the injected polymer solution and the porous media give rise to the mildly viscous adsorption film on the rock surface [31]. Usually, the adsorption effect reduces water permeability by dividing by the actual resistance factor R k (see Equation (1)). RRF (RRF ≥ 1) is a constant parameter called the residual resistance factor, which depends on the rock type.ĉ pmax is the Energies 2020, 13, 5125 5 of 19 maximum adsorbed concentration. The adsorbed concentrationĉ p is expressed by the generalized Langmuir function of the polymer concentration c p (see Equation (2)).
(2) Polymer viscosity effect The polymer increases the viscosity of the aqueous solutions to improve the oil-water mobility ratio, which is an important factor for enhancing oil recovery. Due to the dissolving degree of polymers, a mixing parameter ω is defined to calculate the water effective viscosity µ w,e f f and the polymer effective viscosity µ p,e f f , respectively. Based on the Todd-Longstaff mixing model [32], µ p,e f f is evaluated in Equation (3), where µ f p c p is the viscosity of the full-mixed solution with the c p polymer concentration and µ f mp c pmax describes the viscosity of the c pmax maximum polymer concentration solution. Equation (4) introduces the µ w,e f f calculation method, where µ w,e is the water viscosity of the partially dissolved polymer solution. With the same mixing model, µ w,e can be calculated in Equation (5). µ w is the water-phase viscosity.
(3) Polymer rheology effect Due to the rheology characteristics of non-Newtonian fluids, the viscosity of the polymer solution depends not only on the pressure but also on the shear rate, which differs from the black-oil model. For describing the rheology effect, a shear factor Γ sh is added in mathematical models to modify µ w,e f f and µ p,e f f (see Equation (6)). With the same mechanism model as commercial simulators [27], the shear factor Γ sh is calculated in Equation (7), where f µ,w is a table interpolation function to express the effect of the water's velocity on the sheared water viscosity.
(4) Inaccessible pore volume fraction As the injected polymers in the solution have large sizes, these polymers cannot access some certain smaller pore throats. In a porous media, the volume of the polymers that have not established effective sweeping percolation is called the inaccessible pore volume (IPV) [33]. To describe the flow characteristics of the polymer solutions, S ipv is defined as a fraction of pore volume, which depends on the rock type and polymer properties.
Based on the functions describing the polymer properties, the flow equations are obtained as follows. For polymer component: For water component: For oil component: Saturation equation: Capillary pressure equation: Figure 2 shows the overall solution workflow of the mathematical model with the polymer concentration interpolation model section. The solution in this work is based on the solution of the black-oil model. After the discretization of the partial differential equations, the pressure, saturation and polymer concentration for each grid are solved by the fully-implicit finite difference method. In each iteration step, the interpolation factor Ω is determined by the ratio of the polymer concentration to the maximum polymer concentration. Then, the oil-water relative permeability of the corresponding polymer concentration can be interpolated by Ω. The simulator containing the polymer concentration interpolation model is programmed by Fortran. As the output documents are in the same format as the commercial software, the simulated results could be exhibited with the aid of other software. Figure 3 shows the compared results of the oil recovery and the water cut between our self-developed polymer flooding simulator and the commercial simulator, to demonstrate its reliability [  The simulator containing the polymer concentration interpolation model is programmed by Fortran. As the output documents are in the same format as the commercial software, the simulated results could be exhibited with the aid of other software. Figure 3 shows the compared results of the oil recovery and the water cut between our self-developed polymer flooding simulator and the commercial simulator, to demonstrate its reliability [34]. The simulator containing the polymer concentration interpolation model is programmed by Fortran. As the output documents are in the same format as the commercial software, the simulated results could be exhibited with the aid of other software. Figure 3 shows the compared results of the oil recovery and the water cut between our self-developed polymer flooding simulator and the commercial simulator, to demonstrate its reliability [34].

Solution and Verification
The results of the interpolated relative permeability curves are validated by comparison with the experimental results. In Figure 4, it can be seen that the evaluated relative permeability curves are matched well with the experiment data, which demonstrates the applicability and stability of the interpolation process for polymer concentration variation.    The results of the interpolated relative permeability curves are validated by comparison with the experimental results. In Figure 4, it can be seen that the evaluated relative permeability curves are matched well with the experiment data, which demonstrates the applicability and stability of the interpolation process for polymer concentration variation.

Cases Study
Based on the methodology and model solution, a mechanism model with a five-spot pattern is applied in this work for the case study (see Figure 5). The basic parameters of the mechanism model and the properties of the polymer solution are respectively shown in Table 2 and Figure 6. The production period of the case study includes three parts: depletion for 15 days, waterflooding for 10 years and polymer flooding for 10 years. In each developmental period, the well's performances are changed as well. In depletion, the liquid production of the producer (PROD01 well) is set as 30 m3/day, with a 10 MPa limit of the bottom hole pressure (BHP). In the waterflooding and polymer flooding period, the liquid production of the producer increases to 80 m3/day, and the four injectors are switched to inject liquid at 20 m3/day constantly, wherein the liquid is water or polymer solution

Cases Study
Based on the methodology and model solution, a mechanism model with a five-spot pattern is applied in this work for the case study (see Figure 5). The basic parameters of the mechanism model and the properties of the polymer solution are respectively shown in Table 2 and Figure 6. The production period of the case study includes three parts: depletion for 15 days, waterflooding for 10 years and polymer flooding for 10 years. In each developmental period, the well's performances are changed as well. In depletion, the liquid production of the producer (PROD01 well) is set as 30 m 3 /day, with a 10 MPa limit of the bottom hole pressure (BHP). In the waterflooding and polymer flooding period, the liquid production of the producer increases to 80 m 3 /day, and the four injectors are switched to inject liquid at 20 m 3 /day constantly, wherein the liquid is water or polymer solution (c p = 2.28 kg/m 3 ) for the corresponding period. years and polymer flooding for 10 years. In each developmental period, the well's performances are changed as well. In depletion, the liquid production of the producer (PROD01 well) is set as 30 m3/day, with a 10 MPa limit of the bottom hole pressure (BHP). In the waterflooding and polymer flooding period, the liquid production of the producer increases to 80 m3/day, and the four injectors are switched to inject liquid at 20 m3/day constantly, wherein the liquid is water or polymer solution   To discuss the effect of polymer concentration variations on well production in the polymer flooding period, three cases (Base case, Case A and Case B) with different relative permeability curves are evaluated in the simulator (see Table 3). The simulated results show that the decrease degree of the water cut in the Base case is the lowest, since it ignores the difference in flow behavior between the water phase and polymer solution (see Figure 7a). In increasing the oil production of the producer, Case A has the largest effect of these three cases (see Figure 7b). This is because the displacement efficiency of the water phase in the reservoir has been overestimated by the oil-polymer solution relative permeability curve when the polymer flooding starts to work. The calculated water cut and the oil production of Case B are between those of the Base case and Case A. Compared with Case A, in the sweeping area of the injected polymer the water saturation in Case B increases with the polymer concentration due to the effects of polymer concentration variations on water flow capability (see Figure 8). Although a similar distribution trend is observed in the Base case, the estimated deviation of the polymer solution's mobility may lead to inexact predictions. Thus, Case B could describe the flow behavior of the polymer flooding more accurately.   Case A, in the sweeping area of the injected polymer the water saturation in Case B increases with the polymer concentration due to the effects of polymer concentration variations on water flow capability (see Figure 8). Although a similar distribution trend is observed in the Base case, the estimated deviation of the polymer solution's mobility may lead to inexact predictions. Thus, Case B could describe the flow behavior of the polymer flooding more accurately. Moreover, Case B, Case C and Case D are applied in this work to study the changing degree of the interpolation factor Ω with the polymer concentration (see Table 3 and Figure 9). The respective interpolation factors Ω 1 , Ω 2 , Ω 3 are shown in Table 4 and Figure 10. The calculation results demonstrate that different Ω could influence the timing of controlling the water cut and increasing oil production. Thus, the oil-water relative permeability curves dependent on the polymer concentration improve the applicability of the polymer flooding simulation. Further, the interpolation factor Ω could be set as a flexible parameter in history matching to reflect the interaction between different rock types and different polymer solutions. Moreover, Case B, Case C and Case D are applied in this work to study the changing degree of the interpolation factor W with the polymer concentration (see Table 3 and Figure 9). The respective interpolation factors 1 Table 4 and Figure 10. The calculation results demonstrate that different W could influence the timing of controlling the water cut and increasing oil production. Thus, the oil-water relative permeability curves dependent on the polymer concentration improve the applicability of the polymer flooding simulation. Further, the interpolation factor W could be set as a flexible parameter in history matching to reflect the interaction between different rock types and different polymer solutions.  Table 3. Model setting of the cases study.

Treatment of the Relative Permeability Curve in Polymer Flooding
Base Case (ECLIPSE simulator [27]) Curve at c p = 0.00 kg/m 3 Case A Curve at c p = 2.28 kg/m 3 Case B Curve calculated by the interpolation factor Ω 1 Case C Curve calculated by the interpolation factor Ω 2 Case D Curve calculated by the interpolation factor Ω 3

Field Study
Based on the above analysis of the case study, the oil-water relative permeability curves with polymer concentration variation exhibit a clear difference in the production of the polymer flooding, and the interpolation factor W is set as the adjustable parameter to fit the field data. Here, the production performances of two producers in the Y oilfield are considered to implement the polymer concentration variation in the history matching process. Since the water cut of the main productive zone increases sharply, the Y oilfield employs polymer flooding at 3440 days.
For well S07 and well B12, the water cut matching results of the waterflooding period are within the engineering error. However, the oil-water relative permeability curves for waterflooding could not apply to the history matching of polymer flooding, which is called in this work the simulated results before implementation (see Figures 11a and 12a). To achieve better results, core experimental data of

Field Study
Based on the above analysis of the case study, the oil-water relative permeability curves with polymer concentration variation exhibit a clear difference in the production of the polymer flooding, and the interpolation factor Ω is set as the adjustable parameter to fit the field data. Here, the production performances of two producers in the Y oilfield are considered to implement the polymer concentration variation in the history matching process. Since the water cut of the main productive zone increases sharply, the Y oilfield employs polymer flooding at 3440 days.
For well S07 and well B12, the water cut matching results of the waterflooding period are within the engineering error. However, the oil-water relative permeability curves for waterflooding could not apply to the history matching of polymer flooding, which is called in this work the simulated results before implementation (see Figures 11a and 12a). To achieve better results, core experimental data of relative permeabilities with different polymer concentrations in the verification section are used to obtain the interpolation factor Ω (see the red points in Figures 11b and 12b). With these experimental results, more supposed interpolation factors proportional to the normalized polymer concentration c p /c pmax are interpolated in order to describe the dynamic oil-water relative permeability. Based on the trial water cut matching results, the supposed Ω is adjusted dynamically. Figures 11b and 12b illustrate the comparison results of the Ω adjustment. Correspondingly, the matching results after implementation are closer to the actual field data (see the red lines in Figures 11a and 12a).
results, more supposed interpolation factors proportional to the normalized polymer concentration p pm a x c / c are interpolated in order to describe the dynamic oil-water relative permeability. Based on the trial water cut matching results, the supposed W is adjusted dynamically. Figures 11b and 12b illustrate the comparison results of the W adjustment. Correspondingly, the matching results after implementation are closer to the actual field data (see the red lines in Figures 11a and 12a).

Conclusions
For polymer flooding, the polymer solution causes an alteration in the oil-water phase flow behavior. This alteration will significantly influence the well performances, such as the water cut and oil production rate. In this study, an interpolation method dependent on polymer concentration is applied to dynamically describe the variation in the oil-water relative permeability curves. The following observations are obtained from this study: (1) With the normalized polymer concentration R , a polymer concentration interpolation model for numerical research is proposed to dynamically calculate the corresponding flow mobility of the oil-water phases based on experimental validation.
(2) Compared with that in the commercial simulator, the proposed interpolation method could more effectively enhance the numerical predictions of the oil production rate and the water cut. (3) By combining with the core experimental results, the presented method has been successfully applied in the field history matching of the main productive zone in the Y oilfield, leading to the flexibility improvement of the history matching process.

Conclusions
For polymer flooding, the polymer solution causes an alteration in the oil-water phase flow behavior. This alteration will significantly influence the well performances, such as the water cut and oil production rate. In this study, an interpolation method dependent on polymer concentration is applied to dynamically describe the variation in the oil-water relative permeability curves. The following observations are obtained from this study: (1) With the normalized polymer concentration R, a polymer concentration interpolation model for numerical research is proposed to dynamically calculate the corresponding flow mobility of the oil-water phases based on experimental validation. (2) Compared with that in the commercial simulator, the proposed interpolation method could more effectively enhance the numerical predictions of the oil production rate and the water cut. (3) By combining with the core experimental results, the presented method has been successfully applied in the field history matching of the main productive zone in the Y oilfield, leading to the flexibility improvement of the history matching process. Acknowledgments: The authors thank the anonymous reviewers and all the editors in the process of the manuscript revision.

Conflicts of Interest:
The authors declare no conflict of interest.