Perturbation Solution for One-Dimensional Flow to a Constant-Pressure Boundary in a Stress-Sensitive Reservoir

Most analyses of fluid flow in porous media are conducted under the assumption that the permeability is constant. In some “stress-sensitive” rock formations, however, the variation of permeability with pore fluid pressure is sufficiently large that it needs to be accounted for in the analysis. Accounting for the variation of permeability with pore pressure renders the pressure diffusion equation nonlinear and not amenable to exact analytical solutions. In this paper, the regular perturbation approach is used to develop an approximate solution to the problem of flow to a linear constant-pressure boundary, in a formation whose permeability varies exponentially with pore pressure. The perturbation parameter αD is defined to be the natural logarithm of the ratio of the initial permeability to the permeability at the outflow boundary. The zeroth-order and first-order perturbation solutions are computed, from which the flux at the outflow boundary is found. An effective permeability is then determined such that, when inserted into the analytical solution for the mathematically linear problem, it yields a flux that is exact to at least first order in αD. When compared to numerical solutions of the problem, the result has 5% accuracy out to values of αD of about 2—a much larger range of accuracy than is usually achieved in similar problems. Finally, an explanation is given of why the change of variables proposed by Kikani and Pedrosa, which leads to highly accurate zeroth-order perturbation solutions in radial flow problems, does not yield an accurate result for one-dimensional flow. Approximate solution for flow to a constant-pressure boundary in a porous medium whose permeability varies exponentially with pressure. The predicted flowrate is accurate to within 5% for a wide range of permeability variations. If permeability at boundary is 30% less than initial permeability, flowrate will be 10% less than predicted by constant-permeability model. Approximate solution for flow to a constant-pressure boundary in a porous medium whose permeability varies exponentially with pressure. The predicted flowrate is accurate to within 5% for a wide range of permeability variations. If permeability at boundary is 30% less than initial permeability, flowrate will be 10% less than predicted by constant-permeability model.


Introduction
Most analyses of fluid flow through porous rock formations are conducted under the assumption that the permeability at any given physical location is constant. In many situations, this is a reasonable assumption. However, in some "stress-sensitive" rock formations, the variation of permeability with pore fluid pressure is sufficiently large that it cannot be neglected. Unfortunately, accounting for the variation of permeability with pore fluid pressure renders the governing pressure diffusion equation nonlinear and thus not amenable to exact analytical solutions.
This type of nonlinear pressure diffusion problem is often attacked by using the pseudopressure (Raghavan et al. 1972;Tabatabaie et al. 2017), which is essentially a Kirchhoff transformation. This has the effect of slightly "lessening" the effect of the nonlinearity, by moving the nonlinearity from the second-order space-derivative term, to the first-order time-derivative term. However, it does not fully linearize the governing PDE, and this final linearization is usually accomplished by simply evaluating the resulting nonlinear coefficient at some suitable "average" pressure.
Another approach that has been used to attempt to analytically account for the stress dependence of the permeability as it appears in the pressure diffusion equation is the perturbation approach pioneered by Kikani and Pedrosa (1991). They started by expressing the dimensionless pressure in terms of a new variable, η, as follows: where γ D is a dimensionless parameter that is proportional to the logarithmic derivative of the permeability with respect to pore pressure. (Note that the notation used in Eq. (1), which is that used by Kikani and Pedrosa, is not the same as will be used in the remainder the present paper.) Equation (1) is essentially a Kirchhoff transformation (Vadasz 2010), although this connection has not usually been made when discussing this problem in the petroleum engineering context. Vadasz (2010) used a Kirchhoff transformation to solve a nonlinear diffusion equation in the context of heat conduction and showed that it is a special case of a Cole-Hopf transformation, which actually has some advantages over the Kirchhoff approach. Kikani and Pedrosa (1991) found that, for the problem of radially symmetric flow to a vertical well under conditions of constant production rate, the 0th-order perturbation solution for their auxiliary function η actually provided a very accurate approximation for the pressure, in the late-time regime. Ren and Guo (2018) extended this approach to the problem of radial flow under conditions of constant well-bottom pressure, as well as to a few other problems. Other researchers have also followed this approach, for problems including dual-porosity and other effects, generally stopping at the 0th-order solution (Wang et al. 2015;Huang et al. 2018). However, as will be seen below, the 0th-order KP-type solution does not work well for the problem of flow to a linear constant-pressure boundary. Instead, a "traditional" perturbation approach will be used, to first-order, without using the change of variables embodied in Eq. (1). (1)

Model Development
In this paper, the problem of one-dimensional flow to a planar boundary held at constant pressure will be examined. This model would apply at early times to flow to a hydraulic fracture in an oil reservoir, for example. It would also apply, at early times, to flow from a matrix block to a nearby fracture, in a dual-porosity reservoir. The governing equation for this problem will be briefly outlined, following the derivation given in standard textbooks such as Zimmerman (2018). The conservation of mass principle for one-dimensional flow through a porous medium can be expressed as where ρ [kg/m 3 ] is the fluid density, q [m/s] is the Darcy flux, and [−] is the porosity of the rock. Darcy's law is used to relate the flux to the pressure gradient: where P [Pa] is the pore fluid pressure, μ [Pa s] is the fluid viscosity, and k(P) [m 2 ] is the permeability, which will be assumed to vary with the pore pressure. Inserting Darcy's law into the conservation of mass equation yields Applying the chain rule and product rule, this equation takes the form Invoking the usual assumption that the fluid compressibility cf [1/Pa] and rock pore compressibility c [1/Pa] can be considered constant over the range of pressures of interest, the above equation reduces to A simple order-of-magnitude estimate shows that the ratio of the second term on the left to the first term on the left is c f ΔP , where ΔP is some characteristic pressure change. For liquid flow in a reservoir, c f will be on the order of a GPa, and ΔP will be at most a few tens of MPa, so this second term will be neglected. (An elegant method for treating cases in which this term cannot be neglected, using a Cole-Hopf transformation, has been presented by Marshall (2009)). Hence, the governing equation for this process is where c = c + c f is the total compressibility of the system. In this model, the only source of mathematical nonlinearity is the pressure dependence of the permeability. Equation (7) is the nonlinear partial differential equation that governs the one-dimensional flow of a slightly compressible reservoir fluid in a stress-sensitive reservoir. If flow takes place to a constant-pressure boundary, then the following initial and boundary conditions apply:

Permeability Model
If the pore fluid pressure decreases, the pores will contract, and the permeability of the rock will tend to decrease (Sisavath et al. 2000). For many rock formations, this effect is small enough to ignore. For some "stress-sensitive" formations, the change in permeability will be too large to be neglected. Numerous mathematical models and correlations have been developed to model the reduction in permeability, as a function of effective stress. A commonly used model (Yilmaz et al. 1991;Kikani and Pedrosa 1991), which fits many data sets (Franquet et al. 2004), is that of an exponential dependence of permeability on pore pressure: where k i is the permeability at some reference pressure P i , and the parameter [1/Pa] is often referred to as the permeability modulus. The reference pressure can conveniently be identified with the initial reservoir pressure.

Integral Equation Formulation
McWhorter and Sunada (1990) developed a method for solving a one-dimensional nonlinear diffusion problem, by reformulating the governing partial differential equation as an integral equation, which can be solved numerically to a high degree of accuracy using an iterative approach. Although their formulation was developed in the context of "unsaturated" flow of water in a two-phase air-water system, it can easily be adapted to the present problem, as well as other problems governed by a nonlinear diffusion equation (Bajwa and Blunt 2016). The modified McWhorter-Sunada approach can be used to provide a benchmark for the perturbation solution. A brief derivation of this formulation, specialized to the current pressure-sensitive permeability model, will now be given. Returning to Eq. (7), a Boltzmann variable is first defined as According to the standard procedure for transforming equations into the Boltzmann domain, Eq. (12) takes the form of the following second-order ordinary differential equation: Since the limits t → 0 and x → ∞ both correspond to → ∞ in the Boltzmann domain, the three conditions (8-10) reduce to the following two boundary conditions: A flux in the Boltzmann domain can be defined as follows: Inserting Eq. (16) into Eq. (13), and making use of the chain rule, leads to Differentiating again, making use of Eq. (16), yields (Replacing q( ) with q(P) is permissible, because q and P are each monotonic functions of η.) Integration with respect to P yields where λ is a dummy variable of integration. It follows from Eq. (17) that dq∕dP vanishes at = 0 , which is to say, it vanishes when P = P f . Hence,

Integrating again yields
where ω is another dummy variable of integration. Reversing the order of integration leads to .

Dividing through by q(Pf ) leads to
The flux must vanish when → ∞ , where P = P i , and so it follows from Eq. (23) that

Inserting this expression into Eq. (23) yields
This is an integral equation for the fractional flow function, F(P) . This equation can easily be solved by an iterative procedure, by defining the iterated functions F n (P) , as follows: in which the current estimate F n (P) is inserted in the integrals on the right, to generate a new estimate, F n+1 (P) . The iterations continue until successive functions differ by less than some prescribed amount, such as 10 -5 . A sensible initial guess, which satisfies the bound- The flux in the physical (x, t) domain is now found from Darcy's law and Eq. (16): in which q(x, t) is the flux in the physical domain, and q(P) is the "flux" in the Boltzmann domain, as defined by Eq. (16). Since P = P f when x = 0 and = 0 , the flux at the outflow boundary is given by where q(P) is found from Eq. (24), using the converged value of the fractional flow function F. The numerical results obtained from this solution will be presented and discussed in Sect. 6, where they serve as a benchmark for the results obtained using the perturbation approach.

Perturbation Approach
The classical regular perturbation approach (Hinch 1991) will now be applied to the nonlinear differential equation boundary-value problem defined by Eqs. (7-10). The first step is to insert the exponential permeability function (11) into Eq. (7) and expand out the left-hand side, to yield: where D i = k i ∕ c is the hydraulic diffusivity of the system at the initial pressure, P i . A dimensionless pressure drawdown is now defined as which can be inverted to yield P(x, t) = P i − P D (P i − P f ) . Inserting this expression into Eq. (29) yields and dividing through by −(P i − P f ) , yields The dimensionless combination (P i − P f ) will be denoted by D , where D is a dimensionless stress sensitivity parameter. It follows from Eq. (11) that D = ln(k i ∕k f ) , where k i is the initial permeability before production starts, and k f is the lowered permeability at the outflow boundary due to the drawdown. In terms of this parameter, Eq. (31) can be written as From eqs. (8-10), the initial and boundary conditions for this problem take the form A dimensionless Boltzmann variable can be defined as (Ghez 1988): In the Boltzmann domain, Eq. (32) takes the form of the following second-order ordinary differential equation: The limits t → 0 and x → ∞ both correspond to → ∞ , and so the three conditions (33-35) reduce to the following two boundary conditions: Since D = 0 for a reservoir that is not stress sensitive, the parameter D can serve as the perturbation parameter in this problem. When D = 0 , Eq. (32) reduces to the traditional diffusion equation in the (x, t) domain, and Eq. (37) reduces to the form that the diffusion equation has in the Boltzmann domain, so it is seen that this problem will be a regular perturbation problem (Hinch 1991). Finally, since Eqs. (32) and (37) contain no singular points, the radius of convergence of the perturbation series will be infinite.
The dimensionless pressure is therefore sought in the form of the following series: in which P n D ( ) denotes the nth perturbation function (i.e., n is not a power exponent). This series is inserted into Eq. (37), and all terms, including the exponential, are expanded out in a power series in D , to yield: In order for a power series to identically equal zero, the coefficient of each power of D must vanish. For the ( D ) 0 term, this leads to the following homogeneous ODE for This is nothing other than the unperturbed problem that corresponds to D = 0 , the solution to which, after integrating and satisfying the boundary conditions (38,39), is (Ghez 1988) Note also for future reference that The ( D ) 1 term leads to the following inhomogeneous ODE for P 1 D ( ): In terms of the new variable y( ) = dP 1 D ( )∕d , and after making use of eqs. (43) and (44), this ODE can be written as The general solution of this equation can be written as the sum of the homogeneous solution and a particular solution to the inhomogeneous problem. The homogeneous solution is readily found to be y h ( ) = Ce − 2 , where C is an arbitrary constant. A particular solution to the inhomogeneous equation can be found by using the method of variation of parameters (Coddington 1961), which starts by writing (50) Since the 0th-order solution P 0 D ( ) already satisfies the boundary conditions (38) and (39), it follows that all higher-order functions P n D ( ) must satisfy zero boundary conditions. The condition that P 1 D (0) = 0 shows that D must be 0. The condition that P 1 D (∞) = 0 leads to The definite integral appearing in Eq. (53) can be evaluated (Edward and Murray 1968) to yield from which it follows that The final expression for P 1 D ( ) is found by inserting Eq. (55) into Eq. (52): In principle, the perturbation calculations could be continued, so as to obtain the second-order term, etc. However, it soon becomes apparent that the number of individual terms appearing in P n D ( ) grows very rapidly, and P 2 D ( ) would, for example, contain several dozen terms. Thus, pursuing the perturbation solution to higher orders does not seem to be feasible. Instead, the first-order solution will be used to obtain the flux (exactly) to first order in D , and this result will be extended to larger values of D by invoking an appropriate "effective permeability," as discussed below in Sect. 7.
By inserting Eq. (43) for P 0 D ( ) , and Eq. (56) for P 1 D ( ) , into series (40), the pressure profile, to first order in D , is seen to be given by Pressure profiles, in the Boltzmann domain, are plotted in Fig. 1, for D = 0 and D = 1 . Since the pressures are fixed at both the outflow boundary and in the far field, the additional perturbation function vanishes at both boundaries. As D increases, the pressure gradient at the outflow boundary increases, but the permeability at the outflow boundary decreases by a greater proportion, and so the flux will decrease, as would be expected. This latter fact can be seen by noting that flux is proportional to the area between the pressure drawdown curve and the line P D = 0 , and this area clearly decreases as D increases. The precise decrease in flux due to this stress dependence is calculated in the following section.

First-Order Perturbation Expression for the Flux
The flux of fluid out of the reservoir is found by starting with Darcy's law, and using the chain rule, along with Eqs. (30) and (36), to express the pressure gradient P∕ x in terms of the pressure profile in the Boltzmann domain: The flux at the outflow boundary is found by evaluating Eq. (58) at = 0 , which corresponds to x = 0 . For example, returning to the case of no stress dependence, then k(P) = k i , and Eq. (44) shows that dP D ∕d = −(2∕ √ )e − 2 , and so at the boundary, dP D ∕d = −2∕ √ . Hence, Eq. (58) shows that q( which is the known classical result (Stewart 2011;Zimmerman 2018). The minus sign that appears in the equation for the flux indicates that the fluid is flowing towards the boundary, whereas the x-coordinate axis points away from the boundary.
Returning to the pressure-dependent case, Eq. (58) gives, after making use of Eq. (57) for P D ( ), . Fig. 1 Pressure profiles, in the Boltzmann domain, for the problem of 1D flow to a constant-pressure boundary in a stress-sensitive formation, showing the influence of the dimensionless permeability modulus, D , according to the first-order perturbation solution Equation (44) shows that dP 0 D (0)∕d = −2∕ √ , and Eq. (56) shows that dP 1 D (0)∕d = −2( − 1)∕ 3∕2 , from which it follows that The normalized flux, normalized with respect to the case of no pressure dependence, i.e., the case D = 0 , can be written as Since the bracketed term has already neglected all terms of second-order or higher in D , the entire expression is only reliable to first order. Expanding the right-hand side in a Taylor series in the parameter D , and neglecting the higher-order terms, leads to This expression is exact, to first order. In Fig. 2, it is compared to the numerical values obtained from the McWhorter-Sunada integral equation formulation, Eq. (28). The flux computed from the first-order perturbation solution has 5% accuracy up to about D = 0.8 . Recalling that D = ln(k i ∕k f ) , it follows that D = 0.8 corresponds to k f = 0.45k i . Hence, the first-order solution has 5% accuracy for situations in which the stress sensitivity of the formation has caused the permeability at the outflow boundary to be less than half of the initial permeability. For more extreme cases, which may arise from a greater (60) Fig. 2 Boundary flux, for the problem of 1D flow to a constant-pressure boundary in a stress-sensitive formation, as a function of the dimensionless stress sensitivity parameter, D , according to the numerical solution, the first-order perturbation solution (short dashes), and "effective permeability" approach (long dashes). The flux is normalized against the flux for the case of no stress sensitivity stress sensitivity of the formation, a greater drawdown, or a combination of both factors, the accuracy of the first-order solution degrades further. This can be partially remedied by appealing to the concept of effective permeability, as shown in the following section.

Effective Permeability Concept
For heterogeneous rock masses in which the permeability is constant at each point in space, but varies spatially, the concept of "effective permeability" is often invoked. Roughly speaking, the effective permeability is the value that, when inserted into the solutions that can be obtained for the uniform permeability case, would yield the correct flux for the actual inhomogeneous problem. Exact expressions for the effective permeability can be obtained only for simple cases such as layered rocks, although accurate approximate methods have been devised for more realistic permeability distributions (Renard and de Marsily 1997). In all such scenarios, the effective permeability is essentially some sort of spatial average of the local permeability.
In the present problem, the permeability varies in space due to the fact that the permeability varies with pressure drawdown, and the drawdown varies in space and time. Hence, classical "spatial averaging" upscaling methods for computing effective permeability will not apply. However, a modified version of this concept can still be invoked, as will now be explained.
Based on the form of the permeability variation given by Eq. (11), it follows that the permeability at the outflow boundary is given by k f = k i e − D . The permeability in the far field is of course k i , and so the permeability throughout the formation varies between these two extremes. It seems reasonable and simple to put k eff = k i e − D , for some constant γ that will be expected to lie between 0 and 1. If this expression is inserted into the expression for the flux for the constant-permeability case, the result is All that is known with certainty in this problem is that the first-order flux is given by Eq. (60). Hence, it seems sensible that the factor γ should be chosen so as to make Eq. (63) agree with Eq. (60) to first order in D . Since e − D ∕2 = 1 − D ∕2 + ⋯ , this matching criterion requires that ∕2 = 1∕ , which is to say, = 2∕ . This choice leads to an effective permeability of k eff = k i e −(2∕ ) D , and a boundary flux of which corresponds to a normalized boundary flux of This expression is plotted in Fig. 2, where it is seen to have much better accuracy, over a larger range of values of D , than does the first-order expression given by Eq. (62). In fact, this new flux approximation has 5% accuracy out to about D = 2 . As this value of D corresponds to a situation in which k f = 0.14k i , i.e., a sevenfold reduction in permeability due to the drawdown at the outflow boundary, Eq. (65) is therefore accurate for a very wide range of expected scenarios.

Summary and Discussion
A first-order perturbation solution has been derived for the problem of one-dimensional flow to a planar boundary held at constant pressure, in a pressure-sensitive rock formation in which the permeability varies exponentially with pressure. The solution yields a result that is exact to first order in the parameter D = ln(k i ∕k f ) , where k i is the permeability at the initial reservoir pressure, and k f is the permeability at the outflow boundary. This solution gives a normalized flux (normalized against the flux for the case D = 0 ) of 1 − 0.318 D , which is to say, 1 − 0.318 ln(k i ∕k f ) . Although the problem has been discussed in terms of production of fluid from the reservoir, the result is equally applicable to injection, in which case D will be negative, and the injectivity will be enhanced by the stress sensitivity. The first-order solution was found to have 5% accuracy for values of D up to about 0.8. In order to extend this range of accuracy, an "effective permeability" was then defined so that, when inserted into the analytic solution for the case of no stress sensitivity, it yields a flux that agrees exactly, to first order, with the perturbation solution. This effective permeability was found to be given by The flux obtained using this effective permeability has 5% accuracy up to values of D as large as 2.
Since the "exact" solution to this type of problem can be obtained by numerically solving an integral equation using a rapidly convergent iterative procedure, this raises obvious questions regarding the value of the perturbation solution. One partial answer to this question is that the perturbation method yields a flux that is exact to first order in the parameter D , and any exact result for this type of highly nonlinear problem is certainly welcome. However, the main value of the foregoing exercise is probably as a "proof of concept," to show the feasibility of using perturbation series to solving nonlinear fluid flow problems in porous media. Although the perturbation approach for addressing porous media flow problems was pioneered in 1991 by Kikani and Pedrosa (1991), most researchers who subsequently followed their approach essentially stopped at the 0th-order solution. As shown in detail by Ren and Guo (2018), the 0th-order solution to the problem, obtained after invoking the change of pressure variable suggested by Kikani and Pedrosa (1991), yields very accurate solutions for flow problems in radial geometries, particularly at large times. However, the analogous solution for the linear flow geometry actually has very poor accuracy, as will now be shown.
The Kikani-Pedrosa-type 0th-order solution for the present problem can easily be derived by following the steps taken by Ren and Guo (2018) in their solution for the problem of flow to a vertical well at the center of a bounded circular reservoir under conditions of constant wellbore pressure. In fact, the steps are completely analogous, and so only the key results will be given here. Using the notation of the present paper, the dimensionless pressure variable can be expressed in terms of a new variable, D : (This is essentially a Kirchhoff transformation (Vadasz 2010), although this connection has not usually been made in papers that discuss this problem in a petroleum engineering context). The function D can be expanded in a perturbation series, the 0th-order term of which is the solution to the unperturbed differential equation, subject to the following boundary condition at the outflow boundary: It follows that the KP 0th-order solution for this problem is The flux corresponding to this approximate solution is found by inserting Eq. (69) into the general expression (59): The normalized flux at the boundary, normalized with respect to the case of a reservoir having no pressure sensitivity, therefore takes the form Expansion in a Taylor series yields But it is known from Eq. (62), and from the numerical solution, that the first two terms in the expansion for the normalized flux are 1 − 0.318 D . Hence, the 0th-order KP solution to this problem is not exact to first-order in the stress sensitivity parameter D , and is actually not particularly accurate. Finding a result accurate to first order in D would require computing the D1 term in the perturbation expansion for D .
This result at first seems to conflict with the findings of Ren and Guo (2018), who showed that the 0th-order KP solution for the radial flow problem is extremely accurate. But this accuracy was observed at large times and did not extend to very small values of dimensionless time, as these authors in fact pointed out. Although values of t D = k i t∕ cr 2 w < 1 are not of much practical importance in reservoir engineering or hydrology, the solution developed by Ren and Guo (2018) would yield a flux that diverges (70) 1 3 further below the correct flux as t D → 0 . A hint of this divergence can be seen in their Fig. 2b, if one visually extrapolates the numerical and perturbation fluxes to values below the range of graph, which stops at t D = 1.
Although the radial flow problem exhibits different behavior in the small-time and largetime regimes, the flow problem in a linear geometry has a self-similar solution and therefore does not possess different "early time" and "late-time" regimes (ignoring the eventual influence of any far-field boundary). Moreover, the flux in the radial flow problem must asymptotically approach the flux in the 1D problem, as t D → 0 (Crank 1975, Eq. (58)). The early-time flux in the radial problem must therefore be identical to the flux in the 1D problem, and so, in some sense, the solution to the 1D problem is an "early time" solution. In both the radial flow and linear flow cases, the 0th-order KP-type solution does not yield accurate results for "early times," and hence, the present findings are not contradictory or paradoxical. The flux expression derived in the present paper for the 1D flow problem; on the other hand, it has reasonably good accuracy, even for quite large values of D .