A New Mathematical Model for Pressure Transient Analysis in Stress-Sensitive Reservoirs

For stress-sensitive reservoir, the permeability decreases with the increase of the effective overburden pressure, and pressure transient analysis based on constant rock properties, especially permeability, can lead to significant errors. In this paper, according to the permeability-pressure relationship described by a power function instead of the popularly used exponential function, a new mathematical model for transient fluid flow in stress-sensitive reservoirs is established. The numerical solution is obtained by the fully implicit finite difference method which has been validated by some published analytical solutions before it is used to compute pressure transient responses for stress-sensitive reservoirs. Pressure response curves are plotted and the effects of relevant parameters on both pressure drawdown and buildup responses have been studied.The presentedmodel could provide an alternative method for understanding and predicting the performances for stress-sensitive reservoirs.


Introduction
In practice, the effective overburden pressure of the reservoir will increase with the decrease of the fluid pressure [1].Because the porous media are not always rigid and nondeformable, the increase of the effective overburden pressure may result in plastic deformation or elastic deformation of the reservoir rocks and reduce the permeability, porosity, compression coefficient, and so forth [2,3].To a certain extent, the property of rock is pressure dependent or stress sensitive [4].
The production capacity of oil and gas wells is mainly affected by the permeability in reservoirs [5].In the development of oilfield, the fluid of the reservoir is exploited constantly and the pore pressure in the reservoir decreases gradually.As a result, it could cause the effective overburden pressure to increase, and lead to the compression and deformation of the rocks in the reservoir and a closure of the pores in the rocks.Thus it may cause the permeability to reduce.The changes of permeability inevitably cause the changes of permeable capacity and then production capacity in oil well is affected [6].Therefore, it has been observed that measured flow rates in stress-sensitive reservoirs are sometimes much lower than the predicted ones by transient flow equations based on the assumption of constant rock properties [7,8].
Since the 1950s, the laboratory studies and the mathematical modeling research for stress-sensitive reservoirs have been reported by many investigators.Laboratory experiments have found that the permeability decreases when the effective overburden pressure increases and this phenomenon is even more severe in tight samples [9][10][11][12][13].
In order to describe the transient flow behavior in stresssensitive porous media, several models have been presented.The existing strategies for constructing the stress-sensitive models can be classified into two categories, that is, the pseudopressure approach and the permeability-stress function approach.The pseudopressure approach has been applied to study pressure drawdown, pressure buildup, injection, and falloff testing by several investigators [14][15][16].The main disadvantage of the pseudopressure approach, which may hinder its applications in practice, is that the tabulated properties of rock and fluids versus pressure should be known a priori at each pressure level.In addition, a nonlinear diffusivity term still exists in the diffusivity equation.

Mathematical Problems in Engineering
The permeability-stress function approach, considering the pressure-dependent permeability as a permeability-stress function (e.g., exponential function), is widely used to study the transient flow behavior in porous media through solving a nonlinear diffusivity equation.Based on the exponential function relationship between the permeability and the pore pressure (or the effective overburden pressure), Pedrosa [17] presented the first-order approximate analytical solution for a line-source well producing at a constant rate from an infinite radial flow system using the perturbation technique.Kikani and Pedrosa [18] extended the work of Pedrosa [17] by presenting the second-order approximate analytical solution for the same problem using the same solution methodology and suggested the use of zero-order perturbation solution to investigate the effects of wellbore storage, skin, and outer boundary on the pressure transient response for stress-sensitive reservoirs.Zhang and Ambastha [19] further investigated the effects of stress-sensitive permeability on drawdown and buildup pressure transient behavior using the stepwise permeability model and found that the pressure transient response obtained by Kikani and Pedrosa [18] was not accurate when wellbore storage and skin were taken into consideration, so the numerical method was suggested to study the pressure transient response in stress-sensitive reservoirs.Ambastha and Zhang [20] extended the work of Zhang and Ambastha [19] by presenting three models (i.e., one-parameter model, stepwise permeability model, and two-parameter model) for stress-sensitive reservoirs based on the exponential function (or the modified exponential function) relationship between the permeability and the pore pressure (or the effective overburden pressure).Jelmert and Selseng [21] presented an exponential permeability model through introducing normalized permeability variables and obtained approximate analytical solutions which are computationally simple and readily available.Wu and Pruess [22] presented an integral method for describing the transient flow behavior in stress-sensitive porous media.Zhang and Tong [23] investigated the pressure transient response of the fractal medium in stress-sensitive reservoirs using the self-similarity solution method and the regular perturbation method.Marshall [24] presented a new analytical method for solving the flow equation which contains a squared gradient term and an exponential dependence of the hydraulic diffusivity on pressure.Friedel and Voigt [25] investigated approximate analytical solutions for slightly compressible fluids and real gas flow with constant-rate and constantpressure wellbore boundary conditions using the perturbation technology and Boltzmann transformation.Zhang et al. [26] presented the well test model for stress-sensitive and heterogeneous reservoirs with nonuniform thicknesses and obtained the zero-order approximate analytical solution using the perturbation technology.Zhang et al. [27] extended the work of Zhang et al. [26] to the case of stress-sensitive and heterogeneous dual-porosity reservoirs with nonuniform thicknesses.Qanbari and Clarkson [28] presented a new method to solve the nonlinear partial differential equation for the transient flow of slightly compressible liquid in a stresssensitive formation by introducing a new variable.Ai and Yao [29] presented a new well test model for low-permeability stress-sensitive reservoirs, in which all permeability, porosity, and pore compressibility are pressure dependent.Yi et al. [30] presented a mathematical model of fluid flow in a stress-sensitive reservoir with a horizontal well based on the pressure-dependent characteristics of permeability and porosity.
Although many models for stress-sensitive reservoirs, based on the permeability-stress function approach, have been presented, almost all published models are only confined to the exponential function (or the modified exponential function) relationship between the permeability and the pore pressure (or the effective overburden pressure).Recently, through a lot of actual experimental data, several investigators [12,[31][32][33][34][35][36] have found that the decrease of permeability as a function of the effective overburden pressure can be described by a power function, which usually fits with experimental data better than the exponential function.But surprisingly, although the power function model has served as a good alternative method for describing the relationship between the permeability and the effective overburden pressure, there have been few attempts to introduce the power function model into the seepage flow model for studying the transient flow behavior in stress-sensitive reservoirs.
In this paper, based on the experimental data of actual core samples, we validate that the improved power function model is a good alternative method for describing the permeability-stress relationship in comparison with the oneparameter exponential function model.Then, according to the improved power function model, we present a new mathematical model to study the pressure transient behavior in stress-sensitive reservoirs.The mathematical model is solved by the fully implicit finite difference method which has been validated by some published analytical solutions.Pressure response curves are plotted and the effects of relevant parameters on both pressure drawdown and buildup responses are studied.

Laboratory Measurements of Stress Effect on Permeability.
In order to accurately simulate the stress state of the reservoir rock, the permeability test should be conducted in a three-dimension stress condition.In this experiment, the hydrostatic core holder is usually used to simulate the change of the effective overburden pressure.By keeping the confining pressure constant and changing the fluid pressure in the core through the back pressure valve, the measurement can be conducted at different effective confining pressures to simulate the effective overburden pressure variation.The test steps are as follows.
(1) Determine the initial effective overburden pressure.
The initial effective overburden pressure can be expressed as where  ob can be accurately obtained through the density log or be approximately determined by the following expression: (2) Put the core into the core holder and then adjust the confining pressure and the fluid pressure to set the effective pressure at the initial effective overburden pressure which can be calculated by (1).The permeability at the initial effective overburden pressure can be measured.(3) Keep the confining pressure constant and decrease the fluid pressure in the core.Wait till a stable state is reached in the deformation of the core, and then the permeability can be measured under different effective pressures.(4) Repeat Step (3).The relationship between the permeability and the effective overburden pressure can be obtained.
Five core samples from  field and six core samples from  field were tested, respectively.The permeability-effective overburden pressure curves are shown in Figures 1 and 2, respectively.It is clear that the permeability decreases when the effective overburden pressure increases.

Curve-Fitting Equations for the Permeability-Effective
Overburden Pressure Relation.In order to study fluid flow in stress-dependent porous media, some kinds of mathematical equations are used to describe the relationship between the permeability and the effective overburden pressure.

One-Parameter Exponential Function
Model.The oneparameter exponential function model to represent the relationship between the permeability and the effective overburden pressure is widely used to study the transient flow behavior in porous media.The one-parameter exponential function model is given as [17,18]

Stepwise Permeability Model.
Based on the oneparameter exponential function model, stepwise permeability model to approximately represent the permeability modulus  changing with net confining pressure has been introduced.The stepwise permeability model based on the stepwise permeability modulus is given as [19,20] where (for  = 1, 2, . . ., ) , (5) with  0 =   and (  ) 0 =   .Although the stepwise permeability model seems to be more reasonable than the one-parameter exponential function model for describing the relationship between the permeability and the effective overburden pressure, the stepwise permeability model has not been widely used in practice because the critical pressure is difficult to be determined accurately.

Two-Parameter Exponential Function Model.
Compared with the one-parameter exponential function model, the two-parameter exponential function model to represent stress-sensitive permeability is proposed to allow for the permeability to approach a specified minimum permeability, not necessarily zero, at large effective overburden pressure.The two-parameter exponential function model can be expressed as [20,37] 2.2.4.Power Function Model.In recent years, through a lot of actual experimental data, some scholars have found that a power-function curve-fitting equation can mathematically represent the experimental data for core samples and usually obtain a better correlation between the measured data and calculated values than the exponential-function curve-fitting equation.Therefore, the power function model which is a good alternative method to describe the relationship between the permeability and the effective overburden pressure has received more attention.The power function model is usually written as [12,34] 2.2.5.Improved Power Function Model.Although ( 7) is a good curve-fitting equation used to describe the relationship between the permeability and the effective overburden pressure, ( 7) may hinder the study of the fluid flow in stressdependent porous media, because the value of  depends on the unit of  eff .In order to obtain the dimensionless relationship between the permeability and the effective overburden pressure, the improved power function model is proposed as follows [31][32][33]35]: In what follows, (8) would be used to establish the seepage model for studying the fluid flow in stress-dependent porous media.

The Fitting Comparison between These Kinds of Mathematical Equations.
The experimental data of actual core samples from  field and  field and the ones given by Vairogs et al. [7] are used to analyze the fitting correlation between the measured data and calculated values from the fitting equations, respectively.It should be noted that both the power function model and the improved power function model are based on the power function relation between the permeability and the effective overburden pressure.Therefore, in the following, we focus on the fitting comparison between the one-parameter exponential function model, which has been widely applied in studying the transient flow behavior in porous media, and the improved power function model.
The experimental data of actual core samples from  field and  field are fitted by ( 3) and ( 8), respectively.The fitting parameters (i.e., stress-sensitive coefficient, , and ) and correlation coefficients,  2 , of each core sample are listed in Tables 1 and 2, respectively.As shown in Tables 1 and 2, for core samples from  field and  field, most of correlation coefficients between the measured data and calculated values from the improved power function model are greater than the ones from the one-parameter exponential function model, and the average correlation coefficients between the measured data and calculated values from the improved power function model are also greater than the ones from the one-parameter exponential function model, respectively.
For further analyzing the fitting comparison between one-parameter exponential function model and improved power function model, the measured data reported by Vairogs et al. [7] are used to be fitted by ( 3) and ( 8), respectively.The results including the fitting parameters and correlation coefficients are listed in Table 3, which shows that good agreement between the measured data and calculated values from the improved power function model is obtained for most of core samples and the average correlation coefficient between the measured data and calculated values from the improved power function model is greater than the one from the one-parameter exponential function model.
Consequently, the improved power function model, based on the power function relation between the permeability and the effective overburden pressure, is a good alternative method for describing the permeability-stress relationship in comparison with the one-parameter exponential function model which has been widely used in practice.

3.1.
Assumptions.An isotropic, homogeneous, horizontal, and slab reservoir is bounded by the top and the bottom parallel impermeable planes.The reservoir is filled with a single fluid which is a slightly compressible fluid with constant viscosity.The fluid flow in the reservoir follows Darcy's law, with the influence of gravity force and capillary force being ignored.The initial pressure is assumed to be uniform throughout the reservoir.The formation rock is the sensitivity of permeability to effective stress.Fluid is produced at a constant rate by a finite-radius well with wellbore storage and skin.

Establishment of Mathematical Model.
Based on the improved power function model mentioned above, the mathematical model used to study the transient flow behavior in stress-sensitive porous media is derived in Appendix A. For the convenience of calculation and application, the dimensionless variables, which are defined in Table 4 where all the parameters are explained in the nomenclature, are introduced into the mathematical model.Therefore, the dimensionless mathematical model is as follows.
Dimensionless seepage flow differential equation is as follows: Inner boundary condition for constant-rate production is Outer boundary conditions are the following: lim Based on the mathematical model including ( 9) to (15), several simplified models can be obtained as follows.
(1) When  = 0 and  = 0, the mathematical model is reduced to the conventional radial flow model without considering the effects of the stress sensitivity and the quadratic gradient term, which has been studied for a long time [38,39].
(2) When  = 0 and  ̸ = 0, the mathematical model is reduced to the radial flow model with only considering the effect of the quadratic gradient term, which has attracted attention and has been studied since the 1990s [40][41][42].
(3) When  ̸ = 0 and  = 0, the mathematical model is reduced to the radial flow model with only considering the effect of the stress sensitivity.Because of neglecting the quadratic gradient, in the certain operations, such as hydraulic fracturing, large-drawdown flows, drill-stem test, and large-pressure pulse testing, this simplified model may cause significant error of the predicted pore pressure for stress-sensitive reservoirs.
It should be noted that the present model may not be available for simulating the fluid flows in stress-sensitive reservoirs with lots of very narrow pores (a few nanometer thick pores).That is because the present model is based on Darcy' law, which neglects the microscale effects and is always employed to simulate macroflows in porous media with lots of micrometer pores.For the porous media with lots of nanopores, it is important to consider the microscale effects in simulations.Therefore, the present model cannot be applied to the simulation of fluid flows in porous media with lots of nanopores, while the molecular simulation method can successfully simulate the fluid flows in nanopore [43][44][45].

Solutions of Mathematical Model.
Because the seepage flow differential equation in the mathematical model is a nonlinear differential equation, the analytical solution is hardly obtained.Therefore, in this study, (9) to (15) are solved numerically by finite difference method.The grid blocks in the radial direction are spaced in a geometric fashion.A fully implicit method is used to generate the finite-difference forms of ( 9) to (15), and the Newton-Raphson method is used to solve the resulting system of nonlinear equations.The process of solving the model is shown in Appendix B.
It should be noted that the time step size and space step size are not constants.The time step size and space step size may be increased or decreased according to the convergence and accuracy.
In practice, pressure buildup is widely applied in well test.For pressure buildup response, the "superposition in time" method used for analytical solution of linear differential equation cannot be directly used for stress-sensitive reservoirs, because ( 9) is a nonlinear differential equation.Therefore, in order to obtain the pressure buildup response, the well is directly shut-in, after the specified producing time is reached, by setting the right-hand side of ( 11) and (B.16) to be zero.
It is important to validate the numerical solution before using the numerical simulator to compute the pressure transient response for stress-sensitive reservoirs.In the following, the numerical solution is validated for non-stresssensitive reservoirs by comparing the dimensionless wellbore pressure response obtained by the numerical solution with the one obtained by the well-known Van Everdingen and Hurst's analytical solution [38].The Van Everdingen and Hurst's solution was implemented by Ambastha and Ramey Jr. [39].The analytical solution in the Laplace space should be inverted to the real space using the algorithm proposed by Stehfest [46].outer boundary conditions considering the effects of wellbore storage and skin.The solid lines are analytical results and the scatter points are numerical results.All of them are matched pretty well.Figure 4 shows the case in which the effects of wellbore storage and skin are not considered.The numerical solutions are also in good agreement with the analytic solutions.
In order to obtain the pressure buildup response, the "superposition in time" method is applied in Van Everdingen and Hurst's analytical solution, while the "direct shut-in" method is used in the numerical solution.Figure 5 shows the results of pressure buildup response obtained by the Van Everdingen and Hurst's analytical solution and the numerical solution.A perfect match has been obtained, too.From Figures 3-5, we can conclude that the numerical solution in this study agrees with the Van Everdingen and Hurst's analytical solution very well for both pressure drawdown and buildup response with or without wellbore storage and skin effects for non-stress-sensitive reservoirs.Chakrabarty et al. [40] proposed an analytical solution of the radial flow model considering the effect of the quadratic gradient term which is a simplified model of the proposed stress-sensitive model in this study when  = 0 and  ̸ = 0.The numerical solution is validated by comparing the   6, which also shows an excellent match between the numerical solution and the Chakrabarty et al. 's analytical solution.
Based on the preceding validation efforts, we conclude that the numerical computation method in this study, which has yielded accurate pressure transient response for nonstress-sensitive reservoirs with or without the effect of the quadratic gradient term, can be used to compute the pressure transient response for stress-sensitive reservoirs.

Pressure Transient Characteristics
In this section, we will calculate the dimensionless wellbore pressure (  ) and the derivative (d  /d  ) for a stresssensitive reservoir with the proposed model and numerical computation method.In what follows, the standard loglog typical curves of   and (   ⋅   /  ) versus   /  are obtained, and the effects of relevant parameters on the pressure transient behavior are studied.Basic data used for demonstration in the base case are shown in Table 5.
Figure 7 shows the effect of stress-sensitive coefficient, , on the wellbore pressure transient behavior for an infinite stress-sensitive reservoir.As shown in Figure 7, the entire transient-flow process includes three main flow stages.In early time (stage I called as wellbore storage period), the pressure curve and the pressure derivative curve, which are not affected by stress sensitivity in this flow period, align in a unit slope line.Then the transitional flow period (stage II) and the radial flow period (stage III), in which the stresssensitive coefficient, , has a significant effect on the pressure and derivative curves, can be seen in the typical curves.The positions of pressure and derivative curves ascend with an increasing value of  in the transitional flow period and the radial flow period.
In order to quantify the effect of stress sensitivity on the pressure behavior, we introduce the relative error between the wellbore pressure for non-stress-sensitive reservoirs and the one for stress-sensitive reservoirs with all other parameters kept constant.The relative error is expressed as Figure 8 shows the relative error between the wellbore pressure for non-stress-sensitive reservoirs and the one for stress-sensitive reservoirs with different stress-sensitive coefficients, .As shown in Figure 8, stress sensitivity has no influence on the wellbore pressure transient behavior in early time (stage I called as wellbore storage period).After the wellbore storage period, the relative error increases with time increasing.The magnitude of the relative error is greatly dependent on the stress-sensitive coefficient, .Along with the increase of , the relative error appears higher.In other words, a larger  value could cause a larger deviation of wellbore pressure from  = 0. Figure 11 shows the effect of stress-sensitive coefficient, , on the wellbore pressure buildup for a stress-sensitive reservoir with closed outer boundary.As shown in Figure 11, irrespective of the severity of stress sensitivity, all wellbore pressure buildup responses merge with  = 0 response at late time.It should be noted that the effect of stress sensitivity on the wellbore pressure buildup mainly occurs in early and intermediate time.Before the merger of pressure buildup responses for a stress-sensitive reservoir (i.e.,  ̸ = 0 response) with the responses for a non-stress-sensitive reservoir (i.e.,  = 0 response) in Figure 11, apparent semilog straight lines of slopes higher than that for  = 0 response may be drawn at intermediate time, which may result in underestimated values for initial effective permeability if not considering the effect of stress sensitivity for a stress-sensitive reservoir.The larger the value of  is, the more severe the underestimation of initial effective permeability becomes.

Conclusions
This paper has presented a new mathematical model for studying the pressure transient behavior in stress-sensitive reservoirs based on the improved power function model.The proposed model has been solved by the fully implicit finite difference method.The effects of relevant parameters on both pressure drawdown and buildup responses have been studied.

Mathematical Problems in Engineering
The model presented in this study has provided an alternative method for understanding and predicting the performances for stress-sensitive reservoirs.Several important conclusions can be drawn from this study.
(1) The improved power function model, which is based on the power function relation between the permeability and the effective overburden pressure, provides a good match to the experimental data of actual core samples and could serve as a good alternative method for describing the permeability-stress relationship in comparison with the one-parameter exponential function model.
(2) The numerical computation method proposed in this study, which has been validated by some published analytical solutions, can be used to compute the pressure drawdown and buildup responses for stresssensitive reservoirs.
(3) Irrespective of the severity of stress sensitivity, the pressure curve and the pressure derivative curve always align in a unit slope line in early time.
(4) After the wellbore storage period, stress sensitivity has an important effect on the wellbore pressure transient behavior.The positions of pressure and derivative curves ascend with an increasing value of .The relative error between the wellbore pressure for non-stress-sensitive reservoirs and the one for stress-sensitive reservoirs increases with time and the value of  increasing.
( (6) In all likelihood, a conventional semilog analysis of pressure buildup data from stress-sensitive reservoirs, assumed to be falling in the radial flow period, will underestimate the value of initial effective permeability.The higher the stress sensitivity is, the more severe the underestimation of the initial effective permeability becomes.

A. Mathematical Modeling
The continuity equation for a one-dimensional radial system is given as Darcy's law, under the above assumptions, takes the form The permeability-stress relationship can be expressed by (8), and the fluid compressibility is defined as follows: The pore compressibility is defined as follows: Taking (A.2) to (A.4) and ( 8) into (A.1), the seepage flow differential equation for stress-sensitive reservoirs in a onedimensional radial system is given by the following:

Figure 1 :Figure 2 :
Figure 1: Relationship between permeability and effective overburden pressure for core samples from  field.

Figure 7 :
Figure 7: Effect of stress-sensitive coefficient, , on the wellbore pressure transient behavior for an infinite stress-sensitive reservoir.

Figure 8 :
Figure 8:  The relative error between the wellbore pressure for nonstress-sensitive reservoirs and the one for stress-sensitive reservoirs with different stress-sensitive coefficients, .

Figure 9 :Figure 10 :
Figure 9: Effect of the dimensionless outer boundary radius,   , on the wellbore pressure transient behavior for stress-sensitive reservoirs with closed outer boundary.

Figures 9 Figure 11 :
Figure 11: Effect of stress-sensitive coefficient, , on the wellbore pressure buildup for a stress-sensitive reservoir with closed outer boundary.

)
The outer boundary configurations and the value of   only have influence on the wellbore pressure transient behavior at late time.The start time of the outer boundary reflection is a function of the value of   : with a large   , the outer boundary reflection occurs later.

Table 2 :
The fitting comparison between one-parameter exponential function model and improved power function model for core samples from  field.Sample name   , (10 −15 m 2 ) Improved power function model /  = ( eff / eff ) − One-parameter exponential function model /  = exp[−( eff −  eff )]

Table 3 :
[7] fitting comparison between one-parameter exponential function model and improved power function model for core samples as reported by Vairogs et al.[7].

Table 4 :
Definitions of dimensionless variables.

Table 5 :
Data used in the base case.Skin factor, , dimensionless 2 Density of rock skeleton,   , kg/m 3 2650 Density of formation water,   , kg/m 3 1000 dimensionless wellbore pressure response obtained by the numerical solution with the one obtained by the Chakrabarty et al. 's analytical solution with the effect of the quadratic gradient term.The numerical pressure drawdown solution and the Chakrabarty et al. 's analytical solution for infinite, closed, and constant pressure outer boundary conditions with the effect of the quadratic gradient term are plotted in Figure =   +   .At time  = 0, pressure is distributed uniformly in the reservoir, equal to the initial pressure   .For the convenience of discretizing the mathematical model, let us introduce a new dimensionless space variable, , that is related to the dimensionless radial distance according to  = ln   .Equations (B.2) to (B.8) are discretized by a fully implicit method, and in order to improve the precision of processing, the virtual node will be introduced for the third boundary condition (e.g., virtual node −1 for infinite and constant pressure outer boundary model, and virtual node −1 and virtual node  + 1 for closed outer boundary model).The discrete forms of (B.2) to (B.8) are as follows: