Deconvolution of Wellbore Pressure and Flow Rate

Fikrl J. Kuchuk,
Research Center, SPE, Schlumberger-Doll Research Center, SPE, Schlumberger and LWS Ayestaran,

Richard
Technical

G. Carter,
Services,

NASA Dubai

Langley

FE

1<

Summary. Determination of the inffuence function of a welllresewoir system from ihe deconvolution of weUbore flow rate and pre.ssare is prasented. Deconvolution is fimdamental and is pa&aiady applicable to system identification. A variety of different deconvolution algorithms are presented. The simplest algoritlrn is a dmect method that works we!J for data without measurement noise but that I%& in the me.sence of even small amounts of noise. We show. however. that a modified algorithm that inwoses constraints on the solution set works very welf, even with significant measurementenors.

introduction h reservoir testing, we geeratly know the characteristic features tion integral in connection with amfps of simukam.mrsly mea8of the sywtem tlom its mnstmt-flow-rate and constant-prsssure beured wellbore pressure and flow rate. Although we restrict our havior. hrs it is important to determine the mnstaqt:rate or discussion mostly to determining influence functions for the -pressure behavior of the system for the identification of IIS charconstant-rate case, we do trsat the convolution integral in a generacteristic feature?.. For instance, identification of a one-half slope al manner. In other words, the influence function can also be the on a log-log plot of the pressure data may indicate a vertically fracsolution of the constant-pressure case. tured well, as two parallel straight limes on a Homer graph may k a linear causal system (reservoir), tie relationship behvssn indicate a fractured reservoir. The presence of either wellbore input (the Iime-depndent boundary condition that can be either the storage or flow-rate variations, however, usually masks characterflow rate or pressure) and output (the system response meamred istic system behavior, pardcufady at early tines. For many sysas either the flow rate or pressure) at the welfbore can be described tems, it is desirable to have a wellbore pressure that is free of as a convolution operation. 10We let the quantities measured at the wellbom-stmage and/or variable-flow-rate effects to obtain inforwellbore, above the sandface, be pti=wellbore pressure, mation about the welflreservoir geometry and its parameters. For Q~=mmdative wellbore production, and g~=wellbore flow example, the effects of partiaJ penetration, hydraulic bctures, sarate. lution gas within the vicinity of the weltbore, gas cap, etc., on the The convolution integral is wellbore pressure can k. masked entirely by wdfbore storage, flowrate vwiations, or both. g(t) = f(r)K(t-r)dr Although deconvolution of pressure and flow rate has not been t commonly UWAfor reservoir engineering problems, one can still =K(0)f(t)-f@f(0)+ j fC(7)f(t-r)d~, . (1) fmd a few works on deconvolution (computing intluence function) in the petroleum engineering fiteraturs. Ffruchinson and SikOm, 1 where K(O)=&,. For th~ constant-flow-rate case, g(t) =&#(t)= Jones et al.,2 and Coats et al. 3 presented methcds for determinh,.f(r=O)-PM(Ol> K(t) =APL(t), and jW=AqZJ(O=[q@(r=O) ing the influence function direcOy from field data. %@kJ,, while for the constant-pressure case, sO)=AQ@ = Jargon and van PooUen4 were perhaps fxst to we tie dec&volQ@(t=O)-Q@(f)l, ~(t)=AQj (z), f(t)=~n(t)= ~ (z=o)lution of welfbore-flow-rate and pressure data to compute the q~(t)[, aII~ K(I)= constant-rate behavior (the influence function) of the formation in PW(Ol/P, or 8(~)=Aq~(t)=lq~(t=o) well testing. Bostic et al. 5 used a deconvolution technique to ob- ~j;~ #@)&:,tie t.~ K(OY(O-KW(0) can be replaced by ..taia a constant-rate solution &om a variable-rate and -pressure frisThe functions Apw(t) and Qw(t) or qw(t) are the solutions of the tory. They also extended the deconvolution technique to combine diilkivi~ equation for the ccmstant-tlow-mte or -prsssure case with production and buildup data as a single test. Pasca16 ako used or without wellbore-storage and skin effects. Although it is usualdeconvolution techniques to obtain a constant-rate solution from ly small, the deconvolved pressure will always be affectd by the variable-rate (measured at the surface) and -pressure measurements wellbore volume between the measurement point and the sandfice of a drawdown test. Kutk and Ayestamm7 presented several because the sandface flow rate is different from the wellbore flow deconvolution methods inclding the Laplace transform and curve tit. Thompson et al. 8 and Thompson and Reynoldsg presented a rate, qti, when it is measured by the flowmeter at some wellbore location above the perforations. stable integration procedure for deconvolution. As shown by Coats et al.,3 the general solutions of the difdrsivThk paper focuses on deconvolution methcds. Mathematically, of internal IJOCOIIity WMtion with $e T$ and second w the deconvolution operation can be defined as obtaining solutions ditiom and nonpentic tidal and outer-bmmdaty conditions, satisfy for convolution-type, linear, Volterra integral equations. In resarthe constraints vou testing, it is defined as determining the pressure behavior (iifluence fiction or unit.respmse behavior) of a system from K(t)~OfOr 2~0 . . . . . . . . . . . . . . . . . . . . .. .. . . . . . . . . . ...@) simultaneously measared downlole pressure and flow rate. In other K(t) >OfOrt> O, . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ...(3) words, deconvolution computes the pressure behavior of a well/resemOir system as if the well were prcducing at a constant K(t)SOf Orr 20, . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ...(4) rate. We calf the computed pressure behavior of the system deand K(f) aOforr> O, . . . . . . . . . . . . . . . . . . . . . . . . . . . ...(5) ccmvolved pressure.

.,.,

Convolution

Integral

when the real time is >1 sezcmdand if the diffusivity constant,

Theconvolution integral, which is a special case of the Volterra integral equations, is widsly known for providing techniques for solving time-dependent boundary-value problems. It is also known as the superposition theorem (i.e., Dubamels principle) and has played an iqmrtant part in transient welf-tmt analysis. Jn recent years, there ba8 been more interest in the solution of the convoluCopwhht19908ccieiyal i%roleumEng[mers SPE Formation
EvahMiO!l,

kl@ct, is not very small. Coats et al. 3 amd linear programming with the above constmints to compute K(t) from measured g(t) and
f(r). Here we use these constrains ence function. to compute the system irdJu-

Deconvolutkrn Without Measurement Noise

Deconvolution is the determination of the inffuence function, K, for the constan-ffow-rate case (or the constant-pressure case with 53

March 1990

1400.0

I
, 1300.0 , .+

..
.

~%+ . .

1200.0

~/
i: :/ +

+ demnvolved 1[ . 2-DeconvoIved 3r. Tree, t, hr ,@ lime,t,hr

exact deconvolved

1100.0 ,..2

pressure for the flow rate with 0.1

Fig. lDeconvolved

pressure for the flow rate without erro

the appropriate deftions off and g above) by computing the pressure behavior from any arbitmmy pressure and flow-rate. history for a welfhesewoir system as though the well prcduces at a constant rate. Unfike the conventional and convolution imapretation methods, deconvolution of welfbore pressure and flow rate does not require a prerequisite bMluence funodon. Thus, the computed influence function (deconvolve.d pressure) can Lwused directly for parameter estimation as welf as for system identification. A direct discretimtion soheme for numerical deconvolution5.~8.11 can be written as a lower-triangular linear system, LT. b, witbtbe cm~ncnts of ~repmsenting the vafues of Kand the components of b representing the measured values of g at various dme points. The solution of this system afso can be written dizectly as

n-1

org?(t.)=

K(r)f(tn ;=1 ~,_l

x ] .,

-r)d~,

............

. .

(9)

where Kin Eq. 9 includes skin effwt. Substitution of the simple piecetiso functions for the approximation of K and fever the intervals [r+~ ,2J then yield linear system3 of tie form AX=Z . . . . . . . ... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. (10)

f(tn+, r,+,)]}
and K(~o)=g(tl)/~(tl),

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ..(6I
. . . . . . . . . . . . . . . . . . . . . . . . . . . . ...(7)

where F represents an approximation to Kand ~ contains the pressure measurements g(tn). The form of the A matrix depends on the functions used to approximate Kaodj as well as the scheme used to represent K(I) in terms of the solution vector Z For example, approximzdngf(t) by the piwewise linear timtioosf~t) =j +&i+, -fi)(t-ti)/(ti+l ti) fOr ti~t~ti+l ad K(t) by tie w=me ~nstant functionsK(t) =K(ti) for t;s t= tj+ ~ io Eq. 9 yields a Iower*$@r SYSkm eqtivdent to the deconvolution formufa given by ~S. 6 and 7. In this work, we chmse,f(r) to be the piecewise. finear functions, as before, but use more sophisticated approximating functions for

K(t). Note that if the weflbore storage and skin are zero, then K(f)=~ . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .(11)

here ~=+-i=) =O:N doses Thompson et al. dmcuss the selechon of smtable Oi WILWSin terms of a stable and nonosciflato~ deconvolution solution if the flow-rate data are free. of measurement noise. F@. 6 and 7 reduce to Eq. 23 in Ref. 7 if 0;=% and to F@ 8 in Ref. 5 if~i=l. A variety of deconvolution schemes can be generated by writing Eq. 1 as

.$

g 1300,0 % 3! = 1200.0 c .-

% ~ 1100.0 6

: . H
g(tn)=@,f(tn)+ n-l Zi z ~ K(T)f(tn ~)d~ i-l ,i-l 1500.0,
1400.0

k the early-time approximation for the solution of the *ivity eWatiOn for anY ZYztmnthat is intemalfy bounded by a smooth surface (sphere, cyfinder, t%cture, etc.) with a constant-flow-rate boundary. fftbe wellbore storage is zero with a nonzero skin, then

K(t)=ti+&.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ..(12)

. . . . . . . . . . . . . . (8)

ff both the wellbore storage and skin are nonzero, then K(t)=ti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ..(13)

.% g

1400.0 . .+

~ i 300.0 :. t? .: g 1200.0 g

+ .+

+ .,. +

q ,

+:

exact

exact

0. 1100.0 E , ~2 , ~, Tree, t, hr

+ decunwlved Id

deconwlved

1000.0

, .2

, ~.,

, ~o

Xme, t, hr

Fig. 4-Deconvolved pressure for the flow rate with 0.1% error with finite-element sw of 4At. SPE Formation Evahradon, March 19S0

34

1400.0

1500.01

.-

~ 1100.0 j 0.1

+ 1U*

E /
Ig

$100.0 -.

exact
deconvolved
1(

.
,@

t 000.0I
1[

,@

,IJ.*

Tree, t, hr Fig. 6DeconvoIved pressure for the flowrate with 2% normallydistributed emorwith the monotone-increasing,seoondierlvafive option.

Tree, t, hr ?g. 5Deoonvolved pressure for the flow rate with 0.19 ?rrorwith finite-element size of 8At

is the early-time approximation of the ftite welfbore solution with the wellbore-storage and skin effects. This solution afso approximate.! the case where the wellbore storage is tinite with a zero skin. An appropriate general approximation at early times is thus K(f)= at+~+AP, . . . . . . . . . . . . . . . . . . . . . . . . . . . . ..(14)

for @O,tl]. For afl other time intervals, K is approbated by smcoth piecewise quadratic functions: The solution vector for K, X is given in terms of the unknowns Ap,, K(tI), K(tI), K(t2). .K(tN-l). It is then straightforward, i~somewhat tedious, to express Eq. 9 in the discrete form A 7= b, where A is an NxN+ 1 matrix. Appendix A presents the details of the derivation, along with the procedure for reconstmcting K(t) from the sohmion vector X Fig. 1 demomtrates what happens when we use this discretization on synthetic test data with 100 points and no error in f(t) or

the mesh forf(r) from the tial regions must match with the mesh for K(r) for the initiaf regions when we are integrating K(T)f(tn r). To simpfify the computation, we will increase the mesh size only for K(t). llds is equivalentto simply detinirig groups of adjacmt swondderivative values on the original mesh to be equal to each other. Such equal K(t) values can be represented by a single variable in Z and the corresponding columns of A can be replaced by a siogle column that is the sum of alf the others in the group. BecauseEq. 10 is now overdetennird, we solve the linear least-squares problem by minimizing

1%(~AY)T(~-A7).

. . . . . . . . . . . . . . . . . . . . . . . . . . . (15)

Deconvolution With NOISY Flow.Rate Measurements All measurements, no matter how carefufly obtained, are subject to errors. At the wellbore, errors in pressure measurements usualIy are smaller than errors in ffowmeter measurements. Sometimes errors can be verj large in flowmeter measurements, depending on wellbore geometries and fluid twe. Kuchuk12 ;howed that tbe deco;;olution algorithm, given by RI. 6 or any other direct deconvolution method, does not work if the downbole flow rate contains measurement noises. Let us assume that the dovmhole flowrate with noises can approximated by
~(r)=f(t) +,(t). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ..(16)

g(t). Note that the discretization gives a soIution that is abnost indistinguishable from the qmfytical solution. The wellbore pressure and flow rate used for deconvolution plots (Pigs. 1 thro.gb 8) are computed from a fully penetrated, finite radkd-resenoti model with the reservoir and fluid properties given in Table 1 in Ref. 7. This is a buifduu test case with C. = l,COO, s=1O, and . . =4,4W P. psi [30 MPa]. Note afso that as it stands, Eq. 10 is underdetennined. One way to make the solution uniaue is by defnim K;- ,=0. This mild conoiA qd the last eledition allows us to elii;ate th; last coiti ment of z giving us an NxN finear system that we can solve with standard algorithms. Another approach is to create an overdetermined system by increasing the size of the tinite elements, with a fine mesh for early-time values and a progr.ewively coarser mesh for later-time values. Doing this is somewhat complicated because

For an evenly spaced partition, consider the Riemann SUMfor@. 1 without the skin term:

K;=gn ( fl \
1500.0 .: :1400,0 : , ~oo,o : S!

Ef-i+lKi - i=,

. . . . . . . . . . . . . . . . . . . ..(17)

1500.0 .2 $ 1400.0-

s : z 1300.0k c ..a m ~ 1200.0-

. .

6 1100.0 1u? , ~, Ime, t, hr

exact + dewnvolved I I&

& 51100.0 ; . . 10 t 000.0 102 10 Tree, t, hr

exact - demnvclved .10 10

Fig. 7Deconvolved pressure for the f low rate with 2% normally dlstribtied error with the reduced-monotone option.

Fig. S-Deconvolved pressure for the flow rate with 2% normally distributed error with the nonposiuve secondderivative option. 55

SPE Formation Bvatuadon, March 19S0

TABLE lCPU dumber of Points 1,000

4,400

TIME AND STORAGE REQUIREMENTS Memory CPU Time Requirement, NM (seconds) 39,255 5,811
39,255

Method Monotone Reduced monotone

KSO

45
20

squares problem that can be solved with an active-set, quadraticprograrmoiog method. 13.14 Much care was needed to formulate the probIem so that it was numerically tractable for large data sets (N> 1,033). The elements of 7 were specifically defined as Ap., K(tI), K(rI), K(12). K(!N-l) so that constraints given by Eqs. 2 through 4 could be Press~ a shpk bounds 011the elements of ~
x~20 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ..(2o)

Monotone Reduced monotone

KSO
11,250 Monotone Reduoed monotone

317,300 23,200 317,300 1,263,800 56,248

570 50 850 .. 122

for iS2, Xiso . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ..(21)

hememoryrequirements(NM) for each methcd(i.e., the amount ,fcomputermemoryrequiredto storethe Matrix,4 and associated mrk arrays) are defined as follows. +6P for Krs O and monotoneoptions NM >(M+ 1 P+3(N+ 1) lnd Nu >2P J +3(N+ 1) +6P for the reduced-monotoneoption, /here M is any integer such that M(M+ I)>N and P-M+4, and /here NM is the number of floating point words required for torage. .Ttmlngesuilsrefora vu 111780 flmtlng r a With
point amlerew. .WaJld not CM@!, m VAX VMS bmause of @xcessivemmm.rj ,equlrmenls.

for i> 3, plus the single linear constraint N-1

K(zN)=K(tl)+At

~ K(ti)zO. i=I

. . ..(23)

Substituting

Q.

16 into Eq. 17 yields

n-1

en-i+lKj.

. . . . . . . . . .. . . . . . ... . . . . . ..(18)

(fl+q)

i=l

Each term in Eq. 18 is divided by (ql +C1), which includes the fmt flow-rate measurement and the associated error. ff the error, El, is small, we may have reasonable results from deconvolution. The real problem, however, is the last summation term in Eq. 18, n-l x i=,

6n-i+lK;.

ThiSterm isa summation of the solution, K, times the error, t. It grows as time increases, which makes deconvolution an unstable process. Algorithms for approximating the solution of F.q. 1 woufdb relatively trivial except for the violent effects small amounts of flowrate noise. have on the solution. Figs. 1 through 3 show us what happens when we use Eq. 10 with equally spaced time points on test data with 100 points and errors in f @ of O, 0.1, and 1%, respectively. Wit3 no error, the discretization gives a solution that is almost indistinguishable from the analytical solution. With only 0.1 % noise, though, our solution starts to oscillate wildly. With 1 % noise, the solution is obviously useless. Several other dmcretiza.tion schemes were tied, including Eqs. 6 and 7, but all exhibit the same phenomenon. We considered three techniques for ~g the effects of measurement noise. The first was to increase the size of our tinite elements as previously discussed so that noise will tend to average itself out. Typical resuks of such smoothing are showo h Figs. 4 and 5. We found that the use of larger elements signitlantly helps our solutions when they are compared with the solutions given in Fig. 2. Unfortumtely, for larger errors and large N, there is usually still some oscillation. The second technique was to restrict the solution space by boposing some or all constmimts given by Eqs. 2 through 5. We =fer to our dgoritbm that folfows constraints given by Eqs. 2 through 4 as the ncmpa$itive secondderivative option. Because our apif &. proximation to K is only twice continuously differ~tile, 5 is used also, it should be replaced by

h general, when Eqs. 20 and 21 are satisfied, Eq. 22 also tends to be satisfied in practice. Consequently, we needed only mild safeguards in the program to prevent it from being violated. @. 22 was handled by adding a penal~ term that increased the objective function if the constraint was violati. The penalty comtant was chosen heuristically on tie basis of the element size of A, and provisions were included for automatically increa.siogit ifit proved too smafl. Treating Eq. 22 io this manner allowed us to we specialized, bigMy efftcient algorithms for quadratic programming problems subject only to simple constraints. The monotonically increasing seconddexivative option MD be meated similarly if the components of Y are deiimxl as Ap,, K(tI), K(tI), K(t2) .K(tl).. K(tN)-K(tN_I). Appendix B presents the details Oi tills tomulation. The use of progressively larger linite elements for larger time. values was also highly signiikmt in making our pmbiem computationafly tractable. We used a mesh that kept the ncdes approximately equally spaced on a logarithmic scale and reduced the size of Eq. 15 fiOm NxN tO NxP, where p is rOugMY fi (tie Precise definition of P is given io TabIe 1). Solving this problem thus involves a singfe or orthogonal factorization of Matrix A (A =QR, with Q being NxNor an ortbonormal matrix and R being an upyr tigular matrix) and the solution of a .PXP quadratic program. For very targe N, solving this problem is still fairly timeconmming even tier the savings given from larger titdte elements. For that reason, a third scheme was considered as a modification to the monotone-increasing swondderivative option. This scheme is called reduced monotone. It simply ignores many of the rows of A for large t where things may not change significantly. Table 1 presents CPU time and storage requirements for three different numbers of data points. It ah shows thattbe reduced-monotone @ion reduces the storage requirements considerably. The last technique for deafing with measurement noise is a slight modification of the previous approaches and can be motivated by a detailed e xandnation of our least-squares fomudation. Eq. 15 is mOst aPPrOPMte if the m=SUrement nOise is ~ g(t) (~ P~ss~ UImmIementS). More precisely, the sOIutiOn of Eq. 15 is ~e exWt solution to the overdetermined system AZ =b, whero b @ the smafkst wssible perturbation in the least-squares sense of b (the vector of pressure measurements). It would, of course, be better to tind the srnallost ~rturbation to f(t) (the flow-rate measurements) for which we cao s?lve the problem exactly. This can be fOrumfated as fOllows. L&f be defined as in Eq. 16 and find K(t) and e(t) that minimize I]@) II subject to

g(t)=@~f(t)+@K(r-r)dr
0

..

...

. . ..(23)

K(!) is monotone increasing for t>O,

..

.. ...

.(19)

and the constraints given by Eqs. 20 through 22. Such a reformuIation of he dwonvoltion ~ ~ ~on~~~~-~ation problem Unfortunately, h minimishoufd give the best result we can expxt z,aticm problem givem by Eq. 23 is compwationafly intractable if Nis Large because it is nonlinear. We have therefore made the following simpfitication. Combining Eqs. 16 and 23 for t=tn gives 8(tn)=&s~(~n)+~(tn)]+~ 0 SPE Fommtion Evaluation,March 1990 *V(7)+e(~)]K(t-r)d~, . . ..(24)

which we refer to as the monotone-increasing secondderivative option. fn either case, Eq. 15 becomes a constrained finear least56

we

4350

3750

if L
%..

1.20 11000
9000
8000 7200 s

..

.. . . .... .. . ... . 10000

-. @a
6000 $ 5000 *

1.oo-

:. .: .. .
:

:%

;+% .

0.80-

8 0,60o exact dernn.olved

. Pressure . Flow Rate

, .:

+...

; .2

% . ..

4000 E 3000 2000 . 10$ 1000 0

..OO

k,

. ..

o.40-r 10+

,V2

# , ~2

, ;.1 to

1d

100 Xme,

10

10

I 100

1 1Ot

hr

Fig. 9Downhole pressure and flow rate of a drawdown test for a horizontal well. which can be rewritten

of deconvolved pressure for the presmm lg. 10-Derivative md flow rate given In Ffg. 9. reduced-monotone option. It can be seen from these figures that the deconvohmion with the monotoneXncreasing K constraints gives the best resuft starting from a very early time. The reducedmonotoneoption, wbife not as accurate, may stUf be acceptable if computational spfed is a factor.

2s

Dlscusslo6f
Ducredzatiom other than the ones listed above were also tried. Versions that did not use the square-root term in the timt interval usually were a bit less accurate for a few timesteps. A Iower+rder discretization, which bound ordy the first derivative to be >0, prodnced results that were not 2.ssm.xdb as we desired. A fdgher+rder discredzation, which was based on cubics and which bamdtbe third derivative to be =0, turned out to be bigfdy sensitive to heuristic details (such as the number of .Tcornponents initially fwed and how they were chosen). The direct deconvolution formula, given by Eqs. 6 and 7, is a recursive algoritlun in which K.+* is computed with all previous computed K values (Kn, K.-~, Kn-2, etc.). h other words, the computed value of Kn is not affected by Km+,. The constmineddeconvoludon afgoritbm mmputes the deamvolved pressure in a nonrecursive manner. As in the zolution of a syztem of linear equations, all computed values of K affect each other; i.e., the accura-

Now recall hat K(O)= m for the constant-rate case. Because it is only a slightly cmde approximation, we can discard the term n E 8.s king $''K'(r),(tn-~)d~ ". . . . . . . . . . . . . . . . . . . . . . . . . . ..(2Q t=z ti_* dominated by the term . . . . . . . . . . . . . . . . . . . . ... . . . . . . ..(27)

~tiKr(r),(tn-r)dr,

and can consider e over tie intervaf [In,tn_l] to be constant and eqwal to en. With these simplifications Eq. 23 becomes h g. [J @sfn+ 0 ,

f(T)K (tnr)d~

=enKI.

. . . ....(28)

cyofKn is affected by Kn+l, Kn+2, K.+3

. and Km.l, Kn_l.

Therefore, we can do the following mani@adon. f-et us fix the vahe of K, to be a constant. Then we minimize Z#=lc~, which is equivzfent to minimizing E:= I $., where

an=gn

[1

Ap,fn+ o

f(~)K,(tn~)dr 1

. ..

. .

. . . . ..(29)

ff we include more terms from W. 26, the mathematical analysis

becomes slightly more complicated. It sdU suggests, however, that we can replace the ~ on problem defined by Eq. 23 by the mucbsimpler Eq. 15 if we freeze Omfirstfewvafuezof XTbese values can be determined by solving our discretization once with use of an N stroll enough so that ourf(t) data are still essentially exact. A large number of test runs were made to determine about how much of T to fieezc before solving the full discretization. The best choice seemed to be to fix the fust three elements. Some typical results for tkse different methods are shown in Figs. 6 tiuough 8. The weflbore flow-rate &ta,f(t), bad a 2% normally distributed error induced into it, and 1,~ points were used. Fig. 6 presents cmmputed dewnvolved and actual (analydcafly computed) welbre pressure drop using the monotone-increasing secondderivative option, while Fig. 8 presents computed deccmvolved weflLmre pressure drop using the nonpaitive secondderivadve option. Fig. 7 presents the computed deconvolved pressure drop using the 8PE Formation Evaluation, March 1990

K,. As a result, we recommend that, in practice, a given data set should be partitioned and that deconvolution should be carried out in a sequential mmner. This is particularly necessay for early time. In weU testing, the idlue.ce function of real systems Wiwdfy can be much more complicated than the Miite did system that is used for the above deconvolution examples (Tigs. 1 through 8). Now let us use the amlyticaf solution for a horizontal welf (Example 3 in Ref. 15) and the flow-rate data (real measurements) shown in F]g. 9 to compute wellbore pressure, also shown in Fig. 9. This horizontal weU example is used because it has few diferent flow psricd.z, as shown in Fig. 5 in Ref. 15. Fig. 10 compares tbe derivative of the deconvolved pressure, computed with the con?.traineddeconvolution method with partition for the pressure and flow-rate data set (Fig. 9), with the derivadve of tie amlyticaf solution. As Fig. 10 shows, the derivative of the deconvolved pressure agrees well with the analytical derivative, with the exception of a few points at early time. Note that ewm at early dine, the mend (shape) of the deconvolved-pressure derivative is almost identkaf to that of the a.nalydcaf derivative.
Conclusions

fn this paper, new deconvolution techniques that am numerical wIution.s of the convolution integml are prezented to compute the influence funcdon (the umst.mt-mto or -pmmure solution) of a system from the measured wellbore pressure and flow rate. It is shown 57

that deconvolution of the welbore flow rate and pressure is shaightforward with direct integration methods if the data afe ffee of measurement errors. However, the direct techniques fail if even smaJl amounts of measurement noise are present. Because of the stabdity problem associated with measurement errors, comtrained-miniiation methods are used for the d~onvolution of measured weJJbore pressure and flow-rate data ifermrs are significant in these measurements; particularly in flow rate. The constraints used here satisfy the general solutions of the diffusiviw equation. It is shown that deconvolution with the monotom?increasing second derivative of the intluence function (deconvolved pressure) yields the best results.

2. Jones, S. C., Katz, D. L., and Tek, M. R.: A Generalized Mc-M for Predicdng the Performance of G8sReservoirs Subjectto WateI Drive, paper SPE 428 presented a! tbe 1962 SPE .&mwd Mcaing, LO, AQ. 8eIes, Oct. 7-10. of Aquifer Jatluenct Fumdons Fmm 3. [email protected], K.H. eraL, s~tion Field Data,, JPT (Dec. 19&) 1417-2* Tmns., MN@ 231. 4. Jarson, J.R. and van P.xdlen, H. K.: Whit Response Function From V=Y@-~1. Da~. JpT (Au8. 1963) 965-6% Tr~., .QME. 234. 5. Bosfic, ]. N., Agarwal, R. G., andCarter,R.D.: <CombinedAnalysis of Postfrachuin8 Performance and Pressure Buildup Data for Evaluating an MHF Gas Well, JPT (Oct. 1980) -,-. 1-10 , ~fl --6. Pascal. H.: Advances in Evafuadm Gas W.U Deliverability Using Variable RateTests Under Non-Darcyhow, paper SFE 9841 p-ti at the 1981 SPE/DOE Low Permeability Sympasium, Denver, May
7.7-?Q.

Nomenclature o ~ = ~1~~~~~f ~trix

~ = ma~ix

7. ~u~, F. and Ayestaran, L.: Analysis Simultaneously of Measured pressurendSandface OWRateinTmmimtWeUTesfing, JPT(Feb. a H
1983) 323-34. 8. Thompson, L.G., Jones, J.R., and Reynolds, A.C.: -Analysisof Fres-

b = vector c, = total compressibdity,

CD = f = g =

k =
K =

L =
N = NM =

p = PO =
AP. = P =

psi 1 &Pa-1] dimensionless welJbore storage input to system system output permeabdity, md influence function lower trianguk mshix number of measurements number of memory requirement pressufe original pfessure, psi KPa] PE5SUE drop Ca.Wd by Skill defined in Table 1

sure Buildup Data Influenced by We!JLmrePhase Redishibution,

SPEFE (Oct. 1986) 435-52 Tmm., AJME, 2S3. 9. fhcmpson, L.G. and Fqnolds, A.C.: AMIYSSof VariableJ?a@WeJl-

q = flow rate Q = cumulative flow rate R = upper trisnguk mshix

~ = $&* t = At = F= fJ = W = a = y = 6 = e = p = u = ~ = time timestep vector of (xl, X2. . .xN)T defined by Eq. B-3 deiinti by 3?+ A-12 constant Comtsnt noise in flow-rate measurement constant oil viscosity, cp [Pa.s] defined by Eq. 29 d-y integration vsriabIe

Test Pressure Data UsinS Dubamds Principle, SPEJ?E (Oct. 1986) 453-69. IO. van Everdingem, A.F. and Hurst, W.: Application of the LapPace Transformationto Flow Fmbkms in R&moirs, Trans., AIME (1949) 186, 305-24. 11. Hamming, R.W.: Numerical Methods for Sciemsts.@ Engineers, McGraw-Hill Book Co., New York City (1973) 375-77. 12. Kuchuk, F.J.: New Methods for Estimating Parameters of LowPermeabiliy Reservoirs, paper SPE 16394 presented at the 1987 SPEI130E SvmKQsium Low Permeabilih Resewoirs, Denvm, MaY on 18-19. 13. GiU, P. E., Murray, W., and Wdght, M.U.: Practical Optitition, Academic Press, London (19gl). 14. .skwart, G.W.: buro.fuctim to Mmix Computations,Academic press, London (1973). 15. Kuchuk, F.J. et UL: Wessure Transient Behavior for Ho&ontal WelJs With Gas Cap or Aquifer,, paper SPE 17413 presented at the 1988 SPE California Regional Meztimg, Jmn8 Beach, March 23-25.

Appendix APlecewIse Linear Approximation

To derive Eq. 8, we MUSTevaluate Wch krm of the summation in Eq. 10 in terms of the discredzd forms off(t) and K(t). Now, f(r) is taken to be piecewise linear, so we have for TE[tj-l,fjl, .f(tn-~)=~n-j+[(~-l -~)/At](~n_j--~n_j-t). . . . . . . . . (A-1)

where fj denotes f($) and At is the dtiferemce between successive time points (assumed constant). For each time intewal T@- 1,$1, j> 1, we afso have K(r) =Kj_, and +[(rfi_I)/&](Kj -K}_I) . . .(A-2)

4 = pormi~ Subscripts

D = dimmsionJess i,j,k,n,m = summation indices

.= N = r = w = @= normahzed refe*ece weJJ, wellmre weJlbore flowing

K~=Kj_l+AtK~-l=Kl Furthermore, Kj=Kj-l+(Af/2)(Kj+Kj_

+At

j1 ~ KY. i=*

.(A-3)

I).

..(A4) 8.$.xt-!-d . + (A-5)

superscripts .: perturbation to flow rate = derivative

In the time intewaJ r.GIO,tl], K is approximated Ap,, which can be expressed as K(r)=Ap,+r(Ki or
KI(T)=Ki+2NKfJ-2N3/2 Tl~Kfr.

+2AfK{)-4A#n~K~

T = tIWpOSe Acknowledgments
We thank SchJumberger for permission to pubJish tfds paper. We dso think Jeffrey Joseph of SchJmnberger WeJJ Services for providing helpful discussions,

. . . . . . . . . . . . . . . (A-6)

Hence, for j> 1, from Eqs. A-1 and A-2, !~lK(,~(tn-r)d,= 2 ~~-l+[(T-V-i)/MI(K}-Kj-l)}

References
1. Hutcb!nsnn,T.S. and Sikma, V.I.: A Gem=dized Water-Drive AmJysis,,, JPT (July 1959) 169-77; Tnm.s., AJME, 216. 58

J
-J-1)}dT. . .. . . (A-7)

X{ Yn-j+[($lr)/~l(fn-jf.

SPE Formation Evaluation, March IWO


For j=l, we have ~-~ (K~+UtKff-Ut327-I/2K,)
0.

and A is the NxN+ all a 12

032 a42 ... %2

1 matrix, defined as
Q13OO0...

@rY,*n-r)d-

a2i a72 a23

a33 U43 ... %

0
0 a45 aN5

0
0 o . aN6 ... . .

a31 041 ... )}dr. .. . ... . ~ (A-Q

034 %4 ... aN4

x{~n-j[(~I

T);~l(~.-j~m-j-l

[ where

aNl

fntegratigEqi.

.~-7 and A78 and substiming At +~ 2 ~fn-1 +~fn_2 ~ )

into Eq. 8 yields -l ~

.1.
: GN,N+ 1

o 0 0

aU =f,; alz=fll%

a,~ = 513fl; a23= 1/3(5f2 +fl) + l/2Wl + l/6f1;

1{2WZ+ l/6f1;

g. Apsfn

Ki + j=2 ~ E

~fn-j-l

am fz; an 1/2(f2+2fl);
a31

f3;

U32 ]/2(f3

+5/3f2 + l/3f,)+

4 + ~fn-j

+ ;f..j+l I

Kj+~
)IN(2

yfo+ ~fl
1),*(

ii.+

f. +fn-1 )

% = - l/3(5f3 +fd + 1/2 WZ+1/6fl; a34= l/2Wi + l/6fl ;

Q41f4; U4Z 1/2(f4 +5/3f3 + l/3fZ) + 1/2W3+ l/6fl; = a43= - l/3(5f4+f3)+l12 W3+ 1/6fl; a~= 1/2 W2+ 116fI;

Ki +A#(fn

+fn-l)K~

>12(2fn

+fn_l)K{.

. . . . . . . . (A-9)

a~~= 112WI + 116f1; ...

Using Eq. A-3, we can then obtain At 8n=&fn+y (

2
;f._l+:fn_2

At )
rl+

n-l
z
j=z

~
( ~fn-j-l

ON1fN; ah2 = l12(fN +5i3ftJ_, + l/3fN-2) + 1/2WN-, + l/6fl ; aim = - l/3(5ffi+fN_,) + 1/2WN_, +1/6f1; aN4= 1/2 N-2 + l/6fl; and W
aNj=l/2wN-j+z +l/6fl
fOr J2t. .................(A-16)

+:fn_j+:fn-j+,
)(

K;+At

~
:;

K; +7
) *(2

~fo+;fl
1 )

Appendix BMonotorre.lncreashrg Sec03rd.0erlw3tlveOption Thevector Z for the monotone-increasingsecond-derivativeoption is expressed in terms of differences in successivesecondderivaiive values as Z=[Ap~,AtKi,A&Ki,A# (Kf-Ki). . .@Kfi_l

X K~+At (

~ K~ +; ::)7

fn+fn_l
)

Ki+A#(fn+fn_I)K~

-~@(2fn+fn_1)Ki'

. . . . . . . . . . . . . . . . . . . . . . .. . . ..(A-lO)

K~-2)]T.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ..(B-l)

To make Eq. A-10 computationafly tractable, we must do sonic further manipulation. Reversing the order of the double summation and ccdlecdng terms in Eq. A-10 yields At g=APsfn+~ 5 1 fn+yfn-l+;f._2+ [ n-l ~

are tbe same as Om,sein Matrix A in Eq. A-9. The remaining elements are defined as
a~~ = -(5/3y1 a~j ==- l/3(5f2 ;

The first two columos of MatiA

+f,)+u] ;
U3; a~=U2;

E ~fi+;t.-l j=2 (

a33 = - U3(5f3 +f2)+ UZ; aN=UI;

a43 = -

l/3(5 f4+fJ+

a43=Ul;

+:$-z (5) )21 )1 +y 1):: 2(2

+~fl+~f:

KIA$

~fn+yfn-1

Ki

am= - l/3(5fN+fN-,)+
a13d0Nj=UN-.j+~ f0rj24,

u&~; aN4=uN-~;
. . . . . . . . . . . . . . . . . . . . . . . . . (B-2)

A$ + 2

n-2 ~

fi+

i=2

[(

K: j=2 ~

:&+:t-l

+:3-2

where Um=

mll Z -Wj+~fl j=l 2

.....................

(B-3)

~fo+ifl

K:.

. . . . . . . . . . . . . ... . . . ..(A-1I) a function W

We can simpli& the above expression by deftig

and where Wj is given by Eq. A-12. The ~vector is the same as that given in I@ A-15. The constraints are then xi 20 for i#3, X3 sO, 33.+ 22, and N-1

30 that

K(tN-l)=K(tI)+ 14
1

E
i-z

(KiK~-I)=O.

@4)

w.= ;
j=2

;fi+;&1+;_&2

> . . . . .. . . . . . . . .. CA-12)

where W1=O. Use Eqs. A-12 and A-13 and put everything nmtzixnotationin Eq. A-n to reach ,4==E . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ..(A-13) . .. At2Kl)T.)T, . . . . . . . ..(A-14)

Both Eqs. 22and B-4 can be handled by adding a psnalty term to the optimization problem; hemce, Eq. 15 can SW be solved by the zmne spwialized zfgorithmz used previously that treat onfy zimple constraints on the idvidual mmponents of X S1 Metric Conversion Factors bbl X 1.589873 psi x 6.894757 E01 = m3 E+OO = kpa

where. Z=(Ap,,A2Ki,@Ki ~=(g1,g2...gN)T,

WEFB
. . . . . . . . . . . . . . . . . . . . .. .. . . ..(A-l5) March 1990
original SPE nmnuwbt mceiwd far rav(w Feb. 13,1985. Paper (SPE 13330] accepted la publlc%trcm S@, 19, 1989, Revlsad mamscr[pt wived April 28.1939.

SPEFonnatim.E valuation,

59