1. business and industrial
2. energy
3. renewable energy
4. hydroelectric energy

LAPLACE TRANSFORM DECONVOLUTION AND ITS APPLICATION
TO PERTURBATION SOLUTION OF NON-LINEAR
DIFFUSIVITY EQUATION
by
Mahmood Ahmadi
A thesis submitted to the Faculty and the Board of Trustees of the Colorado School of Mines in partial
fulfillment of the requirements for the degree of Doctor of Philosophy (Petroleum Engineering).
Golden, Colorado
Date: ______________________
Signed: ______________________________
Mahmood Ahmadi
Signed: ______________________________
Dr. Erdal Ozkan
Thesis Advisor
Signed: ______________________________
Dr. Luis Tenorio
Thesis Advisor
Golden, Colorado
Date:_________________________
Signed: ______________________________
Dr. Ramona M. Graves
Professor and Petroleum Engineering
Department Head
ii
ABSTRACT
The primary objective of this dissertation is to extend the conveniences of deconvolution to non-linear problems
of fluid flow in porous media. Unlike conventional approaches, which are based on an approximate linearization of
the problem, here the solution of the non-linear problem is linearized by a perturbation approach, which permits
term-by-term application of deconvolution. Because the proposed perturbation solution is more conveniently
evaluated in the Laplace-transform domain and the standard deconvolution algorithms are in the time-domain, an
efficient deconvolution procedure in the Laplace domain is a prerequisite. Therefore, the main objective of the
dissertation is divided into two sub-objectives: 1) the analysis of variable-rate production data by deconvolution in
the Laplace domain, and 2) the extension of perturbation solution of the nonlinear diffusivity equation governing gas
flow in porous media presented by Barreto (2011) into the Laplace domain.
For the first research objective, a new algorithm is introduced which uses inverse mirroring at the points of
discontinuity and adaptive cubic splines to approximate rate or pressure versus time data. This algorithm accurately
transforms sampled data into Laplace space and eliminates the Numerical inversion instabilities at discontinuities or
boundary points commonly encountered with the piece-wise linear approximations of the data. The approach does
not require modifications of scattered and noisy data or extrapolations of the tabulated data beyond the end values.
Practical use of the algorithm presented in this research has applications in a variety of Pressure Transient
Analysis (PTA) and Rate Transient Analysis (RTA) problems. A renewed interest in this procedure is inspired from
the need to evaluate production performances of wells in unconventional reservoirs. With this approach, we could
significantly reduce the complicating effects of rate variations or shut-ins encountered in well-performance data.
Moreover, the approach has proven to be successful in dealing with the deconvolution of highly scattered and noisy
data.
The second objective of this research focuses on the perturbation solution of the nonlinear gas diffusivity
equation in Laplace domain. This solution accounts for the nonlinearity caused by the dependency of gas properties
(viscosity-compressibility product and gas deviation factor) on pressure. Although pseudo-pressure transformation
introduced by Al-Hussainy et al. (1966) linearizes the diffusivity equation for compressible fluids (gas), the pressure
dependency of gas properties is not completely removed. Barreto (2011) presented a perturbation-based solution
using Green’s functions to deal with the remaining non-lineraities of the gas diffusion equation after pseudopressure transformation. The presented work is an extension of the work of Barreto (2011) into Laplace domain. The
extension of the solution into Laplace domain is an advantage as less effort is required for numerical integration.
Moreover, solutions of different well and reservoir geometries in pressure transient analysis are broadly available in
Laplace domain. Field application of the solution will involve analysis of gas-rate data after deconvolution.
iii
TABLE OF CONTENTS
ABSTRACT .................................................................................................................................................. iii
LIST OF FIGURES ....................................................................................................................................... vi
LIST OF TABLES......................................................................................................................................... ix
ACKNOWLEDGEMENTS ............................................................................................................................ x
CHAPTER 1
INTRODUCTION ................................................................................................................ 1
1.1
Organization of the Thesis .......................................................................................................... 2
1.2
Motivation of Research .............................................................................................................. 3
1.3
Objectives ................................................................................................................................... 4
1.4
Method of Study ......................................................................................................................... 4
CHAPTER 2
2.1
LITERATURE REVIEW ..................................................................................................... 6
Convolution and Deconvolution ................................................................................................. 6
2.1.1 Superposition and Convolution ............................................................................................. 6
2.1.2 Rate Normalization ............................................................................................................... 8
2.1.3 Deconvolution ....................................................................................................................... 8
2.2
Solution of Non-Linear Diffusivity Equation ........................................................................... 10
CHAPTER 3
DEVELOPMENT OF THE CUBIC SPLINE BASED DECONVOLUTION METHOD . 12
3.1
Convolution and Deconvolution ............................................................................................... 12
3.2
Laplace Transformation of Sampled Functions Using Cubic Splines ...................................... 13
3.3
Inverse Mirroring at Boundaries ............................................................................................... 16
3.4
Adaptive Cubic Spline .............................................................................................................. 16
3.5
The Iseger Algorithm................................................................................................................ 19
3.5.1 Verification Examples Using Iseger Algorithm .................................................................. 22
3.5.2 Discontinuous and Piecewise Differentiable Functions ...................................................... 23
3.5.3 Wellbore Pressure Solution for a Combined Drawdown and Buildup................................ 25
3.6
Adaptive Nonparametric Kernel Regression ............................................................................ 25
3.7
Deconvolution of Pressure Responses for a Sequence of Step-Rate Changes .......................... 28
iv
3.8
Field Examples ......................................................................................................................... 32
3.9
Sandface Rate Deconvolution – Muenuier et al., (1985) Example ........................................... 32
3.10
Sandface Rate Deconvolution – Fetkovich and Vienot (1984) Example.................................. 33
3.11
Deconvolution of Variable-Rate Data-Shale-Gas Well ............................................................ 33
CHAPTER 4 LAPLACE TRANSFORMATION SOLUTION TO THE NONLINEAR
DIFFUSIVITY EQUATION............................................................................................................. 35
4.1
Mathematical Model ................................................................................................................. 35
4.2
Asymptotic Expansion (Perturbation) ...................................................................................... 39
4.3
Variable Gas Rate Deconvolution ............................................................................................ 45
4.4
Validation ................................................................................................................................. 45
4.4.1 Solution for Wells in Infinite Slab Reservoirs .................................................................... 45
4.4.2 Solution for Wells in Closed Cylindrical Reservoirs .......................................................... 56
4.5
Discussions ............................................................................................................................... 61
CHAPTER 5
CONCLUSIONS AND RECOMMENDATIONS.............................................................. 66
5.1
Conclusions .............................................................................................................................. 66
5.2
Recommendations .................................................................................................................... 66
NOMENCLATURE ..................................................................................................................................... 68
REFERENCES ............................................................................................................................................. 71
APPENDIX A
EXAMPLE OF CONVOLUTION INTEGRAL ............................................................... 74
v
LIST OF FIGURES
Figure 3.1 Interpolation of discontinuous pressure data using cubic spline. ....................................... 17
Figure 3.2 Inverse Mirroring at discontinuous points (t  50,100, 200) . ............................................ 17
Figure 3.3 Application of inverse mirroring at discontinuous points using cubic spline
interpolation. ............................................................................................................................. 18
Figure 3.4 Inverse mirroring at both ends. .......................................................................................... 20
Figure 3.5 Application of adaptive cubic spline using simulated data. Circles show where the
adaptive cubic spline is required............................................................................................... 20
Figure 3.6 Application of adaptive cubic spline using field data. Circles show where the
adaptive cubic spline is required............................................................................................... 21
Figure 3.7 Numerical Inversion of Unit Step Function by the Iseger algorithm; Effect of the
number of inversion points. ...................................................................................................... 22
Figure 3.8 Effect of nrp in the numerical inversion of a function by the Iseger algorithm. ................ 22
Figure 3.9 Numerical Inversion of a unit step function; comparison of Iseger and Stehfest
algorithms. ................................................................................................................................ 23
Figure 3.10 Numerical Inversion of a function with multiple step changes; comparison of Iseger
and Stehfest algorithms. ........................................................................................................... 24
Figure 3.11 Numerical Inversion of a piecewise differentiable function; comparison of Iseger
and Stehfest Algorithms. .......................................................................................................... 24
Figure 3.12 Numerical Inversion of wellbore pressure change for a combined drawdown and
buildup sequence. ..................................................................................................................... 26
Figure 3.13 Optimum fixed bandwidth obtained from cross-validation method. ............................... 28
Figure 3.14 Local bandwidth factor, i , applied in adaptive Nadaraya-Watson kernel estimator.
Higher lambda means higher dispersion at particular time. ...................................................... 29
Figure 3.15 Nadaraya-Watson kernel estimator using Gaussian kernel regression using fixed and
adaptive bandwidth (   0.5 ) showing improvement over ordinary kernel regression. .......... 29
Figure 3.16 Step-rate sequence for data in Table 2. ............................................................................ 30
Figure 3.17 Pressure changes corresponding to step-rate sequence shown in Figure 3.16. ................ 31
Figure 3.18 Deconvolution of pressure responses in Figure 3.17 for the step-rate sequence in
Figure 3.16. ............................................................................................................................... 31
Figure 3.19 Sandface-rate deconvolution to remove wellbore storage effect; Muenuier et al.
(1985) example. ........................................................................................................................ 32
vi
Figure 3.20 Sandface- Numerical inversion of pressure drop from tabulated data; Fetkovich and
Vienot (1984) example. ............................................................................................................ 33
Figure 3.21 Example of flowing gas rate and corresponding pseudo-pressure. .................................. 34
Figure 3.22 Numerical inversion of pressure drop from tabulated data: Shale gas application for
data shown in Figure 3.21. ........................................................................................................ 34
Figure 4.1 Omega factor (w 
 ct  (  ct )i
) versus dimensionless pressure in an infinite acting
(  ct )i
reservoir. ................................................................................................................................... 52
Figure 4.2 Omega factor (w 
 ct  (  ct )i
) versus dimensionless time in an infinite acting
(  ct )i
reservoir. ................................................................................................................................... 53
Figure 4.3 Dimensionless pseudo-pressures ( mD (0) ) for the data set presented in Table 4.1
(Infinite Acting Reservoir). ...................................................................................................... 53
Figure 4.4 Dimensionless first non-linear term of pseudo-pressure ( mD (1) ) for the data set
presented in Table 4.1 (Infinite Acting Reservoir). .................................................................. 54
Figure 4.5 Dimensionless pseudo-pressure ( mD  mD(0)  mD(1) ) for the data set presented in
Table 4.1 (Infinite Acting Reservoir). ...................................................................................... 54
Figure 4.6 Dimensionless pseudo-pressure ( mD (1) ); comparison of two different rates for the
reservoir and gas ....................................................................................................................... 55
Figure 4.7 Dimensionless pseudo-pressure ( mD (0) ) for data set presented in Table 4.2 (Closed
Cylindrical Reservoir). ............................................................................................................. 62
Figure 4.8 Dimensionless first non-linear term of pseudo-pressure ( mD (1) ) for data set presented
in Table 4.2 (Closed Cylindrical Reservoir). ............................................................................ 62
Figure 4.9 Dimensionless pseudo-pressure ( mD  mD(0)  mD(1) ) for data set presented in Table
4.2 (Closed Cylindrical Reservoir). .......................................................................................... 63
Figure 4.10 Dimensionless pseudo-pressure ( mD (1) ) comparing two different rates for reservoir
and gas properties presented in Table 4.2 (Closed Cylindrical Reservoir). The results
show that, mD (1) is rate dependent. ........................................................................................... 63
Figure 4.11 Comparing Green’s function for both Laplace and time domain at rD  1 . .................... 65
Figure 4.12 Comparing w(tD ' )
mD(0) (tD' )
tD'
Figure A.1 In an infinite acting reservoir
for both Laplace and time domain at rD  1 . .................. 65
mD (tD )
is a decreasing function as time increases. ......... 75
tD
vii
Figure A.2 In an infinite acting reservoir Greens function as a decaying function as time
increases. .................................................................................................................................. 75
tD
Figure A.3 In an infinite acting reservoir

0
mD (tD' )
tD'
GD (tD  tD ' )dtD ' is a decaying function as
time increases. .......................................................................................................................... 76
viii
LIST OF TABLES
Table 3.1 Data for drawdown followed by build up ............................................................................ 25
Table 3.2 Data for deconvolution example ......................................................................................... 30
Table 4.1 Data for infinite acting reservoir. ........................................................................................ 55
Table 4.2 Data for closed cylindrical reservoirs. ................................................................................. 64
ix
ACKNOWLEDGEMENTS
First and foremost, I give thanks to the One above all of us, the omnipresent God, for giving me the strength to
march on and complete my PhD– thank you so much Dear Lord.
I would also like to express my gratitude to my compassionate advisor Dr. Erdal Ozkan for his continuous
support and valuable advice throughout my studies at Colorado School of Mines. I am highly grateful for his
immense gentleness, enthusiasm, endurance, inspiration, and motivation. His exceptional guidance and counsel has
changed my life course several times throughout my PhD research journey. “Thank You” is just not enough - thank
you very, very much.
A special thanks to my co-advisor Dr. Luis Tenorio and my thesis committee member Dr. Paul Martin for their
support, help, and guidance during my research.
I thank my thesis committee members: Dr. Hossein Kazemi, Dr. Vaughan Griffiths, and Dr. John Humphrey for
their help and guidance.
I would like to thank Dr. Mahadevan Ganesh from the Mathematics Department for all the help and insight
throughout my research.
I would also like to thank MI3 Petroleum Engineering for strong support and my colleague Mr. Oscar G.
Gonzalez for his assistance and direction.
The Marathon Center of Excellence for Reservoir Studies (MCERS) at Colorado School of Mines has been
invaluable as well as the camaraderie with my school colleagues: Ali, Ayyoub, Elham, John, Mehdi, Reza, Nasser,
Shirin, Mojtaba, Najeeb, Farshad, and Younki.
I am grateful for Denise Winn-Bower, for helping me with A LOT of administrative issues during my current
study at CSM.
Finally, I want to give a special thanks to my wife Elham for her support and motivation to accomplish this
endeavor and to my great father, mother, and family who have always been my guiding light. Their prayers have
always been my inspiration. I love you all.
x
CHAPTER 1
INTRODUCTION
This dissertation presents the results of a study for a Doctor of Philosophy degree conducted at Marathon Center
of Excellence for Reservoir Studies (MCERS) in the Petroleum Engineering Department at Colorado School of
Mines.
The main objective of the dissertation is to extend the conveniences of the deconvolution to non-linear problems
of fluid flow in porous media.
The usual to approach to extend deconvolution procedures to non-linear problems in oil and gas reservoirs is to
linearize the non-linear diffusion equation in terms of a pseudo-pressure and then to apply the deconvolution. It is
well known that the pseudo-pressure approach does not completely remove the nonlinearity, but, for practical
purposes, the remaining nonlinearity is assumed to be weak and ignored. In this dissertation, the solution of the nonlinear problem is obtained by a perturbation approach, which presents the solution as a series of solutions of linear
problems. This approach permits term-by-term application of deconvolution. One practical problem still remains:
The proposed perturbation solution is more conveniently evaluated in the Laplace-transform domain. The standard
deconvolution algorithms, however, are in the time-domain. Thus, the development of an efficient deconvolution
procedure in the Laplace-transform domain is a prerequisite to accomplish the main objective of this dissertation.
Therefore, the main objective of the dissertation can be expanded into two sub-objectives.
The first sub-objective of the dissertation appertains to the analysis of variable-rate reservoir performance data
in Laplace domain by deconvolution. The approach taken in this work leads to the deconvolution of variable-rate
data in the Laplace domain. Specifically, an approximate function is required to take sampled (tabulated) production
rate and pressure data into Laplace domain. Furthermore, the step-changes in the production rate during shut-in
periods lead to inaccuracy in approximating functions and instability in numerical Laplace inversion algorithms. In
this study, a cubic-spline method with piecewise linear interpolation and boundary mirroring is developed in
Laplace domain to approximate and transform the production rate and bottom-hole pressure into the Laplace
domain. This algorithm accurately transforms sampled (tabulated) data into Laplace domain and eliminates the
numerical inversion instabilities at discontinuous points or boundaries commonly encountered in the piecewise
linear approximations of the data. The developed approach does not require modifications of scattered and noisy
data or extrapolations of the tabulated data beyond the end values.
Rate and pressure measurements of wells usually include some level of noise, and due to the nature of the
deconvolution process (more specifically, deconvolution in Laplace domain), the computed underlying constant-rate
response will display oscillations, which requires some degree of smoothing. To smooth the deconvolved pressure
response, an adaptive approach using a Gaussian and Epanechnikov kernel regression is proposed. The adaptive
kernel regression proposed herein is shown to be more successful than the normal kernel regression.
Since the deconvolution algorithm and the approximating function for tabulated data are in the Laplace domain,
the solution requires a numerical Laplace inversion algorithm. Common Laplace inversion algorithms usually face
accuracy problems in dealing with functions including contributions from step-changes and discontinuities. For this
1
particular problem, a numerical Laplace inversion algorithm developed by Peter Iseger (2006) and introduced to the
Petroleum Engineering field by Al-Ajmi et al. (2008) is used. The application and the accuracy of Iseger’s numerical
Laplace inversion algorithm for functions including step-changes and discontinuities in pressure-transient analysis
have been validated in the work presented by Al-Ajmi et al. (2008). The algorithm presented by Iseger (2006)
removes the restriction of continuity (that is, it tolerates piecewise continuous functions) and provides opportunities
for many practical applications. The Iseger algorithm has been tested for several common conditions requiring the
use of piecewise-differentiable and discontinuous functions including:

the use of tabulated data in the Laplace transform domain

problems involving step-rate changes

build-up tests following a drawdown period

Mini-DST tests with a sequence of drawdown and buildup periods.
Due to the significant contribution to this research, one section of this study will be devoted to Iseger’s numerical
Laplace inversion algorithm and its verification examples. The critical parameters of the Iseger algorithm will also
be reviewed.
The second sub-objective of the dissertation involves the fundamental perturbation solution of the nonlinear gas
diffusivity equation in Laplace domain. The fundamental perturbation solution presented here will be an extension
of the solution developed by Barreto (2011), which is in terms of Green’s functions but accounts for the nonlinearity
caused by the dependency of gas properties on pressure.
The original work of Barreto (2011) was presented in the time domain. The objective of the present work is to
enhance the utility of Barreto’s approach by using the conveniences offered by the properties of Laplace transforms.
The expected benefits from extending Barreto’s solution to Laplace domain include:1) solutions of different well
and reservoir geometries in pressure transient analysis are broadly available in Laplace domain, 2) less effort is
required for numerical integration, and 3) the solution will be utilized in convolution and deconvolution in Laplace
domain as presented herein.
1.1
Organization of the Thesis
This dissertation is presented in five chapters. The preliminary results are also discussed.

Chapter 1, the introduction, describes the research objectives, the methodology, and the motivation of the
study.

Chapter 2 is the literature review of publications and papers along with the proposed research.

Chapter 3 addresses mainly the first sub-objective of the dissertation; the derivation of a cubic spline based
deconvolution method. It includes the following:
o
The mathematical derivation for cubic spline method adapted with piecewise linear interpolation
in Laplace domain used in the deconvolution process.
o
Inverse mirroring at both boundaries to reduce the oscillation of the response function recovered
from the deconvolution process.
2
o
Smoothing of the constant pressure response recovered from the deconvolution process using
adaptive kernel regression.
o
A brief review of Iseger’s algorithm, along with verification examples and critical parameters. For
demonstration purposes, simulated and field data are also examined.

In Chapter 4 a perturbation solution to the nonlinear diffusivity equation is presented in the Laplacetransform domain. Although the non-linearity displayed in the gas diffusivity equation in terms of psudopressure is week, it is used for the demonstration of the solution approach in this chapter because it is the
best-known, non-linear flow problem in porous media.

Chapter 5 reviews the results of the research accomplished in this PhD dissertation and proposes
recommendations for future work.
1.2
Motivation of Research
Interest in the analysis of variable-rate production data (rate and pressure measurements as function of time) has
increased in the oil and gas industry over the last two decades due to the popularity of mini-DST tests and the
increased utilization of unconventional, tight oil and natural-gas resources. In mini-DST tests, a sequence of
production and shut-in periods gives rise to a piecewise-continuous rate and pressure behavior. In the case of
unconventional tight reservoirs, well performance data are usually available in the form of daily or monthly
production and the corresponding tubing-head pressures. However, the standard theory of pressure-transient and
well-performance analysis is based on the assumption of constant production rate or bottom-hole flowing pressure.
Duhamel’s principle expresses the bottom-hole pressures of a variable-rate production case in the form of a
convolution relationship between the variable flow rate and the pressures for the unit constant-rate production (the
influence function). Therefore, due to the immense utility in data analysis, recovering the influence function of the
convolution integral; that is, deconvolution of variable-rate responses to obtain the underlying constant-rate
responses, is of great interest in petroleum engineering.
van Everdingen and Hurst (1949) introduced the application of convolution/deconvolution in Laplace domain for
the solution of common transient-flow problems in porous media. Specifically, they highlighted the conveniences
due to the fact that the convolution of two functions turns into an algebraic product of the functions in the Laplace
transform domain. However, although convolution/deconvolution in Laplace domain provides a convenient means
of generating analytical solutions for many variable-rate problems, its use for the measured (tabulated) data is not
straightforward due to the requirement that discretized (measured and tabulated) rate and pressure data be
transformed into Laplace domain.
Several approaches have been tested regarding the deconvolution of tabulated variable-rate and pressure data in
Laplace domain (Roumboutsos-Stewart (1988) and Onur-Reynolds, (1998)). These algorithms have suffered from
problems in the numerical Laplace inversions of piecewise continuous functions. This was a problem due to the
limitations of the existing numerical inversion algorithms. Another obstacle was the construction of approximating
functions representing sampled (tabulated) data in Laplace domain. The approaches of these algorithms to this
problem required extrapolation of the sampled data beyond the limits of the sampling interval to evaluate the
3
Laplace integral over the semi-infinite positive domain. This usually led to an oscillatory deconvolution result at the
late portion of the tabulated data (usually referred to as tail effect).A regularization algorithm or some type of
regression procedure was also needed to smooth the constant pressure response in all scenarios.
Moreover, complications exist in the application of Duhamel’s principle (convolution/deconvolution) regarding
the flow of real gases in porous media. Application of Duhamel’s principle requires a linear system but the
diffusivity equation of compressible fluids (gas) is non-linear. The conventional approach of using a pseudopressure definition does not completely remove the pressure dependency of the gas properties and the
convolution/deconvolution of gas-well response remains non-rigorous. Therefore, a new approach is required to
extend Duhamel’s principle to variable-rate problems of gas-well performances.
1.3
Objectives
This research aims to improve the applications of convolution/deconvolution in Laplace domain and the
extension of these applications to the non-linear problems of gas flow. The specific objectives for the improvement
of Laplace-domain convolution/deconvolution include the following:
1.
Develop alternatives to existing approaches to transform tabulated data to Laplace domain.
2.
Introduce techniques to reduce oscillatory deconvolution behavior caused by the tail effect or
discontinuities of the input functions.
3.
Consider the impact of the numerical inversion algorithms on the success of Laplace-domain
deconvolution for different methods of transforming tabulated data into Laplace domain.
4.
Examine the possibility of smoothing deconvolved responses without interfering with the physical
signatures of the data.
5.
Rigorously verify the theoretical basis of the ideas used in the above developments.
6.
Demonstrate the use of Laplace-domain deconvolution with simulated and field data and compare
the applications with the existing algorithms.
For the extension of Duhamel’s principle to the non-linear problems of gas flow in porous media, the following
objectives are defined:
1.
Consider the use of the perturbation approach introduced by Barreto (2011) to obtain an
approximate Laplace transformation of the non-linear gas-diffusivity equation.
2.
1.4
Demonstrate the use of the deconvolution algorithm for the analysis of variable gas-rate field data.
Method of Study
The main method of this study is analytical. Existing mathematical-physical techniques are used to formulate the
problems mathematically and to obtain solutions at a desirable form and level of accuracy. Perturbation and infinite
series expansion are used to approximate solutions to non-linear diffusion equation. The non-linear solutions are
obtained in the form of small additions from a series of linearized problems. The solutions of the linearized diffusion
equation are expressed in terms of appropriate Green’s functions. Semi-analytical and numerical techniques are used
to compute solutions and approximating functions are invoked to deal with tabulated data.
4
Because the linearized problems encountered in this research lend themselves to convolution-type solutions,
Laplace transformation is preferred due to the simplification of convolution relationship in the transform domain.
This also brings the necessity to deal with the numerical inversion of the results in the Laplace transform domain.
Various approximating functions are used to represent tabulated data and kernel regression is used to smooth the
solutions.
To validate the developed solutions both simulated and field data are used. The results are then compared to
commercial software results. For computational purposes, Stehfest and Iseger numerical Laplace inversion
algorithms are employed. The computational codes are written in both Fortran 77 and Fortran 90.
5
CHAPTER 2
LITERATURE REVIEW
This chapter is divided into two parts: the first part discusses the convolution and deconvolution of variable-rate
data for well-test analysis and the second part discusses the development of a solution for the non-linear gasdiffusivity equation and application of deconvolution to variable-gas-rate data in the Laplace domain.
2.1
Convolution and Deconvolution
The solution to variable-rate problems has been widely popular in petroleum engineering literature since the
seminal work of van Everdingen and Hurst (1949). In regards to this subject, various approaches and references can
be found in the literature. The approaches available in the literature have been listed and categorized by Ilk (2005).
The three main categories are as follows:
I.
II.
III.
2.1.1
Superposition and Convolution
Rate Normalization and Material Balance Deconvolution
Deconvolution
Superposition and Convolution
The term superposition simply states that, for all linear systems, the net response of the total system at a given
point in space and time is the summation of all the individual parts/stimulus that contributes to the total system.
Petroleum engineers apply superposition in space to construct solutions for bounded reservoirs, complex well
geometries, and multiple-well problems. These applications are also known as the method of images. Superposition
in time is used to construct solutions for pressure buildup and variable-rate problems from constant rate solutions.
The governing equation for fluid flow in porous media is the diffusivity equation. For the flow of liquids, due to
the small and constant fluid compressibility and small pressure gradients, the diffusion equation can be linearized for
most practical purposes and the superposition principle becomes valid. In the case of gas flow, however, fluid
properties are strong functions of pressure and the nonlinear terms in the diffusivity equation cannot be removed by
physically acceptable assumptions. A standard approach is, then, to recast the diffusivity equation in terms of a
pseudo-pressure. The pseudo-pressure definition groups the pressure and pressure-dependent fluid properties in a
new pseudo-variable in which the nonlinearity of the diffusion equation becomes weaker. For some applications,
such as pressure-transient analysis, the remaining nonlinearity of the diffusion equation may be ignored if the
viscosity-compressibility product does not change considerably over the infinite-acting period of the well test. This
assumption, however, has been shown not to hold during boundary-dominated flow of high-flow-rate wells
(Raghavan 1993). One of the major consequences of this condition is the lack of true pseudo-steady state flow in gas
wells, which makes many of the standard performance prediction tools for gas wells questionable. To remove the
nonlinearity completely, a pseudo-time definition has also been proposed, but its validity could not be agreed upon
except for pressure-buildup applications (Raghavan 1993). The details of the use of pseudo-pressure and pseudotime can be found in the works of Al-Hossainy et al. (1966) and Agarwal (1980).
6
In mathematics, particularly in functional analysis, convolution is defined as a mathematical operation on two
functions f (t ) and g (t ) producing a third function that is typically viewed as a modified version of one of the
original functions is translated. In other words, the overlap amount between two functions f (t ) and translated
function of g (t ) can be represented by convolution. The mathematical form of the convolution of two functions f (t )
and g (t ) is defined as:
def t
z (t )  f (t ) * g (t ) 

t

f ( ) g (t   )d  f (t   ) g ( )d
0
(2.1)
0
A discrete form for the right hand side of Eq. (2.1) over a finite interval can be written as:
n
z (t ) 
 f (
i 1 ) g (t
  i 1 ) i
(2.2)
i 1
where z (t ) is the response of the system or the convolution of the functions f (t ) and g (t ) and  is a dummy variable.
This is also called the superposition equation, which was previously defined as the total system response and is
equal to the summation of individual stimulus that contributes to the total system.
In regards to petroleum engineering applications utilizing linearized diffusivity equation, the pressure response
of the reservoir to a sequence of fluid withdrawals at a source location can be obtained by the superposition of the
pressure responses for instantaneous withdrawals of variable volumes of fluid (variable-strength instantaneous
sources) over an interval of time. This operation can be written in the form of the convolution of the rate variations
(input) with the system’s response to an impulse function (the unit instantaneous source), which is also known as
Duhamel’s principle. In petroleum engineering terms, Duhamel’s principle states that the observed wellbore
pressures drop for a variable rate can be written as the convolution of the input rate function and the derivative of the
impulse function or constant-rate pressure response. It is also assumed that the system is in equilibrium initially
 p(r , t  0)  pi  . The convolution integral in dimensionless form is defined in the petroleum engineering literature
as:
tD
pwD (tD ) 

0
dqD ( )
 psD (t D   )d 
d 
tD
q
D ( )
d  psD (tD   )
d
0
d
(2.3)
where pwD is the dimensionless pressure for variable-rate production (system response), qD is the dimensionless
variable production rate, and psD is the dimensionless pressure for the constant-rate solution (influence function),
which is equal to
psD  pD (tD )  s
(2.4)
Here s represents the mechanical skin damage around wellbore. The equivalent equation to Eq.(2.2) or the
discretized form for the right-hand side of Eq.(2.3) in terms of superposition can also be written as:
n
pwD (t D ) 
 (q  q
i
i 1
i 1 ) 
 pD (t
 ti 1 )  s 
(2.5)
7
van Everdingen and Hurst (1949) utilized Duhamel’s principle to obtain a dimensionless wellbore pressure-drop
solution for a continuously varying rate. The method of Odeh and Jones (1965) or the methods of Soliman (1981)
and Stewart et al. (1983) can be derived from van Everdingen and Hurst’s solution. Generally, their approaches use
the exponential integral solution to the constant-rate wellbore response or its semi-log approximation. Fetkovich and
Vienot (1984) and Agarwal (1980) also used Duhamel’s principle to generate solutions for the variable-rate
production problems. All above methods, however, were limited to forward modeling applications since a particular
reservoir model had to be selected in advance for the constant-rate pressure function under the convolution integral.
The analysis of variable-rate production data by the modern regression analysis techniques requires backward
modeling; that is, the constant-rate pressure function under the convolution integral is determined by knowing the
field response.
2.1.2
Rate Normalization
Rate normalization is generally used when the rate is smoothly changing as a function of time. The rate
normalization method was introduced by Gladfelter et al. (1955) for the analysis of wellbore-storage-dominated
pressure buildup responses. Winestock and Colpitts (1965) employed the rate normalization for drawdown test
analysis. The rate normalization method is also used to remove the effect of wellbore storage when the sand face
rates are available. The appeal of the rate normalization method is in its simplicity (it only requires an algebraic
operation of the rate and pressure data). Raghavan (1993) showed that the rate normalization is valid if the
convolution integral can be approximated as follows:
p(t )  q(t ) pc (t )
(2.6)
This approximation is shown to be valid when the rate and pressure are smoothly varying functions of time, which is
a restrictive condition for most production data analysis applications (Raghavan 1993).
2.1.3
Deconvolution
The analysis of pressure drop over time at constant production rate is a primary goal in well-test analysis.
Identifying reservoir characteristics from constant/unit rate pressure response by pressure-transient analysis is a
well-established practice. When variations exist in rate, it is desirable to convert the data to constant-rate production
responses before analysis; in other words, the first step in analysis becomes the identification of the unit-rate
pressure response from the convolution integral. In mathematics, the convolution integral is a linear Voltera integral
equation of the first kind. Deconvolution of variable-rate production, mathematically, is an ill-conditioned inverse
problem (Lamm, 2000), and a small error in the data or in the subsequent calculations results in much larger errors
in the answers. It should be re-emphasized that convolution/superposition and rate normalization requires that a
model to be utilized for the influence function while in the deconvolution of variable-rate production, finding the
influence function is the objective.
The deconvolution approaches of variable-rate production data proposed in the literature can be grouped into two
main categories; time domain methods and spectral methods. The most successful deconvolution algorithm in time
8
domain is the one recently presented by von Schroeter et al. (2004) and Levitan (2003). In the progression of the
time-domain deconvolution algorithms for variable-rate data, instability of the influence function was an issue and
some artificial approaches were designed to make it stable. Coats et al. (1964) employed the prior knowledge of
constant-rate response (non-negativity) before applying constraints in their model. These additional constraints in
the model proposed by Coats (1964) increased the number of unknowns, which leads to an over-determined system.
An over-determined system requires some minimization. Kuchuk (1985) and Baygun (1997) attempted to improve
Coat’s method by applying a least squared approach and by using different constraints. Issues related to the
constraints were solved by the work presented by Schroeter et al. (2004). According to Schroeter et al. (2004),
instead of the rate-normalized pressure derivative itself, they estimated its logarithm, which makes explicit sign
constraints unnecessary. The deconvolution algorithm proposed by Schroeter et al. (2004) is a nonlinear, total least
squares problem, which means both pressure and rate data are influenced by noises. In the presence of a certain level
of noise in the input data, the deconvolved pressure response could still show some instability. Regularization or
eliminating irregular response behavior is also proposed by Schroeter et al. (2004). Various ways for the
regularization of the impulse function may be found in literature. The regularization in the work by von Schroeter et
al. (2004) is based on a measure of overall curvature. Levitan (2004) presented some modifications on Shroeter et
al.’s work proposing some critical recommendations to obtain accurate results.
Spectral methods more specifically, the Laplace transform, became a standard tool in petroleum engineering for
the analysis of transient flow problems after the classical paper of van Everdingen and Hurst (1949). Performing
deconvolution in the Laplace domain simplifies the convolution integral to the product of the transforms of the
functions. Eq. (2.4) in the Laplace domain can be written as follows:

 p wc ( s) 

 p w ( s)
(2.7)

s q( s)
Once the tabulated pressure and the rate data are transformed into Laplace domain, the impulse function for
constant-production rate can be obtained by the simple division of the functions as shown in Eq.(2.7). The results
can be numerically inverted back to the time domain using a numerical Laplace inversion algorithm, such as those
proposed by Stehfest (1970) or Iseger (2006). In spite of the simplicity of using Eq.(2.7), Laplace transformation
suffers from two issues: firstly, Laplace transformation is defined on the entire positive axis of time from zero to
infinity, which requires some extrapolation techniques. Roumboutsos et al. (1988) and Onur and Reynolds (1998)
proposed piecewise linear and polynomial functions for transforming the data into Laplace domain with some
extrapolation techniques; but neither of them was successful to overcome the foregoing issue. Secondly, Laplace
transform of piecewise-continuous functions is well defined; however, discontinuities in the Laplace domain
function cause problems in the main stream numerical inversion algorithms. For example, the numerical inversion
algorithm of Stehfest (1970), which is commonly used in petroleum engineering, cannot handle the discontinuities in
the function to be inverted. In the approach proposed by Ilk (2005), a combination of B-splines and numerical
inversion of the Laplace transform is used. The difference between their methods and other methods in the Laplace
domain is in the fact that B-spline was used to represent the unknown response solution. Similar to the other
methods, a limitation in their method is that the non-negativity of the response function was not ensured.
9
Additionally, their regularization approach seems weaker, which makes their algorithm less tolerant to the noise in
pressure and rate data. A comparative study of recent deconvolution algorithms can be found in the work by Cinar
(2006).
As previously explained, solutions to numerous pressure transient analysis problems in the field of petroleum
engineering are provided in the Laplace transform domain. Usually these solutions do not lend themselves to an
analytical inversion, and numerical inversion is normally the only resort. The common numerical inversion
algorithm presented by Stehfest (1970) is applicable only to continuous functions and fails where discontinuities
emerge. Some other algorithms, such as Bellman (1966), Crump (1976), and Talbot (1979), were tried but their
complexity and inflexibility limited their application. The simple form of the convolution integral in Laplace
transformation and the high interest in deconvolution algorithms encouraged many researchers to look for new
Laplace inversion algorithms. Ilk (2005) applied Gaver-Wynn-Rho’s algorithm presented by Valko and Abate
(2004) in the B-spline deconvolution.
Iseger (2006) presented a new algorithm, which removes the necessity of the function to be continuous. In
addition, the algorithm provides remarkable results in the case of piecewise continuous and piecewise differentiable
functions. In well test analysis, such functions can be tabulated as data in the Laplace transform domain,
deconvolution algorithms, and solutions that include step-rate changes in the mini-DST tests. Adapting cubic spline
with piecewise linear approximation in the Laplace domain and employing inverse boundary mirroring and kernel
regression using a Gaussian or Epanechnikov kernel function is shown in this dissertation to improve the behavior of
constant pressure response obtained by deconvolution.
2.2
Solution of Non-Linear Diffusivity Equation
The general diffusivity equation, which governs fluid flow in porous media is a non-linear equation. For slightly
compressible fluids (liquids) the variation of fluid viscosity and compressibility with pressure can be assumed
negligible. However for compressible fluids (gas) this is not a valid assumption.
The diffusivity equation for the flow of a real gas in porous media is given by:
.(
 ct p p
p
p) 
z
k  z t
(2.8)
Eq.(2.8) is non-linear due to the fact that the gas properties (viscosity, compressibility and real gas deviation factor)
are pressure dependent. The pseudo-pressure transformation introduced by Al-Hussainy et al. (1966) attempts to
linearize the diffusion equation and allows an analytical solution of Eq.(2.8). The pseudo-pressure proposed by AlHussainy et al. (1966) is defined by
p
m( p )  2
  ( p ') z( p ')dp '
p'
(2.9)
pb
where
pb is a reference pressure. Upon applying the pseudo-pressure transformation, the gas diffusivity equation can
be written as
10
2 m 
 ct m
k
t
=
1 m
 (m) t
(2.10)
where the diffusivity term is defined as:
 (m) 
k
(2.11)
 (m)ct (m)
Eq.(2.10) is not completely linear since the diffusivity term given in Eq.(2.11) is a function of pressure. For
practical purposes, the quantity  (m)ct (m) is approximated as constant at the average pressure of the drainage area.
Some other approaches have been tried for the linearization of the gas diffusivity equation. Kale and Mattar
(1980) derived an approximate solution using a perturbation technique for radial flow with constant rate. An exact
solution using a perturbation approach was developed for constant-rate production in an infinite reservoir by Serra
and Reynolds (1990a, 1990b). In their derivation the Boltzmann transformation is employed.
In the case of variable-rate production, Duhamel’s principle can be applied if the governing equation of fluid
flow is linear. For slightly compressible fluids (liquid) this is a routine approach to well test problems. However, for
gas flow, superposition in time cannot be applied, as the gas diffusivity equation is non-linear. If the diffusivity term
 is assumed to be constant, then the theory of superposition can be applied as in the works of Gupta and Andsage
(1967), Samaniego and Cinco-Ley (1991), and recently von Schroeter and Gringarten (2007).
Recently, Barreto (2011) introduced a solution to the diffusivity equation with a pressure-dependent diffusivity
term. He obtained the solution of the nonlinear diffusion equation using a perturbation approach and Green’s
Functions. Barreto’s solution is in the time domain, which requires double and/or triple numerical integrations. In
this research, Barreto’s approach is extended to Laplace domain. There are several advantages to extending
Barreto’s solution to Laplace domain: 1) solutions for different reservoir geometries in petroleum engineering are
broadly available in the Laplace domain (Ozkan and Raghavan, 1991a, b), 2) because the time integral of the
Green’s functions drops off in Laplace domain, less effort will be required for the numerical evaluation of the
solution, 3) the Laplace-domain deconvolution method proposed in this work can be used with the solution of the
nonlinear diffusion equation in Laplace domain. The solution proposed in this work may be applied to
deconvolution of gas well-test data, especially in shale-gas reservoirs.
11
CHAPTER 3
DEVELOPMENT OF THE CUBIC SPLINE BASED DECONVOLUTION METHOD
Despite unequivocal advantages of using sampled well-performance data in the Laplace transform domain, timedomain analyses of pressure and production data have been more popular lately. This is due to the fact that
unresolved problems in the transformation of sampled data into Laplace domain as opposed to the demonstrated
success of the recent real-time deconvolution algorithms. However, the transformation of sampled data into Laplace
domain has a broader range of applications than deconvolution and the limited success of the past approaches to
transforming tabulated data to Laplace domain, such as piece-wise linear approximations, is an algorithmic issue;
not a fundamental defect. Specifically, an adequate algorithm to transform the piecewise-continuous sampled data
into the Laplace space and an appropriate numerical Laplace inversion algorithm capable of processing the
exponential contributions caused by the tabulated data are essential to exploit the potential of Laplace domain
operations.
In this chapter, we introduce a new algorithm, which uses inverse mirroring at the points of discontinuity and
adaptive cubic splines to approximate rate or pressure versus time data. This algorithm accurately transforms
sampled data into Laplace space and eliminates the Numerical inversion instabilities at discontinuities or boundary
points commonly encountered with the piece-wise linear approximations of the data. The approach does not require
modifications of scattered and noisy data or extrapolations of the tabulated data beyond the end values.
Practical use of the algorithm presented in this chapter has applications in a variety of Pressure Transient
Analysis (PTA) and Rate Transient Analysis (RTA) problems. A renewed interest in this procedure has arisen from
the need to evaluate the production performances of wells in unconventional reservoirs. With this approach, we
could significantly reduce the complicating effects of rate variations or shut-ins encountered in well-performance
data. Moreover, the approach has proven to be successful in dealing with the deconvolution of highly scattered and
noisy data. To illustrate the applications, typical field examples including shale-gas wells are presented in this
chapter.
3.1
Convolution and Deconvolution
Laplace transformation of a function, f  t  , defined for all t  0 , is given by

L  f  t   f  s  
e
 st
f  t  dt
(3.1)
0
where, f  s  , is the Laplace transformation of f  t  and
s denotes the Laplace transform parameter. Convolution of
two functions yields the algebraic product of the functions in Laplace domain as follows:
t

0

L  f1  f 2  t  = L   f1   f 2  t    d   f1  s  f 2  s 
12
(3.2)
Eq.(3.2) has important applications in transient flow in porous media associated with variable-rate production. The
solution of pressure drop, p  M , t  , at a point M and time t due to variable production rate, q  t  , is known as
Duhamel’s principle (Duhamel, 1833) and given by the convolution of flow rate, q  t  , and pc  M , t  as follows:
t

p  M , t   q   pc  M , t    d
(3.3)
0
In Eq.(3.3), pc  M , t  is the time derivative of the pressure drop at a point M and time t for a constant (unit)
production rate. In many applications of pressure transient analysis, the constant-rate pressure response at the
wellbore, pc  M , t  , is of interest. To recover the constant-rate response from the measured data, p  M , t  and
q  t  , the deconvolution of Eq.(3.3) is required.
Application of Eq.(3.2) on Eq.(3.3) yields the following algebraic expression for deconvolution in Laplace space:
pc  s  
p  s 
(3.4)
sq  s 
In the field application of the Laplace-domain deconvolution (Eq.(3.4)), the input data, p  M , t  and, q  t  need to
be transformed to the Laplace domain. This requires an approximating function to properly transform sampled
(tabulated) data into Laplace space. In addition, a numerical Laplace inversion algorithm, which handles
discontinuities and singularities, is entailed. All these requirements need to be addressed to develop a successful
deconvolution algorithm in Laplace domain as deconvolution described by Eqs.(3.3) and (3.4) is an ill-condition
problem. Moreover, production-rate data frequently show step-change behavior, which causes instability in the
numerical Laplace inversion.
The objective of this chapter is to contribute to the solution of the following problems:
3.2
Laplace Transformation of Sampled Functions Using Cubic Splines
Interpolation by polynomials of degree n is widely used in practice. For various functions, the higher quality of
the interpolation might be expected with the increasing degree of the polynomials. Unfortunately, this is not always
true and the interpolation may yield oscillatory results by higher-degree of polynomials. Losing the quality of the
interpolation by using higher-degree polynomials has been discussed by Kreyszig (1999). To avoid such
oscillations, spline methods are widely applied. The mathematical idea of the spline is to replace a single highdegree polynomial over the entire interval by several low-degree polynomials (De Boore, 2001). This is expected to
reduce the oscillation of the interpolation. In general, a spline function is a function that consists of piecewise
polynomials joined together with certain smoothing conditions. A spline function also utilizes the nth degree
piecewise polynomials to preserve (n  1)th order derivatives at the data points. A spline function is defined by knots
and the order of the spline. In the theory of splines, the points t0 , t1 ,..., tn , at which the character of the function
changes, are called knots.
13
In this work, the most popular piecewise polynomial, called the natural cubic spline, is utilized. The cubic spline
preserves the first and second derivative continuity at knots. If S (t ) is given over the interval a  t  b with knots
defined by
a  t0  t1  t2  ...  tn  b
(3.5)
then the cubic spline in each subinterval can be written in the following form:
Si (t )  A0i  A1i t  A2i t  A3i t ,
2
i  1, 2,..., n
3
(3.6)
The coefficients Aij need to be determined for each subinterval.
The simplified cubic splines function for each subinterval ti 1  t  ti with values yi 1 and yi at ti 1 and ti respectively
can be address in the form of (Atkinson, 1985)
(ti  t )3 M i 1  (t  ti 1 )3 M i
Si (t ) 
+
6(ti  ti 1 )
(ti  t ) yi 1  (t  ti 1 ) yi
(ti  ti 1 )

(3.7)
1
6
(ti  ti 1 )  (ti  t ) M i 1  (t  ti 1 ) M i 
i  2, 3,..., n
In Eq.(3.7), M represents the second derivative and ti 1  t  ti . Eq. (3.7) can also be rearranged as follows:
Si (t ) 
M i  M i 1
3
t +
6(ti  ti 1 )
 3t 2i 1 M i  3t 2i M i 1
+
6(ti  ti 1 )

3ti M i 1  3ti 1 M i
+
6(ti  ti 1 )
yi  yi 1
ti  ti 1
+
1
6
t
2
 ti M i 1  ti M i  ti 1M i 1  ti 1M i 
 t 3i M i 1  t 3i 1 M i ti yi 1  ti 1 yi


6(ti  ti 1 )
ti  ti 1


 1 2
2
  t i M i 1  ti ti 1 M i  ti 1ti M i 1  t i 1 M i
 6








i  2, 3,..., n

t

(3.8)
ti -1  t  ti
with
Applying Eq.(3.1) on cubic spline over the interval t1  t  tn for tabulated data, such as pressure or production rate,
yields

L S (t )  S ( s ) 

 S (t )e
0


dt  S (t )e
0

tn
t1
 st
 st


dt  S (t )e
 st

dt  S (t )e
t1
tn
 st
tn

 st
dt
(3.9)
t1

t1
For the set of tabulated data from t1 to tn both terms

dt = S (t )e
S (t )e  st dt and
0
 S (t ) e
 st
dt vanish because the cubic spline
tn
becomes zero outside the interval. Application of Eq.(3.9) to individual terms of Eq.(3.8) yields:
The Laplace transform of the first term:
14
tn
 M  M i 1 3 
M i  M i 1 3  st
L i
t =
t e dt
 6(ti  ti 1 )  t 6(ti  ti 1 )

0
t1

M1  M 0
 6(t  t )
1
t0
t2
3  st
t e
dt 
0
M 2  M1
 6(t
t1
2
 t1 )
tn
3  st
t e
dt ...
M n  M n 1
 6(t
tn1
  t 3i 3t 2i 6ti 6   st
  2  3  4  e
s
s
s 
M i  M i 1   s

2
6(ti  ti 1 )   t 3
3t i 1 6ti 1 6   st
i 1
   s  2  3  4  e
s
s
s 
 
n
i
n


i 2
i 1
3  st
t e
 tn 1 )







dt
(3.10)
i  2, 3,..., n
with ti -1  t  ti
The Laplace transform of the second term:
 3ti M i 1  3ti 1 M i
L
6(ti  ti 1 )

t1


3t1 M 0  3t0 M 1
t0
6(t1  t0 )
2

t =

tn

3ti M i 1  3ti 1 M i
t0
tn
2  st
t e
2  st
t e
6(ti  ti 1 )
dt ...

tn 1
3tn M n 1  3tn 1 M n
6(tn  tn 1 )
  t 2i 2ti 2   st
  2  3  e
s
s 
3ti M i 1  3ti 1 M i   s

2
6(ti  ti 1 )
  t i 1 2ti 1 2   st


 2  3 e


s
s 
  s


i 2
2  st
t e
dt
(3.11)
i
n
dt
i 1







i  2, 3,..., n
with
ti -1  t  ti
The Laplace transforms to the third and the fourth terms, respectively:
  3t 2i 1 M i  3t 2i M i 1 yi  yi 1
+


6(ti  ti 1 )
ti  ti 1
L 

 1
 +  ti M i 1  ti M i  ti 1M i 1  ti 1M i 

 6
n


i 1
 
 
 t

 
 
 
 3t 2i 1 M i  3t 2i M i 1 yi  yi 1
+

6(ti  ti 1 )
ti  ti 1

 1
 +  ti M i 1  ti M i  ti 1M i 1  ti 1M i 
 6
i  2, 3,..., n
with


 ti 1

   2
  s s


ti -1  t  ti
15
  st  ti 1 1

e 

 s s2
i
  st 

e


i 1
(3.12)
  t 3i M i 1  t 3i 1 M i ti yi 1  ti 1 yi


  6(ti  ti 1 )
ti  ti 1
L
 1 2
2
   6 t i M i 1  ti ti 1 M i  ti 1ti M i 1  t i 1 M i

n


i2



 



 t 3i M i 1  t 3i 1 M i ti yi 1  ti 1 yi
+

ti  ti 1
 6(ti  ti 1 )
 1 2
2
  t i M i 1  ti ti 1 M i  ti 1ti M i 1  t i 1M i
 6

i  2, 3,..., n


 1
 st
 st
 
e e
 s



i
i 1

(3.13)



with ti -1  t  ti
In all of above equations, L denotes Laplace operator and s denotes the Laplace transform parameter.
3.3
Inverse Mirroring at Boundaries
Discontinuous points in sampled data, such as step changes in production rate or build up in the pressure data,
may cause oscillations in the approximation functions obtained (for example, by cubic spline). Transforming the
data into the Laplace domain with such oscillations in the approximation function increases the error in the
deconvolved constant pressure response. Figure 3.1 shows an example of the oscillatory behavior around
discontinuous points caused by the application of cubic-spline approximation of the data. In addition, sampled
functions are normally available over a finite interval while Laplace transformation requires the function be defined
over the positive semi-infinite domain. To take the sampled data to Laplace domain, extrapolation from zero to the
first sampling point and from the last sampling point to infinity is required. If the behavior of the tabulated function
beyond the endpoints is known [e.g. constant wellbore storage, radial flow, pseudo steady state, etc. (Al-Ajmi et al.,
2008)], then the Laplace transform of the function can be generated by the procedures suggested by Roumboutsos et
al. (1988) and Onur and Reynolds (1998). As previously noted, however, due to the property of the cubic spline, the
Laplace integration over the regions of extrapolation vanishes. To remove the remaining oscillations of the
approximating function at discontinuity points, we propose the use of inverse mirroring at these points. In this
approach, the function is extended beyond the points of discontinuity by using its inverse mirror image and the
cubic-spline interpolation is applied to each extended function obtained by individual inverse mirroring. This
reduces the oscillations of the function known as the tail effect. Figure 3.2 shows an example of inverse mirroring at
discontinuous points (t  50,100, 200) . Figure 3.3 shows the interpolation results from the application of the inverse
mirroring and cubic spline.
3.4
Adaptive Cubic Spline
Inverse mirroring at discontinuous points can be also used to extend the sampled data at both ends of the table
where the behavior of the function is unknown. In this case, the two discontinuous points are the two ends of the
data set. Figure 3.4 shows an example of inverse mirroring at the data boundaries of an arbitrary function. The
inverse mirroring at the first data point may create negative values.
transforming the extended function into the Laplace domain.
16
These data points are rejected while
Figure 3.1 Interpolation of discontinuous pressure data using cubic spline.
Figure 3.2 Inverse Mirroring at discontinuous points (t  50,100, 200) .
17
Figure 3.3 Application of inverse mirroring at discontinuous points using cubic spline interpolation.
Field data may not be as smooth as was shown in Figure 3.3. The oscillatory nature of the field data may cause
the same effect as discontinuities on the cubic spline interpolation. The application of inverse mirroring for the
entire set of data in these cases would make the procedure impractical. Using inverse mirroring only at the two ends
of the data set as shown in Figure 3.4 may be a partial solution for this problem. Inverse mirroring at the first data
point is used to fill the gap from zero to first data point. However, It is used in the last data point to shift the tail
effect may occur during Laplace inversion to the right.
An alternative solution for data with discontinuity at any point between the first and last data points is to use an
adaptive cubic spline approach. In the adaptive cubic spline approach, piecewise linear approximations are
substituted for the function in intervals where the cubic spline approximation causes large oscillations. The
mathematical expression of the Laplace transformation of the piecewise linear approximations used in selected
segments are as follows:
The equation of a straight line over a segment can be written as:
y  mx  b
(3.14)
where m is the slope of the line passing through knots over a particular segment, and b represents the intercept of the
line.
Slope
m
m
for a particular line between knots can be found from the following relation
yi  yi 1
(3.15)
ti  ti 1
Similarly, the intercept of the line can also be found as follow
18
b  yi 1  mti 1  yi 1 
yi  yi 1
ti  ti 1
ti 1
(3.16)
The Laplace transformation of a piecewise-linear function over an interval from ti 1 to ti is then written as follow:
ti


L  y (t )  y ( s ) 
ti
ti 1
 ti
 m 
 s


y (t )e  st dt 
ti
ti 1
1 
e
s2 
 sti
 ti 1


(mt  b)e  st dt 
 s
ti
mte  st dt 
ti 1

1 
e
s2 
 sti 1
 be
 st
dt
ti 1
(3.17)
 b   sti  sti1 
  s e  e


where m and b are given in Eq. (3.15) and Eq.(3.16), respectively. A cutoff value is used to decide where the switch
between the cubic-spline and piecewise-linear approximations should occur. The cutoff value is defined as the
absolute value obtained from the difference between the cubic spline approximation and the linear approximation in
each individual segment in time domain. For practical purposes, the adaptive cubic spline approach combined with
inverse mirroring at data-set boundaries has proven to be an excellent combination for deconvolution purposes in
Laplace domain. Figure 3.5 and Figure 3.6 illustrate the application of the adaptive cubic spline with linear
approximation for both simulated and field examples.
3.5
The Iseger Algorithm
As previously noted, the Stehfest (1970) algorithm is a standard numerical Laplace inversion algorithm used in
petroleum engineering. The Stehfest algorithm, however, cannot handle discontinuous functions. As previously
discussed, in the application of the deconvolution method in transient flow problems, step-wise behavior of
production rate and other discontinuities make the Stehfest algorithm inapplicable. A newer algorithm presented by
Iseger (2006) removes the restriction on discontinuities. The Iseger algorithm is based on Poisson’s summation
formula in the form of Fourier series. In this algorithm, Poisson’s summation relates an infinite sum of Laplace
transform values to Z-transforms of the function’s values. The infinite sum is approximated by a finite sum based on
the Gaussian quadrature rule, and the time domain values of the function are computed by a Fourier Transform
algorithm (Al-Ajmi et al., 2008). The practical application of Iseger’s algorithm in transient-flow problems was
introduced by Al-Ajmi et al. (2008). In this work, the Iseger algorithm is used for deconvolution applications.
The appeal of the Iseger algorithm is in its ability to compute inverse Laplace transforms of functions with
discontinuities, singularities, and local non-smoothness, even if the points of discontinuity and singularity are not
known a-priori. The implementation of the algorithm may be slightly more involved than the Stehfest (1970)
algorithm, but it is comparable to the implementation of the other algorithms with similar capabilities.
A critical parameter used in the Iseger algorithm is  computed from the following relationship:

T
M
(3.18)
where T is the period for which the inversions are computed and
M
is the number of points at which the inverse
Laplace transforms are computed. To obtain accurate results, Al-Ajmi et al. (2008) found that 1  should be an
( )
integer éë 1 D ÎN ùû and the following condition must be satisfied:
19
Figure 3.4 Inverse mirroring at both ends.
Figure 3.5 Application of adaptive cubic spline using simulated data. Circles show where the adaptive cubic
spline is required.
20
Figure 3.6 Application of adaptive cubic spline using field data. Circles show where the adaptive cubic spline is
required.
M  10M ; that is, M  10T
(3.19)
To demonstrate the use of these conditions, we consider the inversion of the unit step function for a period of 20 hr
as shown in Figure 3.7. By Eq. (3.19), the minimum number of inversion points required for this case is M  200 .
If 220 points are used, the condition in Eq. (3.19) is satisfied and 1   11  N . Figure 3.7 shows how we obtain
an accurate inversion of the unit step function with 220 data points. If we use 103 data points, for example, the
condition in Eq. (3.19) is not satisfied and if we use 223 data points, the condition in Eq. (3.19) is satisfied but
1   N . For both cases, inversions are not sufficiently accurate. We also note from our experience, that if the
characteristics of the function change in short intervals (e.g., rate and pressure variations due to short drawdown and
buildup periods in a mini-DST test), a large
M
may be required for accurate inversions.
The Iseger algorithm uses M 2  nrp  M over sampling points to compute inversions at
M
data points. It was
recommended by Iseger (2006) to use nrp  8 especially for well-behaved functions. This particular choice of nrp
was dictated by the Fast Fourier Transform originally used by Iseger, which required M and nrp to be a power of 2.
For our applications, however, we did not obtain stable inversions with nrp  8 . We have found that using nrp  3
with Discrete Fourier Transforms generated more stable inversions for our applications. Figure 3.8 shows the effect
of nrp on numerical inversion of a typical pressure function encountered in fluid flow problems by using the Iseger
algorithm.
21
2
Unit Step Function (One Step Change)
220 data
103 data
1.5
233 data
1
0.5
0
-0.5
-1
0
5
10
15
20
Time, hr
Figure 3.7 Numerical Inversion of Unit Step Function by the Iseger algorithm; Effect of the number of inversion
points.
3.5.1
Verification Examples Using Iseger Algorithm
In this section, we present examples to verify the success of the Iseger algorithm in the inversion of functions
that are piecewise differentiable or discontinuous. We compare the results to the Stehfest algorithm to delineate the
differences from standard algorithms.
1.4
1.2
Function
1
0.8
0.6
0.4
nrp=3
nrp=7
0.2
0
0
5
10
15
20
Time, hr
Figure 3.8 Effect of nrp in the numerical inversion of a function by the Iseger algorithm.
22
3.5.2
Discontinuous and Piecewise Differentiable Functions
Figure 3.9 shows the inversion of a unit step function (Heaviside function) with the pole at t  10 hr using the
Iseger and Stehfest algorithms. For the Iseger algorithm, inversion was performed with M  220 and nrp  3 . Three
values of the parameter, N  6, 8, and 12, were used for the Stehfest inversions ( N controls the number of
functional evaluations in the Stehfest algorithm and, theoretically, higher N values yield better inversions). Figure
3.9 shows that the Iseger algorithm recovers the true step-wise character whereas the Stehfest algorithm smears the
function around the point of discontinuity ( t  10 hr).
In Figure 3.10, we consider a function with multiple step changes. This function corresponds to the rate sequence
used in a synthetic deconvolution example shown below. For the results in Figure 3.10, we used M  512 and
nrp  3 in the Iseger algorithm and N  6, 8, and 12 in the Stehfest algorithm. The smearing of the function at the
points of discontinuity with the Stehfest algorithm completely destroys the local and global characteristics of the
original function. The success of the Iseger algorithm in this example is remarkable.
To show that special algorithms, such as the Iseger algorithm, may be required even when the function is
continuous but only piecewise differentiable, we present the results in Figure 3.11. To generate the function in this
example, we used the function in Figure 3.10 as the rate sequence and generated the corresponding pressure
changes. The Laplace transformation of this tabulated pressure-change function was obtained by the RoumboutsosStewart (1988) algorithm (the Onur-Reynolds, 1998, algorithm generated the same results) and inverted back by the
Iseger and Stehfest algorithms. As shown in Figure 3.11, while the Stehfest algorithm returns a smooth,
continuously differentiable function, the Iseger algorithm accomplishes an accurate inversion that preserves all the
characteristics of the piecewise differentiable function.
2.0
Unit Step Function
1.5
Original Function
Iseger
1.0
Stehfest, N=6
Stehfest, N=8
Stehfest, N=12
0.5
0.0
-0.5
0
5
10
15
20
Time, hr
Figure 3.9 Numerical Inversion of a unit step function; comparison of Iseger and Stehfest algorithms.
23
Original Rate
Iseger
Stehfest, N=6
Stehfest, N=8
Stehfest, N=12
q, rbbl/d
300
200
100
0
0
4
8
12
16
Time, hr
Figure 3.10 Numerical Inversion of a function with multiple step changes; comparison of Iseger and Stehfest
algorithms.
300
Original Function
Iseger
Stehfest
p, psi
200
100
0
0
4
8
12
16
Time, hr
Figure 3.11 Numerical Inversion of a piecewise differentiable function; comparison of Iseger and Stehfest
Algorithms.
24
3.5.3
Wellbore Pressure Solution for a Combined Drawdown and Buildup
This is a classic example where the Stehfest algorithm does not yield good results. For the particular example
here, we assume that the well is produced at a constant-rate until time t p  8 hr and then shut in. This rate sequence
introduces a discontinuity into the problem. Pertinent data for this example is given in Table 3.1.
Table 3.1 Data for drawdown followed by build up
q
B
kh
C
rw
rbbl/d
rbbl/stb
md.ft
bbl/psi
ft
1
1000
0.01
736
η
Sf
0.3
md-psi/cp
0
29300
Correa and Ramey (1986) provided the following solution in the Laplace domain for the wellbore pressure
change for this problem:
pwf





 st

K0 s  rw  S f s  rw K1 s  rw 1  e p
141.2qB

141.2C  
kh

s  s  rw K1 s  rw 
s K 0 s  rw  S f s  rw K1

kh








s  rw  

(3.20)

The term 1  e  stp in Eq. (3.20) introduces the effect of the step-rate-change at t p  8 hr into the solution. In Eq.
(3.20), S f represents skin factor and s denotes Laplace transform parameter.
Figure 3.12 shows the results of the numerical inversion of the solution in Eq. (3.20) for the data in Table 1. As
expected, the Stehfest algorithm cannot accurately compute the behavior of the solution at the point of rate
discontinuity ( t p  8 hr). Chen and Raghavan (1996) suggested a practical approach to handle this problem by using
the Stehfest algorithm. Their approach is based on the superposition solution of the pressure buildup problem and
requires the inversion of two drawdown solutions with a time shift. This approach, however, is not convenient for all
applications that involve step-rate changes. As shown in Figure 3.12, inversion by the Iseger algorithm yields
accurate results and does not require any special treatment of the solution (we used M  1024 and nrp  3 for this
example).
3.6
Adaptive Nonparametric Kernel Regression
As previously mentioned, deconvolution is an ill-conditioned problem. Small errors in the input data (in this
case, pressure and production rate) create large variations in the deconvolved constant rate pressure response
(impulse function). The oscillations in the constant-rate response retrieved from the deconvolution process are
magnified in the calculation of its time derivative and cause difficulty in the pressure-transient interpretation. To
alleviate this problem, several approaches such as smoothing, filtering, and regularization [Ilk et al., (2005) & von.
Schroeter et al. (2004)] have been proposed.
In this work, to smooth the deconvolved data, we apply the adaptive Nadaraya-Watson (1964) kernel regression
25
900
Iseger
Stehfest
800
700
p, psi
600
500
400
300
200
100
0
0
2
4
6
8
10
12
14
16
18
Time, hr
Figure 3.12 Numerical Inversion of wellbore pressure change for a combined drawdown and buildup sequence.
using a kernel function. The kernel regression is a non-parametric technique with the objective of finding a
relationship between a dependent and independent variables. For the purposes of this work, several kernels (e.g.,
Gaussian, Epanechnikov, Tricube, and Triweight) were tested. We have found that the Gaussian and Epanechnikov
kernels yielded better results than the others. In general, the Gaussian kernel is unbounded and has a smooth
behavior compared to the other kernels; however, if the number of the data points is large (in the order of
thousands), the Epanechnikov kernel becomes faster than the Gaussian kernel.
The Nadaraya-Watson (1964) kernel estimator of the regression function is obtained as:
n
Y K (
x  xi
i
g ( x | h) 
i 1
n
 K(
x  xi
i 1
)
h
h
  x  
(3.21)
)
where K is a kernel function with a bandwidth h . The bandwidth, h , of the Nadaraya-Watson kernel estimator
controls the smoothing level of the estimation (the larger the h , the smoother the result). The bandwidth, h , plays a
very important role in the performance of the kernel estimator. Various methods of choosing bandwidth h are
available. In this work we used the common method of cross-validation. The cross-validation method is often
preferred because it is easily computable and applicable for any regression model; but it is time consuming.
The following expression is the general Gaussian function;
 ( x | , ) 
1
 2
2

 (x  ) 


2
2





exp 
(3.22)
26
In Eq.(3.22) and more specifically in Kernel regression,  is known as the centering parameter and  as the
dispersion smoothing parameter (  h) .
The general form of Epanechnikov kernel function is also presented as following;
 ( x) 
3
4
1  x   x  1
2
(3.23)
The deconvolved pressure response obtained from Eq. (3.4) usually shows a smooth behavior at early times and
some oscillations are observed at late times (the tail effect). If a constant smoothing factor is used in the kernel
regression for the entire data set, while sufficient smoothing is achieved at late times, the early time data is over
smoothed. The over smoothing at an early time may mask some important features, such as dual-porosity system
behavior in pressure transient analysis, which misleads the data interpretation. (In some cases, an appropriate
constant smoothing factor may not be found to smooth the entire data.) To overcome the problem, an adaptive
kernel regression approach is recommended as follows:
x  xi
n
 Y K ( h( x ) )
i
g ( x | h( x)) 
i 1
n
i

K(
i 1
x  xi
h( xi )
  x  
(3.24)
)
Abramson (1982) proposed a method, which uses a value of h( xi ) proportional to


f ( xi )1/ 2 where f ( xi ) is
a prior kernel
density estimator with a fixed bandwidth. Silverman (1986) gave a three-step algorithm for the Abramson-type
estimator. In the first step, a prior kernel density estimator

f ( xi ) with
a fixed bandwidth is obtained. In the second
step, the local bandwidth factor i is defined as


i   f ( xi ) / J 



(3.25)


1
where J  1  f ( xi ) and f ( x) 
n
n
 x  xi 
.
h 
 K 
nh
i 1
i 1
In Eq. (3.25), 0    1 is the sensitivity parameter. In the last step, the modified adaptive Nadaraya-Watson kernel
estimator with the modified local bandwidth factor i can be written as
n
g ( x | h) 

i 1
n
Yi
i

i 1
K(
1
i
K(
x  xi
hi
x  xi
hi
)
  x  
(3.26)
)
In Eq. (3.26), the adaptive bandwidth is taken as h( xi )  hi and h is obtained from the cross-validation method. The
adaptive kernel estimation is equivalent to the kernel estimation with fixed bandwidth when   0 . When   1 , then
the adaptive kernel estimation is equivalent to the nearest neighbor estimation. Abramson and Silverman
emphasized that taking   0.5 leads to good results.
As an example, we have applied smoothing with fixed and adaptive bandwidth Nadaraya-Watson kernel
27
estimator using Gaussian kernel regression to a noisy set of data. We applied cross-validation to obtain the fixed
value of h in Eq. (3.21).Figure 3.13 shows the optimum fixed bandwidth obtained from cross-validation and Figure
3.14 illustrates the local bandwidth factor i . The original and the smoothed data are shown in Figure 3.15. Both
fixed and adaptive bandwidth approaches successfully smooth the data. As expected, the adaptive bandwidth
approach provides smoother results than the fixed bandwidth approach.
3.7
Deconvolution of Pressure Responses for a Sequence of Step-Rate Changes
To demonstrate the use of the adaptive cubic spline in deconvolution applications, we consider the step-rate
sequence and corresponding bottomhole pressure in Figure 3.16 and Figure 3.17, respectively. The pertinent data for
this example are given in Table 2. For this variable-rate (step-rate changes) example, we used the Laplace domain
deconvolution given in Eq. (3.4) to generate the constant rate responses (with wellbore storage) shown in Figure
3.18. The Laplace transformation of the tabulated pressure responses, pw (t ) , were generated by the adaptive cubic
spline algorithm combined with an analytical solution for tabulated production rate. The Iseger numerical Laplace
inversion algorithm is utilized in this example.
Figure 3.13 Optimum fixed bandwidth obtained from cross-validation method.
28
Figure 3.14 Local bandwidth factor, i , applied in adaptive Nadaraya-Watson kernel estimator. Higher lambda
means higher dispersion at particular time.
Figure 3.15 Nadaraya-Watson kernel estimator using Gaussian kernel regression using fixed and adaptive
bandwidth (   0.5 ) showing improvement over ordinary kernel regression.
29
Table 3.2 Data for deconvolution example
rw
pi
μ

ct
h
kh
ft
psi
cp
fraction
psi-1
ft
md.ft
5000
1
0.1
3.00E-06
0.3
30
1000
Figure 3.16 Step-rate sequence for data in Table 2.
30
Sf
0.0
C
B
bbl/psi
rbbl/stb
0.01
1
Figure 3.17 Pressure changes corresponding to step-rate sequence shown in Figure 3.16.
Figure 3.18 Deconvolution of pressure responses in Figure 3.17 for the step-rate sequence in Figure 3.16.
31
3.8
Field Examples
Three field examples of the adaptive cubic-spline deconvolution will be presented herein. The first two are the
classical examples of Meunier et al. (1985) and Fetkovich and Vienot (1984) which deal with variable sandface flow
rate due to wellbore storage effect. The third example is a variable-rate production case in an unconventional
reservoir.
3.9
Sandface Rate Deconvolution – Muenuier et al., (1985) Example
The first standard example is an 8-hour buildup test with after flow effects considered by Muenuier et al. (1985).
In this example, the underlying dual-porosity reservoir behavior is masked by wellbore storage and the objective of
deconvolution is to reveal the underlying reservoir behavior. Figure 3.19 shows the original data and the
deconvolved pressure response. Adaptive cubic spline, inverse mirroring, adaptive kernel regression, and Iseger
numerical inversion algorithm are employed in this example. The results from adaptive cubic spline deconvolution
and adaptive kernel regression display the valley in derivative responses, which is characteristic of dual-porosity
systems (the depth of the valley is used to estimate the storativity ratio of the dual-porosity medium). The valley in
the derivative responses was masked in the original data by wellbore storage. The results of both the fixed and the
adaptive bandwidth utilized in kernel regression are shown in Figure 3.19. The fixed bandwidth value used in kernel
regression was obtained from cross-validation to avoid over smoothing. A discussion on how over smoothing can
affect the storativity estimation in this example is given by Al-Ajmi et al. (2008).
ΔP and dΔP/dlnt,Psi
100
10
1
0.1
0.01
0.001
Original Pressure
Original Pressure Derivative
Deconvolved Pressure- No Smoothing
Deconvolved Pressure Derivative-No Smoothing
Deconvolved Pressure-Constant Bandwidth Smoothing
Deconvolved Pressure Derivative-Constant Bandwidth Smoothing
Deconvolved Pressure- Adaptive Bandwidth Smoothing
Deconvolved Pressure Derivative-Adaptive Bandwidth Smoothing
0.01
0.1
1
10
Time, hr
Figure 3.19 Sandface-rate deconvolution to remove wellbore storage effect; Muenuier et al. (1985) example.
32
3.10
Sandface Rate Deconvolution – Fetkovich and Vienot (1984) Example
This example is a buildup test with after flow effect. The data for this example are from Phillips’ Oil Well Nr. 1,
which is hydraulically fractured. Figure 3.20 shows the results of adaptive cubic spline deconvolution using kernel
regression. In this example, the Epanechnikov kernel with constant bandwidth is applied. The Iseger algorithm is
also used for the numerical inversion from Laplace space to time domain. Deconvolution results are satisfactory
without smoothing in this example. Smoothing provides some improvement on derivative data but the early-time
data seems to be distorted by the constant-bandwidth smoothing.
3.11
Deconvolution of Variable-Rate Data-Shale-Gas Well
The final example is the application of adaptive cubic spline deconvolution in an unconventional, shale-gas
reservoir. Figure 3.21 shows the variable production rate and corresponding pseudo-pressures for this example. In
Figure 3.22, the deconvolved pressure and corresponding pressure derivative is shown. To remove the oscillation
and the effect of the noise, kernel regression utilizing Epanechnikov kernel with constant bandwidth obtained from
cross-validation is applied.
For comparison purposes, the time domain deconvolution method (von. Schroeter et al., 2004) using commercial
software (Ecrin 4.20) is also provided in Figure 3.22. We have found that the deconvolution results for this example
were strongly dependent on the input parameters of the time-domain deconvolution. In addition, to obtain a similar
trend from adaptive cubic spline deconvolution, over smoothing was required.
Figure 3.20 Sandface- Numerical inversion of pressure drop from tabulated data; Fetkovich and Vienot (1984)
example.
33
Figure 3.21 Example of flowing gas rate and corresponding pseudo-pressure.
Figure 3.22 Numerical inversion of pressure drop from tabulated data: Shale gas application for data shown in
Figure 3.21.
34
CHAPTER 4
LAPLACE TRANSFORMATION SOLUTION TO THE NONLINEAR DIFFUSIVITY EQUATION
Non-linear fluid-flow problems arise in porous media due to various conditions including pressure-dependent
rock and fluid properties, such as the flow of real gases (pressure-dependent viscosity and compressibility of real
gases) and stress-dependent porosity and permeability of the porous medium. In some of these cases, depending on
the application ranges, the non-linearity may be reduced to a weak form and may be ignored while, in others,
alternate approaches are required. In this chapter, flow of real gases in porous media is considered for the
demonstration of the solution procedure, as it is probably the best-known non-linear fluid flow problem in porous
media. The remaining non-linearity of the real-gas diffusion equation after pseudo-pressure transformation is
relatively weak; therefore, the results of this chapter should not be evaluated in terms of the significance of the
magnitude of the non-linear contributions demonstrated by the solutions. Here, the objective is to demonstrate the
procedure on a well-known problem and leave the application of the procedure to other non-linear problems as a
choice for the reader.
The diffusivity equation, which describes fluid flow in porous media, is used in the formulation of the problem
in this Chapter. It is well known that the diffusivity equation is derived from three fundamental physical principles:
a) conservation of mass, b) equation of motion, and c) equation of state (EOS).
The following assumptions are considered in the derivation of the diffusivity equation in this Chapter:
4.1

Homogeneous and isotropic 2D reservoir of infinite extent in

Uniform thickness, permeability, and porosity

Darcy flow

Permeability independent of pressure

Isothermal gas flow

No gravity

No compositional changes
x and y directions
Mathematical Model
The continuity equation in porous media is given by;
.   v  

     q( x, y, z, t )
t
(4.1)
where  is fluid density, v is velocity,  is porosity of porous media, and q is fluid volume per unit volume. The
equation of flux (Darcy’s Law) implies that the velocity is proportional to the negative of the gradient of the
potential
v
k


(4.2)
35
where
p


pb
dp '
(4.3)

p
pb is the pressure at a datum (reference pressure). In Eq.(4.3) that the potential,  

pb
dp '

, known as work done
by the flow. By combining Eq.(4.2) with Eq.(4.3), the flux law can be written as follows;
k
p
 ( p)
v
(4.4)
Substituting Eq.(4.4) into Eq.(4.1) yields

 
k
.   
p        q( x, y, z, t )
 ( p)  t

(4.5)
The equation of states for real gases is given by

pM
Z ( p) RT
(4.6)
where, M is the molecular weight of the gas, Z is the gas compressibility factor, R is the ideal gas constant, and T
is absolute temperature.
Let’s expand the second term in the left hand side of Eq.(4.5);



      
t
t
t
(4.7)
 1  1  
  


  t  t 
(4.8)
 1  p 1  p 
  


  p t  p t 
(4.9)
The compressibility of the gas and the compressibility of the rock as a function of pressure are given by;
cg 
1 d
 dp
(4.10)
cr 
1 d
 dp
(4.11)
Substituting Eq.(4.10) and Eq.(4.11) into Eq.(4.9) yields

p
     cg  cr
t
t


(4.12)
Let’s define
cg  cr  ct
(4.13)
Then Eq.(4.12) can be written in the form
36

p
     ct
t
t
(4.14)
By substituting Eqs.(4.14) and (4.6) into Eq.(4.5)
.(
k
pM
 ( p) Z ( p) RT
p)   ct
p
pM
pM

q( x, y, z, t )
t Z ( p) RT Z ( p) RT
(4.15)
Cancelling out the molecular weight from both sides, yields:
.(
 ( p) ct p
p
p
p
p) 

q( x, y, z, t )
Z ( p)  ( p)
k
t Z ( p)  ( p) Z ( p)k
(4.16)
The initial and boundary conditions are defined as
p( x, y, z, t  0)  pi
(4.17)
lim p( x, y, z, t )  pi
(4.18)
x, y, z  
Al-Hussainy and Ramey (1966) proposed the concept of real gas pseudo-pressure to obtain a linear equation for gas
flow in porous media. The real gas pseudo-pressure is given by
p
m( p)  2
p'
 Z ( p ) ( p )dp
'
'
(4.19)
'
pb
Defining the pseudo-pressure difference by
pi
m( p)  m( pi )  m( p)  2
p'
 Z ( p ) ( p )
'
pb
'
p
'
dp -2
p'
 Z ( p ) ( p )
'
'
pi
'
dp =2
pb
p'
 Z ( p )( p )dp
'
'
(4.20)
p
and substituting Eq. (4.20) into Eq. (4.16) yields
2 (m( p)) 
 ( p) ct m
k
t

2p
q( x, y, z, t )
Z ( p )k
(4.21)
The corresponding initial and boundary conditions in terms of pseudo-pressure are written as
m( p)( x, y, z, t  0)  0
lim
(4.22)
m(p)( x, y, z, t )  0
(4.23)
x, y, z  
If we define a diffusivity term by
 ( p) 
k
 ( p) ct
(4.24)
then Eq.(4.21) can be written in the compact form as
2 (m( p)) 
1 m
2p

q( x, y, z, t )
 ( p) t
Z ( p) k
(4.25)
37
Equation (4.25) is a non-linear and non-homogeneous partial differential equation. The source term,
2p
q( x, y, z, t ) , in the right hand side of the Eq.(4.25) should be treated separately. If we use the mass balance
Z ( p)k
(4.26)
 q ( x, y, z, t )  sc qsc ( x, y, z, t )
and replace  by  
PM
, Eq.(4.26) can be re-formulated in the form given by
Z ( p) RT
P M
PM
q( x, y, z, t )  sc
qsc ( x, y, z, t )
Z ( P) RT
Z sc RTsc
(4.27)
Substituting Eq.(4.27) into Eq.(4.25) yields
2 (m( P)) 
m
1
 ( p) t

2 Psc qsc ( x, y, z )T
Tsc k
(4.28)
Let’s define the following dimensionless parameters:
tD 
kt
 (  ct )i L
xD 
2
x
L
 i
, yD 
t
(4.29)
L2
y
L
, zD 
z
L
(4.30)
where L is defined as a characteristic length. Equation (4.28) can be expressed in terms of the dimensionless
parameters given by Eqs. (4.29) and (4.30) as follows
~

m 2 psc q sc ( x, y, z ) T
 (m( P))  2 i

Tsc k
L  ( P) tD
(4.31)
2
If we define the dimensionless pseudo-pressure, mD , by
mD 
khTsc
m( p)
2qref pscT
(4.32)
and also define the term
k

 (  ct )i
 ct
RD  i =
=
k
 ( p)
(  ct )i
 ct
(4.33)
Equation (4.31) can be written as
~
mD h q sc ( x, y, z ) L2
 mD  RD

tD
qref
(4.34)
2
Finally, the initial and boundary conditions are expressed in dimensionless form as follows:
38
mD ( xD , yD , zD , tD  0)  0
(4.35)
lim mD ( xD , yD , zD , tD )  0
(4.36)
xD , yD , zD  
If we define
~
h q sc ( x, y, z ) L2
g ( xD , y D , z D ) 
qref
(4.37)
and
w
 ct  (  ct )i
(  ct )i
(4.38)
combining Eqs.(4.33) and (4.38) yields
RD  1  w
(4.39)
Substituting Eq.(4.39) into Eq.(4.34) yields the following compact form:
2 mD  (1  w)
mD
 g ( xD , yD , zD )
tD
(4.40)
or
2 mD 
mD
m
 w D  g ( xD , yD , zD )
tD
tD
The nonlinear term w
(4.41)
mD
with g ( xD , yD , zD ) are the source terms in Eq.(4.41).
t D
In Eq. (4.39) where RD  1 it means that there are no changes in the viscosity-compressibility product which is
equivalent to w  0 . In the above equation, the factor w(mD ) is a function of dimensionless pseudo-pressure and the
term w
mD
is implicit and it can be computed by an iterative method or by an asymptotic expansion of mD
t D
(perturbation method) presented by Barreto (2011), and Peres et al. (1990). In the iterative approach, the iteration
starts from the solution of a slightly compressible fluid that is also the solution for the case where the diffusivity
term  is constant at its initial value. The results presented by Barreto (2011) demonstrated that in an asymptotic
expansion (perturbation) method, the solution truncated at the first order term showed excellent results. In this work,
we follow Barreto’s (2011) findings and only evaluate the first order term of the perturbation.
4.2
Asymptotic Expansion (Perturbation)
By definition, the theory of perturbation is the study of the effects of small disturbances. Based on the
fundamental theorem of perturbation, the dimensionless solution given in Eq.(4.41) can be written in the form of
39
mD  m(0)   m(1)   2 m(2)  ... 
D
D
D


k
mD( k )
(4.42)
k 0
which represents a power series in  where the coefficients mD ( k ) represents the order of the solution. The solution is
obtained by taking the limit as   0 . If  is sufficiently small and if the coefficients m(0) , m(1) , m(2) ,.... are
D
D
D
independent of  , then
m(0)  m(1)  m(2)  ...  m( N )
D
D
D
(4.43)
D
Application of perturbation theory to Eq.(4.40) disturbed by the introduction of variable  multiplied by w
factor which is responsible for non-linearity yields
2 mD  (1   w)
mD
 g ( xD , yD , zD )
tD
(4.44)
where   0 yields the partial differential equation becomes linear.
2 mD  2 m(0)   2 m(1)   2 2 m(2)  O( 3 ) 
D
D
D

 
2 k
mD( k )
(4.45)
k 0

 m(0)

m(1)
m(2)
mD
m ( k )
2
D
D
D

(1   w)
 (1   w)


 O( 3 )   (1   w)  k D
 tD

tD
tD
tD
tD
k 0



(4.46)
The left hand side of Eq.(4.40) is obtained by subtracting Eq.(4.45) from Eq.(4.46)
2 mD  (1   w)


mD
mD( k )

2 k mD( k )  (1   w)  k
 g ( xD , yD , zD )
tD k 0
tD
k 0


(4.47)
Expanding Eq.(4.47) and rearranging
 m(0)

m(1) 2 m(2)
D
D

 m   m    m  ...  (1   w)
 D 
 ...   g ( xD , yD , zD )
D
D
D
 tD

tD
tD


2 (0)
2 (1)
2 2 (2)
(4.48)
Furthermore, Eq.(4.47) can be rearranged and collected in the form

m(0)
  2 m(0)  D  g ( xD , yD , z D )
D

tD


m(2)
m(1)
2  2 (2)
D
  m 
w D
D

tD
t D



m(1)
m(0)
     2 m(1)  D  w D
D


t D
t D



  ....  0






(4.49)
subject to the following initial and boundary conditions


k
mD( k ) ( xD , yD , zD , tD  0)  0
(4.50)
k 0
40


k
mD ( k ) ( xD , yD , zD , t D )  0
(4.51)
k 0
xD , y D , z D  
Barreto (2011) confirmed that if the perturbation theory applies, the solution truncated at the first order of
perturbation is close to the slightly compressible solution. In other words, the m(1) , m(2) ,.... solutions are small and
D
D
the mD solution is close to m(0) . Then, the changes in compressibility-viscosity product or w can be approximately
D
estimated from m(0) while solving for m(1) . Accordingly, w in the calculation of m(2) can be calculated from the
D
D
D
solution obtained from the previous step m(1) .
D
Referring to the perturbation theory noted previously (Eqs.(4.42) and (4.43)), the individual terms in Eq.(4.49)
can be treated separately as follows:
The order 0 term is:
2 m(0) 
D
m(0)
D
tD
 g ( xD , y D , z D )  0
(4.52)
subject to initial and boundary conditions;
m(0)
( xD , yD , zD , tD  0)
D
(4.53)
m(0)
( xD , yD , zD , t D )  0
D
(4.54)
lim xD , yD , zD  
The order 1 term is:
2 m(1)

D
m(1)
D
tD
 w(0)
m(0)
D
tD
0
(4.55)
subject to initial and boundary conditions;
m(1)
( xD , yD , zD , tD  0)
D
(4.56)
m(1)
( xD , yD , zD , t D )  0
D
(4.57)
lim xD , yD , zD  
The order 2 term is:
2 m(2)

D
m(2)
D
tD
 w(1)
m(1)
D
tD
0
(4.58)
subject to initial and boundary conditions;
41
m(2)
( xD , yD , zD , tD  0)
D
(4.59)
m(2)
( xD , yD , zD , t D )  0
D
(4.60)
lim xD , yD , zD  
The order
k term can also be written as
2 m(Dk ) 
m(Dk )
tD
 w( k 1)
m(Dk 1)
tD
0
(4.61)
subject to initial and boundary conditions;
m(Dk ) ( xD , yD , zD , tD  0)
(4.62)
m(Dk ) ( xD , yD , zD , t D )  0
(4.63)
lim xD , yD , zD  
All of the preceding linear partial differential equations can be solved by Green’s Functions.
The source terms of partial differential equation of order zero to order k are as follows:
f (0)  g ( xD , yD , zD )
f
(k )
( k 1)
w
m( k 1)
D
tD
(4.64)
k  1, 2,...
(4.65)
where the factor w is defined as
k 1
m
w( k 1)  w
D
(i )
k  1, 2,...
k  1, 2,...
(4.66)
i 0
The general solution of the non-homogeneous partial differential equation with homogeneous boundary conditions
can be shown in the form of the Green’s function. The Green’s function
GD ( xD , xD' , yD , yD' , zD , zD' , tD , tD' )
(4.67)
is introduced as a solution due to a concentrated source at xD  xD ', yD  yD ', zD  zD ' acting instantaneously only
at tD  tD' subject to
G
(k )
D

G (0)
D
tD
  ( xD  xD ' ) ( yD  yD ' ) ( zD  zD ' ) (t D  t D ' )
(4.68)
with the initial and boundary conditions given by
G( k ) ( xD , xD' , yD , yD' , zD , zD' , tD  0, tD' )  0
(4.69)
D
and
42
G ( k ) ( xD , xD' , yD , yD ' , zD , zD ' , tD , tD ' )  0
(4.70)
D
lim xD , yD , zD  
where  ( xD  xD ' ) is the delta Dirac function of the appropriate dimension.
The solution of non-homogeneous partial differential equations of order zero to k in terms of Green’s function
can be then written as
tD
m
(k )
D
 f

(k )
( xD' , yD' , zD' , tD' )G ( k ) ( xD , xD' , yD , yD' , z D , z D' , t D , t D' )dxD' dyD' dz D' dt D'
D
(4.71)
0 D
It is noted that for a given reservoir geometry and its boundary conditions, the Green’s function is the same.
Therefore, in partial differential equations of order zero to k , the same Green’s function can be used.
G( k ) ( xD , xD' , yD , yD' , zD , zD' , tD , tD' )  GD ( xD , xD' , yD , yD' , zD , zD' , t D , t D' )
D
(4.72)
Substituting Eq.(4.72) into Eq.(4.71), m( k ) for orders zero to k can be rewritten as follows:
D
tD
m(0) 
D
  g(x
'
D
, yD ' , zD ' , tD ' )GD ( xD , xD ' , yD , yD ' , zD , zD ' , tD , tD ' )dxD ' dyD ' dzD ' dt D '
(4.73)
0 D
tD
m
(1)
D
 w

(0)
0 D
tD
m(2) 
D

w(1)
0 D
m(0)
D
t 'D
m(1)
D
t 'D
GD ( xD , xD' , yD , yD ' , zD , z D ' , t D , t D ' )dxD ' dyD ' dz D ' dtD '
(4.74)
GD ( xD , xD ' , yD , yD ' , z D , z D ' , t D , t D ' )dxD ' dyD ' dz D ' dt D '
(4.75)
and
tD
m
(k )
D

 w
( k 1)
0 D
m( k 1)
D
t 'D
GD ( xD , xD ' , yD , yD ' , z D , z D ' , t D , t D ' )dxD 'dyD 'dz D 'dt D '
k  1, 2,...
(4.76)
To perform the double integration numerically, Barreto (2011) used a numerical integration package called CUBA
(Hahn-2005a, 2005b). The generated solution was obtained directly in the time domain. In this research the Laplace
transformation of the solution is presented. The Laplace transformation of the solution in Eq. (4.71) is
 (k )
mD
t D



 L
f ( k ) ( xD ' , yD ' , z D ' , t D ' )GD ( xD , xD ' , yD , yD ' , z D , z D ' , t D , t D ' )dxD 'dy D 'dz D 'dt D ' 
 0 D


t D



= L  f ( k ) ( xD ' , yD ' , z D ' , t D ' )GD ( xD , xD ' , yD , yD ' , z D , z D ' , t D , t D ' )dxD ' dyD ' dz D 'dt D ' 
0
D 

 
where
43
(4.77)
f (0)  g ( xD , yD , z D )
f
(k )
( k 1)
( xD , y D , z D )  w
'
'
'
m( k 1)
(4.78)
k  1, 2,...
D
t 'D
The integration over time in Eq. (4.77) is a convolution integral. Utilizing the Laplace transformation property for
the convolution of two functions given by Eq.(3.1), Eq.(4.77) can be written as a product of the functions in the
Laplace domain; that is,
 (k )
 (k )

mD = f

( xD' , yD' , zD' ) G D ( xD , xD' , yD , yD' , z D , z D' )dxD' dyD' dz D'
(4.79)
D
Note, however, that the use of the Laplace transformation of the solution with sampled data requires a proper
approximating function for the first term presented in Eq.(4.79). As previously explained, this was the motivation for
developing the cubic splines approximation for the Laplace transformation of sampled data under current research.
Our experience has shown that the cubic spline should be used when the Iseger (2006) numerical Laplace
inversion algorithm is utilized. The Stehfest (1970) numerical Laplace inversion algorithm did not yield accurate
results when the cubic spline interpolation was used for the approximation function demonstrated in Eq. (4.79).
Since the solutions provided in this chapter are in dimensionless form and the Iseger algorithm has a linear sampling
limitation, it is extremely time consuming to solve for large values of time compared with the Stehfest algorithm. To
be able to use the Stehfest algorithm for the developed solutions in the Laplace domain, the approximation function
introduced by Onur and Reynolds (1998) or Roumboutsos-Stewart (1988) provides a better option than the cubic
spline.
The Onur-Reynolds (1998) algorithm for a sampled function in the Laplace domain is
f s 
0
s 0 1
N 1


i 1
   0  1    0  1, st1   

fi e sti  e sti1
s2

f1e st1 f N e st N

s
s
(4.80)

    1, st N  
s  1 
Onur and Reynolds (1998) suggest that the constants,  0 and  0 , be determined by fitting a least-square straight
 
line of the form ln f t    0 ln t  ln 0 through the sampled data in the interval t1  t  t1 where ln t1 t1  Ll with



0.1  Ll  1 . The constants,  and  , should be determined by a similar procedure for the interval t N  t  t N that



includes a minimum of five data points and satisfies ln t N t N  Lr with 0.1  Lr  1 .
The Roumboutsos-Stewart (1988) algorithm for the Laplace transform of a sampled function is given by
f s  


f 0 f N stN
f  1  e  st1

e
 0 2

s
s
s
N 1

i 1

f i e  sti  e  sti 1
s
2
 f  e
s
N
2
stN
(4.81)
f i are the values of the function at points t i and f i is the linear slope of the data in the interval ti  t  ti1 .
Ozkan and Raghavan (1990) provided the Green’s function for different reservoir geometries and different
boundary conditions of fluid flow problems in porous media. The numerical integration in Laplace domain is
required to invert the solution back to the time domain using the Stehfest algorithm.
44
4.3
Variable Gas Rate Deconvolution
The research study presented by this work was focused on the solution of non-linear gas diffusivity equation
using perturbation method. Application of Duhamel’s principle requires a linear system but the diffusivity equation
of compressible fluids (gas) is non-linear. Al-Hussainy and Ramey (1996) introduced the concept of pseudopressure to linearize the gas diffusivity equation, however, as previously explained, it does not completely remove
the pressure dependency of the gas properties and the convolution/deconvolution of gas-well responses remains to
be non-rigorous. Using the pseudo-pressure approach, Levitan and Wilson (2010) applied time domain
deconvolution on gas reservoir with significant pressure depletion where pressure dependency of gas properties
becomes important. In the approach presented in this work, the pseudo-pressure concept is used to partially linearize
the problem. The remaining non-linearity is handled by the perturbation method described above and
convolution/deconvolution of gas well responses is applied more rigorously.
4.4
Validation
Two examples were used to validate the ideas and mathematical formulations presented above. The first example
is the solution for wells in an infinite acting reservoir with impermeable boundaries at the top and bottom. The
second example is the solution for wells in closed cylindrical reservoirs with impermeable boundaries at top,
bottom, and rD  reD .
4.4.1
Solution for Wells in Infinite Slab Reservoirs

q
The solution for an instantaneous point source of strength
acting at t   in an infinite, homogeneous, and
 ct
anisotropic reservoir is presented by Nisle (1958) is as follows:

p( M , M ', t   ) 
q
8 ct  x y z  (t   )
3/2


( M  M ')2 /  

exp 


4(t   )


(4.82)
In Eq. (4.82) the location of the observation point and the source are indicated as M and M ' , respectively. For a
three-dimensional Cartesian coordinate system,
(M  M ')2 /   ( x  x ')2 /  x  ( y  y ')2 /  y  ( z  z ') 2 /  z
(4.83)
where  x ,  y ,  z describe the diffusivity constants in x , y , and z directions defined as follows:
 
k
 ct
where   x, y, or z
(4.84)

q ( )
The pressure drop due to a continuous point source of strength
over time interval 0    t is
 ct
45
1
p(M , M ', t ) 
 ct
t


q( ) S ( M , M ', t   )d
(4.85)
0
In Eq. (4.85) S (M , M ', t   ) is the instantaneous point source in an infinite reservoir corresponding to unit strength


 q( )  1 and is given by,
  ct



S ( M , M ', t   ) 
1
8  x y z  (t   )
3/2


( M  M ')2 /  
exp  


4(t   )


(4.86)
Similarly, the pressure drop for instantaneous sources distributed over a line or any other source geometry in an
infinite acting reservoir can be obtained as;
p( M , M ', t ) 
1
 ct


q( M w ') S ( M , M w ', t   )dM w '
(4.87)
w

where q( M w ') indicates the instantaneous withdrawal volume of fluids per unit dimension of the source geometry
and M w is the point on the source  w . According to Eq. (4.87), the pressure drop for an infinite line source located
at x ' , y ' , and
  z   may
p( M , M ', t   ) 

be obtained as follows (Ozkan 2003):
1
8 ct  x y z  (t   )
3/2
 ( z  z ') 2 /  z
 q( x, y, z ) exp  

4(t   )




 ( x  x ')2 /  x  ( y  y ') 2 /  y
exp  

4(t   )


dz '





(4.88)

q( )
Assuming uniform flux along the line source with unit strength
 1 , then instantaneous point source for an
 ct
infinite line source in an infinite reservoir can be written as follows:
S ( x, x ', y, y ', t   ) 
 ( x  x ')2 /  x  ( y  y ')2 /  y
exp  

4(t   )
4  x y (t   )

1




(4.89)
In Cylindrical coordinate system (r , , z) , we define,
d 2  ( x  x ')2  ( y  y ')2  r 2  r 2 ' 2rr 'cos(   ')
(4.90)
For isotropic conditions (  x   y  r ), substituting Eq. (4.90) into Eq. (4.89) yields
S (r , r ', t   ) 
 r 2  r 2 ' 2rr 'cos(   ') 
1
exp  


4r (t   )
4r (t   )


In dimensionless form, Equation (4.91) becomes
46
(4.91)
S D (r , r ',, t D  t D ') 
 r 2  rD 2 ' 2rD rD 'cos(   ') 
1
exp   D


4 (tD  tD ')
4(tD  tD ')


(4.92)
The definition of t D was presented in Eq. (4.29) and the dimensionless radial distance rD is given by
rD 
r
rw
(4.93)
Recalling Eq. (4.71), the perturbation solution of the non-linear diffusion equation truncated at the first order
perturbations is given by
tD
mD ( xD , yD , zD , tD ) 
'
'
'
'
(0)
'
'
'
' 

'
'
'
'
(0)
'
'
'
' mD ( xD , yD , z D , t D )
 g ( xD , y D , z D , t D )  w ( x D , y D , z D , t D )

tD '


0 D

(4.94)
 GD ( xD , xD' , yD , yD ' , zD , zD ' , tD , tD ' )dxD ' dyD ' dz D ' dt D '
In Cylindrical coordinates, Equation (4.94) can be written as follows:
2 tD 
mD (rD , , tD ) 
   f (r
D
'
, ', tD ' )GD (rD , rD ' , , , t D , t D ' )rD ' drD ' dtD ' d '
(4.95)
0 0 0
The source term is given, in Cylindrical coordinates, by (Cole et.al 2010):
f ( xD ' , yD ' , z D ' , t D ' )  f (rD ' ,  ', t D ' ) 
f (rD ' , t D ' )
2 rD '

mD (0) (rD ' , t D ' ) 
'
'
 g (rD , t D ) 

t D '


'


2 rD




(4.96)
Substituting Eq. (4.96) into Eq. (4.95) and replacing the line source solution obtained from Eq.(4.92) yield

2 t D  
mD (rD ,  , t D ) 

  
0 0 0

g (rD ' , t D ' )  w(rD ' , t D ' )
2 rD '


mD (0) (rD ' , t D ' ) 

t D '




(4.97)
 r 2  rD 2 ' 2rD rD 'cos(   ') 
1
'
'
'
exp   D
 rD ' drD dt D d

4 (t D  t D ')
4(
t

t
')
D
D


Interchanging the order of integration, Eq.(4.97) becomes,
(0)
'
' 

'
'
'
' mD (rD , t D )
 g (rD , tD )  w(rD , t D )

tD '


tD  

2 rD '


mD (rD , , tD ) 



0 0
2
 r  rD ' 
 r r 'cos(   ') 
1
'
'
'
exp   D
  exp  D D
rD ' drD dt D d

4 (t D  t D ')
4(
t

t
')
2(
t

t
')
D
D
D
D



 0
2
2

47
(4.98)
The integral over  ' can be simplified as (Cole et.al 2010):
2
 rD rD 'cos(   ')  '
 rD rD ' 
d  2 I 0 

2(t D  t D ')

 2(t D  t D ') 
 exp 
0
(4.99)
Substituting Eq. (4.99) into Eq. (4.98), it simplifies to
tD 
(0)

(rD ' , t D ' ) 
'
'
'
' mD
 g (rD , t D )  w(rD , t D )

t D '


0 

mD (rD , t D ) 
0
 r 2  rD 2 '   rD rD ' 
1

exp   D
I
dr ' dt '
 4(t  t ')  0  2(t  t ')  D D
4 (t D  t D ')
D
D
D
D




(4.100)
Let’s define,
 r 2  rD 2 '   rD rD ' 
1
exp   D
I
 4(t  t ')  0  2(t  t ') 
4 (t D  t D ')
D
D
D
D


 
S D (rD , rD ', t D , t D ') 
(4.101)
Then, taking Laplace transformation of Eq. (4.100) and using the Laplace transform property of convolution integral
presented in Eq. (3.2), Eq.(4.100) can be written as


m D (rD , t D ) 

 h( r
'
D

, t D ' ) S D (rD , rD ', t D , t D ')drD '
(4.102)
0
where


mD(0) (rD' , t D' ) 
h(rD' , t D' )  L  g (rD' , t D' )  w(rD' , t D' )

tD '


(4.103)

and S D (rD , rD ', tD , tD ') is the Laplace transform of the source function in Eq. (4.101) which is given by the following
expressions (Carslaw and Jaeger 1954):
 1
I (r ' s ) K 0 (rD s )

 2 0 D
S D (rD , rD ', s )  
 1 I (r s ) K (r ' s )
0 D
0 D

 2
rD  rD '

(4.104)
rD  rD '
Recalling Eq.(4.102), the total pressure drop in dimensionless form can be written as,







m D (rD , t D )  m D 0 (rD , t D )  m D1 (rD , t D )  h(rD ' , t D ' ) S D ( rD , rD ', s) drD '
0


mD (0) (rD ' , t D ' )  
'
= L  g ( rD ' , t D ' )  w( rD ' , t D ' )
 S D ( rD , rD ', s) drD
'

t


D

0 


=
 L  g (r
D
'

, t D )  S D ( rD , rD ', s )drD

'
'
0


mD (0) (rD ' , t D ' )  
'
 L  w( rD ' , t D ' )
 S D ( rD , rD ', s) drD
'
t D


0 

where
48
(4.105)


(0)
mD
(rD , t D ) =


L  g (rD ' , t D ' )  S D (rD , rD ', s)drD '


(4.106)
0
and


(1)
mD
(rD , t D ) 

mD (0) (rD ' , t D ' )  
'
L  w(rD ' , t D ' )
 S D (rD , rD ', s )drD
'

t


D


0

(4.107)
Similar to Eq. (4.32), let’s define the dimensionless pseudo-pressure for radial geometry as follows:
mD 
 khTsc
qref pscT
m( p)
(4.108)
Let us Recall Eq. (4.37) for radial geometry:
~
h q sc ( x, y, z ) L2
g ( xD , yD , z D )  2
qref
(4.109)
Defining

~
q ( x, y , z , t )
q sc ( x, y, z, t )  sc
h
(4.110)
we can write

q ( x, y, z, t ) L2
g ( xD , yD , z D )  2 sc
qref
(4.111)
Transforming the coordinates from Cartesian to Cylindrical, Eq. (4.111) can be written as:

q (r , , t ) L2
g (rD , , t D )  2 sc
qref
(4.112)

If we define q sc (r , t ) as the source density of the surface of a cylinder ( 2 rh ), then


q (r , t )
q sc (r , , t )  sc
2 r
(4.113)
Substituting Eq. (4.113) into Eq.(4.112) yields

q sc (r , t ) 2 

L
q sc (r , t ) L2 q sc (r , t ) L
2

r
g (rD , , t D )  2


qref
rqref
rD qref
(4.114)
Also because (Cole et.al 2010)
g (rD , , tD ) 
g (rD , tD )
2 rD
(4.115)
we obtain
49


q (r , t ) L
q (r , t ) L
g (rD , tD )  2 rD g (rD , , t D )= 2 rD sc
 2 sc
rD qref
qref
(4.116)
This solution corresponds to a distributed source in the space 0  r   . For the case of a single cylindrical source at
r  r0 ,

q sc (r , t )  qsc (t ) (r  r0 )
(4.117)
In Eq. (4.117), qsc (t ) represents volumetric rate of fluid withdrawal from the porous media at time t .
Substituting Eq. (4.117) into Eq. (4.116) yields
g (rD , tD ) 
2 qsc (t ) L
 (r  r0 )  2 qD (tD ) (rD  rD0 )
qref
(4.118)
where qD (tD ) represents the dimensionless production rate from the cylindrical source at rD 0 and t D defined by
qD (tD ) 
qsc (t )
qref
(4.119)
and  (rD  rD0 ) is the Dirac delta function. In Eq.(4.118), the following property has been applied
 (cx) 
 ( x)
(4.120)
c

Recalling Eq. (4.106) for m D0 (rD , tD ) and substituting Eq. (4.118), we obtain


(0)
mD

(rD , s)=

L 2 qD (t D ) (rD  rD 0 )  S D (rD , rD ', s)drD '
(4.121)
0
Also applying the sifting property of the Dirac delta function given by

  (r
D
 rD 0 ) S D (rD , rD ', s)drD '  S D (rD , rD 0 , s)
(4.122)
0
Eq. (4.121) can be written as follows:


(0)
mD
(rD , s)=


L  2 qD (t D )  S D (rD , rD 0 , s)drD '
(4.123)
0
Considering uniform flux along the line source and using qref  qsc or qD (tD )  1 , we have

mD
(0)
(rD , s) =
2 
S D (rD , rD 0 , s)
s
For the case where the source is at
(4.124)
rD0  0 (Line source solution), substituting Eq. (4.124) for

S D (rD , rD0  0, tD , tD ') and using rD '  rD , Eq. (4.124) is simplified as follows:
50

(0)
mD
(rD , s ) =
2
s
1
I 0 (rD ' s ) K 0 (rD s )
rD '  0
2
(4.125)
1
= K 0 ( rD s )
s
As we are interested in the solution at the wellbore, rD  1 (or r  rw ), we obtain the following expression from Eq.
(4.125):

mD
(0)
(rD  1, s) =
1
K0 ( s )
s
(4.126)
Equation (4.125) is the general solution for a well in an infinite acting reservoir in Laplace domain. Once
(0)

mD
(rD , s) is known then the first order of perturbation can be numerically evaluated.

(1)

Recalling Eq. (4.107) for m D (rD , s) and substituting S D (rD , rD0 , s) from Eq. (4.104), we have

1

m (0) (r ' , t ' ) 
 K0 (rD s ) L  w(rD' , t D' ) D D D  I 0 (rD ' s )drD'
 2
tD'


0 
 (1)

m D (rD , s)  


(0)
'
' 

1
'
' mD (rD , t D )
'
I
(
r
s
)
L
w
(
r
,
t
)

 K 0 (rD ' s )drD
0
D
D
D
 2
'

t


D


0 



rD  rD '
(4.127)
rD  rD '
(1)
Our interest would be the solution of m D (rD , s) at the wellbore, rD  1 (or r  rw ); that is,

1

m (0) (r ' , t ' ) 
 K0 ( s ) L  w(rD' , t D' ) D D D  I 0 (rD ' s )drD'
 2
tD '


0 
 (1)

m D (rD  1, s)  


(0)
'
' 

1
'
' mD (rD , t D )
'
I
(
s
)
L
w
(
r
,
t
)

 K 0 (rD ' s )drD
0
D
D
 2
'

t


D


0 


rD  rD '
(4.128)
rD  rD '
Evaluating Equation (4.128) requires numerical integration over the interval of 0  rD   . The results are, then,
numerically inverted back to time domain (we have used the Stehfest numerical inversion algorithm for this
purpose). Also, it is required to evaluate w(rD' , tD' )
mD (0) (rD' , tD' )
tD '
at each rD ' and then either utilize Onur and
Reynolds (1998) (Eq. (4.80)) or Roumboutsos-Stewart (1988) (Eq. (4.81)) algorithm to take the tabulated data as a
function of t D ' into Laplace domain.
Figure 4.1 and Figure 4.2 show w as a function of dimensionless pressure and dimensionless time for an infinite
acting reservoir and gas properties presented in Table 4.1.
51
Figure 4.3 through Figure 4.5 show the results for the example data set provided in Table 4.1. Figure 4.3
illustrates the solutions for mD (0) and Figure 4.4 shows the solution of mD (1) . Notice that the values of mD (1) shown in
Figure 4.4 have a decline trend as time increases.
This point will be explained later in this chapter under the discussion section and an example will be given in
Appendix A. This behavior can also be physically explained. Initially, a large pressure drop is created in a small
vicinity of the wellbore where the pressure-dependency of the viscosity and compressibility has a significant
influence. As the radius of influence increases in time, the pressure behavior is dominated by the viscositycompressibility product farther away from the wellbore where the pressure drop and the change in viscositycompressibility product are smaller.
Figure 4.5 shows the mD solution ( mD  mD(0)  mD(1) ) obtained from this work. The effect of the flow rate on
the solution of mD (1) has also been investigated. Peres et al. (1190) and Barreto (2011) showed that the first order
dimensionless pseudo pressure ( mD (1) ) is rate dependent. For the same reservoir and gas properties shown in Table
4.1, two different flow rates are tested in Figure 4.6: 1,500 MscfD and 10,000 MscfD. As expected, higher
production rate yields larger values of ( mD (1) ). Therefore, higher production rates magnify the effect of nonlinear
terms in the solution.
Figure 4.1 Omega factor (w 
 ct  (  ct )i
) versus dimensionless pressure in an infinite acting reservoir.
(  ct )i
52
Figure 4.2 Omega factor (w 
 ct  (  ct )i
) versus dimensionless time in an infinite acting reservoir.
(  ct )i
Figure 4.3 Dimensionless pseudo-pressures ( mD (0) ) for the data set presented in Table 4.1 (Infinite Acting
Reservoir).
53
Figure 4.4 Dimensionless first non-linear term of pseudo-pressure ( mD (1) ) for the data set presented in Table 4.1
(Infinite Acting Reservoir).
Figure 4.5 Dimensionless pseudo-pressure ( mD  mD(0)  mD(1) ) for the data set presented in Table 4.1 (Infinite
Acting Reservoir).
54
Figure 4.6 Dimensionless pseudo-pressure ( mD (1) ); comparison of two different rates for the reservoir and gas
properties presented in Table 4.1 (Infinite Acting Reservoir).
Table 4.1 Data for infinite acting reservoir.
Pi
q
Temp
k
h
rw
Φ
psia
Mscf/day
F
md
ft
ft
fraction
1500
140
1
140
0.3
0.1
5000
55
SG
Cf
1/psi
0.7
3E-6
4.4.2
Solution for Wells in Closed Cylindrical Reservoirs
Ozkan and Raghavan (1991a and 1191b) presented Laplace domain solutions for a point source and a fully
penetrating well in closed cylindrical reservoirs for different boundary conditions.
The pressure drop solution in Laplace domain presented by Ozkan and Raghavan (1191a, and 1991b) for a point
source in a closed cylindrical reservoir with impermeable boundaries at zD  0 , zD  hD , and rD  reD is given by
the following expression,

 pp 

2


q 
 K0
2 k hD s 

cos n
n 1

Ik 

k 






 I 
sRD 
srD
k
k 


Ik
 
srD K k
I k

sreD

sreD
 cos k    

zD
z  
n 2 2
cos n D  K 0  s  2 RD 

hD
hD  
hD

 



n 2 2  
n 2 2

I k  s  2 rD  K k  s  2 reD 




2 2


h
h
D
D
n
 
 cos k      
s  2 rD  







hD
n 2 2

I k  s  2 reD 



h
 
D



(4.129)
where
RD2  rD2  rD '2  2rD rD 'cos(   ')
(4.130)
 n  u  n2 2 / hD 2
(4.131)
 2n1  u  (2n  1)2  2 / hD 2
(4.132)
and
In all above equations
k
/L
kz
hD  h
(4.133)
Here, L denotes the characteristics length. The corresponding line-source solution in a closed cylindrical reservoir
[source is located at an arbitrary point (rD ', ') ] is given by

p

q 
K0
2 ks 



 I 
sRD 
k 
k
srD

Ik

 
srD K k
I k

sreD

sreD
 cos k    


(4.134)
If we write the solution for a vertical well at the coordinate center (rD '  0) , then
lim I k
rD 0

0
srD  
1

k 0
(4.135)
k 0
and
56

p

q 
K0
2 ks 



sRD 
I0

 
srD K1
I1

sreD


sreD 



(4.136)

In our terminology, Eq. (4.136) corresponds to the solution of m D

(0)
(rD , s) . To formulate the first term of
(1)
perturbation, m D (rD , s) , we start from the Laplace domain solution for a point source in a cylindrical reservoir
presented in Eq. (4.129). Let us define,
_
_
G ( xD  xD ', yD  yD ', z D  z D ', s) 
 p p ( xD , y D , z D , s )
__
(4.137)
~
q ( xD , y D , z D )
and denote the dimensionless Green’s function by
_
_
GD ( xD  xD ', yD  yD ', z D  z D ', s) 
p D ( xD , yD , zD , s)
__
(4.138)
~
q D ( xD , y D , z D )
If we also assume uniform flux, then the flow rate is given in Laplace domain by
__
~
~
q( xD , yD , zD )
q ( xD , y D , z D ) 
s
(4.139)
where
__
~
~
q( xD , yD , zD )
q D ( xD , y D , z D ) 
qref s
(4.140)
Transforming the coordinates from Cartesian to cylindrical, Eq. (4.139) and Eq. (4.140) can be written as follows:
__
~
~
q( xD , yD , zD )
q (rD , , z D )  2
s
(4.141)
and
~
__
~
q ( xD , y D , z D )
qD ( xD , yD , z D )  2
qref s
(4.142)
Then Eq. (4.137) becomes
_
_
G (rD , rD ',  , s) 
 pp
(4.143)
~
q
2
s
or
57




K
sR

Ik


0
D
4 2 k hD 
k 



z
z  
n 2 2
2 cos n D cos n D  K 0  s  2 RD 

hD
hD  
hD
n 1

 

_
G (rD , rD ',  , s) 


n
Ik  s  2

2 2

hD
n
s  2 rD  


hD

I k 


2 2


Ik 

k 




srD

Ik

 
srD K k
I k

sreD

sreD
 cos k    
The dimensionless point source solution is obtained by multiplying both side of Eq. (4.144) by

2 kh _
1 
G (rD , rD ',  , s) 
 K0
q
2 q 


2 cos n
n 1


Ik 

k 




(4.144)

 

n

rD  K k  s  2 reD 
 


h
D
 
 cos k      



n 2 2 
s  2 reD 


hD
 

2 2

  
sRD 
k 
Ik
srD

Ik

 
srD K k
I k

sreD

sreD
2 kh
:
q
 cos k    
zD
z  
n 2 2 
cos n D  K 0  s  2 RD 

hD
hD  
hD

 


n 2 2  
n 2 2 

I k  s  2 rD  K k  s  2 reD 




2 2


h
h
D
D
n
 
 cos k      
s  2 rD  





hD

n 2 2 

I k  s  2 reD 



hD
 


(4.145)
_
Also, the dimensionless Green’s function, G(rD , rD ', , s) is obtained to be
_
2 kh _
G (rD , rD ',  , s )  GD (rD , rD ',  , s ) 
q

qref


K0
q 

2

qref
1

2 cos n
n 1


Ik 

k 





  
sRD 
k 
Ik
srD

Ik

 
srD K k
I k

sreD

sreD
 cos k    


zD
z 
n 2 2
cos n D  K 0  s  2 RD 


hD
hD 
hD







n 2 2 
n 2 2

I k  s  2 rD  K k  s  2 reD 




2 2


h
h
D
D
n 


 cos k      
s  2 rD  






hD

n 2 2

I k  s  2 reD 



hD


 

58
(4.146)
The solution for a line, surface, or volume source in a closed cylindrical reservoir can be obtained as:
_
_
S D (rD , rD ', s)   G D (rD , rwD ', s)dM w
(4.147)
w
where  w is the length, surface, or volume of the source. Also, the constant production rate from the length, area, or
~
q
.
w
the volume of the source,  w , is denoted by q so that q 
Applying of Eq. (4.147) on Eq. (4.146), the line source solution may be obtained as follows:
_
S D (rD , rD ', , s) 

K
q  0
2

qref
h

K
=
~  0
q 
2
qref
1



  
sRD 
k 
Ik
srD

 I 
srD

srD
sRD 
k 
k


Ik
Ik

I k

 
srD K k

sreD
 
srD K k
I k

sreD


sreD
sreD
 cos k    


 cos k    
(4.148)


Recalling Eq. (4.135), we have
_
S D (rD , rD ', , s) 

K
~  0
2 q D 
1

sRD  I 0


I0

 
srD K1
I1

sreD


sreD 



(4.149)
where
RD2  rD2  rD '2  2rD rD 'cos(   ')
(4.150)
Now, using the following theorem for Bessel function, K0
K0

sRD

 
  Ik
 k 


 
 I
k

 k 

 
srD K k

sRD


srD ' cos k (   ')
rD  rD '
(4.151)

 

srD cos k (   ')
srD ' K k
rD  rD '
and recalling Eq. (4.135), Eq. (4.151) can be written as,
K0

sRD

I0




I

 0


 
srD K 0
 
srD ' K 0
srD '
srD


rD  rD '
(4.152)
rD  rD '
Substituting Eq. (4.152) into Eq. (4.149), the line source solution can be obtained as follows:
59

 1
 ~
 2 q D
_

S D (rD , rD ', , s)  

 1
 ~
 2 q
D


I
 0


 
 
srD '  I 0
srD K 0
srD
I0


 
srD K1
I1

sreD


sreD 



rD  rD '
(4.153)

I
 0


 
 
srD  I 0
srD ' K 0
srD
I0


 
srD K1
I1

sreD


sreD 



rD  rD '
The Laplace domain solution for the first term of perturbation can now be written as

(0)

(rD ' , t D ' ) 
'
' mD
 w(rD , t D )

t D '
 S (r ',  , s)r ' dr ' d '
L
D
D
'

 D D
2

r
D
0 



2 
(1)

m D (rD , s) 
0
=
=
1
2
1
2
2 
(4.154)

m (0) (rD ' , t D ' )  
'
L  w( rD ' , t D ' ) D
 S D ( rD ',  , s) drD d '
'

t


D

0 

0


m (0) (rD ' , t D ' ) 
L  w(rD ' , t D ' ) D

t D '


0 

2

S
D ( rD
',  , s)d ' drD '
0
where

1
~
 qD
2


S D (rD ', , s)d '  

0
1
~
q
 D


I
0


I
0

  
srD K0
 sr  K  sr  
   I sr
  
srD '  I 0
srD
I0
D
1
1

  sr   I  sr 
srD ' K0
D
0
I0
D
eD
rD  rD '
eD
 sr  K  sr 
I  sr 
D
1
1
eD
eD
(4.155)




rD  rD '
Substituting Eq. (4.155) into Eq. (4.154), we have


I
 1  I sr K sr '  I sr 0
0
D
0
D
0
D
 ~ 
 2 q D 
 
(0)
'
' 


'
' mD (rD , t D )
 drD '
 L  w(rD , tD )
tD '

 0 
(1)



m D (rD , s)  



I0
 1 
 ~  I 0 srD ' K 0 srD  I 0 srD
 2 q D 

(0)
  
(r ' , t ' ) 
'
' m
 L  w(rD , tD ) D D' D  drD '
tD

 0 

 
 


 
srD K1
I1

sreD


sreD 



rD  rD '


 
 

(4.156)

 
srD K1
I1


60
sreD


sreD 



rD  rD '

~
~
(1)
Considering qref  q or q D (tD )  1 and evaluating m D (rD , s) at the wellbore ( rD  1 ), Eq.(4.156) can be written
as

 
I
 1  I s K sr '  I s 0
0
0
D
0
 2 
 
 
(0)
'
' 


'
' mD (rD , t D )

L
w
(
r
,
t
)
 D D
 drD '

tD '

 0 

 (1)

m D (rD  1, s)  

 
I0
1 
 2  I 0 srD ' K 0 s  I 0 s
 

(0)
  
(r ' , t ' ) 
'
' m
 L  w(rD , tD ) D D' D  drD '
tD

 0 
  
  
 
srD K1
I1

sreD


sreD 



rD  rD '


    

(4.157)
 
srD K1
I1

sreD


sreD 



rD  rD '


(1)
The solution to the first term of perturbation, m D (rD , s) , is evaluated numerically and inverted back into time
domain using a standard numerical inversion algorithm. The process of handling the integrations is similar to that
explained for the infinite acting reservoir system. The results obtained from this study ( mD (0) , mD (1) and mD ) are
shown in Figure 4.7 through Figure 4.10 for the data set presented in Table 4.2. Figure 4.7 illustrates the results for
mD (0) . Figure 4.8 depicts the mD (1) solution from this study. The difference here is in the boundary dominated flow
part, Once the effect of the boundary is felt, to maintain constant production rate, larger pressure drop is required
and this has a significant effect on mD (1) .
Figure 4.9 shows the mD solution ( mD  mD(0)  mD(1) ) obtained from study. The effect of the production rate on
mD (1) is also shown in Figure 4.10 for the reservoir and gas properties presented in Table 4.2. Two different rates
are tested: 2,000 MscfD and 10,000 MscfD. As expected, mD (1) is larger for the higher production rate indicating
more significant contribution of the non-linear terms to pressure responses.
4.5
Discussions
Barreto (2011) and Kale and Mattar (1980) showed that the solution for the first perturbation term ( mD (1) ) reach
a plateau in an infinite acting reservoir as time increases. However, our results discussed above illustrated a decline
as time increased. Comparing our solution with Kale and Mattar (1980), it has been found that, in their solution the
second order term and higher order differentials have been neglected and their solution reduced to an ordinary
differential equation. This assumption removes the Green’s function effect from their solution, which explains the
differences from the results of this study.
61
Figure 4.7 Dimensionless pseudo-pressure ( mD (0) ) for data set presented in Table 4.2 (Closed Cylindrical
Reservoir).
Figure 4.8 Dimensionless first non-linear term of pseudo-pressure ( mD (1) ) for data set presented in Table 4.2
(Closed Cylindrical Reservoir).
62
Figure 4.9 Dimensionless pseudo-pressure ( mD  mD(0)  mD(1) ) for data set presented in Table 4.2 (Closed
Cylindrical Reservoir).
Figure 4.10 Dimensionless pseudo-pressure ( mD (1) ) comparing two different rates for reservoir and gas
properties presented in Table 4.2 (Closed Cylindrical Reservoir). The results show that, mD (1) is rate dependent.
63
Table 4.2 Data for closed cylindrical reservoirs.
Pi
q
Temp
k
h
rw
re
psia
Mscf/day
F
md
ft
ft
ft
140
1
140
0.3
800
5000
2000
Φ
SG
fraction
0.1
Cf
1/psi
0.7
3E-6
Similar to Kale and Mattar (1980), Barreto (2011) also showed that mD (1) should reach a plateau as time
increases. The solution provided by Barreto (2011) has the same fundamental approach and assumptions as in this
study, except that it is in time domain. The difference between our results and those of Barreto is expected to be in
the different numerical evaluation procedures used. As shown in Appendix A, it can be verified theoretically that
mD (1) should be a decreasing function of time increases. To verify our results further, we have used several other
approaches also.
Figure 4.11 is a comparison between the Green’s function in time domain and Laplace domain for rD  1 . The
results verify the accuracy of our numerical inversion of the Green’s function from the Laplace domain to time
domain. Similarly, in Figure 4.12, the Laplace domain and time domain results for w(tD' )
mD(0) (t D' )
t D'
are compared
for rD  1 . Except for some points of numerical instabilities of the Laplace inversion, the two results match very
well. Figure 4.11 and Figure 4.12 show that both the Green’s function and w(tD' )
mD(0) (t D' )
t D'
are positive and
deceasing functions of time. Therefore, when these two functions indicate a declining trend for different rD , then the


expression w(rD ', tD' )
0
mD(0) (rD ', t D' )
t D'
drD ' must decline in time also.
64
Figure 4.11 Comparing Green’s function for both Laplace and time domain at rD  1 .
Figure 4.12 Comparing w(tD ' )
mD(0) (tD' )
tD'
65
for both Laplace and time domain at rD  1 .
CHAPTER 5
CONCLUSIONS AND RECOMMENDATIONS
5.1
Conclusions
There were two primary objectives carried out under this research study. The first objective was to improve the
deconvolution of variable-rate data in the Laplace domain. The main hypothesis of the first objective was that the
limited success of Laplace-domain deconvolution and transformation of tabulated data into Laplace domain was due
to algorithmic issues. The conclusions of the first objective can be summarized as follows;

A new algorithm, which used inverse mirroring at the points of discontinuity and adaptive cubic spline
method to approximate rate or pressure versus time data is introduced.

One simulated and three field examples demonstrated the success of the algorithms presented in
this work in deconvolution of variable-rate production data.

Nadaraya-Watson kernel estimator with the application of adaptive bandwidth using
Epanechnikov and Gaussian kernel also demonstrated a successful smoothing approach.

This research showed that the smoothing parameters and the choice of the regression algorithm
could significantly alter the behavior of the deconvolved pressure and its derivative.

The algorithms presented in this dissertation have broader application in a variety of PTA and RTA
problems.
The second objective of this research study focused on the perturbation solution of the nonlinear diffusivity
equation for flow in porous media. Although the procedure is applicable to other non-linear flow problems also, the
particular solution developed in this dissertation accounts for the nonlinearity caused by the dependency of gas
properties (viscosity-compressibility product changes) with pressure. This work is an extension of the original work
of Barreto (2011) into Laplace domain. The analytical solution developed in this study is illustrated by two
examples (infinite acting and closed cylindrical reservoirs). The computed results indicated some differences from
those of Barreto’s (2011) and Kale and Mattar (1980). We have verified that theoretically that mD (1) should be a
decreasing function as time increases.
5.2
Recommendations
The following are the areas, which have potential to improve and could extend the use of the work presented
herein:
1.
Numerical Laplace inversion algorithm presented by Iseger applied in this work is more accurate
compare to Stehfest algorithm when the function to be inverted includes discontinuities. However,
the Iseger algorithm is significantly slower than the Stehfest algorithm since it requires linear
sampling. Modification of the Iseger algorithm for logarithmic sampling may improve its speed
without compromising the accuracy.
66
2.
Cubic splines combined with piecewise linear approximation presented in this work demonstrated
a very good combination to take sampled data into Laplace domain. However, for more accuracy,
this can be extended for different combinations such as second order cubic splines or second order
polynomials combined with cubic splines.
3.
In the solution for non-linear gas diffusivity equation presented in this work, Stehfest algorithm is
utilized. Using cubic splines with Iseger numerical inversion algorithm is not practical due to its
slow performance. Once the Iseger algorithm is modified for logarithmic sampling, it is
recommended to replace the Stehfest algorithm by the Iseger algorithm to remove any instabilities
due to the Stehfest algorithm.
67
NOMENCLATURE
Bo
Oil formation volume factor, rbbl/stb
b
Intercept of a straight line used in Eq.(3.16)
cf
Formation compressibility factor
C
Wellbore storage (bbl/psi)
ct
Total compressibility, psi-1
fi
Value of f (t ) at time ti
f t 
Arbitrary function
f s 
Laplace transform of f (t )
f i
Linear slope of the data in the interval ti  t  ti 1
h
Constant kernel bandwidth parameter and reservoir thickness
g ( x | h)
Nadaraya-Watson kernel estimator used in Eq.(3.21)
GD
Dimensionless Green’s function
h( xi )
Adaptive bandwidth
I0
Bessel function of first kind order zero
I1
Bessel function of first kind order one
J
Arithmetic average of data used in Eq.(3.25)
K0
Bessel function of second kind order zero
K1
Bessel function of second kind order one
K
Kernel function
L
Characteristic length, ft
L
Laplace transform operator
L-1
Inverse Laplace transform
M
Number of points at which the inverse Laplace transforms are computed, Eq. (3.18)
mD
Dimensionless pseudo-pressure
m( p)
Pseudo-pressure
m
Slope of a straight line used in Eq. (3.14)
N
Parameter used in the Stehfest algorithm
Pi
Initial reservoir pressure, psi
PwD
Dimensionless bottomhole pressure
68
PD
Dimensionless pressure
pM , t 
Pressure drop at a point M and time t , psi
pc M , t 
Time derivative of pressure drop at point M and time t for constant unit rate, psi
pw
Wellbore pressure drop, psi
pw
Laplace transform of wellbore pressure drop, psi
Pwc
Constant rate pressure response, psi
q
Production rate, rbbl/d, stb/d
q
Production rate in Laplace space, rbbl/d, stb/d
qD
Dimensionless rate
rw
Wellbore radius, ft
re
Reservoir boundary radius, ft
s
Laplace-transform parameter
Si
Cubic spline over a segment
SG
Gas specific gravity
T
Period for which the inversions are computed, Eq. (3.18)
t
Time, hr
tD
Dimensionless time
Z ( p)
Gas deviation factor
69
Greek
α
Constant used in Eq. (4.80)
α0
Constant used in Eq. (4.80)
B
Formation volume factor, rbbl/stb
Β
Constant used in Eq. (4.80)
β0
Constant used in Eq. (4.80)

Parameter used in the Iseger algorithm

Porosity of reservoir rock, fraction
μ
Viscosity, cp
η
Diffusivity constant, md-psi/cp
i
Bandwidth adaptive parameter

Standard deviation and dispersion parameter

Gamma function

Perturbation parameter
70
REFERENCES
Abate, J. and Valkó, P. P. 2004. Multi-Precision Laplace Transform Inversion. Computers Math.Applic.48: 629-636.
Abramson, I. 1982. On bandwidth variation in kernel estimates-a square root law. The Annals of Statistics, Vol. 10.
Page 1217-1223.
Agarwal, R.G. 1980. A New Method to Account for Producing Time Effects When Drawdown Type Curves Are
Used to Analyze Pressure Buildup and Other Test Data," Paper SPE 9289 presented at the SPE Annual
Technical Convention and Exhibition, Dallas, TX, 21-24Sep.
Al-Ajmi, N., Ahmadi, M., Ozkan, E., and Kazemi, H.2008. Numerical Inversion of Laplace Transforms in the
Solution of Transient Flow Problems with Discontinuities. Paper SPE 116255 at SPE Annual Technical
Conference and Exhibition. Denver, CO, 21-24 Sep.
Al-Hussainy, R., Ramey Jr., H. J. and Crawford, P.B. 1996. The flow of real gases through media, Journal of
Petroleum Technology, SPE Paper 1243 A, vol. 18, pp.624-636.
Atkinson, K. 1985. Elementary numerical analysis,3rd edition. John Wiley & Sons, INC.
De Boore, C. 2001. A Practical Guide to Splines. Applied Mathematical Sciences,Vol 27.
Barreto, J., A.B. 2011. Non Linear Gas Diffusivity Equation Solution by Green’s Functions. PhD Thesis (in
Portuguese). Universidade Estadual Norte Fluminense, Macae’, Brazil.
Baygun, B., Kuchuk, F.J., and Arikan, O. 1997. Deconvolution Under Normalized Autocorrelation Constraints,
SPEJ (September) 246.
Bellman, R., Kalaba, R. E., and Lockett, J. A. 1966. Numerical Inversion of Laplace Transform: Applications to
Biology, Economics, Engineering, and Physics. American Elsevier Publishing Co. Inc., New York, NY.
Chen, C. C. and Raghavan, R.1996. An Approach to Handle Discontinuities by the Stehfest Algorithm. SPEJ (Dec.):
363.
Cinar, M., Ilk, D., Onur, M., Valko, P.P, and Blasingame, T.A. 2006. A Comparative Study of Recent Robust
deconvolution Algorithms for Well-Test and Production-Data Analysis. Paper SPE 102575 presented at the SPE
Annual Technical Conference and Exhibition, San Antonio, Tx, 24-27 Sep.
Coats, K. H., Rapoport, L. A., McCord, J. R., and Drews, W. P. 1964. Determination of Aquifer Influence Functions
from Field Data. JPT (Dec.): 1417-1424.
Correa, A. C. and Ramey, H. J. Jr. 1986. Combined Effects of Shut-in and Production: Solution with a New Inner
Boundary Condition. Paper SPE 15579 presented at the SPE Annual Technical Conference and Exhibition, New
Orleans, LA.5-8 Oct.
Crump, K. S. 1976. Numerical Inversion of Laplace Transforms Using a Fourier Series Approximation. J. ACM
233: 89-96.
Duhamel, J. M. C. 1833. Mémoire sur la méthode générale relative au mouvement de la chaleur dans les corps
solides polongé dans les milieux dont la température varie avec le temps. J. de Ec. Polyt. (Paris) 14: 20-77.
Fetkovich, M. J. and Vienot, M. E. 1984. Rate Normalization of Buildup Pressure by Using Afterflow Data. JPT
(Dec.): 2211-2224.
Gaver, D. P. 1966. Observing Stochastic Process and Approximate Transform Inversions. Oper. Res. 14, 3: 459.
71
Gladfelter, R.E., Tracy, G.W., and Wilsey, L.E. 1955. Selecting Wells Which Will Respond to ProductionStimulation Treatment, Drill. and Prod. Prac.,117-129.
Gupta, K. C., and Andsage, R. L. 1967. Application of variable rate analysis technique to gas wells, Fall Meeting of
the Society of Petroleum Engineers of AIME, New Orleans, Louisiana, SPE Paper 1836.
Hahn, T., 2005a. Cuba - a library for multidimensional numerical integration, Computer Physics Communications,
no. 168, pp. 78-95.
Ilk, D., Valkó, P. P., and Blasingame, T. A. 2005. Deconvolution of Variable-Rate Reservoir Performance Data
Using B-Splines. Paper SPE 95571 presented at the SPE Annual Technical Conference and Exhibition. Dallas,
TX. 9-12 Oct.
Ilk, D. 2005. Deconvolution of Variable Reservoir Performance Data Using B-Splines. Master of Science Thesis.
Texas A&M University.
Iseger, P. D. 2006. Numerical Transform Inversion Using Gaussian Quadrature. Probability in the Engineering and
Informational Sciences, 20. 1-44.
Kale, D. and Mattar, L., 1980. Solution of a non-linear gas flow equation by the perturbation technique, JCPT,
PETSOC 80-04-06, pp. 63-67.
Krayszig, E. 1999. Advanced Engineering Mathematics. John Wiley & Sons, INC.
Kucuk, F.J. and Ayestaran, L. 1985. Analysis of Simultaneously Measured Pressure and Sandface Flow Rate in
Transient Well Testing, JPT (February), 323.
Lamm, P. K. 2000. A Survey of Regularization Methods for First-Kind Volterra Equations. Surveys on Solution
Methods for Inverse Problems. ed. D. Colton, H. W. Engl, A. Louis, J. R. McLaughlin, W. Rundell. Springer,
Vienna and New York: 53-82.
Levitan, M. M. 2003. Practical Application of Pressure-Rate Deconvolution to Analysis of Real Well Tests. Paper
SPE 84290 presented at the SPE Annual Technical Conference and Exhibition, Denver, CO. 5-8 Oct.
Levitan, M.M., Crawford, G.E., Hardwick, A. 2004. Practical Considerations for Pressure-Rate Deconvolution of
Well Test Data. Paper SPE 90680 presented at the Annual Technical Conference and Exhibition, Houston,
Texas, 26-29 Sep.
Levitan, M.M., Wilson, M.R. 2010. Deconvolution of Pressure and Rate Data From Gas Reservoirs With Significant
Pressure Depletion. Paper SPE 134261 presented at the Annual Technical Conference and Exhibition, Florence,
Italy, 20-22 Sep.
Mendes, L.C.C., Tygel, M., and Correa, A.C.F. 1989. A deconvolution algorithm for analysis of variable-rate well
test pressure. Paper SPE 19815 presented at the SPE Annual Technical Conference and Exhibition, San
Antonio, TX. 8-11 Oct.
Meunier, D., Wittmann, M. J., and Stewart, G. 1985. Interpretation of Pressure Buildup Test Using In-Situ
Measurement of Afterflow. JPT (Jan.): 143-152.
Nadaraya, E. A. 1965. On nonparametric estimates of density functions and regression curves. Theory of probability
and its application. Vol. 10, Issue 1, page 186-190.
Nisle, R.G. 1958. The effect of partial penetration on pressure build up in oil wells. SPEJ 971-G (March 1958).
Odeh, A.S. and Jones, L.G. 1965. Pressure Drawdown Analysis, Variable-Rate Case, JPT (Aug): 960-964; Trans.,
AIME 234.
72
Onur, M. and Reynolds, A. C. 1998. Numerical Laplace Transformation of Sampled Data for Well-Test Analysis.
SPE REE (June): 268-277.
Ozkan, E. 2003. Applied Mathematics of Fluid Flow in Porous Media. Lecture Notes, Laboratory Exploration of
Petroleum Engineering.
Ozkan, E. and Raghavan, R. 1997. Some Strategies to Apply the Stehfest Algorithm for a Tabulated Set of
Numbers. SPE Journal. 2 (Sept.): 363-372.
Ozkan, E. and Raghavan, R. 1991a. New Solutions for Well-Test-Analysis Problems: Part 1 - Analytical
Considerations. SPEFE 42 (1): 359-368. SPE-18615-PA. “DOI: 10.2118/18615-PA.”
Ozkan, E. and Raghavan, R. 1991b. New Solutions for Well-Test-Analysis Problems: Part 2 - Computational
Considerations and Applications. SPEFE 42 (1): 369-378. SPE-18616-PA. “DOI: 10.2118/18616-PA.”
Peres. A. M. M., Serra. K. V., and Reynolds. A. C. 1990. Toward a unified theory of well testing for nonlinearradial-flow problems with application to interference tests, SPE Formation Evaluation, vol. 5, pp. 151_160,
June, SPE Paper 18113.
Raghavan, R.1993. Well Test Analysis, Prentice-Hall, New Jersey.
Roumboutsos, A. and Stewart, G. 1988. A Direct Deconvolution or Convolution Algorithm for Well-Test Analysis.
Paper SPE 18157 presented at the SPE Annual Technical Conference and Exhibition. Houston, TX. 2-5 Oct.
Samaniego, F., and Cinco-Ley, H. 1991. Transient pressure analysis for variable rate testing of gas wells, Low
Permeability Reservoirs Symposium, Denver, Colorado, SPE Paper 21831.
Silverman, B.W. 1986. Density estimation for statistics and data analysis. Published in monographs of statistics and
applied probability, London: Chapman and Hall.
Soliman, M.Y. 1981. New Technique for Analysis of Variable Rate or Slug Test. Paper SPE 10083presented at the
SPE Annual Technical Convention and Exhibition, San Antonio, TX.05-07 Oct.
Stehfest, H. 1970. Numerical Inversion of Laplace Transforms. Communications, ACM13 (1): 47–49.
Talbot, A. 1979.The Accurate Numerical Inversion of Laplace Transforms, J. Inst. Maths. Applics., 23: 97-120.
van Everdingen, A. F. and Hurst, W. 1953. The Application of the Laplace Transformation to Flow Problems in
Reservoirs, Trans. AIME 198: 171-176.
Von Schroeter, T., Hollaender, F., and Gringarten, A. C. 2004. Deconvolution of Well-Test Data as a Nonlinear
Total Least-Squares Problem. SPEJ (Sept.): 375-390.
von Schroeter, T. and Gringarten, A. C., 2007. Superposition principle and reciprocity for pressure transient analysis
of data from interfering wells, SPE Journal, SPE Paper 110465.
Winestock, A.G. and Colpitts, G.P. 1965. Advances in Estimating Gas Well Deliverability, JCPT (July-Sept.) 111119.
Watson, G.S. 1964. Smooth regression analysis. Sankhyā: The Indian Journal of Statistics, Series A (1961-2002),
Page 359.
73
APPENDIX A
EXAMPLE OF CONVOLUTION INTEGRAL
Here, a simple example is used to demonstrate that the convolution integral is not necessarily increasing as time
increases. The dimensionless pressure solution for an infinite acting reservoir in time domain is given by
mD (rD  1, t D ) = 
1
1
Ei(
)
2
4t D
(A.1)
The counterpart of this solution in Laplace domain can be written as follows;

m D (rD  1, s) =
1
K0 ( s )
s
Figure A.1 depicts the pressure derivative
(A.2)
mD (tD )
mD (tD )
in an infinite acting reservoir. The function
is a
tD
tD
positive and declines as time progresses.
Figure A.2 illustrates the source function for an infinite acting reservoir as a function of dimensionless time. The
source function in time domain is given by
S D (t D , t D ') 

 

1
1
1
exp  
 I0 

4 (tD  tD ')
 2(tD  tD ')   2(tD  tD ') 
(A.3)
Equivalently, the source function in Laplace domain can be written as follows:

S D (s) 
1
I 0 ( s ) K0 ( s )
2
(A.4)
As time increases, the source function declines due to the exponential term. The fact that both the pressure
derivative,
tD

0
tD

0
mD (tD' )
tD'
mD (tD' )
tD
'
mD (tD )
, and the source function are positive, does not indicate that their convolution integral,
tD
GD (tD  tD ' )dtD ' , must increase as t D increases. Figure A.3 confirms that the convolution integral,
GD (tD  tD ' )dtD ' , is not necessarily an increasing function; it might be a decreasing as t D increases because
the convolution integral of two functions is not a simple sum of the products of the two functions.
74
Figure A.1 In an infinite acting reservoir
mD (tD )
is a decreasing function as time increases.
tD
Figure A.2 In an infinite acting reservoir Greens function as a decaying function as time increases.
75
tD
Figure A.3 In an infinite acting reservoir

0
mD (tD' )
tD
'
GD (tD  tD ' )dtD ' is a decaying function as time increases.
76