Regular Article
Addressing nonlinear transient diffusion in porous media through transformations
^{1}
The Raghavan Group, PO Box 52756, Tulsa, OK 74152, USA
^{2}
Kappa Engineering, 11767 Katy Freeway, #500, Houston, TX 77079, USA
^{*} Corresponding author: raghavan.raj@gmail.com
Received:
25
May
2021
Accepted:
12
November
2021
The nonlinear differential equation describing flow of a constant compressibility liquid in a porous medium is examined in terms of the Kirchhoff and ColeHopf transformations. A quantitative measure of the applicability of representing flow by a slightly compressible liquid – which leads to a linear differential equation, the Theis equation – is identified. The classical Theis problem and the finitewellradius problem in a system that is infinite in its areal extent are used as prototypes to address concepts discussed. This choice is dictated by the ubiquity of solutions that depend on these archetypal examples for examining transient diffusion. Notwithstanding that the Kirchhoff and ColeHopf transformations arrive at a linear differential equation, for the specific purposes of this work – the estimation of the hydraulic properties of rocks, the Kirchhoff transformation is much more advantageous in a number of ways; these are documented. Insights into the structure of the nonlinear solution are provided. The results of this work should prove useful in many contexts of mathematical physics though developed in the framework of applications pertaining to the earth sciences.
© R. Raghavan & C. Chen, published by IFP Energies nouvelles, 2021
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Nomenclature
B : Formation volume factor [L^{3}/L^{3}]
c _{D} : Dimensionless compressibility
C _{D} : Dimensionless storage constant
L : Width of linear system [L]
m(p): Pseudopressure, see Equation (59) or Equation (60) [M^{2}/L^{4}/T2] or [M/L/T^{2}]
p′: Logarithmic derivative [M/L/T^{2}]
p*(x_{D}, y_{D}, t_{D}): Pressure transformation; see Equation (18)
Subscripts
Superscript
1 Introduction
The transient diffusion equation that describes the flow of slightlycompressibility 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 wellworn 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 (BarrosGalvis 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 ColeHopf transformation midway for the Neumann boundarycondition as they wish to employ the OzkanRaghavan (1991) solutions for sources and sinks. Source/sink solutions, however, are difficult to address through the ColeHopf transformation (see below). Interestingly, this transformation is also useful to study multiphaseflow 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 AlHussainy 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 ColeHopf transformation is expressed in terms of an integral. Although AlHussainy 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 porevolume 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 OzkanRaghavan (1991) tables.
An appealing feature of the ColeHopf 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^{1} 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 OzkanRaghavan (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 ColeHopf transformation has on the boundary conditions; specifically, that it precludes application to an important problem, namely the linesource 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 ColeHopf 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
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
or
with the subscript 0 representing the reference state. If we were to assume that ρ/ρ_{0} ≈ 1 + c(p −p_{0}) – the slightly compressible fluid – then the second terms on the righthand 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
We require a solution to equation (4) be subject to the following conditions. For a quiescent system where equilibrium prevails everywhere we require
further, as the system is infinite in its areal extent we also require
and
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
and
where α_{p} depends on the mode of production; that is,
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
with
and
2 Analytical solutions
As discussed in Raghavan (1993) and other places the ColeHopf transformation involves defining a new dependent variable, p^{⋆}. That is, we consider the expression of the form
where ℶ is a constant. With the use of equation (18), the ColeHopf transformation, equation (14) leads to the standard diffusivity equation
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
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 ColeHopf 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
with
and
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
The inversion of equation (25) yields the standard result for constantpressureproduction in a semiinfinite system, namely,
Simplification yields
and for small values of the exponent, the expression
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
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:
and we seek a solution of equation (30) subject to
and
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 linesource well because the solution for , the Laplace transformation for p^{⋆}(r_{D}, t_{D}) is
Here K_{1}(·) 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
where K_{0}(·) 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 linesource 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 OzkanRaghavan 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 linesource well is
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 righthand side of equation (35) yields,
where f(t) is
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
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,
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 longterm 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 longterm expressions may be determined through the Integral, I(p, q, x), defined by
For large values of x,
where ) is the Riemann Zeta function, with
We use p = 1, q = c_{D} with the expression
with to obtain
and now
Rather than using the Inversion Integral, I(p, q, x), we may rewrite the expression in equation (35) as
and determine the longterm measure of the expression in equation (47) for by considering s → 0 for K_{0}(x) and K_{1}(x) which are given, respectively, by
and
On using the first term of the two expansions, we may write the righthand side of equation (47) as
According to Ritchie and Sakakura (1956)
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
and
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
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
where −Ei(−u) is the Exponential Integral with the restriction that exp{c[p(r) − p_{0}]} ≈ 1 + ϵ, or c[p(r) − p_{0}] ≈ ϵ.
Incidentally, the analogue of the solution given in equation (56) for a linesource well that completely penetrates the formation given in equation (18) of Braeuning et al. (1998) which is obtained by the sourcesink approach, in our terminology yields the solution at the wellbore to be
2.5 Formulation in terms of pseudopressure; the Kirchhoff transform
We define the pseudopressure, m(p), through a new variable, defined by
or by
where B is the formation volume factor. The compressibility of the liquid in terms of B is
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
On substituting the expression for the equation of state we may show that the lefthand side of equation (62) will yield the lefthand side of equation (56); thus, on reexpressing this solution in terms of pressure, p(r_{D}, t_{D}), we obtain equation (57).
Equation (57) applies to a linesource 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
3 Computational results
The goals of the computations noted below are threefold. 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 finitewellradius 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 ColeHopf 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 ColeHopf 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 singlephase 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^{−2}, 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^{−2} 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 AlHussainy 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.
Fig. 1 Pressure responses at a well in terms of the pressure derivative, , illustrating the influence of dimensionless compressibility, c_{D}. 
The advantage of the pseudopressure transformation is demonstrated in Figure 2. Solutions for three values of compressibility, c_{D}, obtained through the ColeHopf 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 finitewellradius solution obtained through the expression
Fig. 2 Advantages of the pseudopressure concept to incorporate seconddegree terms. The circles, squares and diamonds are the ColeHopf solutions (dashed lines) when expressed in terms of pseudopressures. Alignment with the classical finitewellradius solution shown as a unbroken line is excellent. Norms based on the concept of a slightlycompressible fluid apply when pressures are evaluated through the pseudopressure concept or Kirchhoff transformation. 
The circles, squares and diamonds that essentially form a single curve and fall in alignment with the finitewellradius solution represent solutions of the dashed lines, the ColeHopf 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 seconddegree terms; most importantly, when done so, the semilogstraight 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 slightlycompressible 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 singlephase 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 OzkanRaghavan tabulations for other well configurations and is given further credence in the solutions given in Figure 3 which we assess next.
Fig. 3 Correlation of pressure distributions along the lines indicated in Mueller and Witherspoon (1965) for a linesource well. For values of r_{D} > 350 pressure responses may be correlated in terms of for the time range shown for the value of c_{D} considered. At longer times, solutions have to be expressed in terms of and r_{D}. 
Mueller and Witherspoon (1965) note that for r_{D} ≥ 20, the linesource solution of Theis (1935) and the finitewell 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 multiplewell 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 MuellerWitherspoon scale for c_{D} = 10^{−1}. A careful examination indicates that solutions are approximately correlatable in terms of for r_{D} ≥ 10^{2}; a dependence on r_{D} does exist for r_{D} > 10^{2} 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^{−1}, computations indicate the limit of r_{D} of 10^{2} 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_{0})] 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 MuellerWitherspoon (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 wellworn, 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 ColeHopf 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 slightlycompressible fluid holds if c_{D} ≤ 10^{−2}. This bound addresses conditions under which classical solutions apply. Second, using the linesourcewell problem that is intractable through the ColeHopf transformation, the advantages that accrue through the Kirchhoff transformation are explicitly demonstrated. Although solutions have been discussed in the context of the linesource solution, the observations we have made are of a more general import for our observations are not restricted to the linesource idealization and apply to all well completion schemes and models (see, e.g., the OzkanRaghavan 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 OzkanRaghavan (1991) tabulations; furthermore, no theoretical framework needs to be developed as all rules and expressions that are available for a slightlycompressible 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) twoterm 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
 AlHussainy R., Ramey H.J. Jr, Crawford P.B. (1966) The flow of real gases through porous media, J. Pet. Tech. 18, 624–636. [Google Scholar]
 BarrosGalvis N., Samaniego V.F., CincoLey H. (2018) Fluid dynamics in naturally fractured tectonic reservoirs, J. Petrol. Explor. Prod. Technol. 8, 1–6. https://doi.org/10.1007/s1320201703208. [Google Scholar]
 Bateman H. (1915) Some recent researches on the motion of fluids, Mon. Weather Rev. 43, 163–170. [Google Scholar]
 Bertini L., Giacomin G. (1997) Stochastic Burgers and KPZ equations from particle systems, Comm. Math. Phys. 183, 571–607. [Google Scholar]
 Braeuning S., Jelmert T.A., Vik S.A. (1998) The effect of the quadratic gradient term on variablerate welltests, J. Pet. Sci. Eng. 21, 203–222. [Google Scholar]
 Burgers J.M. (1940) Application of a model system to illustrate some points of the statistical theory of free turbulence, Proc. Nederl. Akad. Wetensch 43, 2–12. [Google Scholar]
 Burnell J.G., McNabb A., Weir G.J., Young R. (1989) Twophase boundary layer formation in a semiinfinite porous slab, Trans. Porous Med. 4, 395–420. [Google Scholar]
 Carslaw H.S., Jaeger J.C. (1959) Conduction of heat in solids, 2nd edn., Clarendon Press, Oxford, 510p. [Google Scholar]
 Chakrabarty C., Farouq Ali S.M., Tortike W.S. (1993) Analytical solutions for radial pressure distribution including the effects of the quadraticgradient term, Water Resour. Res. 29, 1171–1177. [Google Scholar]
 Chen C., Raghavan R. (1995) Modeling a fractured well in a composite reservoir, SPE Form. Eval. 10, 241–246. https://doi.org/10.2118/28393PA. [Google Scholar]
 Chin L.Y., Raghavan R., Thomas L.K. (2000) Fully coupled geomechanics and fluidflow analysis of wells with stressdependent permeability, SPE J. 5, 32–45. https://doi.org/10.2118/589680PA. [Google Scholar]
 Cole J.D. (1951) On a quasilinear parabolic equation occurring in aerodynamics, Quart. Appl. Math. 9, 225–236. [Google Scholar]
 Duhamel J.M.C. (1833) Mémoire sur la méthode générale relative au mouvement de la chaleur dans les corps solides plongés dans les milieux dont la température varie avec le temps, J. de Éc. Polyt. Paris 14, 20–77. [Google Scholar]
 Earlougher R.C. Jr (1977) Advances in well test analysis, SPE Monograph Series 5, 264. [Google Scholar]
 Finjord J. (1986) Curling up the slope: Effects of the quadratic gradient term in the infiniteacting period for two dimensional reservoir flow, SPE16451MS, Society of Petroleum Engineers. [Google Scholar]
 Finjord J., Aadnoy B.S. (1989) Effects of the quadratic gradient term in steadystate and semisteadystate solutions for reservoir pressure, SPE Form. Eval. 4, 413–417. [Google Scholar]
 Forsyth A.R. (1906) Theory of differential equations, Vol. VI,Cambridge University Press, p. 102, Ex. 3. [Google Scholar]
 Friedrichs K.O. (1948) Formation and decay of shock waves, Comm. Pure Appl. Math 1, 211–245. [Google Scholar]
 Hopf E. (1950) The partial differential equation u_{t} + uu_{x} =u_{xx}, Comm. Pure Appl. Math. 3, 201–230. [Google Scholar]
 Horner D.R. (1951) Pressure buildup in wells, in: Brill E.J. (ed.), Proceedings of the Third World Petroleum Congress, Vol. II, Leiden, pp. 503–521. [Google Scholar]
 Jacob E.E. (1947) Drawdown test to determine effective radius of artesian well, Trans. Am. Soc. Civ. Eng. 112, 1047–1070. [Google Scholar]
 Jaeger J.C., Carslaw H. (1943) Heat flow in the region bounded internally by a circular cylinder, Proc. R. Soc. Edinb. A: Math. Phys. Sci. 61, 223–228. https://doi.org/10.1017/S0080454100006233. [Google Scholar]
 Jelmert T.A., Vik S.A. (1996) Analytic solution to the nonlinear diffusion equation for fluids of constant compressibility, J. Pet. Sci. Eng. 14, 231–233. [Google Scholar]
 Kikani J., Pedrosa O.A. Jr (1991) Perturbation analysis of stresssensitive reservoirs, SPE Form. Eval. 6, 379–386. [Google Scholar]
 Kirchhoff G. (1894) Vorlesungen über mathematische Physik, Theorie der Wärme, Vol. 4, Barth, Leipzig. [Google Scholar]
 Leibenzon L.S. (1953) Underground hydraulics of water, oil, and gas, Vol. 2, Izd. AN SSSR (Publ. Acad. Sci.), Moscow, pp. 163–209. [Google Scholar]
 Matthews C.S., Russell D.G. (1967) Pressure buildup and flow tests in wells, SPE Monograph Series 1, 163. [Google Scholar]
 Marshall S.L. (2009) Nonlinear pressure diffusion in flow of compressible liquids through porous media, Transp. Porous Med. 77, 431–446. [Google Scholar]
 Miura R.M. (1968) Kortewegde Vries equation and generalizations. I. A remarkable explicit nonlinear transformation, J. Math. Phys. 9, 1202–1204. [Google Scholar]
 Mueller T.D., Witherspoon P.A. (1965) Pressure interference effects within reservoirs and aquifers, J. Pet. Tech. 17, 471–474. [Google Scholar]
 Muskat M. (1934) The flow of compressible fluids through porous media and some problems in heat conduction, Physics 5, 71–94. [Google Scholar]
 Odeh A.S., Babu D.K. (1988) Comparison of solutions of the nonlinear and linearized diffusion equations, SPE Reserv. Eng. 3, 1202–1206. [Google Scholar]
 Ozkan E., Raghavan R. (1991) Some new solutions to solve problems in well test analysis: I – Analytical considerations, SPE Form. Eval. 63, 359–368. [Google Scholar]
 Phillips W.R.C., Mahon P.J. (2011) On approximations to a class of Jaeger integrals, Proc. R. Soc. A. 467, 3570–3589. https://doi.org/10.1098/rspa.2011.0301. [Google Scholar]
 Raghavan R. (1993) Well test analysis, Prentice Hall, Englewoods Cliffs, NJ.p. 40, Ex. 5. [Google Scholar]
 Raghavan R. (1976) Well test analysis: Wells producing by solution gas drive, Soc. Pet. Eng. J. 16, 196–208. https://doi.org/10.2118/5588PA. [Google Scholar]
 Raghavan R., Scorer J.D.T., Miller F.G. (1972) An investigation by numerical methods of the effect of pressuredependent rock and fluid properties on well flow tests, SPE J. 12, 267–275. [Google Scholar]
 Raghavan R., Ozkan E. (1994) A method for computing unsteady flows in porous media, Pitman Research Notes in Mathematics Series, Vol. 318, Longman Scientific & Technical, Harlow, UK. [Google Scholar]
 Ren J., Guo P. (2017) Nonlinear flow model of multiple fractured horizontal wells with stimulated reservoir volume including the quadratic gradient term, J. Hydrol. 554, 155–172. https://doi.org/10.1016/j.jhydrol.2017.09.005. [Google Scholar]
 Ritchie R.H., Sakakura A.Y. (1956) Asymptotic expansions of solutions of the heat conduction equation in internally bounded cylindrical geometry, J. Appl. Phys. 27, 1453–1459. [Google Scholar]
 Singh S.R., Sagar B. (1980) On Jacob’s approximation in flow through porous media, Water Resour. Res. 16, 377–380. https://doi.org/10.1029/WR016i002p00377. [Google Scholar]
 Stehfest H. (1970a) Algorithm 368: Numerical inversion of Laplace transforms [D5], Commun. ACM 13, 47–49. [Google Scholar]
 Stehfest H. (1970b) Remark on algorithm 368: Numerical inversion of Laplace transforms, Commun. ACM 13, 624. [Google Scholar]
 Theis C.V. (1935) The relationship between the lowering of the piezometric surface and the rate and duration of discharge using groundwater storage, EoS Trans. AGU 16, 519–524. [Google Scholar]
 van Everdingen A.F., Hurst W. (1949) The application of the LaPlace transformation to flow problems in reservoirs, Trans. AIME 186, 305–324. [Google Scholar]
 Vadasz P. (2010) Analytical solution to nonlinear thermal diffusion: Kirchhoff versus ColeHopf transformations, J. Heat Transfer 132, 121302. https://doi.org/10.1115/1.4002325. [Google Scholar]
 Van Dusen M.S. (1930) Note on the theory of heat conduction, Bur. Stand. J. Res. 4, 753–756. https://doi.org/10.6028/jres.004.050. [Google Scholar]
 Wang Y., Dusseault M.B. (1991) The effect of quadratic gradient terms on the borehole solution in poroelastic media, Water Resour. Res. 27, 3215–3223. https://doi.org/10.1029/91WR01552. [Google Scholar]
 Yeh H.D., Wang C.T. (2007) Largetime solutions for groundwater flow problems using the relationship of small p versus large t, Water Resour. Res. 43, W06502. https://doi.org/10.1029/2006WR005472. [Google Scholar]
All Figures
Fig. 1 Pressure responses at a well in terms of the pressure derivative, , illustrating the influence of dimensionless compressibility, c_{D}. 

In the text 
Fig. 2 Advantages of the pseudopressure concept to incorporate seconddegree terms. The circles, squares and diamonds are the ColeHopf solutions (dashed lines) when expressed in terms of pseudopressures. Alignment with the classical finitewellradius solution shown as a unbroken line is excellent. Norms based on the concept of a slightlycompressible fluid apply when pressures are evaluated through the pseudopressure concept or Kirchhoff transformation. 

In the text 
Fig. 3 Correlation of pressure distributions along the lines indicated in Mueller and Witherspoon (1965) for a linesource well. For values of r_{D} > 350 pressure responses may be correlated in terms of for the time range shown for the value of c_{D} considered. At longer times, solutions have to be expressed in terms of and r_{D}. 

In the text 