Well Test Analysis of Naturally Fractured Vuggy Reservoirs with an Analytical Triple Porosity – Double Permeability Model and a Global Optimization Method

The aim of this work is to study the automatic characterization of Naturally Fractured Vuggy Reservoirs via well test analysis, using a triple porosity-dual permeability model. The inter-porosity flow parameters, the storativity ratios, as well as the permeability ratio, the wellbore storage effect, the skin and the total permeability will be identified as parameters of the model. In this work, we will perform the well test interpretation in Laplace space, using numerical algorithms to transfer the discrete real data given in fully dimensional time to Laplace space. The well test interpretation problem in Laplace space has been posed as a nonlinear least squares optimization problem with box constraints and a linear inequality constraint, which is usually solved using local Newton type methods with a trust region. However, local methods as the one used in our work called TRON or the well-known Levenberg-Marquardt method, are often not able to find an optimal solution with a good fit of the data. Also well test analysis with the triple porosity-double permeability model, like most inverse problems, can yield multiple solutions with good match to the data. To deal with these specific characteristics, we will use a global optimization algorithm called the Tunneling Method (TM). In the design of the algorithm, we take into account issues of the problem like the fact that the parameter estimation has to be done with high precision, the presence of noise in the measurements and the need to solve the problem computationally fast. We demonstrate that the use of the TM in this study, showed to be an efficient and robust alternative to solve the well test characterization, as several optimal solutions, with very good match to the data were obtained.


* Corresponding author
Re´sume´-Analyse des puits d'essai de re´servoirs vacuolaires naturellement fracture´s avec un mode`le de triple porosite´-double perme´abilite´et une me´thode d'optimisation globale -L'objectif de ce travail est l'e´tude de la caracte´risation automatique de re´servoirs de pe´trole fracture´s ve´siculaires via l'analyse des puits d'essai, avec un mode`le de triple porosite´et de perme´abilite´double. Les parame`tres qui doivent eˆtre identifie´s comme ceux du mode`le, sont l'e´coulement entre les trois media, les ratios de stockage, la porosite´, le ratio de perme´abilite´, l'effet de stockage, la perme´abilite´de peau et totale. Dans ce travail, nous avons effectueĺ 'interpre´tation de l'essai dans l'espace de Laplace, en utilisant des algorithmes nume´riques pour transfe´rer les donne´es obtenues dans un espace de temps complet en un espace de Laplace. Le proble`me de l'interpre´tation des tests dans l'espace de Laplace est pose´comme un proble`me d'optimisation des moindres carre´s non line´aires avec boıˆte contrainte et ine´galite´line´aire contrainte, ce qui est ge´ne´ralement re´solu en utilisant des me´thodes de type Newton avec une re´gion de confiance. Cependant, les me´thodes locales comme celle utilise´e dans notre travail appele´e TRON, ou la me´thode bien connue de Levenberg-Marquardt, ne sont pas souvent en mesure de trouver une solution optimale avec un bon ajustement des donne´es. É galement la caracte´risation automatique des re´servoirs du pe´trole fracture´s vacuolaires via l'analyse des puits d'essai avec le mode`le de porosite´-perme´abilite´a`double triple, peut, comme la plupart des proble`mes inverses, donner des solutions multiples avec une bonne corre´lation des donne´es. Pour faire face a`ces spe´cificite´s, nous avons utilise´une approche avec la me´thode d'optimisation globale appele´e Tunneling Method (TM). Dans l'adaptation de l'algorithme, nous prenons en compte des proble`mes comme le fait que le l'estimation des parame`tres doit eˆtre formule´e avec une grande pre´cision, comme la pre´sence de bruit dans les mesures et comme la ne´cessite´de re´soudre le proble`me de calcul rapidement. Nous de´montrons dans cette e´tude que l'utilisation de TM est une alternative efficace et robuste pour re´soudre la caracte´risation des puits d'essai, du fait que plusieurs solutions avec un bon ajustement aux donne´es aient e´te´obtenues.

INTRODUCTION
During the last decades, well testing has advanced to one of the most powerful tools for determining complex reservoir characteristics. The causes of this progress are the development of high accuracy pressure measurement devices, powerful computers and new interpretation methods. More accurate data and more computer power for analysis of these data are available. New models that incorporate more insight on the pressure behavior in the reservoir are also established.
Automated well-test analysis is today a commonly used tool for formation evaluation. This procedure employs numerical methods to solve the inverse problem of estimating reservoir parameters from the analysis of pressure and rate data given an appropriate flow model. The inverse problem is usually stated in the form of a Nonlinear Least Squares (NLS) minimization problem or a least absolute value problem, see for instance [1][2][3][4][5][6]. This procedure is typical for the nonlinear regression field.
The study of flow processes in naturally fractured vuggy reservoirs has recently received considerable attention because a number of such reservoirs have been found worldwide with significant oil and gas production and reserves. Triple-porosity models for such reservoirs have been proposed by different authors [7][8][9][10][11][12][13] as extensions to the double-porosity model of Warren and Root [14]. Among these models the one presented by Camacho et al. in [10] is the only one that considers primary flow in both networks of fractures and vugs, and it is a triple porosity-dual permeability model. There is evidence that in some vuggy carbonate rocks the vugs are touching forming a network of directly connected vugs and creating flow pathways unrestricted by matrix resistance, see for instance [15][16][17][18]. The triple porosity-dual permeability model of Camacho et al. [10] depends on two storativity ratios, one for the fractures and one for the vugs, and three inter-porosity flow parameters which take into account fluid transfer between matrix, fractures, and vugs. It also depends on the permeability ratio between fractures and vugs.
The triple porosity-dual permeability model consists on a system of partial differential equations that can be solved analytically to obtain the pressures in the three considered media (fractures, matrix, vugs) and also the wellbore pressure. The analytical solutions are obtained in Laplace space with the aid of the Laplace transformation. The Stehfest method [19] is then applied to return the solutions to real time. The nonlinear regression problem for the model of Camacho et al. [10] amounts to the search of the model parameters that performs the best fit, of the pressure curve. To check the goodness of the fit, we also compute the pressure logarithmic derivative curve. As stated above the search for the best fit parameters is formulated as a NLS minimization problem.
Because any minimization method will require the solution of the model many times, the use of Stehfest algorithm to return each solution in Laplace space to real time may cause the optimization to be unacceptably slow, and due to the ill-posedness of this inverse procedure, it can produce pressure curves showing oscillatory behavior. To avoid this inconvenience, some authors proposed to conduct the parameter identification process in Laplace space, see [20][21][22]. We will then perform the parameter identification in Laplace space, using the method proposed by Bourgeois and Horne [21] for the transformation of pressure data to Laplace space.
The inverse parameter estimation problem can be illposed in the Hadamard sense [23] because several parameters sets may fit the data sufficiently well to a certain precision. Furthermore in some flux models, different values of the parameters of the model could produce an identical analytical solution. This implies that for the solution of the corresponding minimization problem, robust global optimization methods, able to get several global solutions with good match should be employed. It is known that typical gradient-based Newton type methods, for example the commonly employed Levenberg-Marquardt method [24,25], are local methods in the sense that they are only able to find one local solution. Gradient based local methods may be even forced to stop even at non optimal points when the objective function has sharp or flat valleys.
In this work, we will employ the Tunneling Method (TM) of Go´mez and Levy [26], Levy and Montalvo [27] and Levy and Go´mez [28] for the solution of the NLS fitting problem. The TM is an optimization method capable to find several global optimal solutions for objective functions with several local minima. Also, it is able to get the optimal solution with the required precision (which is an important issue in our application) even in the presence of flat or sharp valleys. The TM version we are using in this work is gradient based and therefore the derivatives of the objective function with respect to the parameters are needed. The development of the parameter identification process in Laplace space allows the development of an analytic gradient.
We will proceed as follows: in Section 1, we describe the triple porosity-double permeability model and present its solution for certain initial and boundary conditions. In Section 2, we introduce the well test analysis problem in the form of a nonlinear least squares optimization problem and the procedure for transformation of real pressure data to Laplace space. In Section 3, we describe the TM for global optimization. In Section 4, we present and discuss some results for synthetic examples and real field cases exhibiting one or multiple minima with good match. Finally, some conclusions are provided.

THE TRIPLE POROSITY MODEL
In [10], Camacho et al. introduced a model for the simulation of high secondary porosity, mainly vuggy porosity, in naturally fractured reservoirs. They proposed a pressure-transient triple porosity model based on the pseudo-steady inter-porosity flow approximation, i.e. they considered the fluid transfer between matrix, vugs, and fractures systems directly proportional to the difference in the volume average macroscopic pressure of matrix, vugs, and fractures. Besides considering the interaction between these 3 media, the model includes the possibility of having primary flow through the system of vugs in addition of having flow through the fractures. It is a triple porosity-dual permeability model that includes as special cases triple porosity-single permeability models as the ones presented in Abdassah and Ershaghi [7], Liu et al. [8], and Wu et al. [9,11].
In developing the model following assumptions are made: -rock properties are constant in each medium; -the reservoir is of uniform thickness with impermeable lower and upper boundaries; -the fluid is slightly compressible with constant viscosity and the fluid flow is single-phase and isothermal; -fluid flow into the wellbore is radial and both fractures and vugs feed the well. Using dimensionless variables the differential equation for the pressure on the fractures is the following: For the pressure on the matrix blocks, assuming that the permeability of the matrix is negligible compared to that of the fractures and vugs system, the equation is: The pressure on the vugs is governed by the following equation: The dimensionless variables are given by: where p Dj is the pressure in the medium j using index j ¼ f for fractures, j ¼ v for vugs or j ¼ m for matrix and p 0 is the value of the initial pressure considered as equal in all media.
The definitions of the model parameters are given by the following relations considering that k T ¼ k f þ k v is the total permeability: where: The parameter k ij is the inter-porosity flow shape factor between medium i and medium j, j is the permeability ratio and x j are the storativity ratios. Note that in the definitions of k mf and k mv , we have used k m because, in the absence of capillary forces, we expect that under primary production conditions, the fluid goes from matrix to vugs and fracture networks. For the case of the triple porosity-single permeability model, i.e. when there is only primary flow through the fractures or through the vugs network, the vugs and fractures permeabilities, respectively, in the above definitions are set equal to zero except in the numerator of k vf . For these cases, j ¼ 1 and 0, respectively. Thus, the parameter j takes values between zero and one.
Considering the presence of mechanical skin for both fractures and vugs and the wellbore storage effect, new parameters s f ; s v and C D are introduced to the model. The parameters s f ; s v are related to the mechanical skin for fractures and vugs, respectively, while C D is the dimensionless wellbore storage coefficient.
Taking into account the skin, the following boundary and initial conditions must be satisfied: Considering the wellbore storage effect, the dimensionless wellbore pressure in Laplace space in the well has the form: where u is the Laplace variable and: The analytical solution of Equations (1-3) for the boundary and initial conditions in (6a-d) can be obtained in Laplace space, through the use of the Laplace Transform. Considering the first of the equalities of (8), we obtain: where x ¼ ðj; x f ; x v ; k mf ; k mv ; k vf ; s f ; s v Þ is the parameter's vector that describes the model's dependence on some of the parameters. The equations for K 1 ðx; uÞ, K 2 ðx; uÞ, K 11 ðx; uÞ, K 12 ðx; uÞ, K 21 ðx; uÞ and K 22 ðx; uÞ are given in Appendix A.

Special Ill-Posed Cases of the Triple Porosity Model where Different Values of the Parameters Result in an Identical Analytical Solution
The parameter j in (1) and (3) is assumed to be in the open interval 0; 1 ð Þ and that reflects the issue of primary flow in the networks of fractures and of vugs. When the primary flow is restricted to either the network of fractures or the network of vugs, the model of Camacho et al. [10] can also be used to describe the transient pressure behavior in the reservoir. In these cases, the parameter j is set to the value of one for primary flow in the network of fractures or to zero for primary flow in the network of vugs. This approach is followed in the papers of Camacho et al. [10], Liu et al. [8], and Wu et al. [9] for primary flow in the network of fractures.
When j ¼ 1, Equation (3) becomes a differential equation similar to Equation (2). Correspondingly when j ¼ 0, Equation (1) is also a differential equation similar to Equation (2). Let us then set j ¼ 1. Since the flow into the wellbore is only through the network of fractures the boundary conditions are only needed for p Df . We consider the following boundary conditions: that correspond to an infinite reservoir and the constant flow rate condition respectively. The initial conditions are the same as (6d). The solution p wD , considering wellbore storage C D and mechanical skin s (see Eq. A-17 of Camacho et al. [10]), is: where: Note that in this case, since j ¼ 1, only one skin is considered and the wellbore storage coefficient appears in (12), the vector of parameters x is then of the form The function g x; u ð Þ present in (13) and (14) is a rational function of u given by: where: From Equations (12)(13)(14)(15)(16)(17)(18)(19)(20), it is evident that the model dependence on parameters x f , x v , k mf , k mv and k vf is only given in the function g x; u ð Þ through the coefficients of this rational function. In this situation, the model stated by Equations (1-3) where j ¼ 1, boundary conditions (10,11), initial conditions (6d) and solution (12) is illposed. The ill-posedness in this case derives from the fact that for different sets of parameters the function g x; u ð Þ remain the same. In other words, different sets of parameters when substituted in Equations (16)(17)(18)(19)(20) result in the same set of coefficients a 2 , a 1 , a 0 , b 2 and b 1 . Therefore the model dependence on the parameters is not unique. As an example, we give in Table 1 two different sets of parameters x f , x v , k mf , k mv and k vf that, when substituted in (16)(17)(18)(19)(20), always result in the following set of coefficients: and therefore producing an identical g x; u ð Þ. The procedure to obtain the values shown in Table 1 is described in Appendix B. There can exist up to four different sets of parameters resulting in the same set of coefficients. The non-uniqueness on the dependence of the parameters of the triple porosity model, when primary flow is only taking place in one of the networks of fractures or vugs, allows us to assert that in this case the use of the model for formation evaluation is not very well suited. To determine all the parameters involved in the solution from pressure data only, information about x f and x v from well logs and cores is required to eliminate the non-uniqueness.
This behavior is not evident when j 2 0; 1 ð Þ as we haven't been able so far to determine the analytical well-posedness of the triple porosity model in this case and further research is needed. Nevertheless when j 2 0; 1 ð Þ other sources of non-uniqueness arose in the inverse parameter estimation problem in the form of a NLS minimization problem as it will be presented in sections below. For example, noisy measurements or the inherent limitations of double precision arithmetic can contribute to the existence of multiple minima of the NLS minimization problem.

CHARACTERIZATION IN LAPLACE SPACE
To reduce the number of parameters of the model, we assume that the mechanical skin in fractures s f and vugs s v are the same and equal to s. Therefore the aim of this work is to find the dimensionless parameters k T , j, x f , x v , k mf , k mv , k vf , C D and s that, within a given accuracy, better fit the dimensionless wellbore pressure as given by Equations (7) and (9) to pressure data obtained from a well test and its log-derivative curve.
Assume we have a set of wellbore pressure measurements D t ¼ ðt i ; p i Þ f g i¼1;...;N obtained from a pressure log during a well test. The measurements are fully dimensional data, therefore to perform the well test parameter identification process in Laplace space we need to transform the fully dimensional data to Laplace space. This transformation is accomplished according to the algorithm presented by Bourgeois and Horne in [21]. Also according to Equation (8) on page 17 of [21], we have the relation: where Áp w ðu f Þ is the Laplace transform of the fully dimensional wellbore pressure drop, u f denotes the fully dimensional Laplace variable and p wD ðuÞ is given by Equation (17). Let us denote: we assume that the value of a w is known and from (21), we get: From Equations (9) and (7) it is evident that p wD ðuÞ depends on parameters j, x f , x v , k mf , k mv , k vf , C D and s. Then (22) defines also the dependence on k T : The final relation between the Laplace transform of the fully dimensional wellbore pressure drop and the triple porositydouble permeability model solution is complete.
We denote with: the vector of unknown parameters and with: the model solution. We also compute the Laplace transform of the fully dimensional pressure drop data ..;M of points in the Laplace space according to the algorithm of Bourgeois and Horne. Through this process we get a data setD u ¼ ðu i ;p i Þ f g i¼1;...;M . The best fit x Ã is then given as the solution to the nonlinear least squares optimization problem: The inequality constraint derives from the fact that we set the correspondence: and in the triple porosity model, the parameters x f and x v fulfill the relation x f þ x v < 1. Also x min and x max are lower and upper bounds for the parameters in x. The minimization problem (25) is then a NLS problem with box constraints and one linear inequality constraint. As we will describe next, we will use gradient-based optimization methods to solve (25). Note that the optimization problem is stated in Laplace space, therefore the time consuming Laplace transform inversion of the triple porosity-double permeability model solution by the Stehfest method at each step of the optimization method is completely avoided.

THE TUNNELING GLOBAL OPTIMIZATION METHOD
The problem defined by (25) can be considered as an inverse problem, indeed according to Oliver et al. [29]. Inverse problems that deal with inexact data, as pressure measurements, are almost never well-posed. In our case, the ill-posedness derives from the fact that we can only expect to find the exact solution of (25) for synthetic data. The exact solution x Ã of (25) is the one that realizes Eðx Ã Þ ¼ 0 for every u in dimensionless Laplace space. Obviously, we cannot guarantee the fulfillment of this condition since we are forced to deal with a finite number M of data points, the inherent noise of pressure measurements and limited arithmetic precision. Therefore, we expect (25) to have several optimal solutions in the sense that they fit the data sufficiently well and satisfy the optimality conditions. We will call these solutions as global solutions and denote them by x Ã G i . On the other hand, the optimization problem (25) can also present local minima. A local minimum of (25) is any x Ã L i that fulfills the constraints and the optimality conditions but such that they does not fit the data sufficiently well. Any deterministic optimization method employed for the solution of (25) will proceed iteratively starting from an initial point and try to construct a sequence x k that converges to one of the x Ã G i . Deterministic optimization methods can be classified as global or local methods. Global methods are designed to stop when the value of Eðx k Þ is sufficiently close to the value of one of the global optima x Ã G i and therefore, Eðx k Þ e for e sufficiently small. Local methods can only guarantee the convergence of x k to a local minima and therefore can get stuck in a x k such that Eðx k Þ > e.
The considerations above made explicit the need for a robust global optimization method, able to get several global optima for the solution of (25). Despite several works involving the use of global optimization methods (see for example Buitrago et al. . Actually, the Levenberg-Marquardt method can be considered as a standard for the solution of nonlinear least squares problems and in general it significantly outperforms gradient descent and conjugate gradient methods for medium sized problems, see for instance Madsen et al. [33]. Nevertheless, as similar gradient-based methods, Levenberg-Marquardt is a local method. In order to get several solutions, the algorithm should be run several times and it may converge to the same local solution many times. The shape of the objective function in (25) is particularly complex as it may have very flat or sharp valleys. Gradient based local methods may be forced in this case to stop even at non optimal points. As it will be shown in the numerical results section this happens in a large number of cases.
In this work, we will use the TM of Go´mez and Levy [26], Levy et al. [34], Levy and Go´mez [28] and Levy and Montalvo [27] for the solution of (25). The TM has been designed to find the global optimum of general non-convex smooth functions by sequential improvement of local optima. It is able to find a sequence of minima, using only one initial guess (one run), that converges to the global solution. It has been proved to be robust for problems where multiple global minima at the same level of the objective function value are expected, see for instance Go´mez et al. [35], Nichita et al. [36] and Nichita and Go´mez [37]. There also exist parallel implementations of the TM, see Go´mez et al. [38]. The TM is also capable to deal with objective functions presenting very flat or sharp valleys. Resuming, the TM owns all the needed characteristics to report the global optima that characterizes the very ill-posed problem (25).
The TM is designed to solve problems like: FðxÞ where X is a subset of R n that could be defined through box constraints, and non-linear equality and inequality constraints, although the code we are using in this work, only solves box constraints and a linear inequality. The objective function F, F : R N ! R, is supposed to be differentiable.
The basic idea of the TM is to tunnel from one valley of the objective function to another finding a sequence of local minima x Ã 1 ; x Ã 2 ; . . . ; x G with decreasing function values. That is: where x G is the global minimum of FðxÞ. The following sequence of phases is repeated until the global minimum is found. Minimization phase: starting from any initial value x 0 k find a local minimum x Ã k with the gradient rF x Ã k À Á ¼ 0, using a local gradient based optimization method. We use in this phase the trust region Newton method for bound-constrained problems TRON of Lin and More [ 39]. TRON uses a gradient projection method to generate a Cauchy step, a preconditioned conjugate gradient method with an incomplete Cholesky factorization to generate a direction, and a projected search with trust region, to compute the step.
Tunnelization phase: from x Ã k find a point in another valley x tu k with lower or equal value of the objective function of the last local minimum, that is Fðx tu k Þ Fðx Ã k Þ.

S. Go´mez et al. / Well Test Analysis of Naturally Fractured Vuggy Reservoirs with an Analytical Triple Porosity -Double Permeability Model and a Global Optimization Method
In the tunnelization phase the following inequality problem is solved: If a solution x tu is found, it would certainly be in another valley of the objective function. The above inequality problem is solved by placing a pole at x Ã to destroy the local minimum already found and create a transformed problem using the tunneling function T . Examples of tunneling functions are: -classical tunneling function from Levy and Montalvo [27]: -exponential tunneling function from Barron and Go´mez [40]: where k Ã is the strength of the pole. Solving the inequality problem consists in finding x tu such that: T c ðx tu Þ 0 or T e ðx tu Þ 0; x tu 6 ¼ x Ã Descent directions can be taken to solve this inequality problem since T c x ð Þ and T e x ð Þ are continuous for x 6 ¼ x Ã . Therefore, it is possible to use the same local minimization gradient-based TRON method employed in the minimization phase, with appropriate stopping conditions, to find a point x tu with F x tu ð Þ F x Ã ð Þ. With x tu found, we set x tu ! x 0 0 and a new minimization phase is started. If no x tu is found after starting from a predetermined number of initial points of the tunneling phase, the last minimum of F is considered the putative global minimum. Checking global optimality is the most computationally expensive part of the algorithm.
In order to avoid going back to minima already found, it is necessary to keep " turned-on" the poles used to destroy the minima found before, the tunneling function then becomes: In the search for the solution of the tunneling phase, it may happen that the iteration may find a critical point of the tunneling function, and a descent direction cannot be found. Then a mobile pole has to be used, to be able to continue with the iteration. The new expression of the tunneling function would now be: The same strategy is used, when the objective function has flat or sharp valleys, and a descent step cannot be performed. This use of mobile poles, allows the method to arrive to global minima (or local) with the necessary precision.
The tolerances are very important and they highly influence the performance of the tunneling phase. A discussion on tolerances can be found in Castellanos and Go´mez [41] and specific tolerances required to search for multiple global minima at the same level are given by Nichita et al. [36]. The general stopping conditions and the most important features of the TM as placing mobile poles and pole strength among others are described in detail by Go´mez et al. in [35].
Another very useful feature of the TM is its ability to move away from points where local optimization methods like Levenberg-Marquadt or TRON get trapped. This happens very often in parameter estimation inverse problems, like the one we consider in this work, or in problems where the shape of the objective function has valleys very sharp or very flat. The method employs mobile poles as described in [35]. The overall efficiency of the TM for different kind of applied optimization problems as well as for standard global optimization benchmark problems is considered in References [35,38,42,43].

NUMERICAL RESULTS
The results of the automated well test analysis for the nine parameters identification problem defined in (25) will be presented in this section. We begin with the results obtained for synthetic generated data to show the robustness of the TM. All test were done on a PC with a Intel Core 2 Duo E7400 (3M Cache, 2.80 GHz, 1066 MHz FSB) CPU and 4 GB RAM.

Synthetic Generated Data
In order to implement a comprehensible test, we choose a set of significant values of the parameters and solved a NLS problem similar to (25) for synthetic data generated by all reasonable combinations of these values. The identification problem with synthetic data is not exactly the same as in (25) because for synthetic data it is only necessary to deal with eight parameters as it will be explained below.
Let us define the following sets A x ¼ 0:01; 0:5; 0:9 f g , A k ¼ 10 À9 ; 10 À6 ; 10 À3 È É , A j ¼ 0:01; 0:5; 0:99 f g , A C D ¼ 1 000; 500 000 f gand A s ¼ À3; 50 f g, then from all a in the Cartesian product set: we take the ones that represent meaningful reservoir characteristics. This results in 3 240 different a i 2 A and for all of these we generate synthetic dimensionless wellbore pressure curves in Laplace space according to Equation (7). It should be noted that for synthetic data generated according to (7) and (9) the value of the parameter k T is not needed. The parameter k T appears when transforming from real time to Laplace space and conversely through Equation (22). That is not the situation when real field data is considered and where the formulation of the NLS problem is the one presented in (25) taking into account Equation (24). Therefore, the corresponding NLS problem involves only eight parameters but is identical to (25) in every other aspect. Resuming, we have 3 240 different a i 2 A of the form: The dimensionless wellbore pressure curves in Laplace space corresponding to each a i are generated for 600 values of u logarithmically spaced in interval 10 À8 ; wherep ðiÞ k ¼ p wD ðu k ; a i Þ according to (7). Assuming all of these data sets as measurements we solved the associated NLS problems with the TM to identify the values of the parameters in a i . That is we solved the NLS problems: for all i ¼ 1; . . . ; 3 240.
To measure the robustness and efficiency of the TM, we first make a comparison with a local method to study the ability to get a global minimum, with good match to the data. For this, we also solved the NLS problems (27) with the TRON local method. The objective function tolerance was set as 10 À10 for the two optimization methods and the maximum allowed number of iterations is high enough for not being the cause for interruption of the computations. Under these conditions the robustness of the TM method for the total of 3 240 synthetic cases is shown in Table 2.
It is evident that the TM outperforms the TRON method. The efficiency of the method can be seen through the CPU time taken by the TM to solve the whole set of synthetic cases. This is shown in Table 3.   In Table 4 and Figure 1, the results obtained by the TM for one of the 3 240 synthetic cases without noise are shown. Note that p w denotes the difference between p 0 and p j , and p w the corresponding pressure derivative function. As we can see in Table 4, the values of x f and x v are interchanged between Min 1 and Min 2. These behavior can be attributed to the change of the value of j that is approximately 1 in Min 1 and approximately 0 in Min 2.

Synthetic Generated Data with Noise
To further validate the capabilities of the TM, we selected a subset of the 3 240 synthetic data sets and    added noise to the corresponding synthetic dimensionless wellbore pressure curves. We first truncate the pressure data to two digits of precision since this is approximately the precision of the measurement equipment. Then, we add Gaussian white noise of 1 percent of the pressure change. The so obtained data set is denoted byp ðiÞ k . We obtain a NLS problem similar to (27) but where the data setp ðiÞ k is replaced by the noisy data setp ðiÞ k . It is common practice in the oil engineering community to filter a noisy data set before performing a well test or other kind of analysis depending on the data set. Filtering real data set obtained from pressure gauges involves three steps: number of data points reduction, outlier removal and denoising. The first two steps are not needed in the synthetic noisy data set case. Therefore, we only conduct the denoising step following the three approaches cubic splines, characteristic functions (corresponding to the most common flow geometries in well testing) and wavelets as described in Go´mez et al. [44] and call the obtained data as filtered data. The best results are obtained using the cubic splines denosing method and the results reported in this work are always obtained using this method. After denoising, we proceed to the solution of the NLS problem with the TM. The results  Tables 5 to 7. As is shown in Table 5, the TM is able to find four different global minima with good match to the filtered data of case 623. The objective function tolerance in this example has been set to e ¼ 10 À7 . The fourth minimum, called Min 4 in Table 5 gives the best approximation to the exact parameter values. Nevertheless the remaining three minima also provide good match to the filtered data (Fig. 2). The use of a local method like TRON or Levenberg-Marquardt can only guarantee to find one of these global minima and not necessarily the one that reproduces the exact values.
In Table 6, the minima found for case 2 141 are shown. This example is representative of the difficulties in the interpretation of reservoirs with triple porosity. The derivative curve (Fig. 3) is very flat and can be interpreted as double or single porosity. Nevertheless in the values of the exact parameters can be seen that despite x f is small, k vf is large (at the upper bound) reflecting a possible triple porosity behavior. The values of Min 3 are the ones in better agreement with the exact parameter values.
In Table 7, the exact values and the ones of the only minimum obtained by the TM are shown. Here, the agreement between identified parameters and exact values is very good except for k mf .

Real Field Cases
The capabilities of the TM were proved in the preceding section for synthetic generated data. We also conducted a set of tests for real field data that we will call as Field 1, Field 2 and Field 3. These fields represent different wells at different reservoirs. The data belong to well pressure transient test conducted on reservoirs that are known to present a connected network of vugs in addition to the network of fractures. Therefore, the data from these well pressure transient tests are suitable to be modeled with the triple porosity-double permeability model presented in this work.
Before using the real measurements in the formulation of the NLS problem, the real data should be filtered. As mentioned above the filtering involves three steps: number of data points reduction, outlier removal and denoising. Each of the data sets Field 1, Field 2 and Field 3 comprise more than ten thousand points and therefore these three steps were necessary. For a description of the filtering approaches followed in this work see [44] and [45].
From the filtered data and known field characteristics, we compute the corresponding dimensionless pressure data curve in Laplace domain according to the method presented in Bourgeois and Horne [21]. The dimensionless pressure data curve is computed for a set of 600 logarithmically spaced data points u i in interval 10 À8 ; 10 3 Â Ã . The value of the parameter k T is included through the use of Equation (22). This pressure data curve is used as measurements in the NLS problem (25).  The results obtained with the automated well test analysis performed with the TM for the solution of (25) are shown in Table 8 and Figures 5 to 7. The objective function tolerance for field cases has been set to e ¼ 10 À7 .
For data sets Field 1 and Field 2, the TM only found one set of optimal parameters with the required objective function precision. For data set Field 3, the TM found two optimal solutions. The correspondence of the pressure and pressure derivative curves in real time with the filtered data is shown in Figures 5 to 7.
To conclude this section, we point out that the computing times of the TM for all the examples, presented in the tables and figures range from three seconds to eight minutes. For example, the computing times of the TM for the real cases Field 1 to Field 3 are 3:16, 1:19 and 1:42 minutes respectively. These computing times are typical for all the tests conducted and this shows the efficiency of TM for the kind of problems treated here.

CONCLUSIONS
In this work, we employed the triple porosity-double permeability model of Camacho et al. [10] considering primary flow in both networks of fractures and vugs. The model was selected because we wanted to perform automated well test analysis of real reservoirs that exhibit this type of behavior.
In this work, we show that for the triple porositydouble permeability model, in the case when j ¼ 0 or j ¼ 1, several different sets of parameters may exist with one identical pressure response. This different sets can therefore be interpreted as multiple interpretations of the well test analysis. When 0 < j < 1, this situation can not be concluded so far. Plot of the real pressure drop data and the corresponding pressure drop derivative vs real time of Field 1. Also plotted are the filtered data curves and the curves obtained with the optimal parameters computed by the TM.  Plot of the real pressure drop data and the corresponding pressure drop derivative vs real time of Field 2. Also plotted are the filtered data curves and the curves obtained with the optimal parameters computed by the TM. Plot of the real pressure drop data and the corresponding pressure drop derivative vs real time of Field 3. Also plotted are the filtered data curves and the curves obtained with the two sets of optimal parameters computed by the TM.
The automated well test analysis, posed in the form of a NLS, is known to be a very ill posed problem that may have multiple optimal solutions with good match to the data. Therefore for the solution of this kind of problems, stable global optimization methods are needed. We selected the TM of Levy and Montalvo [27] and Levy and Go´mez [28] for its abilities to find a set of global optimal solutions. The TM is also able to deal with objective functions presenting very flat or sharp valleys which is the case in automated well test analysis. Local optimizations methods may fail dealing with this kind of problem, see Table 2.
The global optimization NLS problem is formulated and solved in Laplace domain as the model solutions are know analytically in this domain. This approach allows to compute analytical gradients of the objective function and also avoids the time consuming task of returning to real time at each step of the computations with the Stehfest method. Considerable time savings of about 60 percent, computing time are obtained.
We conducted an extensive set of tests with synthetically generated data, both with and without noise. The synthetic data was generated by choosing a representative set of parameter's values. We obtained 3 240 representative cases where eight parameters were to be identified. The results obtained by the TM for this set of synthetic data, without noise, show the robustness of the method. In 100 percent of the cases the TM was able to find at least one global optimal solution whereas the local method TRON was only successful in 59 percent of the cases. The computing times of the TM in both synthetic and real field cases are extremely small. The computational efficiency of the TM has been demonstrated also by measuring the CPU time for by both, the synthetic and the real field cases.
We also conducted the automated well test analysis of three real field cases. We properly filtered the real data set and transformed it to Laplace domain. As for synthetic generated data, the results obtained in this work for real field cases reveal two main facts: one is the existence of multiple interpretations with good match to the data and the other is the necessity of other sources of information to be able to decide which one of these interpretations is the correct characterization when using pressure data only. The existence of multiple interpretations imposes the use of global optimization methods with the characteristics of the TM. An example of additional information could be approximate values of x f and x v from well logs and core data analysis. In a current work in progress, we follow the approach of considering x f and x v as given and the results obtained so far are very encouraging in the sense that the TM finds only one solution.
k vf ¼ b 1 À c 21 k mf À Á c 12 À a 1 À c 11 k mf À Á c 22 c 12 c 23 À c 13 c 22 Substituting (B6) and (B7) in (B3), we obtain a quadratic equation for determining k mf depending on a 0 , a 1 and b 1 . Again, for certain values of a 0 , a 1 and b 1 , two solutions are possible.
It is therefore possible to obtain up to four different parameters sets x f ; x v ; k mf ; k mv ; k vf È É that when substituted on (16) and (19)