# Variable-permeability well-testing models and pressure response in low-permeability reservoirs with non-Darcy flow.

1. Introduction

In 1973, the consolidation problem with initial gradient was investigated, and the approximate solution of the moving boundary model with threshold pressure gradient was obtained (Schmidt and Westman, 1973; Pascal, 1973). Then several modified expressions of Darcy's law were developed (Huang, 1998), based on the flow rate-pressure gradient curves of core flood experiments. Feng and Ge (1985, 1988) investigated the low-velocity non-Darcy flow mechanisms in single-porosity and dual-porosity reservoirs and derived the well testing solution with wellbore storage and skin effect based on regular constant-production model of single phase flow. By employing the moving boundary method of the thermal conduction, Wu (1990) presented the second order integral solution of Bingham fluid flow with threshold pressure gradient. Based on the Gringarten type curves, Cheng and Xu (1996) further introduced the solution and type curves of apparent wellbore radius model. Many other researchers have studied on well testing problem and got many modified mathematical models, solutions and interpretation methods (Li and Liu, 1997; Warren, 1993; Hegeman et al., 1993; Liao and Lee, 1993; Luo and Wang, 2011; Guo et al., 2005; Liu et al., 2014).

All the researchers discussed above considered the existence of threshold pressure gradient in low-permeability reservoirs. However, the practical experiments indicate that the flow curves are non-linear at very low velocity, and the existing low-velocity non-Darcy model is only an approximation, and it cannot accurately represent the flow behavior in low permeability reservoirs (Zhou et al., 2002). So developing a well testing interpretation model and the analytical method, by actual flow, will improve the accuracy of the well testing analytical interpretations significantly. This paper presents the concept of variable-permeability in the well testing of low-permeability reservoirs; meanwhile, wellbore storage effect and skin effect are taken into consideration.

2. One-dimensional variable-permeability well-testing model

2.1. Variable-permeability effect

The non-Darcy flow process in low permeability reservoir is shown in Figure 1. The linear section (above point [J.sub.n]) obeys Darcy's law, and there are mainly three methods to describe the curved section (Huang, 1998): (1) simplifying into a straight line passing through the origin with a slope different from the straight line above point [J.sub.n], which is simple but cannot reflect the flow accurately especially in low permeability reservoirs; (2) power law expression, which is accurate but difficult for mathematical simulation; (3) converting into linear relationship with threshold pressure gradient, reflecting the threshold pressure but resulting in lower velocity in big pores under low pressure gradient condition. The third method, threshold pressure gradient model, is used in most cases in recent years; however, there exist some deviations based on this model.

This paper develops a variable-permeability model by dividing the curved section from [J.sub.1] to [J.sub.n] into many linear sections with different slopes and assuming that Darcy's law can be applied in these linear sections. In this way, the simulation can describe the percolation more accurately, which is close to power law method but easy to calculate. According to the definition of permeability, the slope of the curve at each point stands for the permeability at the corresponding pressure gradient. The permeability rises with the increase of the pressure gradient under a low pressure gradient condition while it maintains constant when the pressure gradient is higher than a critical value. This phenomenon is called the variable-permeability effect.

[FIGURE 1 OMITTED]

In Figure 1, the slope at any point stands for the permeability at the correlated pressure gradient. In this way, the relationship between permeability and pressure gradient can be obtained, as shown in Figure 2. Thus, permeability is converted from a physical meaning into a mathematical meaning, and the flow simulation can be solved by the mathematical method. The processing procedures are as follows:

(i) Through the initial experimental data, plot the flow velocity versus pressure gradient curve V-[DELTA]P, then plot the K-[DELTA]P curve presenting the relationship between permeability and pressure gradient;

(ii) Divide the curve into several linear sections, and consider the permeability of each linear section as a constant;

(iii) The finite difference algorithm is employed to discretize the non-linear flow differential equations, and then set up a non-linear equation system in which the pressure is the unknown, and the permeability is a function of pressure gradient, K= K (dp/dr). An alternating iterative method is used to solve the equations: during the iteration process, initialize the permeability first and make it linear, obtain the pressure solution using iteration method to solve the equations. Get the pressure gradient in each section, and then obtain the permeability of each section according to the plot of permeability versus pressure gradient. Substitute this permeability into the non-linear flow equations and make it linear again. Circulate the alternating iterative process until the pressure tends to be stable and satisfies the accuracy requirement. This result is the pressure at the correlated time step, and then go to calculate the next time step.

[FIGURE 2 OMITTED]

2.2. Well testing model

2.2.1 Assumptions

(i) A well is located in the center of an infinite, horizontal, homogeneous, isopachous and low-permeability reservoir with constant production rate;

(ii) Single phase compressible fluid flow in the low-permeability reservoir. The flow does not obey Darcy's law, and variable-permeability effect is taken into consideration;

(iii) Wellbore storage effect and skin effect are considered;

(iv) The effect of gravity and capillary pressure is negligible;

(v) The total compressibility is constant.

2.2.2 Mathematical model

According to the assumptions discussed above, the condition equations are as follows:

v = - K (dp/dr)/[mu] [dp/dr] (1)

State equations:

[empty set] = [empty set] +[C.sub.f] (P-[P.sub.0]) (2)

[rho] = [[rho].sub.0] [1 + [C.sub.L] (P-[P.sub.0])1 (3)

Mass conservation equation:

[partial derivative] ([rho] [empty set])/[partial derivative]t + div ([rho]v) = 0 (4)

Substitute Equations 1-3 into Equation 4:

[1/r] [[partial derivative]/[partial derivative]r] (K (dp/dr) x r [partial derivative]P/[partial derivative]r) =[empty set][micro][C.sub.t] [partial derivative]p/[partial derivative]t (5)

By considering the effects of wellbore storage effect and sk in effect, the correlated initial condition, and boundary conditions are:

Initial condition:

P(r,t)[|.sub.t=0] = [p.sub.e] (6)

External boundary condition:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (7)

Inner boundary conditions:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (8)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (9)

Kinematic equation:

Equations 5-9 compose the general non-Darcy flow well-testing model in the low-permeability reservoir with variable-permeability effect under radial coordinate.

2.3. Numerical simulation

2.3.1 Model discretization

For one-dimensional radial flow, the pressure gradient is higher near the wellbore. So, non-uniform grid, with denser grids near the wellbore, is employed for precise calculation, which can solve the non-convergence problem and decrease iteration time. Set r = [r.sub.w][e.sup.x], time step [DELTA]t = T/ (m-1), and

The flow equation was discretized as:

[a.sub.i][p.sup.j.sub.i-1] + [b.sub.i][p.sup.j.sub.i] + [c.sub.i][p.sup.j.sub.i-1] =[d.sub.i], i = 1,2, n-2 (10)

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Inner boundary condition was discretized as:

[b.sub.0][P.sup.n+1.sub.wf] + [c.sub.0][P.sup.n+1.sub.1] = [d.sub.0] (11)

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (12)

External boundary condition was discretized as:

[a.sub.n-1][P.sup.n+1.sub.wf] + [b.sub.n-1][P.sup.n+1.sub.1] = [d.sub.n-1] (13)

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

implicit difference method is employed to discretize Equations 5-9.

Forward elimination and backward substitution algorithm were employed to solve the non-linear equation system composed by Equation 10-13, and pressure distribution and bottom hole pressure (BHP) in reservoirs with time were then obtained.

2.3.2 Iterated misconvergence solution

During numerical simulation, the approximate solution may not converge and may oscillate near the exact solution. The main reason for the iterated misconvergence is that the non-linear flow curve is not divided into proper sections. At [t.sub.n], the pressure gradient is [DELTA][p.sup.i] ([t.sub.n]) , and the correlated permeability obtained is [K.sup.i+1] ([t.sub.n]), then [p.sup.i+1] ([t.sub.n]) and [DELTA][p.sup.i+1] ([t.sub.n]) can be calculated. After that, calculate the new permeability [K.sup.1+2] ([t.sub.n]), similarly [K.sup.1+3] ([t.sub.n]), [K.sup.1+4] ([t.sub.n]), ..., [K.sup.n+1] ([t.sub.n]) can be calculated. However, some problems are likely to appear:

[DELTA][P.sup.i] ([t.sub.n]) = [DELTA][P.sup.i+2] ([t.sub.n]) [not equal to] [DELTA][P.sup.i+1] ([t.sub.n])

[K.sup.i+1] ([t.sub.n]) = [K.sup.i+3] ([t.sub.n]) [not equal to] [K.sup.i+2] ([t.sub.n])

The oscillation of the solution makes it difficult to solve the variable-permeability model, and the geometric mean of permeability is used to solve this problem. Calculate the permeability [K.sup.i+1] ([t.sub.n]) using the pressure gradient [DELTA][p.sup.i+1] ([t.sub.n]) obtained after the ith iteration. ([t.sub.n]) for next iteration is:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Then keep doing this iteration process until max [absolute value of [K.sup.m+1] ([t.sub.n])-K ([t.sub.n])] [less than or equal to] [epsilon].

2.4. Well Testing Analysis and Interpretation

The relationship between flow velocity and the pressure gradient is presented based on three groups of experimental data to investigate the influence of non-Darcy effect (Fig. 3a), where the non-Darcy effect increases gradually in curve 1, 2, 3. Figure 3b shows the relationship between the corresponding permeability and pressure gradient. The model discussed above is used to interpret the field test data (Table 1), and the type curves are further obtained, as shown in Figure 4.

[FIGURE 3 OMITTED]

[FIGURE 4 OMITTED]

In the early section of Figure 4, the type curves in log-log scale coincide as a 45o line, indicating the influence of wellbore storage effect. In the transition section, the pressure derivative curves drop down after reaching the peak values. The peak value depends on the value of [C.sub.D][e.sup.2S]. Larger value of [C.sub.D][e.sup.2S] results in larger peak value and later appearance of the peak. In the later section, the type curves upturn to various degrees. The longer duration of the variable-permeability effect, the higher the type curves upturn; meanwhile, the pressure derivative curve achieves horizontal for Darcy flow since there is no variable-permeability effect.

Compare the following three kinds of type curves: (i) the Gringarten type curves of Darcy flow (Gringarten et al., 1979); (ii) the type curves with threshold gradient model (Cheng and Xu, 1996); (iii) the type curves based on the variable-permeability model of this paper, using the data in curve 1 in Figure 3a. The comparison of the type curves of these three different models is shown in Figure 5. Horizontal sections appear in the pressure derivative curves of infinite and homogeneous reservoirs (type i), and this is the feature of typical radial flow. Compared with the type curve based on gradient threshold model (type ii), the type curves based on variable-permeability model (type iii) result in gentle upturn in the later section, and its calculation method is more stable than that based on moving boundary model with threshold pressure gradient.

[FIGURE 5 OMITTED]

3. Two-dimensional variable-permeability well-testing model

To take the boundary effect into consideration, we discuss the solution of the percolation problem in two-dimensional reservoirs.

3.1. Well testing model and numerical simulation

Taking boundary conditions into account, and the two-dimensional well testing model with variable-permeability effect is as follows:

Control equation:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (14)

Initial condition:

P (x,y,t)[|.sub.t=0] = [P.sub.e] (15)

Inner boundary conditions:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (16)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (17)

External boundary condition with one impermeable boundary:

[partial derivative]p (x, y, t)/[partial derivative]|[sub.y=L] = [p.sub.e] (18)

Other boundary conditions corresponding to other forms can be obtained successively:

Flow control in Equation14 was discretized as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (19)

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

[q.sub.i,j]=q for well point grids, and [q.sub.i,j]=0 for non-well-point grids.

Consider the grids within wells as sink & source items. Since the pressure gradient near the well is large, linearize the non-boundary grids. The inner boundary condition was discretized as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (20)

where

[r.sub.e]= 0. 14 [square root of [ ([DELTA]x).sup.2] + [ ([DELTA]y).sup.2]] for rectangle grids, and [r.sub.e] = 0.208[DELTA]x for square grids.

External boundary conditions were discretized as:

[P.sub.1,j]=[P.sub.m,j]=[P.sub.e], (j = 1, 2, ...,n) (21)

[P.sub.i,1]=[P.sub.i,n]=[P.sub.e], (i = 1, 2, ...,m) (22)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (23)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (24)

A non-linear system is formed by Equations 19-22. Gauss-Seidel method was employed to solve the above equations, and pressure distribution and bottom hole pressure in reservoirs with time were further obtained.

3.2. Type curves of well testing models with boundaries

Considering different boundary types, we obtained the type curves based on the model discussed above, as shown in Figures 6-8.

Figure 6 demonstrates the influence of non-Darcy effect on type curves in reservoirs with one impermeable boundary and shows that the pressure derivative curve will become horizontal with a value of 0.5 before reaching the impermeable boundary for Darcy flow. After reaching the impermeable boundary, the pressure derivative curve initially upturns and then become horizontal again with a value of 1. For non-Darcy flow, the pressure derivative curves upturn during radial flow period. The severer the non-Darcy effect results in the larger amplitude of the 'concave', the earlier the derivative curves upturn, the longer duration of the upturned radial non-Darcy section, and the weaker the effect of the boundary (nonDarcy effect covers the boundary effect).

[FIGURE 6 OMITTED]

The impact of impermeable boundaries with different shapes and distances on type curves are shown in Figure 7, including parallel faults, orthogonal faults, three faults, and closed boundaries (the dimensionless distances from the well to the faults all equal 500). The four Darcy flow type curves at the bottom of Figure 7 reveal that how fast the pressure derivative curves upturn in the later section is related to the number of the faults. If the distances from the well to the faults are equal, the more the boundaries are, the earlier the curves upturn and the larger the amplitudes are. The pressure derivative curves of closed boundary reservoirs rise with a slope unity line in the later section. If there are two faults near the well, the orthogonal faults show the second radial flow plateau and the pressure derivative curve of the parallel faults rise with a 0.5 slope line. The non-Darcy flow well testing derivative curves show that pressure derivative curves in the radial flow section upturn due to non-Darcy effect, lying above the 0.5 line of Darcy flow. In the later section, the pressure derivative curves upturn further due to the impact of the faults, later than the Darcy flow. More faults result in larger amplitude upturn. If there are closed boundaries near the well, the derivative curves of Darcy flow and non-Darcy flow coincide with each other and show a slope unity line. The variable-permeability well testing curves of four different faults reveal that more faults lead to the earlier upturn of the non-Darcy well testing curves. The non-Darcy pressure derivative curve of orthogonal faults upturns less obviously. Therefore, to detect the distance from the well to the faults, the testing time should be appropriately extended when there are fewer faults near the well.

[FIGURE 7 OMITTED]

Figure 8 demonstrates the influence of non-Darcy effect on type curves in a reservoir with a constant pressure boundary. In the later section, the pressure derivative curves initially upturn and then drop down after the second hump. It drops down faster and faster until the decline rates reach that of Darcy flow. The larger the amplitudes of the concave parts are (the more severe the non-Darcy effect is), the higher the second humps are and the larger the rate of decline of the derivative is at the end of the curves.

[FIGURE 8 OMITTED]

4. Conclusions

This paper employs the concept of variable-permeability effect and develops one-dimensional and two-dimensional non-Darcy well testing models with variable-permeability and boundaries in low permeability reservoirs. Forward elimination and backward substitution algorithm are used to address the difference equation system of the models. The geometric mean of permeability was employed to solve the oscillation of the solution, and the reliable solutions are obtained. The general non-Darcy flow well-testing model based on variable-permeability effect and its calculation method are more stable and more practical than those based on threshold pressure gradient.

The type curves of well testing with variable-permeability effect are further investigated. In the transition section, the pressure derivative curves drop down after reaching the peak values, and the peak value depends on the value of [C.sub.D][e.sup.2S]. A larger value of [C.sub.D][e.sup.2S] results in larger peak value, sharper dropping down rate and later appearance of the peak. In the later section, the type curves upturn. The longer duration of the variable-permeability effect, the higher the type curves upturn.

For the well-testing model with one impermeable boundary, the radial flow section upturns again after first upturning under the influence of boundary. The pressure derivative curves upturn along with the increasing of the number of boundaries under the influence of non-Darcy effect. When the impermeable boundary is far away from the well, the curves are close to the non-Darcy type curves of the infinite reservoir and the boundary effect is not evident. The more the boundaries are, and the less the distance from the well to the boundary is, the earlier the non-Darcy curves upturn and the large the amplitudes are. For constant pressure boundary, the fast the derivative curves drop down after the second hump.

Acknowledgements

This work was sponsored by National Science and Technology Major Projects (No. 2016ZX05009-004, No. 2016ZX05013002-005) and National Natural Science Foundation of China (No. 51304223).

References

Cheng, S. Q. and Xu, L. X. (1996). Type curve matching of well test data for non-Darcy flow at low velocity (in Chinese). Petroleum Exploration & Development, 23 (4), 50-53.

Darcy, H. (1856). Les Fontaines Publiques de la Ville de Dijon. Victor Dalmont, 304-311.

Das, A.K. (1997). Generalized Darcy's law including source effect. Journal of Canadian Petroleum Technology, 36 (6), 57-59. DOI: http:// dx.doi.org/10.2118/97-06-06

Feng, W. G. and Ge, J.L. (1985). Non-Darcy flow at low velocity through single medium and double media (in Chinese). Petroleum Exploration & Development, 1, 56-62.

Feng, W.G., and Ge, J.L. (1988). Influence of after flow and skin effect for the non-Darcy flow at low velocity through single medium (in Chinese). Petroleum Geology & Oilfield Development in Daqing, 7 (2), 45-50.

Gringarten, A.C., Bourdet, D.P., Landel, P.A. and Kniazeff, VJ. (1979). A comparison between different skin and wellbore storage type-curves for early-time transient analysis. SPE Annual Technical Conference and Exhibition, Society of Petroleum Engineers, Las Vegas. DOI: http://dx.doi.org/10.2118/8205-MS

Guo, D.L., Zeng, X.H., Zhao, J.Z. and Liu, C. Q. (2005). Model and method of well test analysis for wells with vertical fracture. Applied Mathematics and Mechanics (English Edition), 26 (5), 571-578. DOI 10.1007/BF02466330

Hegeman, P.S., Hallford, D.L. and Joseph, J.A. (1993). Well-test analysis with changing wellbore storage. SPE Formation Evaluation, 8 (3), 201-207.

Huang, Y.Z. (1998). The mechanism of low permeability reservoir seepage (in Chinese). Petroleum Industry Press, Beijing, 80-88.

Li, F.H., and Liu, C.Q. (1997). Pressure transient analysis for unsteady porous flow with start-up pressure derivative (in Chinese). Well Testing, 6 (1), 1-11.

Liao, Y.Z, and Lee, W.J. (1993). New solutions for wells with finite-conductivity fractures including wellbore storage and fracture-face skin. Eastern Regional Conference & Exhibition, Society of Petroleum Engineers, Pittsburgh. DOI: http://dx.doi.org/10.2118/26912-MS

Liu, Y.W., Ou-Yang, W.P., Zhao, P.H., Lu, Q. and Fang, H.J. (2014). Numerical well test for well with finite conductivity vertical fracture in coalbed. Applied Mathematics and Mechanics (English Edition), 35 (6), 729-740. DOI: 10.1007/s10483-014-1825-6

Luo, E.H. and Wang, X.D. (2011). A study on transient flow under threshold pressure gradient in dual-pore media with low permeability (in Chinese). China Offshore Oil and Gas, 23 (5), 318-321.

Pascal, H. (1973). Consolidation with threshold gradients. International Journal for Numerical and Analytical Methods in Geomechanics, 5, 1201-1215.

Pascal, H. (1981). Nonsteady flow through porous media in the presence of a threshold gradient. Acta Mechanica, 39, 207-224.

Prada, A. and Civan, F. (1999). Modification of Darcy's law for the threshold pressure gradient. Journal of Petroleum Science and Engineering, 22 (4), 237-240. DOI: 10.1016/S0920-4105 (98)00083-7

Schmidt, J.D. and Westman, R.A. (1973). Consolidation of porous media with non-Darcy flow. Journal of the Engineering Mechanics Division, 99, 1201-1215.

Warren, G.M. (1993). Numerical solutions for pressure transient analysis. SPE Gas Technology Symposium, Society of Petroleum Engineers, Calgary. DOI: http://dx.doi.org/10.2118/26177-MS

Wu, Y.S. (1990). Theoretical studies on non-Newtonian and Newtonian fluid flow through porous media. Ph. D. dissertation, University of California, Berkeley, 174-208.

Yan, Q.L., He, Q.X. and Wei, L.G. (1990). A laboratory study on Percolation characteristics of single phase flow in low-permeability reservoirs (in Chinese). Journal of Xi'an Petroleum Institute, 5 (2), 1-6.

Zhou, Y.Y., Peng S.M. and Li Y. (2002). Fully implicit simulation model for low-velocity non-Darcy flow (in Chinese). Petroleum Exploration & Development, 29 (2), 90-93.

Naichao Feng (1), Shiqing Cheng (1) *, Weixing Lan (2), Guoquan Mu (3), Yao Peng (4), Haiyang Yu (1)

(1.) MOE Key Laboratory of Petroleum Engineering China University of Petroleum Beijing Beijing 102249, China

(2.) E&D Research Institute of Liao He Oil Field Company, Panjin, Liaoning 124010, China

(3.) E&D Institute of Changqing Oil Company, Xi'an, Shanxi 710021, China

(4.) Department of Petroleum Engineering, University of Texas at Austin, Austin 78712, USA

* Corresponding author Shiqing Cheng. E-mail: chengsq973@163.com

Record

Accepted for publication: 05/02/2016

doi: http://dx.doi.org/10.15446/esrj.v20n1.54144
```Table 1. Field test data and reservoir parameters

q ([m.sup.3] / d)    [P.sub.e] (MPa)   [C.sub.t] ([Mpa.sup.-1])

5.0                 15.0             0.0014

q ([m.sup.3] / d)   [empty set]   [micro] (mPa x s)   B   [r.sub.w] (m)

5.0                 0.12          1.0                1.0   0.108

q ([m.sup.3] / d)    s   K ([micro][m.sup.2])

5.0                 0   1.0x[10.sup.3]
```