© R. Raghavan & C. Chen, published by IFP Energies nouvelles, 2021

Nomenclature

B : Formation volume factor [L³/L³]

c : Compressibility L T²/M

c _D : Dimensionless compressibility

C _D : Dimensionless storage constant

h : Thickness [L]

k : Permeability [L²]

L : Width of linear system [L]

m(p): Pseudopressure, see Equation (59) or Equation (60) [M²/L⁴/T2] or [M/L/T²]

p : Pressure [M/L/T²]

p′: Logarithmic derivative [M/L/T²]

p*(x_D, y_D, t_D): Pressure transformation; see Equation (18)

q : Flux [L³/T]

s : Skin factor

t : Time [T]

α _p : Normalizing constant

Γ(x): Gamma function

η : Diffusivity [L²/T]

ρ : Density [M/L³]

μ : Viscosity [M/L/T]

ζ(n): Riemann zeta function

ϕ : Porosity [L³/L³]

Subscripts

D: Dimensionless

t: Total

w: Wellbore

Superscript

–: Laplace transform

1 Introduction

The transient diffusion equation that describes the flow of slightly-compressibility fluids in porous media developed in Theis (1935), van Everdingen and Hurst (1949) and other places has been the linchpin for the evaluation of pressure responses for over eight decades with the assumption that the compressibility of the liquid is small with respect to pressure gradient – second degree terms are negligibly small. The assumption becomes important as we focus on lower permeability rocks. The consequences of this assumption may be examined through the Cole (1951)–Hopf (1950) or Kirchhoff (1894) transformations. Here we consider anew the transformations in the context of transient diffusion; the aspects addressed here are yet to be examined. We consider anew the matter of linearity; although it may be considered to be a well-worn topic, a quantitative measure to assess the suitability of the classical solution is yet to be available. Furthermore, in the process of arriving at such a measure which may be quite useful to have at one’s fingertips, at least in our opinion, a few glaring omissions were noted. These deficits are rectified. The Cole (1951)–Hopf (1950) transformation is central to many areas of mathematical physics from the study of turbulence, growth of interfaces, to understanding of shock waves and flow in porous media (Barros-Galvis et al., 2018; Bateman, 1915; Bertini and Giacomin, 1997; Burgers, 1940; Finjord, 1986; Finjord and Aadnoy, 1989; Friedrichs, 1948; Marshall, 2009; Miura, 1968). In essence this remarkable transformation turns a nonlinear equation into a linear equation. Although familiar to many who address nonlinear problems in the areas noted, it is also used quite often with no particular attribution (Braeuning et al., 1998; Chakrabarty et al., 1993; Jelmert and Vik, 1996; Odeh and Babu, 1988; Ren and Guo, 2017; Singh and Sagar, 1980; Wang and Dusseault, 1991). The citations above address a broad range of specific problems pertaining to fluid extraction and injection in both hydrocarbon production and groundwater literatures; but no overarching guidelines as to conditions under which classical solutions are adequate are provided as presented here. Because they are specific, not all authors provide solutions of general import; for example, expressions for pressure distributions are not addressed. Furthermore, as will be discussed in detail below, some works abandon the Cole-Hopf transformation midway for the Neumann boundary-condition as they wish to employ the Ozkan-Raghavan (1991) solutions for sources and sinks. Source/sink solutions, however, are difficult to address through the Cole-Hopf transformation (see below). Interestingly, this transformation is also useful to study multiphase-flow problems (see Burnell et al., 1989), it is not unusual to find an illustration of this transformation to linearize a nonlinear equation (Forsyth, 1906; Raghavan, 1993).

The Kirchhoff transformation is in many ways a better option; for the problem under consideration the transformation also converts the nonlinear differential equation to a linear one. Introduced by Kirchhoff (1894) for steady flow it was later extended by Van Dusen (1930) to unsteady flow, the problem of interest in this work. For flows in porous media the work of Al-Hussainy et al. (1966) is often cited for introducing this transformation; they do, however, note that such a transformation is proposed in Leibenzon (1953). That solution unlike the Cole-Hopf transformation is expressed in terms of an integral. Although Al-Hussainy et al. (1966) are usually quoted and recognized in terms of transient pressure tests, one other significant contribution of that work is to recognize that inclusion of nonlinear terms indicates that the life of the resource increases by an order of magnitude. Incorporation of second degree terms results in similar differences for liquid flow (see Raghavan et al., 1972) who addressed the combined influences of variations in liquid compressibility, viscosity and also rock properties like permeability and pore-volume compressibility with pressure and also multiphase flow (Raghavan, 1976). The goals of the two methods are identical, nevertheless for all problems of interest we recommend the Kirchhoff transformation as it pertains to the evaluation of the properties of rocks; the basis for our preference is outlined. We do recognize that our recommendation is at variance with other works (see, e.g., Vadasz, 2010) who notes that the Kirchhoff transformation is inconvenient for obtaining results directly in terms of pressure. Therefore, we consider afresh both transformations in context of transient diffusion in porous media. Specifically, we address methods to arrive at general conclusions wherein the nonlinear and linear solutions agree in a general, quantitative way rather than the vague, qualitative assertion that pressure gradients need to be small with respect to compressibility for the linear diffusion equations to apply as asserted in van Everdingen and Hurst (1949), Matthews and Russell (1967) and Earlougher (1977) and in all other texts. Such a quantitative measure is quite useful to have for the range of scales addressed nowadays is quite large. The measure given, which to our knowledge has been unavailable until now, is advantageous in that it applies to all the well configuration combinations in the Ozkan-Raghavan (1991) tables.

An appealing feature of the Cole-Hopf transformation is that it is particularly useful in that it provides solutions directly in terms of pressure; but disadvantages do exist. For example, the Neumann boundary condition transforms to the Robin boundary condition¹ and as a result one may not use the Ozkan and Raghavan (1991) tabulations directly to determine the needed solutions (see, e.g., Braeuning et al., 1998; Jelmert and Vik, 1996; Ren and Guo, 2017). This obstacle is overcome by application of Duhamel’s theorem which may have to be applied multiple times should pressure distributions be desired. Nevertheless if the goal is to use the Ozkan-Raghavan (1991) tabulations one may, of course, resort to Muskat (1934) and work in terms of density ab initio or use the Kirchhoff transformation as indicated here. Thus for completeness we also discuss the significant influence the Cole-Hopf transformation has on the boundary conditions; specifically, that it precludes application to an important problem, namely the line-source solution (or any other source/sink solution) which yields the Exponential Integral solution, the nucleus of many methods of evaluation such as the Horner (1951) and Jacob (1947) techniques. In this context we also explore procedures to address nonlinearities as they pertain to the Theis (1935) solution through procedures in Muskat (1934) and the Kirchhoff transformation. In concluding this section, we should note that working with the Cole-Hopf transformation is essentially equivalent to working with density.

1.2 The partial differential equations, associated conditions

As is usual we consider a uniform, homogeneous porous rock of permeability, k, porosity, ϕ, and thickness, h, that contains a fluid of density, ρ, compressibility, c, and viscosity, μ; we use the symbol η to denote the diffusivity of the system. The differential equation representing transient diffusion in terms of the pressure distribution, p(x, y, t), to be solved is

(1)

where the symbols x and y represent the coordinates, and the symbol, t, is time. The equation of state representative of the fluid is given by

(2)

or

(3)

with the subscript 0 representing the reference state. If we were to assume that ρ/ρ₀ ≈ 1 + c(p −p₀) – the slightly compressible fluid – then the second terms on the right-hand side of equation (1) may be dropped to arrive at a linear equation, the standard option as first shown explicitly in van Everdingen and Hurst (1949).

We presume p_i to be the initial pressure which is identical in all parts of the system at t = 0 that is infinite in its areal extent. Thus on defining Δp = p_i − p(x, y, t) we may write

(4)

We require a solution to equation (4) be subject to the following conditions. For a quiescent system where equilibrium prevails everywhere we require

(5)

further, as the system is infinite in its areal extent we also require

(6)

and

(7)

We assume constant properties and consider both Neumann and Dirichlet boundary conditions; consequently, we express certain dimensionless terms differently. Dimensionless pressure, p_D(x_D, y_D, t_D), dimensionless time, t_D, dimensionless distances, (x_D, y_D), and dimensionless compressibility, c_D, are defined, respectively, as

(8)

(9)

(10)

and

(11)

where α_p depends on the mode of production; that is,

(12)

(13)

In the above expressions, ℓ is the reference length; the expressions given here apply for situations where p_i > p(x, y, t) or where p_i < p(x, y, t); that is, production or injection, respectively. The symbol c_D for dimensionless compressibility is yet to appear in the literature and should not be misconceived to be the symbol for the dimensionless storage constant, namely C_D. Additional boundary conditions are defined when specific well conditions are addressed; as we noted earlier, we consider two representative problems.

As properties are assumed to be constant, equations (4)–(7) may be written, respectively, as

(14)

with

(15)

(16)

and

(17)

2 Analytical solutions

As discussed in Raghavan (1993) and other places the Cole-Hopf transformation involves defining a new dependent variable, p^⋆. That is, we consider the expression of the form

(18)

where ℶ is a constant. With the use of equation (18), the Cole-Hopf transformation, equation (14) leads to the standard diffusivity equation

(19)

provided that ℶc_D = −1. Incidentally, we note that many options exist to define the argument of ln(·). For example, Braeuning et al. (1998) define

(20)

We are now in a position to consider problems of interest to us. Use of equation (18) does lead to complications because the transformation will require the modification of the boundary conditions.

In concluding this section, we simply note that p^⋆ is essentially the density defined in a particular way. This point is rather important as it has consequences to this discussion. As already mentioned the matter of linearity may be simply addressed along the lines of Muskat (1934); one does not need to resort to the Cole-Hopf transformation if this were the option one ultimately decides to follow.

2.1 Flow in linear coordinates; production at constant pressure, p_wf; Dirichlet condition

We consider flow in a linear system that is quiescent at a pressure, p_i, initially with the pressure held constant at a value p_wf at the end x = 0 where the other end extends to infinity; that is, we consider the region, x ≥ 0. The relevant equations to be solved are

(21)

with

(22)

(23)

and

(24)

Results are normalized in terms of the width of the system.

In terms of the Laplace transformation, s, the solution of the above system of equations is

(25)

The inversion of equation (25) yields the standard result for constant-pressure-production in a semi-infinite system, namely,

(26)

Simplification yields

(27)

and for small values of the exponent, the expression

(28)

is obtainable.

2.2 Axisymmetric flow to a finite well radius; constant rate, q; Neumann condition

We consider axisymmetric flow to a well of radius, r_w, in a system that is infinite in its areal extent. The well rate, q, is assumed to be constant; that is, the gradient at the wellbore is given by

(29)

where r is the radial coordinate. The problem that needs to be considered in terms of p^⋆(r_D, t_D) is now given by:

(30)

and we seek a solution of equation (30) subject to

(31)

(32)

and

(33)

with c_D now defined by the expression in equation (12). On comparing equation (29) with equation (33), we see that the boundary condition at the wellbore changes from that of a Neumann specification to a Robin type; essentially from one that specifies the flux to one which is of the radiative type. This change does not permit us to obtain the solution for a line-source well because the solution for , the Laplace transformation for p^⋆(r_D, t_D) is

(34)

Here K₁(·) is the modified Bessel function of the second kind of order 0, and A is an arbitrary constant. Using equation (33), we may determine A and thus

(35)

where K₀(·) is the modified Bessel function of the second kind of order 1. As we are considering diffusion with radiation, the solution, equation (35), may be found in Carslaw and Jaeger (1959) and many other places. Now it is clear that the singularity that would exist in the denominator of the corresponding expression of for a line-source well precludes addressing this boundary condition.

The point made here regarding the change in the nature of the boundary condition is not always explicitly recognized. For example, the solutions of Ozkan and Raghavan (1991) have been used to determine p^⋆ directly; that is, without accounting for change in the nature of the boundary condition, equation (33). We make note of this issue because from a perfunctory consideration of Braeuning et al. (1998) it may appear that they have been successful in applying the Ozkan-Raghavan solutions to the problem under consideration. But their work considers only the wellbore. To obtain pressure distributions one then uses the sandface rate and applies Duhamel’s (1833) theorem to estimate pressure distributions. The above observation also applies to the work of Ren and Guo (2017) who follow the lead of Braeuning et al. (1998).

Incidentally, the analogue of equation (33) for a line-source well is

(36)

For classical diffusion p_D(r_D → 0, t_D) is undefined and we should expect the same for the present case. More on this later.

For now, however, we proceed with equation (35). Inversion of the right-hand side of equation (35) yields,

(37)

where f(t) is

(38)

The symbols J_n(x) and Y_n(x) are, respectively, the Bessel function of the first kind of order n[(n = 0) and (n = 1)] and the Neumann Bessel function of the second kind of order n, again, [(n = 0) and (n = 1)]. Or, p^⋆(r_D, t_D) may be determined as given in Jaeger and Carslaw (1943) by

(39)

Both forms given are useful in different contexts.

As we discuss below a simple expression of p^⋆(r_D, t_D) is unavailable. If an expression for p^⋆(r_D, t_D) were available, then the pressure distribution may be readily determined as in the Dirichlet problem, for example,

(40)

Computations of Jaeger Integrals of the type shown above are still of interest (see Phillips and Mahon, 2011), although methods similar to those proposed by Stehfest (1970a, b) appear to be the choice in many cases.

2.2.1 Asymptotic solutions

A number of asymptotic solutions are given in the literature; by their nature, they differ in specific details. Here we touch on a few options, principally to address issues of practical import such as pressure buildup and multirate tests. Again, because of the nature of the governing differential equations, they would have to be addressed through p^⋆. For now, we focus on long-term behaviors although the observations of Finjord (1986) are not in agreement with this perspective. We assume r_D = 1; that is, we consider the wellbore response.

Jaeger and Carslaw (1943) show that long-term expressions may be determined through the Integral, I(p, q, x), defined by

(41)

For large values of x,

(42)

where ) is the Riemann Zeta function, with

(43)

We use p = 1, q = c_D with the expression

(44)

with to obtain

(45)

and now

(46)

Rather than using the Inversion Integral, I(p, q, x), we may rewrite the expression in equation (35) as

(47)

and determine the long-term measure of the expression in equation (47) for by considering s → 0 for K₀(x) and K₁(x) which are given, respectively, by

(48)

and

(49)

On using the first term of the two expansions, we may write the right-hand side of equation (47) as

(50)

According to Ritchie and Sakakura (1956)

(51)

where in equation (51), the column vector, (−1, p)^T, represents the binomial coefficient, and Γ(·) represents the Gamma function. Expressions for the coefficients of the first term are given in Ritchie and Sakakura (1956) and more recently in Yeh and Wang (2007) leading to the expression in equation (45).

2.3. Response caused by a Line Source

This solution permits us to address tests commonly used to consider connectivity (interference tests, hydraulic tomography) in the porous rock. Options available include the use of density as in Muskat (1934) or pseudopressures as in Raghavan (1993). It is also particularly useful in solution techniques such as the method of sources and sinks (Carslaw and Jaeger, 1959; Raghavan and Ozkan, 1994).

2.4 Formulation in terms of density, ρ(r_D, t_D)

We work with density, ρ(r_D, t_D), and q_m the mass production rate. On solving the following system of equations

(52)

(53)

(54)

and

(55)

where Δρ(r_D, t_D) = ρ_i − ρ(r_D, t_D), the development used by Theis (1935) with no modifications, results in the expression for the distribution for density

(56)

Now on expressing the solution in terms of dimensionless pressure, p_D(r_D, t_D), through the definition of density in terms of pressure, we obtain

(57)

where −Ei(−u) is the Exponential Integral with the restriction that exp{c[p(r) − p₀]} ≈ 1 + ϵ, or c[p(r) − p₀] ≈ ϵ.

Incidentally, the analogue of the solution given in equation (56) for a line-source well that completely penetrates the formation given in equation (18) of Braeuning et al. (1998) which is obtained by the source-sink approach, in our terminology yields the solution at the wellbore to be

(58)

2.5 Formulation in terms of pseudopressure; the Kirchhoff transform

We define the pseudopressure, m(p), through a new variable, defined by

(59)

or by

(60)

where B is the formation volume factor. The compressibility of the liquid in terms of B is

(61)

We may replace the integrand in equation (59) by ρ/ρ_sc in which case m(p) would have the same dimensions as pressure, p. Also, there are advantages to using the second definition of m(p) in terms of practical considerations. A formal solution in terms of pseudopressure yields

(62)

On substituting the expression for the equation of state we may show that the left-hand side of equation (62) will yield the left-hand side of equation (56); thus, on re-expressing this solution in terms of pressure, p(r_D, t_D), we obtain equation (57).

Equation (57) applies to a line-source well; in the previous case no such solution was possible. As ln(1 + x) ≈ x as x → 0, we have for small enough c_D, as c_D → 0 the Theis (1935) solution, namely

(63)

3 Computational results

The goals of the computations noted below are three-fold. First, we establish quantitatively the bounds for which the classical solutions, the linchpins for determining the properties of reservoir rocks for close to a century, apply; this limit carries over to all well configurations as, ultimately, radial or pseudoradial flow will prevail. Second, responses to the finite-well-radius solution are compared through both transformations. The basis to apply the Kirchhoff transformation is established analytically. Third, because the Theis solution plays an important role in evaluating interference or tomographic tests, the Mueller and Witherspoon (1965) criteria are explored and new guidelines are provided to consider situations should nonlinear terms be important. Insights into the structure of the nonlinear solutions are provided. All solutions presented here are yet to be available in the literature.

The following steps were taken to ensure the accuracy of our solutions. First, we verified that the Cole-Hopf solutions do approach the classical solutions in the limit; that is, as c_D → 0. Second, concurrence with the two asymptotic solutions, equations (45) and (50), is obtained for 0.01 ≤ c_D ≤ 3.0. Third, the Cole-Hopf solutions are reexpressed in terms of pseudopressures to show concurrence with the solutions for a slightly compressible fluid. Additional details are given in the discussion below. It appears to us that this result also establishes the validity of our recommendation to use the pseudopressure procedure when the compressibility is high even though a single-phase liquid is produced.

Figure 1 examines the influence of dimensionless compressibility, c_D, on responses at a well in terms of the logarithmic derivative, , as the slopes of the pressure response curves are central to the evaluation of rock properties as well as the well condition, namely the skin factor, s. For values of c_D ≤ 10⁻², the conventional value of 1/2 is obtained. For larger c_D, as time increases, the slopes depend on both dimensionless time, t_D, and dimensionless compressibility, c_D. Of particular importance, of course, is the fact that the classical semilog straight line that forms the nucleus of many methods of analysis is no longer applicable if c_D is large enough, and that values of c_D > 10⁻² would negate many of the basic techniques developed over the years to estimate rock properties through transient tests. Except for the upper bound given here that is yet to be available, the result that the slopes depend on c_D if it were large enough, may be intuitively expected from the works of Al-Hussainy et al. (1966) and Raghavan et al. (1972). In concluding this phase of our discussion, we emphasize that the criterion presented is based on the slope of the curves and not the value of dimensionless pressure.

The advantage of the pseudopressure transformation is demonstrated in Figure 2. Solutions for three values of compressibility, c_D, obtained through the Cole-Hopf transformation shown as dashed lines are considered; these values of c_D are large enough for the conventional semilogarithmic approximation to not exist. The solutions assume r_D to be 1. The top unbroken line is the finite-well-radius solution obtained through the expression

(64)

The circles, squares and diamonds that essentially form a single curve and fall in alignment with the finite-well-radius solution represent solutions of the dashed lines, the Cole-Hopf responses, when represented in terms of the dimensionless well pseudopressure, m_wD. This result unequivocally establishes the power of the pseudopressure concept to address the matter of second-degree terms; most importantly, when done so, the semilog-straight line exists. Furthermore, the information in Figure 2 suggests that for evaluation of well performance and for the evaluation of rock properties through transient diffusion it would be appropriate to reexpress measured pressure responses obtained through pseudopressures when second degree terms are important as the classical norms for a slightly-compressible fluid developed in the literature would hold and one does not have to appraise the role of the magnitude of the value of compressibility with respect to pressure gradients in evaluating rock properties. As we noted earlier, such a result establishes the validity of our recommendation concerning the pseudopressure approach when the compressibility is high even though a single-phase liquid is produced. The results establish the accuracy, viability and concurrence between the two methods explored in this work. This observation carries over to all Ozkan-Raghavan tabulations for other well configurations and is given further credence in the solutions given in Figure 3 which we assess next.

Mueller and Witherspoon (1965) note that for r_D ≥ 20, the line-source solution of Theis (1935) and the finite-well radius solution are indistinguishable if ; in addition if r_D were 1, namely the wellbore, then the two solutions are within one percent if . The first observation permits use of the Theis solution for evaluating multiple-well or interference tests and the second permits us to use the Theis solution at the flowing well if times are long enough. Figure 3 considers responses in terms of the Mueller-Witherspoon scale for c_D = 10⁻¹. A careful examination indicates that solutions are approximately correlatable in terms of for r_D ≥ 10²; a dependence on r_D does exist for r_D > 10² but it is negligibly small; this observation, however, applies to the time range considered in Figure 3. For larger times, the differences in the solutions are significant; the slopes depend on r_D. For c_D > 10⁻¹, computations indicate the limit of r_D of 10² suggested above should be increased. After examining results up to c_D = 1, we suggest that, for all practical purposes, the influence of distance may be ignored if r_D were ≥350 and . Most important for our needs, however, is the fact that the curves do not merge and become independent of r_D if times are long enough. Naturally, again, it is needless to say that through inspection of equation (62) the use of the pseudopressure concept resolves concerns should they exist.

For completeness, we note that both transformations may be used for situations beyond those addressed here. If we were to restrict our discussion to analytical options, then situations where porosity and permeability are exponentially dependent on pressure of the form exp[β_k(p − p₀)] where β_k is a constant have been addressed (Kikani and Pedrosa, 1991; Marshall, 2009). For testing purposes, to evaluate rock properties all the observations made thus far hold as these additions merely change the meaning of c_D. One point that needs to be considered is that the constant β_k must be small; otherwise limiting conditions may be approached in unrealistically short times. Experience with assuming an exponential dependence for permeability in situations such as the Ekofisk oil field in the North Sea and others indicate that such an assumption is of limited utility (Chin et al., 2000). Finally, we should note the Mueller-Witherspoon (1965) criterion should not be relaxed should interference tests be evaluated through pseudopressures.

4 Discussion and concluding remarks

One impetus for this work is the evaluation of systems that produce highly compressible fluids such as shales. Although, as mentioned in the Introduction, the topic may appear to be well-worn, the ideas explored here are transformational in many senses. The role of compressibility on transient diffusion has been evaluated and a quantitative measure to ensure the applicability of analyses in many contexts is given. Although, the role of compressibility through the Cole-Hopf transformation has been evaluated numerous times as we have cited, the manner in which we have considered the problem and the observations we have made are yet to be noted. Furthermore practical options to address limitations noted by means of the Kirchhoff transformation are discussed. Two concrete contributions should be noted. First, it is shown, for the first time, that the classical assumption of a slightly-compressible fluid holds if c_D ≤ 10⁻². This bound addresses conditions under which classical solutions apply. Second, using the line-source-well problem that is intractable through the Cole-Hopf transformation, the advantages that accrue through the Kirchhoff transformation are explicitly demonstrated. Although solutions have been discussed in the context of the line-source solution, the observations we have made are of a more general import for our observations are not restricted to the line-source idealization and apply to all well completion schemes and models (see, e.g., the Ozkan-Raghavan tabulations) used to evaluate flow in porous rocks; the same holds true when fluid density is used as the dependent variable. Additionally, there is no change in the nature of the boundary condition and is thus convenient for considering other well configurations especially in 2D. In using density, ρ(x, t), or pseudopressure, m[p(x, t)], as the dependent variable we may directly use the Ozkan-Raghavan (1991) tabulations; furthermore, no theoretical framework needs to be developed as all rules and expressions that are available for a slightly-compressible fluid (including those for multirate tests) apply. Additionally, if pseudopressures are used one may go beyond the constant compressibility idealization and may also incorporate changes in viscosity (see Raghavan et al., 1972). In passing, we note that we have compared responses predicted by equation (58) with the rigorous solution. The Braeuning et al. (1998) approximation is adequate if c_D ≤ 0.1. For larger values of c_D reasonable agreement in slopes is obtained at large times (within 5%), however, during the transition period differences are significant. For other situations, the validity of the Braeuning et al. (1998) two-term approximation remains an open question. For problems of the type considered by Ren and Guo (2017), as noted here, the analysis in terms of pseudopressures is the best option; in addition, the formulation of Chen and Raghavan (1995) is recommended.

---

References

All Figures