#### Abstract

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 heavy-oil 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.

#### 1. 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–3]. Heavy oil is a complex fluid that exhibits high viscosity and high density resulting in various difficulties during exploitation [4–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–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–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–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–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 heavy-oil 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 heavy-oil 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.

#### 2. 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.

##### 2.2. 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

##### 3.1. Basic Assumptions

The key of variable TPG numerical simulation is the specific characterization of TPG in the mathematical model of the flow model. The basic assumptions are as follows before establishing the mathematical model: (1)The model considers the dissolution of gas in oil and water, and the water component is only present in the water phase(2)The flow in the reservoir is isothermal(3)The model considers the TPG of the oil phase, but the gas phase and water phase still satisfy Darcy’s law(4)The reservoir rock is slightly compressible and anisotropic, the reservoir fluid can be compressed, and the capillary force and gravity are considered(5)The oil phase and the gas phase are instantaneously phase-balanced(6)The apparent viscosity and the flow correction coefficient of the oil phase are introduced to describe the viscosity anomaly

##### 3.2. Equation of Motion

Theoretical research shows that the flow curve of heavy oil has the characteristics that the arc is not obvious and the transition is fast, so the PTPG flow model is more in line with the flow of the oil phase, while the water and gas phases still meet the Darcy flow model. The equation of motion of the three components can be expressed as

##### 3.3. Flow Differential Equation

###### 3.3.1. Flow Equation

Substituting the equation of motion into the continuity equation to obtain the flow equation, transforming it into the volume conservation form under the ground standard condition, and adding the solution conditions to form a complete mathematical model, one can have (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):

#### 4. Simulator Development

##### 4.1. Numerical Model Solution

The nonlinear characteristics of the solution of variable TPG numerical simulation are very obvious. The conductivity in the direction between grids is used as an example:

In the black oil model, only 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(3)Before the pressure and saturation of the next iteration step are solved, the nonlinear coefficients are calculated from the reservoir physical property parameters calculated in the previous iteration step(4)The Jacobi matrix is assembled and calculated on the basis of nonlinear coefficient calculation, and then, the Jacobi matrix is solved by an efficient algorithm to obtain the value of pressure and saturation(5)Check the convergence conditions and material balance; if it meets the requirements, enter the next period of calculation

##### 4.2. Simulator Development

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. Operation Method of the Simulator

The self-developed simulator is simple in operation and easy to use. It can be realized only by adding the heavy oil keyword when simulating the nonlinear flow of heavy oil in water drive. The specific methods used are as follows:

###### 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 HEAVYEXP 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 and . The application example is shown in Table 1.

###### 4.3.2. 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 pressure gradient for each time step are output. The application example is shown in Table 1.

###### 4.3.3. 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 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.

#### 5. 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.

##### 5.1. 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.

**(a) Saturation field in the third year**

**(b) Pressure field in the third year**

**(c) Saturation field in the sixth year**

**(d) Pressure field in the sixth year**

##### 5.2. 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.

#### 6. 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 self-developed 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 employed to approximately simulate the TPG. In the STATIC case, based on the actual field experimental data, the TPG of the three layers are, respectively, 0.00298 MPa/m, 0.00258 MPa/m, and 0.00236 MPa/m, so the threshold pressures between partition 1 and partition 2, partition 3 and partition 4, and partition 5 and partition 6 are 0.0596 MPa, 0.0516 MPa, and 0.0472 MPa, respectively.

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.

**(a) Saturation field in the third year**

**(b) Pressure field in the third year**

**(c) Saturation field in the sixth year**

**(d) Pressure field in the sixth year**

##### 6.2. 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.

##### 6.3. Field Application

The oil field is a complete shallow and low-amplitude oil reservoir covered with anticline, bottom water, edge water, and sufficient natural energy. The reservoir thickness ranges from 9.1 m to 35.2 m, with an average of 22.8 m. The target area is a medium- and high-porosity and high-permeability reservoir, with an average porosity of 23.9% and an average permeability of .

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 self-developed 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.

**(a) NgII3**

**(b) NgII4**

**(c) NgII6**

**(d) NgII7**

#### 7. 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(3)Based on the definition of the effective displacement pressure gradient, the inefficient displacement area caused by the variable TPG could be observed intuitively. Also, the accuracy of the field history matching can be greatly improved via the aid of the presented flow model

#### Nomenclature

: | Flow velocity of oil, gas, and water (m/s) |

: | Relative permeability of oil, gas, and water (dimensionless) |

: | Density of oil, gas, and water (kg/m^{3}) |

: | Pressure of oil, gas, and water (MPa) |

: | Flow potential of oil, gas, and water (MPa) |

: | Viscosity of oil, gas, and water (mPa·s) |

: | Depth of the reservoir (km), gravity acceleration constant (m/s^{2}) |

: | Threshold pressure gradient (MPa/m) |

: | Flow correction factor (dimensionless), apparent viscosity of oil phase (mPa·s) |

: | Flow experiment fitting parameters (f) |

: | Fluid density of oil, water, and gas under standard conditions (kg/m^{3}) |

: | Oil, gas, and water inflow or production in each unit of time and volume of rock (kg/s) |

: | Saturation of oil, gas, and water (f) |

: | Volume factor of oil, gas, and water (dimensionless) |

: | Dissolved gas-oil ratio, dissolved gas-water ratio (dimensionless) |

: | Porosity (f), time (s) |

: | Capillary force between oil and water, capillary force between oil and gas (MPa) |

: | Well point coordinates |

: | Bottom hole flowing pressure (MPa) |

: | Number of iterations (f) |

: | Grid block volume (m^{3}) |

: | The difference between the th iteration and the th iteration of and at time steps |

: | th time step, time step (s) |

: | The equivalent supply radius of the grid where the well is located (m), shaft radius (m) |

: | Reservoir thickness (m), distance from the production well (m) |

: | Stable flow production of -phase fluid between the well and grid in a certain time step (m^{3}) |

: | Stable flow velocity of -phase fluid between the well and grid in a certain time step (m^{3}/s) |

: | Bottom hole flow potential of -phase fluid (MPa) |

: | The average flow potential of -phase fluid in the equivalent supply radius (MPa) |

: | Volume conversion factor () |

: | Accumulate to step |

: | Number of grids in , , and 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.

#### Acknowledgments

This research was funded by the State Key Laboratory of Offshore Oil Exploitation and CNOOC Research Institute Ltd. (CCL2018RCPS0034RON).