Numerical Investigations of the Effect of Nonlinear Quadratic Pressure Gradient Term on a Moving Boundary Problem of Radial Flow in Low-Permeable Reservoirs with Threshold Pressure Gradient

The existence of a TPG can generate a relatively high pressure gradient in the process of fluid flow in porousmedia in low-permeable reservoirs, and neglecting the QPGTs in the governing equations, by assuming a small pressure gradient for such a problem, can cause a significant error in predicting the formation pressure. Based on these concerns, in consideration of the QPGT, a moving boundary model of radial flow in low-permeable reservoirs with the TPG for the case of a constant flow rate at the inner boundary is constructed. Due to strong nonlinearity of the mathematical model, a numerical method is presented: the system of partial differential equations for the moving boundary problem is first transformed equivalently into a closed system of partial differential equations with fixed boundary conditions by a spatial coordinate transformation method; and then a stable, fully implicit finite difference method is used to obtain its numerical solution. Numerical result analysis shows that the mathematical models of radial flow in low-permeable reservoirs with TPGmust take theQPGT into account in their governing equations, which ismore important than those of Darcy’s flow; the sensitive effects of the QPGT for the radial flow model do not change with an increase of the dimensionless TPG.


Introduction
Due to a continuously decreasing crude oil output from conventional reservoirs and a high international gasoline price in recent years, unconventional reservoirs such as lowpermeable reservoirs and shale oil and gas reservoirs have become urgent development resources in the petroleum industry.Consequently, considerable attention has been paid to the relevant research on the kinematic principles for the fluid flow behavior in these unconventional reservoirs [1][2][3][4][5] at present.Abundant experimental and theoretical analyses [6][7][8][9][10][11][12][13] have demonstrated that the fluid flow, in low-permeable porous media, does not obey the classical Darcy's law: the seepage velocity is not proportional to the formation pressure gradient, and there exists a threshold pressure gradient (TPG) ; the fluid flow happens only if the formation pressure gradient is larger than TPG.
Much research on these relevant moving boundary models has been conducted [14][15][16][17][18][19][20][21][22].The computed formation pressure distributions corresponding to these moving boundary problems of the fluid flow in the porous media with TPG show big difference from the ones based on Darcy's law (see Figure 1): the formation pressure gradient is much steeper, it decreases until up to zero at a certain value of a dimensionless distance from a well, that is, the position of a moving boundary, and the pressure distribution curve shows a property of compact support [21]; whereas, for Darcy's flow problem, the formation pressure drop can propagate to any infinite distance transiently according to the exact analytical solution [20], and the formation pressure gradient is much more smooth.Actually, the pressure distribution difference between Darcy's flow and fluid flow in the porous media with TPG can be explained through the angle  between the dimensionless formation pressure curve and the  dimensionless distance at the place of moving boundary, as shown in Figure 1: for the case of the existence of TPG, the tangent of  is equal to the derivative of pressure with respect to distance at the place of moving boundary, that is, the TPG; and, as  decreases gradually and is equal to zero, it becomes a limit case, which corresponds to Darcy's flow; that is, TPG is equal to zero.And if a constant value of TPG is given, the tangent of  will not change transiently.Therefore, it can be concluded that the formation pressure gradient is relatively much higher for the fluid flow in porous media with large TPG.The main physical reason is that the presence of TPG can make the formation pressure drop propagate more slowly, which causes a larger pressure gradient in a relatively shorter pressure disturbed distance from a production well.
It is well known that the rock porosity and fluid density are dependent on the formation pressure, and their formula are commonly expressed in the nonlinear forms of exponential functions [23].Therefore, in the mathematical modeling, the deduced governing equation, by substituting these equations of state and a kinematic equation into a mass balance equation, always contains a nonlinear quadratic pressure gradient term (QPGT), where the coefficient of this term is proportional to the fluid compressibility.In general, the nonlinear QPGT is usually neglected by assuming small fluid compressibility or a small formation pressure gradient [24].The error by neglecting the nonlinear term may be acceptable for most routine engineering applications in the development of conventional reservoirs.However, it has been realized that the linearization is not applicable for large values of time [25,26], and certain operations such as hydraulic fracturing, high injection or production rates of wellbore, and large pressure pulse testing can generate a high formation pressure gradient.
For conventional models of fluid flow in porous media based on Darcy's law, Odeh and Babu [27] studied the analytical solutions of the nonlinear model of slightly compressible fluid flow through porous media with consideration of QPGT by a Laplace transform and pointed out that the nonlinear solutions showed the difference in the pressure change for injection and pumping conditions.Wang and Dusseault [28] presented an analytical solution for pore pressure coupled with deformation in a porous medium by taking the quadratic gradient term into account, and the deviations from the existing solutions were identified in cases of high pressure gradients.Chakrabarty et al. [29,30] gave the analytical dimensionless pressure solutions for radial flow systems by a Laplace transform and concluded that the linearized analysis by neglecting QPGT could lead to serious errors in some cases such as high injection rates.Braeuning et al. [31] studied the effect of QPGT on the variable-rate well-tests, and concluded that the linearization error depended on the magnitude of wellbore damage, pseudoskin, and a nonlinear flow parameter.Tong et al. [25,26] presented the exact analytical solutions of nonlinear transient flow models and dual-porosity models including QPGTs by the generalized Weber transform and Hankel transform; it was concluded that the effect of the QPGT should be considered in large time well testing in reservoir engineering.Li et al. [32] took the effects of both QPGT and wellbore storage into account to build a mathematical flow model in a fractal multilayer reservoir and presented its analytical solution by a Laplace transform.Dewei et al. [33] obtained the analytical solution of a mathematical model of transient seepage flow including QPGT by a Laplace transform, and they demonstrated that QPGT could not be ignored for the unsteady flow model with consideration of a wellbore skin effect.Bai et al. [24] constructed a nonlinear dual-porosity model incorporating QPGT in fracture space and obtained its analytical solution using a Hankel transform, and their analysis showed that the presented nonlinear model appeared to be suitable for simulating naturally fractured reservoirs, subjected to a high injection or production rate or significant fracture compressibility.Nie et al. [34,35] studied nonlinear flow models with a QPGT for both a double-porosity reservoir and a tripleporosity reservoir and obtained a solution through a variable substitution for linearization; it was found that the influence of QPGT was very distinct, and the parameter values of nonlinear model explanation law were much accurate than those of linear model explanation.Yao et al. [36] established the mathematical model of transient flow in a double-porosity fractal reservoir by considering QPGT and solved the model by a Laplace transform; it was demonstrated that the relative errors, caused by ignoring QPGT in the constructed model, may amount to 10%.Guo and Nie [37] discussed the origins of nonlinear flow issues in underground formations by using quadratic pressure gradients to deduce diffusion equations and concluded that the nonlinear models could describe fluid flow in underground formations realistically and accurately.
As far as we know, for the moving boundary problems of fluid flow in porous media with TPG, which are based on modified Darcy's law [8] and have been widely involved in engineering applications for low-permeable reservoirs, the effect of QPGT has not yet been taken into account in their governing equations [14][15][16][17][18][19][20][21][22].However, this nonlinear QPGT may play an important role in the temporal and spatial pressure profiles due to a relatively high formation pressure gradient, generated in fluid flow in porous media with TPG.What is more, in modern times, with the development of advanced analysis methods and improved resolutions of pressure measurement machines [29], it is necessary to study quantitatively the effect of QPGT on these moving boundary problems.
Hence, based on these concerns, in this paper, QPGT is incorporated in the governing equation for a basic moving boundary problem of radial flow in an infinite reservoir with TPG for the case of a constant flow rate at the inner boundary.Owing to the existence of the moving boundary and the nonlinearity of partial differential equations, it is really difficult to obtain its exact analytical solution.Here, we adopted a spatial coordinate transform based finite difference method to solve the nonlinear moving boundary model for the radial flow case.The numerical method has been strictly verified for solving such moving boundary problems with good accuracy by a one-dimensional flow case in the recently published literature [21].Furthermore, the effect of this QPGT on the numerical solutions of the moving boundary problem can be discussed and analyzed quantitatively, by using the numerical results with respect to different values of a dimensionless TPG.

Mathematical Modeling by
Considering QPGT The problem considered involves radial flow in an infinite reservoir with TPG for the case of a constant flow rate at the inner boundary; the porous medium is homogeneous, isotropic, and isothermal; the single-phase horizontal flow does not have any gravity effect; for a basic problem, skin effect and wellbore storage are not considered here; finally, the fluid and porous medium are slightly compressible.The fluid density is as follows [20]: where  is the fluid density, kg⋅m −3 ;   is the initial fluid density, kg⋅m −3 ;   is the initial pressure, MPa;  is the pressure, MPa; and   is the compression coefficient of the fluid, MPa −1 .The porosity of the porous medium is as follows [20]: where  is the porosity of the porous medium, fraction;   is the initial porosity, fraction; and   is the compression coefficient of the porosity, MPa −1 .
The modified Darcy's law for the fluid flow in the porous medium with TPG is as follows [8]: where  is the permeability of the porous medium, 10 −3 ⋅m 2 ;  is the fluid viscosity, mPa⋅s; r is the radial distance, m;  is the seepage velocity, m⋅s −1 ; and  is the TPG, MPa⋅m −1 .The continuous equation for the radial flow in the porous medium is as follows [15]: where  is time, h and  is the moving boundary, m.Substituting ( 1)-( 3) into (4), the governing (mass balance) equation for the radial flow in the low-permeable reservoir, considering the nonlinear QPGT, can be deduced as follows (see Appendix): where   is the total compression coefficient, MPa −1 .The initial condition is as follows: The inner boundary conditions with constant flow rate are where  is the constant flow rate, m 3 ⋅d −1 ; h is the reservoir thickness, m;   is the wellbore radius, m; and B is the volume factor, dimensionless.
The moving boundary conditions are: The moving boundary conditions (8) physically mean that the seepage flow only happens at the area near the wellbore inside the moving boundary, where the formation pressure gradient is larger than TPG; however, outside the moving boundary, the formation pressure gradient is smaller than TPG (equal to zero), so there is no seepage flow behavior, and the formation pressure also keeps the initial pressure; in addition, on the moving boundary, the pressure gradient is just equal to TPG; with the time increasing, the moving boundary also moves outside gradually.The existence of moving boundary is the main difference between the models of fluid flow in porous media with TPG and the classical Darcy's flow models.
Equations ( 5)-( 8) together form the moving boundary problem of radial flow in the infinite reservoir with TPG for the case of a constant flow rate at the inner boundary, in consideration of QPGT.
Define the following dimensionless variables: where   is the dimensionless radial distance,   is the dimensionless time,   is the dimensionless pressure,   is the dimensionless compressibility,   is the dimensionless TPG, and  is the dimensionless moving boundary.Consider Equations ( 10)-( 15) form the dimensionless mathematical model, considering QTPG, for the radial flow in the infinite reservoir with TPG for the case of a constant flow rate at the inner boundary.
From (15), we have Differentiating two sides of ( 16), with respect to   , we have Substituting ( 14) into (17) yields Letting   = (  ) on both sides of (10) and then substituting (14) yield By ( 18) and ( 19), the velocity of the moving boundary can be written as follows: From (20), it can be seen that the existence of the QPGT (the dimensionless compressibility   is not equal to zero) can slow down the velocity of the moving boundary of the radial flow in low-permeable reservoir with TPG.

Numerical Method
Due to the existence of the moving boundary in the dimensionless mathematical model, the flow region is not fixed and expands outward continuously with time increasing.In order to overcome this difficulty in the space discretization for the transient flow region with the moving boundary, a spatial coordinate transformation method is used as follows [21,38]: Through (21), the dynamic flow region for the moving boundary problem [0, (  )] can be transformed to a fixed region [0, 1], and the dimensionless pressure   (  ,   ) can be transformed as a new unknown function  of two variables  and   , that is, (,   ) equivalently; correspondingly, the differential variables can be transformed, respectively, as follows: Substituting ( 22)-( 24) into (10) yields Substituting ( 22)-( 24) into ( 13)-( 15) and ( 20), respectively, yields =1 = 0, From ( 26), we have Substituting ( 29) and ( 30) into (25) From ( 21), ( 11) can be transformed into From ( 26) and ( 27), the following equation can be obtained [21]: Equations ( 31)-( 33) and ( 28) together form a closed system of partial differential equations with the fixed boundary conditions with respect to   (,   ), which is equivalent to the dimensionless mathematical model, that is, ( 10)- (15).The new partial differential equation (31) shows strong nonlinearity, indirectly indicating the strong nonlinearity of the original, untransformed moving boundary problem incorporating the QPGT.Here, a stable, fully implicit finite difference method [21,39] is used to obtain its numerical solutions.The first derivative is replaced by a first-order forward difference; the second derivative is replaced by the second-order central difference [21], and then the difference equation for (31) can be written as follows: where  denotes the total number of spatial grid subintervals with the same length; Δ is the length of a grid subinterval, which is equal to 1/N; i denotes the index of the spatial grid Mathematical Problems in Engineering from the well;  denotes the index of a time step; and Δ  denotes the time step size [21].From ( 28), we have Then, from ( 34) and ( 35), the difference equation corresponding to the ( − 1)th spatial grid can be expressed as follows: The difference equation of ( 33) is as follows: From (32), the initial conditions are obtained as follows: 0  = 0,  = 0, 1, 2, . . .,  − 1.
The difference equation of ( 21) is The difference equation of ( 30) is Substituting ( 40) into (39) yields By (41), numerical solutions of (,   ) can be transformed as the ones of   (  ,   ) in the process of numerical solutions at every time step [21].From (41), it can be further concluded that the presented spatial coordinate transformation lets the time-dependent space discretization in the spatial coordinate of   be transformed as a time-independent space discretization in the new spatial coordinate of y and then makes the numerical solutions by the finite difference method more applicable and simpler [21].2-4 show the effect of the QPGT on numerical solutions of the model, with respect to the dimensionless formation pressure distribution, dimensionless transient wellbore pressure, and dimensionless transient distance of moving boundary, respectively, under different values of the dimensionless TPG.From Figures 2-4, it can be clearly seen that with an increase of the dimensionless TPG, the effect of the QPGT on the mathematical model solutions become more and more obvious; the reason can be explained as follows: the bigger the TPG, the steeper the formation pressure gradient, which generates the solution deviations, resulting from neglecting the QPGT in the governing equation, more seriously.Moreover, from Figures 2-4, it can also be indicated that the dimensionless transient wellbore pressure for the case of neglecting the QPGT is larger than that for the case of considering the QPGT; and the dimensionless transient distance of moving boundary for the case of neglecting the QPGT is also larger than that for the case of considering the QPGT.The QPGT can slow down the velocity of the moving boundary in low-permeable reservoir with TPG.

Effect of QPGT under Different Values of TPG. Figures
Tables 1 and 2 show the selected data of calculated relative errors   from the already computed dimensionless formation pressure distribution and dimensionless transient wellbore  pressure between the two cases considering and not considering the QPGTs, under different values of dimensionless TPG, which correspond to Figures 2 and 3, respectively.The relative error   is equal to the absolute error divided by the magnitude of the exact value, that is, the solutions of the models considering the QPGT.In order to clearly figure out the change behavior of relative errors, the curves regarding the relative errors by the data from Tables 1 and 2 are also plotted in Figures 5 and 6, respectively.
From Figure 5, it can be clearly indicated that, for any case with the same value of the dimensionless distance, the larger the dimensionless TPG, the larger the relative error of the dimensionless formation pressure; furthermore, the larger dimensionless TPG can make the relative errors grow more quickly with an increase of dimensionless distance.For the cases with the values of the dimensionless TPG   = 0.034, 0.067, and 0.1, the relative errors have largely exceeded 5% in the whole region.
From Figure 6, it can be clearly indicated that, for any case with the same value of dimensionless time, the larger the dimensionless TPG, the larger the relative error of the dimensionless transient wellbore pressure; and the larger dimensionless TPG can make the relative errors grow more    quickly with an increase of dimensionless time.For the cases with the values of the dimensionless TPG   = 0.067 and 0.1, the relative errors have largely exceeded 5% in the whole time.
From Figures 5 and 6, it can also be concluded that the relative errors corresponding to the case with the value of the dimensionless TPG   = 0.034, which is relatively close to Darcy's flow, stay a lower level for the change of relative errors of the dimensionless transient wellbore pressure with the increase of dimensionless time, whereas the change of relative errors of the dimensionless formation pressure distribution with the increase of dimensionless distance is still very considerable.
In conclusion, the QPGT has a great effect on the mathematical model solutions of radial flow in low-permeable reservoirs with TPG, especially for a large value of the dimensionless TPG; the larger the dimensionless distance and dimensionless time, the bigger the effect of the QPGT on the dimensionless formation pressure and the dimensionless transient wellbore pressure.In summary, the mathematical models of radial flow in low-permeable reservoirs with TPG should take the QPGT into account in their governing equations, which is more important than those of Darcy's flow.

Conclusions
The existence of a TPG can lead to relatively steep formation pressure gradients, and, thus, it is not appropriate to neglect the QPGT for the fluid flow in porous media with the TPG.
In consideration of the QPGT in the governing equation, a numerical method is presented for investigation of the effect of QPGT on a moving boundary problem of radial flow in low-permeable reservoirs with TPG for the case of a constant flow rate at the inner boundary.Numerical result analysis shows that the QPGT plays an important role in the model solutions with respect to the formation pressure distribution, transient wellbore pressure, and transient distance of the moving boundary.The mathematical models of radial flow in low-permeable reservoirs with TPG should take the QPGT into account in their governing equations, which is more important than those of Darcy's flow.Besides, the sensitive effects of the QPGT on the numerical solutions for the radial flow model do not change with an increase of the dimensionless TPG.The presented research supports theoretical foundations for further improving technologies of well testing and numerical simulation in developing lowpermeable reservoirs.

Figure 1 :
Figure 1: Comparison of computed dimensionless formation pressure distribution curves.

)Figure 2 :
Figure 2: The effect of QPGT on dimensionless formation pressure distribution under different values of dimensionless TPG.

Figure 3 :
Figure 3: The effect of QPGT on dimensionless transient wellbore pressure under different values of dimensionless TPG.

Figure 4 :Figure 5 :
Figure 4: The effect of QPGT on dimensionless transient distance of moving boundary under different values of dimensionless TPG.

Figure 6 :
Figure 6: Relative error curves for computing dimensionless transient wellbore pressure by neglecting quadratic pressure gradient term.

)Figure 7 :
Figure 7: The sensitive effect of QPGT on dimensionless formation pressure distribution.

Figure 8 :
Figure 8: The sensitive effect of QPGT on dimensionless transient wellbore pressure.
to cancel the variables /  and  yields

Table 1 :
Comparison data for the computed dimensionless formation pressure distribution and relative errors.Computed dimensionless formation pressure distribution   (  , 1 × 10 5 )

Table 2 :
Comparison data for the computed dimensionless transient wellbore pressure and relative errors.