New Regularization Algorithms for Solving the Deconvolution Problem in Well Test Data Interpretation*

Two new regularization algorithms for solving the first-kind Volterra integral equation, which describes the pressure-rate deconvolution problem in well test data interpretation, are developed in this paper. The main features of the problem are the strong nonuniform scale of the solution and large errors (up to 15%) in the input data. In both algorithms, the solution is represented as decomposition on special basic functions, which satisfy given a priori information on solution, and this idea allow us significantly to improve the quality of approximate solution and simplify solving the minimization problem. The theoretical details of the algorithms , as well as the results of numerical experiments for proving robustness of the algorithms, are presented .


Introduction
The pressure-rate deconvolution problem transpires for model-identification processes in pressure-transient well test data interpretation for solving the first-kind Volterra integral equation of convolution type [1] that is given as ]. [0, ] [0, : This equation is a consequence of Duhamel's superposition principle, which holds for a wellbore/reservoir system at certain conditions, namely, when the system is considered as a linear one.The principle states that the measured pressure drop p is the initial reservoir pressure (pressure in undisturbed reservoir), ) (t p is pressure measured in the wellbore, in particular, on the surface or at bottomhole, and T is the total duration of the well test.Functions ) (t g and ) (t p u are unknown and needed to be estimated by a deconvolution algorithm as an inverse problem.
Functions, ( ) u p t , ( ) g t , and such as ( ) = ( ) / ln u tg t dp t d t contain important information about the wellbore/reservoir system.Normally, they are plotted as a function of time as semilog or log-log graphs [2].For instance, an initial part of a log-log plot of ) (t p u versus time yields the wellbore storage coefficient.Reservoir permeability and wellbore skin damage can be calculated from a log-log plot of ) (t tg versus time or a semilog plot of ) (t p u during the middle time period.Identification of certain geological features such as faults, fractures, reservoir boundaries, etc., is possible from middle and late time periods.
Peculiarities of the deconvolution problem in well test interpretation are in the following: 1) real input data can be discontinuous and oscillating (see data 2); 2) a solution ) (t g is rather multi-scale (see data 1); 3) these specific features of the problem cause serious difficulties for reconstruction of a desired solution and its interpretation.
As is well known [3], reservoir response ) (t g (the solution of the diffusivity equation) satisfies the following a priori given constraints: ( 1) ( ) / 0, = 0,1,2, It is worth noting that during development of the algorithms using a priori constraints (2), one usually employs not more than three following conditions [4]: Let us briefly review a few methods presented in the petroleum literature for solving a pressure-rate deconvolution problem in linear setting (1).Linear programming methods containing a priori constraints (3) were developed by Coats et al. [3] and Gajdica et al. [5].Different rules for discretization of the convolution integral with different modifications, which resulted in explicit recursive formulas, were proposed in a number of works.Hutchinson and Sikora [6] used a trial and error procedure with correcting multipliers, Katz et al. [7] used a linear recursive formula and tried considering a priori information (3), Jargon and van Poollen [8] added logarithmic interpolation to a recursive formula, Thompson and Reynolds [9] proposed three methods based on different discretization rules, Bostic et al. [10] and Kuchuk [11] simply used pure recursive formulas without any modifications.In the work by Kucuk and Ayestaran [12] two methodologies for deconvolution were proposed, namely, Laplace transform deconvolution and curve-fit approximations.A linear least-squares approach with a priori constraints (3) was proposed by Kuchuk et al. [4].Baygun et al. [13] also used the linear least-squares method with a special type of inequalities called the normalized autocorrelation constraints.Spectral methods, based on a convolution theorem, which in the case of Equation ( 1) holds approximately for both Laplace and Fourier transforms, were developed by Roumboutsos and Stewart [14], Mendes et al. [15], Bourgeois and Horne [16], Onur and Reynolds [17].
Unfortunately, almost all methods using recursive formulas failed on pressure data with noise up to 1%.Linear programming methods and linear least-squares approaches worked insufficiently when errors in flow rate data were 1 to 2%.Spectral methods are the simplest, but they mostly failed with noisy data.
In the early 2000s, von Schroeter et al. [18,19] presented a completely new methodology for the pressurerate deconvolution inverse problem.First, they transformed the linear Equation ( 1) to a nonlinear one This replacement makes the condition 0 ) (  t g to hold automatically.Second, the total nonlinear leastsquares method was used with the Tikhonov regularization in which a constraint was imposed on the curvature of the solution.This methodology not only estimates the solution ) (t g but also allows for correcting the flow rates ) (t q and the initial reservoir pressure 0 p if they contain large errors.But Levitan [20] showed that the pressure-rate deconvolution algorithm fails when inconsistent input data are used.He proposed a modified methodology.In a subsequent paper by Levitan et al. [21], valuable suggestions were made for practical applications of the pressure-rate deconvolution.Pimonov et al. [22,23] further modified the deconvolution algorithm of von Schroeter et al. [18,19] and Levitan [20].The modification was done by introducing a more general objective function in the Levenberg-Marquardt minimization procedure.This algorithm was successfully applied by Pimonov et al. [22,23] to deconvolution of drillstem test data and wireline formation tester data.Recently, a totally different approach for solving the nonlinear Equation (1) was proposed by Ageev et al. [24], where the authors used an iteratively regularized modified Levenberg-Marquardt method.
In parallel to von Schroeter et al. [18,19] and Levitan [20] methodologies, research has been continued to solve the linear convolution integral equation (1) directly.For instance, a projection method based on B-splines representation of the solution was proposed by Ilk et al. [25].The method appears to be quite robust, although it does not take into account errors in flow rates and the initial reservoir pressure.Andrecut and Madni [26] and Andrecut [27] studied three methods: 1) Tikhonov regularization, 2) Krylov subspace projection method, and 3) stochastic Monte-Carlo.The authors stated that the stochastic Monte-Carlo method was more robust than the other two.Overall disadvantages of these three methods come from the fact that flow rate errors are not accounted for.
Finally, Cinar et al. [28] compared efficiency and robustness of pressure-rate deconvolution algorithms of von Schroeter et al. [18,19], Levitan [20], and Ilk et al. [25].First, they stated that the regularization approach used by Ilk et al. [25] seems to be weaker; hence, this makes the algorithm less tolerant to the errors in pressure and rate data compared to the von Schroeter et al. [18], [19] and Levitan [20] algorithms.Second, these last two algorithms have certain efficiency and robustness.
This work is a continuation of a study done by Vasin et al. [29,30].We now propose two algorithms for solving the linear convolution integral Equation (1) and pre-sent a theoretical analysis of the algorithms for pressurerate deconvolution.
The first algorithm proposed is based on a quasi-solutions method, i.e., an initial problem (1) is reduced to a minimization problem   2 2 min : on a compact set Q (see ( 6)), which is a closure of the set Q .It is very important that a solution of the minimization problem is found as an expansion by a special basis, which elements satisfy a priori constraints.
It is worth noting that an attempt to develop a deconvolution algorithm based on discrepancy minimization has been done by Kuchuk et al. [4].However, an unsuitable minimization procedure was used, and it led to unstable deconvolved results.The minimization algorithm proposed in our work appears to be more effective.At the same time, the problem with an initial approximation for the solution is successfully resolved.Due to using a special basis, we have possibility to choose always 0 = ) (t g as an initial point (guess) for iterative algorithms used in the numerical minimization.Usually special procedures are required to find an appropriate starting point for iterative processes of minimization (see papers mentioned above).
The second algorithm uses the Tikhonov regularization method with constraints (3), i.e., an approximate solution is calculated by regularized conditional minimization where a stabilizing functional  is nonlinear (nonquadratic), which is opposite to traditional scheme (see [31]).This feature of the functional  significantly improves quality of an approximate solution.The functional is defined explicitly further in Section 4, contained detailed description of the second algorithm.To deconvolve data with a large measurement errors in flow rates , ) (t q a modified version of the algorithm is also developed.In the modification, flow rates at the estimated solution ) (t g are additionally corrected using quadratic Tikhonov regularization.
Our paper is organized as follows: investigation of the algorithm based on the quasi-solutions method is given in Section 2; stability of discrete approximations based on collocation method is proven in Section 3; Section 4 is devoted to the study of the Tikhonov regularization algorithm; Section 5 presents the details of the scheme for flow rates ) (t q correction; synthetic well test examples are presented in Section 6 to illustrate the methodology.

Quasi-Solutions Method and Its Features
Before studying the quasi-solutions method, several assumptions are made on the input data.Let ( ) q t be a continuous function and (0) 0.
q  Let us assume that an integral operator A in (1 On this pair of spaces (as well as on the spaces of continuous functions), operator A is a compact one (completely continuous).Thus, A does not have a continuous inverse operator , A so the problem is essentially ill-posed and the solution is strongly sensitive to errors in the input data.
Further, we will consider the set Q of a priori constraints, which is wider than the initial set , Q namely, (in contrast to the set , Q which is nonclosed), so a basis with non-smooth (e.g., piecewise linear) functions could be used.
Real data ( ) q t and ( ) with the following estimations of measurement errors: , .
Equation ( 1) might not have a solution for a given right because the actual range of A is nonclosed.Thus, in a general case under the solution of problem (1), we will understand that it is a quasi-solution (see [32], p. 36); i.e., a function ˆ, g which realizes minimum of discrepancy: on a compact set Q (see lemma 2.2).Equation ( 8) is solvable for any function because the set Q is a compact and the operator A is continuous.
Solvability conditions for (1) and ( 8) are given in the following lemma.
Proof.By differentiating both parts of Equation (1), we obtain the Volterra equation of the second kind which has a unique solution due to generalized principle of contracting mappings [33].As long as any solution of Equation ( 1) is a solution of Equation ( 9), then (1) cannot have several solutions; therefore, ( ) 0 Ker A  .Since an objective functional of ( 8) is strictly convex at these conditions, a quasi-solution is also unique.
Further, only the case for {0} = ) ( ker A is considered.The general situation is discussed at the end of the section (see Theorem 2.2).
L T Proof is similar to one for the set Q of monotone functions (see [34]).
Let us consider problem (8) with approximate data that satisfy relations (7): where A and the approximation conditions (7) hold.Then problem (10) is uniquely solvable and the following convergence holds: where g ˆ is the solution of problem (8).
In general case, when , ) ( ker A ) ( ker  A are not trivial sets, the problems (8), (10) are non-uniquely solvable.Then the following theorem holds.

Ĝ
be the sets of their corresponding solutions. Then, The proofs of Therorems 2.1, 2.2 can be obtain on the basis of the general theorem on the quasi-solution method for operator equations [32].
Thus, quasi-solutions are stable with respect to errors in the input data and } { ,  g generates a regularized family of approximate solutions.In turn, the quasi-solutions method generates a regularizing algorithm [31].

Construction and Study of a Basis
For a discrete approximation of the equation 0 ( ) ( ) = ( ), 0 , a collocation method is used.Let In this case, integral Equation ( 12) can be approximated by a system of algebraic linear equations T t i  i t are colloca- tion points (time points of pressure ) (t p measurements in the well), , G : which is generated by 1  n elements of the basis, and n R is n-dimensional Euclidean space with the norm , = ).
The basic functions are built with connection to the grid on which the solution is searched by the following rule: are non-negative, non-increasing, and convex functions; i.e., they satisfy the same a priori constraints as the desired solution.This property of the basis is very important, since allows us to get a stable approximation of the desired solution for Equation (1) with real well tests.
For , k we have a formula for calculation of the coefficients where the functions ) ( k g are calculated by recurrent relations

Convergence of Discrete Approximations
To complete the investigation of discretization for problem (10), we approximate the set of a priori constraints Q by a sequence of the sets where the mapping n A is defined by Formula (13), .
The properties (1-3) are established in the same way as for the case of the integral Fredholm operator (see [39]).Since , ), which strongly converges to .g Because Reference to the work by Vasin [35] completes the proof of the theorem.
It is easy to see that problem ( 18) is equivalent to the problem  19) is resolved by the conditional gradient method (see [41]).

Tikhonov Regularization Method with
Non-Linear Stabilizer

Convergence of Regularized Solutions
Numerical processing synthetic data with relatively large errors in flow rates ( ) q t and pressures ( ) p t have shown that the quasi-solutions method did not give a smooth solution, which would be appropriate for interpretation.Therefore, it is necessary to add regularization by a non-linear (non-quadratic) stabilizer  and a priori constraints: where i.e., in contrast to the set , Q we additionally introduced the constraint from below on the function .) (t g From solution of the initial problem (1) it is seen that this constraint could be considered as true and we do not lose generality.
Further, we explain the reasons for the current type of the stabilizer .] [g  It is well known [18,19] that at transition from a linear Equation ( 1) to a nonlinear statement ( 4), the replacements of the variable are used.Now, if in a nonlinear deconvolution statement the stabilizing functional is taken in the form and an inverse replacement is done, then the functional is obtained.Eliminating less significant terms with the first derivative, we obtain stabilizer (21).
It is worth noting that using a traditional norm in Sobolev space ] [0, Proof.Let us denote as ] [g  the objective functional of problem (20).Let n g be a minimizing sequence and .
is a weakly compact one, so there exists a weakly converging sub- Due to the embedding theorems (see [42,43]), existing of such weakly converging subsequence implies uniform convergence This is a consequence of the following relations: The first term in the right-hand-side of the relation tends to zero due to both boundedness of weakly converging sequence and uniform convergence (24); the second term tends to zero due to weak convergence of For this function, the following inequalities hold: and they give an estimation  [42,43], we obtain (23).Remark.After discretization of problem (20) on the base of the scheme described in secton 3, the conjugate gradients method (see [41]) is used.obtained by the Tikhonov regularization algorithm, which makes interpretation of the solution on the log-log derivative plot almost impossible.In this case, an additional correction of the function ) (t q is suggested using the following procedure:

Correction of the Flow Rate Function
1) solve problem (20) 3) use the Tikhonov method (see [31]) for an approximate calculation of ) (t q from Equation (25) in the form where the regularization parameter     is chosen from the relation     ; : q is a noisy flow rate,  is its relative error lev- el, and is the solution of (26); 4) resolve problem (20) using the calculated ; ) ( ~t q 5) repeat the steps (1-4) several times; the calculated ) (t g is considered as an approximate solution of an initial problem (1).
This procedure could be implemented also for the quasi-solutions algorithm.Numerical experiments with the second algorithm (Tikhonov regularization) showed that after several iterations on g and , q the function ) (t q was corrected and the quality of the solution was improved.

Numerical Experiments
The algorithms were tested on two synthetic data (see Figures 1 and 6 For the second set of synthetic data (Figure 6), experiments were performed at noisy ) (t p with a 1% error and exact , ) (t q as well as at ) (t p with a 1% error and ) (t q with a 15% error.In Figures 2, 3, 4, and 5 results of calculations for synthetic data 1 (Figure 1) with pressure errors to 0.5% and 5% ( ) (t q is exact) are presented on a log-log scale.Despite the fact that these data contain small amounts of information because 0 = ) (t q over a large part of the time period, the solution is reconstructed with appropriate accuracy. Exact and numerical solutions for synthetic data 2 with 1% errors in pressure and exact ) (t q as well as          The relative error of ) (t g is 10.7% and 4.45% for ( ) with 1% errors in pressure and 15% maximum error in flow rates are shown in Figures 7, 8, 9 and 10.It is obvious that the solution is reconstructed with sufficient accuracy when there are no errors in flow rates; however, significant distortions are visible at 15% error in flow rates.It is necessary to correct flow rates according to the procedure described in Section 5.After 5 to 15 iterations of this correction, the quality of the solution is improved significantly (compare Figure 10 and Figure 11).It is worth noting that the Tikhonov method gives a more smooth solution for the calculations than with the quasisolutions method, which makes easier interpretation.

Conclusions
Two regularization methods for solving the deconvolution problem are proposed and all stages of the algo-rithms are theoretically established.The effectiveness of the methods is achieved by using all of the a priori information about the solution, choosing the non-quadratic stabilizing functional, and using a special basis for approximation of the solution.
Numerical experiments performed on synthetic data with error levels to 5% for pressures ) (t p and to 15% for flow rates ) (t q show effectiveness of the algorithms we have developed. The advantage of the algorithms is that in iterative minimization processes the initial approximation can be easily chosen, namely as a zero function.The quasi-solutions method allows for avoiding the cumbersome procedure for selecting regularization parameters.Moreover, the special basis for establishing approximate solutions makes an a priori constrain (6) to be true, and this property significantly improves the quality of ap-proximate solution.It is worth noting that in our work, a theoretical foundation for all stages of the developed algorithms is performed, which contrasts works by others on the pressure-rate deconvolution problem, known to us and referenced in the literature.
the algorithm from theorem 3. A discrete analogue of the minimization problem(10) has the form

Figure 1 .
Figure 1.Synthetic data 1 for pressures and flow rates.

Figure 2 .
Figure 2. Solution by the quasi-solutions method without flow rates correction at 0.5% error in ) (t p and exact ) (t q .

Figure 3 .
Figure 3. Solution by the Tikhonov regularization method without flow rate correction at 0.5% error in ) (t p and exact ) (t q .

Figure 4 . 5 .
Figure 4. Solution by the quasi-solutions method without flow rate correction at 5% error in ) (t p and exact ) (t q .

Figure 6 .
Figure 6.Synthetic data 2 for pressures and flow rates.

Figure 7 .
Figure 7. Solution by the quasi-solutions method without flow rate correction at 1% error in ) (t p and exact ) (t q .

Figure 8 .
Figure 8. Solution by the Tikhonov regularization method without flow rate correction at 1% error in ) (t p and exact ) (t q .

Figure 9 .
Figure 9. Solution by the quasi-solutions method without flow rate correction at 1% error in ) (t p and 15% error in ) (t q .

Figure 10 .
Figure 10.Solution by the Tikhonov regularization without flow rate correction at 1% error in ) (t p and 15% error in ) (t q .

Figure 11 .
Figure 11.Solution by the Tikhonov regularization with flow rate correction at 1% error in ) (t p and 15% error in ) (t q .
for noisy ,