Numerical Simulation Study on Waterflooding Heavy Oil Based on Variable Threshold Pressure Gradient

The heavy-oil flow in porous media is characterized by non-Darcy law with variable threshold pressure gradient (TPG) due to the large fluid viscosity. However, available analytical and numerical models hardly consider this effect, which can lead to erroneous results. This paper is aimed at presenting an innovative approach and establishing a numerical simulator to analyze the heavyoil flow behavior with waterflooding. The apparent viscosity of the oil phase and flow correction coefficient characterized by the TPG were applied to describe the viscosity anomaly of heavy oil. Considering the formation heterogeneity, the TPG was processed into a variable related to mobility and the directionality. The discretization and linearization of the mathematical model were conducted to establish a fully implicit numerical model; the TPG value on each grid node was obtained through oil phase mobility interpolation, and then, the Jacobi matrix was reassembled and calculated to solve pressure and saturation equations. The corresponding simulator was thus developed. The pre-/postprocessing module of the simulator is connected to ECLIPSE; then, an efficient algorithm is introduced to realize a fast solution. Results show that considering the TPG will not only reduce the waterflooding area but also reduce the oil displacement efficiency because of aggravating the nonpiston phenomenon and interlayer conflict. The numerical simulation study on the TPG of heavy oil provides theoretical and technical support for the rational development and adjustment of water-driven heavy oil.


Introduction
Due to the tremendous growth of oil demand and the depletion of low-viscosity oil reservoir, heavy oil plays an increasingly great role in petroleum supply for the world and more and more attention is turned to heavy-oil development [1][2][3]. Heavy oil is a complex fluid that exhibits high viscosity and high density resulting in various difficulties during exploitation [4][5][6]. As the viscosity of heavy oil increases, its non-Newtonian characteristics become more obvious leading to the threshold pressure gradient (TPG) in the low-pressure gradient area [7][8][9]. Therefore, it is necessary to consider its impact on the flow characteristics in the numerical simulation [10].
Different from the TPG in tight oil reservoirs caused by ultralow permeability, the TPG in heavy-oil reservoirs is mainly attributed to the abnormal viscosity [11][12][13]. A large number of mechanism studies have shown that the TPG described in heavy-oil simulation could be summarized into the following three methods: pseudo threshold pressure gradient (PTPG) model, piecewise nonlinear model, and variable permeability model. The PTPG model is the most used model which processes the tangent of the arc segment to treat the curved segment of the flow curve equivalently, while the model reduces the fluid movable range [14][15][16]. To improve the description accuracy of the relationship between the flux and the pressure gradient, the piecewise nonlinear model was presented, such as the exponential function [17] and the power function [18]. The model expressions of these functions are complicated and involve the judgment of the critical point between the linear flow part and the nonlinear flow part, so the practical application of the model has certain difficulties [19][20][21]. From another perspective, Liu et al. proposed the variable permeability model to equivalently characterize the TPG effect as the rock permeability changes dependent on the variable pore pressure, which is between the traditional TPG methods and Darcy flow model [22]. In the variable permeability model, the TPG could be dynamically evaluated with the basis of formation pressures. Additionally, for the heavy-oil non-Newtonian property, some fractal models of the initial pressure gradient in porous media were proposed [23,24].
At present, numerous scholars have researched various nonlinear flow models to describe the TPG behavior. Although the TPG has been considered in numerical simulation, few studies have been involved in the effect of the heavyoil fluidity on the TPG, and the ignorance may cause misleading developmental plans. Moreover, all the existing models could not conveniently obtain the convergent results by the implicit pressure and explicit saturation (IMPES) solution method.
In this paper, based on the analysis of the actual experimental data, a variable TPG flow model considering heavyoil mobility is presented to modify the oil phase flux. With the presented variable TPG flow model, the full-implicit numerical model was established and the corresponding self-developed simulator was compiled. Then, after the validation compared with the commercial simulator ECLIPSE, the analysis for the different cases is discussed such as Darcy flow, static TPG, and variable TPG. Finally, the simulator is applied to the studied heavy-oil blocks to explore the fluid distribution and percolation characteristics.

The Variable TPG Flow Model considering Heavy-Oil Mobility
The Darcy flow model cannot objectively describe the dynamic characteristics of heavy-oil reservoirs in waterflooding development and exaggerate the flow capacity due to the ignorance of the TPG. In this paper, a TPG flow model considering heavy-oil mobility is presented based on the core displacement experimental data.
2.1. The Critical Percolation Viscosity. In order to improve the accuracy of the presented model in numerical simulation, the critical percolation viscosities are regarded as the most important parameter in the experiments. Through the regression analysis of the experimental results, the critical percolation viscosities under different permeability levels are obtained, as shown in Figure 1. The critical viscosity increases with the core permeabilities, and the increasing degree is gradually slowing down. Therefore, the TPG is related to the flow capability, and its effect on the seepage behavior needs to be considered when the viscosity of the studied oil block is higher than the corresponding critical viscosity.

Flow
Model. Figure 2 illustrates the relationship between the experimental TPG and the fluidity of the two oil blocks in the studied reservoir, which could be regressed to a better power relationship based on the fitting curves.
With the above analysis, the TPG must first be processed into a fluidity-related variable. Firstly, the relationship model of TPG and fluidity should be determined based on the experimental data fitting and substituted into the oil phase equation of motion. The expression of the TPG of heavy oil is listed as follows: To consider the directionality of the flow correction coefficient based on introducing the oil phase flow correction coefficient and close to the actual situation, the oil phase flow correction factor is defined as follows: 3. Mathematical Model (1) Auxiliary equation (2) Initial conditions (3) Boundary conditions Inner boundary condition (fixed well production): Inner boundary condition (fixed bottom hole pressure): Outer boundary conditions (closed):  Figure 3: Flowchart of the numerical simulation calculation.
In the black oil model, only K rli±1/2,j,k is a strong nonlinear term. The variable start pressure gradient flow model introduces the calculation of nonlinear flow correction coefficients, which strengthens the nonlinear characteristics of conductivity calculation. It put forward higher requirements for the speed and stability of numerical calculation. The full implicit method is used to solve the model in order to meet the calculation needs. The specific steps are as follows: (1) The mathematical model of waterflooding heavy-oil flow is discretized by difference, and the differential equations are transformed into nonlinear algebraic equations (2) The system of nonlinear algebraic equations is transformed into a system of linear coefficient equations with unknown pressure and saturation  Based on the numerical simulation method of variable TPG, the corresponding numerical simulator can be compiled, so as to make a deep study on the variation law of the flow field in the heavy-oil field and make up for the limitations of commercial numerical simulation software. The flowchart of the numerical simulation method of variable TPG is drawn as shown in Figure 3. As the current numerical simulation method is difficult to describe the flow characteristics of heavy oil objectively, this paper uses the Fortran 90 language which is suitable for large-scale numerical calculation to develop a variable TPG simulator. The self-developed simulator can simulate the nonlinear flow of heavy oil, so it is named NONLINEAR. The simulator consists of three parts: preprocessing is used to build the numerical model, which docks with the preprocessing of the ECLIPSE. The keyword of the ECLIPSE is interpreted, and the keyword of heavy oil (HEAVYOIL, HEAVYEXP) is added to the data file. The part of the model calculation is used to calculate the numerical model of preprocessing. Given the strong nonlinear problem in the nonlinear flow of heavy oil, the fully implicit method is used to solve the pressure and saturation, and the advanced algorithm can be introduced, which makes the software have the advantages of stable calculation and good convergence. The postprocessing part is connected with ECLIPSE to realize the compilation of data files (EGRID, INIT, and RST) required by the ECLIPSE. The calculation results of this software include the TPG field (DPX, DPY, and DPZ) and effective displacement pressure gradient field (EDPX, EDPY, and EDPZ) which can be checked by the postprocessing of ECLIPSE. 4.3.1. Preprocessing Part. In this part, the key recognition of ECLIPSE can be directly added to the data file generated by ECLIPSE, that is, adding the keyword HEAVY OIL to the RUNSPEC part of the data file and adding the keyword HEA-VYEXP to the PVT part, to realize the nonlinear flow simulation of waterflooding heavy oil. The first item under the HEAVYEXP keyword represents on/off status, and the second and third items represent the values of the experimental fitting parameters a and b. The application example is shown in Table 1.

Postprocessing Part.
In this part, based on the completion of writing the data files required by ECLIPSE, the postprocessing of ECLIPSE can be used to view the simulation results. By adding keywords that output the TPG (STARTDPX, STARTDPY, and STARTDPZ) and effectively displacement pressure gradient (EDPX, EDPY, and EDPZ) below the RPTRST keyword, the TPG and effective displacement pres-sure gradient for each time step are output. The application example is shown in Table 1.

Model Computing
Section. Similar to the ECLIPSE software, this part uses a completely implicit method to solve, and in order to ensure the speed and convergence of the  7 Geofluids nonlinear calculation, a fast solver was developed. The users can copy the data file with the nonlinear flow keyword to the simulator calculation interface and enter the data file path and then click the "Enter" button to start the calculation.

Simulator Validation
To verify the rationality and accuracy of the newly developed simulator, an anti-five-spot black oil pattern consistent with the characteristics of the studied heavy-oil field is established to achieve the symmetry test and the simulation results are compared with the commercial simulator ECLIPSE. The main parameters of the model are shown in Table 2.

Symmetry Validation.
Symmetry validation is a kind of software test method for reservoir numerical simulation, which is to establish a symmetrical well pattern in the homogeneous model, set up the same parameters of each well, and observe whether the pressure and saturation fields are symmetrical. As shown in Figure 4, the saturation field and pressure field are well symmetrical, which proves that the calculation result of the simulator is reliable and stable.

Well Performance Validation.
To test the stability and the accuracy of the self-developed simulator, the evaluated well performance in the commercial simulation software is applied for comparison. As shown in Figures 5-7, by comparing the oil production rate, water cut, and pressure, it can be found that the errors of the two simulators are extremely small, which proves the rationality and accuracy of the self-developed simulator.

Results and Discussion
6.1. Variable and Static TPG. In order to research the variable effect of the TPG, the mathematical relationship between the TPG and fluidity of the Y3 block is substituted into the selfdeveloped simulator, which is called the VARIABLE case. For comparison, the results from the commercial simulator ECLIPSE are used in this work, which is called the STATIC case. Due to the lack of the variable TPG function in ECLIPSE, usually, the THPRES keyword, which represents the minimum additional flow resistance to be overcome when the fluid flows between adjacent equilibrium zones, is After validating the symmetry tests of the pressure distribution and the saturation distribution (see Figure 8), the comparison of the well performances between the VARIABLE case and the STATIC case demonstrates that the formation pressure in the STATIC case remains better resulting in the higher oil production rate and the lower water cut (see Figures 9-11). Therefore, the relative equivalent static method for the TPG process causes erroneous evaluation, and the presented flow model considering the variable TPG is more realistic for the waterflooding developments. Also, processing the prepared data for TPG becomes much simpler.

Variable TPG and Darcy
Flow. Based on the simulator development, the cases, respectively, considering the variable TPG and the classical Darcy model are simulated to compare the production indexes. As shown in Figures 12 and 13, the variable TPG is equivalent to adding an equivalent resistance to the displacement process leading to the deterioration of the final oil recovery. Additionally, the existence of TPG aggravates the fingering phenomenon of the water phase resulting in the sharp increase in water cut. Therefore, the variable TPG lowers the waterflooding developmental efficiency.  According to the history matching data, it can be seen in Figures 14 and 15 that the oil production rate using the commercial simulator is always higher than the actual field data and the corresponding water cut is always lower. However, considering the variable TPG in the self-developed simulator, the history matching for the well performance improves obviously. Therefore, ignoring the variable TPG will have a negative impact on the fitting results.
To further intuitively distinguish the main waterflooding percolation channels and the heavy-oil retention area, the effective displacing pressure gradient is led into the selfdeveloped simulator, which is defined as the difference between the displacing pressure gradient and the variable TPG. The calculated corresponding distribution is shown in Figure 16. The gray areas illustrate the ineffective displacement area, provide an important guidance for the studied reservoir to research the heavy-oil flow characteristic, and optimize the adjustment developmental programs.

Conclusions
For heavy-oil waterflooding development, the TPG caused by the high viscosity could significantly affect the residual oil distribution and the performances of producers, such as oil production rate, water cut, and producing pressures. In this study, the relationship between the oil flow capability and the TPG is used to dynamically describe the heavy-oil flow behavior. The conclusions from this study are listed as follows: (1) A variable TPG model considering heavy-oil mobility is proposed on the basis of the experimental data analysis (2) With the validation of the self-developed simulator, the proposed flow model is employed in the simulation to effectively improve the description of the heavy-oil flow characteristic and significantly enhance the prediction of the well production index  10 Geofluids P o , P g , P w : Pressure of oil, gas, and water (MPa) Φ o , Φ g , Φ w : Flow potential of oil, gas, and water (MPa) μ o , μ g , μ w : Viscosity of oil, gas, and water (mPa·s) D, g: Depth of the reservoir (km), gravity acceleration constant (m/s 2 ) G: Threshold pressure gradient (MPa/m) f , μ oa : Flow correction factor (dimensionless), apparent viscosity of oil phase (mPa·s) a, b: Flow experiment fitting parameters (f) ρ osc , ρ gsc , ρ wsc : Fluid density of oil, water, and gas under standard conditions (kg/m 3 ) q o , q g , q w : Oil, gas, and water inflow or production in each unit of time and volume of rock (kg/s) S o , S g , S w : Saturation of oil, gas, and water (f) B o , B g , B w : Volume factor of oil, gas, and water (dimensionless) R so , R sw : Dissolved gas-oil ratio, dissolved gas-water ratio (dimensionless) Φ, t: Porosity (f), time (s) P cow , P cog : Capillary force between oil and water, capillary force between oil and gas (MPa) x w , y w , z w : Well point coordinates P wf : Bottom hole flowing pressure (MPa) l: Number of iterations (f) V b : Grid block volume (m 3 ) δ: The difference between the l + 1th iteration and the lth iteration of P o and S w at n + 1 time steps n, Δt: Nth time step, time step (s) r e , r w : The equivalent supply radius of the grid where the well is located (m), shaft radius (m) h, r: Reservoir thickness (m), distance from the production well (m) Q l : Stable flow production of l-phase fluid between the well and grid in a certain time step (m 3 ) v l : Stable flow velocity of l-phase fluid between the well and grid in a certain time step (m 3 /s) Ф lwf : Bottom hole flow potential of l-phase fluid (MPa) Φ l : The average flow potential of l-phase fluid in the equivalent supply radius (MPa) α c : Volume conversion factor (α c = 5:614583) m: Accumulate to step m n x , n y , n z : Number of grids in x, y, and z directions.

Data Availability
In this paper, the data in Figures 1 and 2 are obtained by the experimental method. The data in Figures 6-8 and 10-14 are obtained by comparing the simulator developed by ourselves with ECLIPSE. Figures 15 and 16 are the actual mine data. These data are provided in the supporting file.

Conflicts of Interest
The authors declare that they have no conflicts of interest.