Determination of Geostatistical Parameters Using Well Test Data

Résumé — Détermination de paramètres géostatistiques par l’utilisation des données d’essais de puits — Des études de réservoirs de plus en plus nombreuses font appel à une description probabiliste du réservoir. Ceci permet de quantiﬁer les incertitudes sur les prévisions de production, à condition de disposer d’un modèle géostatistique convenablement paramétré. Les paramètres les plus courants sont la (ou les) longueur(s) de corrélation et la variance de la perméabilité. Ces paramètres sont déterminés à l’aide de considérations géologiques, auxquelles on ajoute éventuellement des mesures pétrophysiques sur carottes à l’échelle du laboratoire. Il serait toutefois intéressant de disposer d’une technique de mesure directe de ces paramètres, en utilisant des données investigant des échelles plus importantes du réservoir considéré. Dans cet article, nous examinons l’apport des essais de puits, qui ont l’avantage d’être disponibles au début de l’exploitation du gisement, et dont on peut montrer qu’ils “moyennent” la perméabilité du réservoir sur une zone correspondant à leur rayon d’investigation. La méthode proposée Abstract — Determination of Geostatistical Parameters Using Well Test Data — In this paper we describe a new method to obtain estimations of the geostatistical parameters (GPs) such as the correlation length, l c and the permeability variance, σ 2 ln from well test data. In practical studies, the GPs are estimated using geological and petrophysical data, but often, these data are too scarce to give precise results. The proposed method uses the Bayesian inversion theory, in conjunction with a fast evaluation of well tests that implies upscaling techniques. The method was tested using synthetic well-test data performed on some training images, and estimations of the underlying correlation length, l c and permeability variance, σ 2 ln were recovered. These estimations give a correct order of magnitude of the actual values, but as noticed in similar methods, the uncertainties are high. Once the GPs are estimated, other well established techniques can be used to get well-test matched reservoir images consistent with the geostatistical model. We will see that excellent well test data are needed, and that the method could be improved using multiple well test data.


Mathematical notations δ (r)
Dirac delta-function Weighting kernel for ring permeabilities G 2D and G Represent the 2D well test simulator and his fast approximation h Depth of the reservoir (m) J Jacobian matrix of G (mD −1 ) k j j=1...N r Ring permeabilities (mD) k e ff Effective permeability (mD) k g Geometric mean of the permeabilities (mD) k w Wellbore permeability (mD) k (t) Estimation of the apparent permeability at time t (mD) k app (t) Apparent permeability at time t (mD) L (Θ) Negative of the log-likelihood function

INTRODUCTION
The use of geostatistical models to improve production forecasts is becoming more and more popular among reservoir engineers. In order to reduce the forecasting uncertainties, an important challenge is to history match production data, using the information provided by geostatistics in a consistent manner. This problem was investigated by many authors, using various optimization techniques and the Bayesian paradigm [1][2][3][4][5][6][7][8]. In general, the underlying geostatistical model, as well as the geostatistical parameters (GPs) such as the correlation length, l c and the log-permeability variance, σ 2 ln are supposed to be fixed. Most often, in practical studies, the GPs are estimated using geological information and statistical analysis of the available petrophysical data [9,10]. In a lot of cases, these data are too scarce to give precise results, and more particularly if the typical interwell spacing is very large compared to the correlation length. This problem received some attention in the oil-oriented technical literature. Yadavalli et al. [11,12] using Sagar [6] and Feitosa' s [13] approach were able to get estimations of the GPs from a single well test. Nevertheless, in their approach, they derived the variogram of a two-dimensional permeability field from the computation of the variogram of a radial one-dimensional equivalent permeability field that leads to the correct well test. Their main assumption is that the two variograms are equivalent. We will show that this is not the case in general and that the relation between the two variograms is rather complex (Eq. (25)). It is therefore not obvious to give an interpretation of their results. Yortsos [14] proposed an alternative technique that works when a lot of independent well test data are available.
The determination of the GPs from head measurement data received considerable attention in the hydrogeology literature (see [15][16][17] and references therein). The authors, in these studies, used a Bayesian framework that formalized the history matching process using least square optimization techniques. In the present paper, we will adapt this general framework to the particular case of well test data. The main idea is to evaluate the joint probability density of having any given reservoir image and any given set of dynamic data measurements, for fixed GPs. The GPs estimation is done using the maximum likelihood principle which states that among the GPs, the probability density to get the actual dynamic data must be maximum. The method was tested using synthetic well-test data on some training images, and estimations of the underlying correlation length, l c and permeability variance, σ 2 ln were found. As noticed in the litterature dealing with the geostatistical approach of inverse problems [18,17,19,20], we will find that the estimates give the good order of magnitude of the actual values but the uncertainties are high. Once the GPs are estimated, one can thus obtain well-test matched reservoir images consistent with the geostatistical model using standard inversion methods.

WELL TESTS AND APPARENT PERMEABILITY
Well tests are commonly used in reservoir engineering for better understanding the structure of the reservoir. Typically the investigated scale is hectometric, and important properties such as the average permeability, the reservoir boundary positions, skin and wellbore storage effects can be estimated from a well test interpretation.
Thanks to the improvement of pressure gauges, modern well test interpretation techniques make an intensive use of the pressure derivative curve. This tool, as well as the associated concept of instantaneous apparent permeability are very powerful techniques to interpret well tests in a geostatistical context. In a first part, we present our main assumptions and introduce the apparent permeability as well as its main properties.

The Diffusivity Equation
Hereinafter, we will work with the following assumptions: -the reservoir is 2D and of infinite dimension -the flow is single-phase -the porosity, φ and the viscosity µ are constant -the fluids and the rock are considered slightly compressible and the total compressibility c t is constant -skin and wellbore storage effects are negligible Within this framework, considering that a well located at the origin is pumping at a constant rate Q from t = 0, the pressure variation P (r, t) and the permeability map, k (r) are related by the diffusivity Equation [21][22][23]: h denotes the reservoir thickness and the second term of the right hand side is a source term. It is the product of a spatial Dirac function δ (r) by a Heaviside function Y (t) (which is equal to 0 if t < 0 and 1 if t ≥ 0). This approximation (known as the line source approximation) holds if the well radius is infinitesimally small with respect to the typical scale of the reservoir.
To get a well defined problem, we must fix some boundary conditions. Initially, the pressure variation is uniformly equal to 0 (i.e.∀r, P (r, t = 0) = 0) and the pressure variation vanishes far from the well (i.e.∀t, P (r → ∞, t = 0) = 0).

The Apparent Permeability
This diffusivity Equation can be solved analytically for a constant permeability, k 0 [22,23]. Typically, the product k 0 h can be derived from the bottomhole pressure P 0 (t) = P 0 (r → 0, t): This relation is the basis of the derivative technique in well test interpretation. Additional physical phenomena such as skin, wellbore storage, and practical constraints transform well test interpretation into a full speciality. An important technical literature exists on the subject, with intensive use of type-curve analysis and numerical optimization methods [21,22,24].
Equation (2) can be generalized in the heterogeneous case defining the apparent permeability concept: Of course, for a homogeneous reservoir, k app (t) is independent of time and is equal to the constant permeability, k 0 . We need also an estimation of the radius of investigation which corresponds to the range investigated by the well test: Here, the dimensionless constant C, the order of magnitude of wich is unity, depends on which definition is used [24]. For the sake of simplicity, we will set C = 1 throughout the paper.
In a geostatistical context, the relation between k app (t) and the permeability field k (r) is very complex and a gridded reservoir simulator is needed. Hereinafter we will consider that k (r) is a random variable (1) , and that the actual permeability field is only a particular realization of k (r).
In this statistical framework, some important asymptotic results can be derived using perturbation methods [26] that can be extended in a systematic manner using a Feynman (1) Gelhar [25], in the introduction of his book ( § 1.3), mentions some key points for working with random variables and stochastic processes. This approach can formally handle the uncertainties due to the lack of measurements of porous media characteristic parameters.
diagram approach [27]. For short times, the apparent permeability tends to the wellbore area permeability, k w : lim t→0 k app (t) = k w (5) and for long times, the apparent permeability tends to the effective permeability k eff corresponding to the upscaled reservoir permeability : It can be shown that the typical time to reach this last limit is such that R app (t) = l c , when the radius of investigation is in the order of magnitude of the permeability correlation length. This means that the well test is performing a homogenization of the reservoir permeability. This is the physical basis of the proposed method where the geostatistical parameters will be estimated from well tests.

FAST EVALUATION OF WELL TESTS
Since we will use an optimization method to evaluate the geostatistical parameters and, therefore, an important amount of simulations will be required, a fast approximation of well tests is needed. The approximation we will use was first proposed by Oliver [4] for low permeability contrasts where a first order expansion of the diffusivity Equation (1) was considered. This approach was further extended by Sagar [6] and Feitosa [13] for high permeability contrasts. They found that the apparent permeability can be considered as an averaging process over volumes of increasing size with respect to time. More precisely, the approximation of the apparent permeabilityk (t) is a weighted sum of a "homogenized permeability field" that contains N r ring permeabilities The weights G ij depend on a kernel function G (r, t) (2) : represent the zero, one and two order modified Bessel functions Definition of ring permeabilities and ring radii for a homogenized field containing 6 rings. Permeability weighting kernel √ zG (z) vs. z = r 2 /t.
We have plotted in Figure 2 the function √ zG (z) vs. z = r 2 /t which appears in Equation (9). The maximum of this function is for r = 0.92R app (t) and 99% of the whole contribution comes from the area contained between the two radii r min = 0.12R app (t) and r max = 2.34R app (t). Furthermore, we have a normalization condition: and, assuming that r 1 = 0 and r N r = ∞: The next step is to relate the ring permeabilities k j to the original 2D permeability map. Feitosa [13] used different averaging formulae (arithmetic, harmonic and geometric) to compute k j from a 2D lognormal permeability field. He compared the resulting estimations with some reference results obtained using a complete finite difference well test simulation. He found that the best match of the apparent permeability is obtained when the geometric mean is used. This result is consistent with the general properties of the apparent permeability given by Equation (6). For short times, independently of the averaging rule,k (t) will tend to the wellbore permeability. For long times, using (6),k (t) has to converge to the effective permeability which is precisely the permeability geometric mean in the 2D lognormal case [28]. Since, in the following, we will only consider lognormal permeability fields, the ring permeabilities are the geometric mean of the fine scale permeabilities contained in the ring: where K l denotes the Cartesian permeabilities, Ω j are the indices of the gridblocks overlapping the ring j and ω l is a factor that represents the fraction of the gridblock's area contained in the ring.

GEOSTATISTICAL PARAMETER ESTIMATION USING WELL TEST DATA
Kitanidis and Vomvoris [19] estimate geostatistical parameters using steady-state pressure and permeability measurements for a 1D flow problem. Our present goal is to generalize the method for transient flow. Now, the data considered are measurements of apparent permeability values. Nevertheless, the same formalism can be applied. The interested reader can refer to different papers [19,18,17] and to the references therein.
Let us now present the geostatistical approach and some notations. We will consider that we have N t apparent permeability measurements and a permeability map described by the value of the permeability at the centers of N m gridblocks.
We define an observation space and a parameter space where: -Y is a N m parameter vector containing the parameters of the problems, which are, in our case, the permeability logarithm map.
The parameter mean Y = Y is considered as known.
is a known function of a few geostatistical parameters contained in a vector Θ (variance, correlation length, nugget effect, etc). The variogram shape is assumed to be known. -d is a N t vector containing all the observations which, in this case, are the inverse of the apparent permeability values for a set of time steps.
Calling G 2D the transfer function (the simulator) that maps the parameter space into the observation space, we obtain the following relation: v is an random error vector. We suppose that v = 0 (the estimation is non-biased) and that the error covariance matrix R 2D = vv t does not depend on Θ.
We denote p (d, y|Θ), the probability density function PDF of the model parameters and of the observations, p (d|y, Θ), the PDF of the observations assuming that the model parameters are known and, p (y|Θ), the PDF of the model parameters assuming that the geostatistical parameters are known. Using Bayes theorem, we obtain We define the marginal PDF of the observations as the sum over all possible configurations of the model parameters of p(d, y|Θ): The sum is taken over N m dimensions. Kitanidis [17] shows that geostatistical parameter estimation can be derived from two different methods: -the Maximum A Posteriori (MAP) method: Θ estimation is obtained by maximizing p (d, Y|Θ) with respect to both Y and Θ. -the Geostatistical Approach (GA) method: the estimation of Θ is obtained by maximizing p (d|Θ) with respect to Θ. This method was widely used in the literature [19,18,20,29,30]. Therefore two different estimations are obtained: Θ MAP and Θ GA . Kitanidis [17], by studying asymptotic results, showed that Θ MAP is a biased estimation whereas Θ GA is not. Moreover, he demonstrates, using an illustrative example, that the bias is highly dependent on the model discretization.
Therefore, we shall focus on the GA method. We assume that the two PDFs p (d|Y, Θ) and p (Y|Θ) are both Gaussian: The second assumption is true since we consider lognormal permeability fields. In order to make notations easier, we assume implicitly the dependence of C on Θ. Equation (18) becomes: with The GA approach involves the maximization of (21) to estimate the GPs. Here we have to point out that the sum is taken over N m dimensions which implies typically 10 6 gridblocks. In practice, it is impossible to evaluate such an integral. Furthermore, the GA method involves an optimization step, and consequently many flow simulations which are very CPU time consuming. Therefore, we will have to simplify the problem: we will use the fast evaluation of the apparent permeability described in the previous section. Considering that the estimation is close to the simulation results, Gautier [31] shows that we have: y is the vector containing N r radial permeabilities defined by (12). Its mean and covariance matrix are derived from the 2D permeability field: ρ represents the correlation function. The computation of the matrix C is expensive since it involves a large number of products. Nevertheless, if matrices are stored, it has to be done only once for a configuration of ring and for a particular variogram shape. From (7), the observations and the parameters are now related by: G is the matrix containing the weights of the kernel function G ij defined by (8). The simplified problem involves the maximization of p (d|Θ) that can be written in the new parameter space: with Usually, instead of maximizing the observation PDF with respect to the GPs, we minimize the negative log-likelihood function L (Θ), which is defined by: If the relation G (y) between the model parameters and the observation is linear, the summation in (27) can be handled analytically using Gaussian integration formula. When the relation is non-linear, we use an approximate method known as the Laplace method or in hydrology as the quasi linear geostatistical approximation [18]: it assumes that I (y) is a peaked function, and that the near peak area represents the main contribution to the integral. The method can be divided into four steps: 1. Choose a first guess for the GPs, Θ k=0 2. First optimization.
Θ k is fixed. Findỹ which maximizes I (y) using a Gauss-Newton method: 3. Computation of the negative log-likelihood for Θ k . Compute L (Θ k ) analytically by linearizing I (y) nearỹ (Using Laplace Method) 4. Second optimization.
Computation of the gradients (∂L (Θ) /∂Θ) and of the Hessian matrix Ξ evaluated at Θ k : Go back to step 2 until convergence is reached. More mathematical details can be found in Appendix A.

Well Test Example
In this part, we will briefly describe some numerical results. Reference well test simulations are performed using SIMTEST, a 3D finite volume simulator developed at the IFP [32]. First the diffusivity Equation (1) is solved. Then, using the bottomhole pressure measurements, we can compute the apparent permeability for a set of times.
Since time increases geometrically the pressure logderivative is computed using [22]: Figure 3 displays a zoom of a permeability field realization around the well (50 gridblocks in each direction). Map of a particular realization of a permeability field (in logscale) with σ 2 ln = 1, l c = 6 gridblocks, k w = k g = 100 mD.  1 Properties of the grid, of the well and of the fluid and rock.

Well properties and boundary conditions
Well radius r w = 8 · 10 −2 m Well rate Q = 1.157 · 10 −3 m 3 /s BC on the reservoir border no flux

Fluid and rock properties
Porosity Total compressibility c t = 10 −4 bar The grid and the flow properties are found in the Table 1.
The permeability field has been generated by a sequential Gaussian simulation provided by the GSLIB library [9]. It has a lognormal distribution with a variance, σ 2 ln = 1 and a correlation length, l c = 6 gridblocks. The well permeability and its geometric mean are k w = k g = 100 mD. A well test has been performed for 100 days. Figure 4 represents the apparent permeability with respect to time. We can see that the well test can be divided in four main periods: 1. The well test is only sensitive to the wellbore area and the apparent permeability is equal to the well permeability. 2. The well test is in a transient period and the fluctuations are more important. Actually, the well test performs a permeability average over volumes increasing with time.     This averaging property pointed out by Oliver [4] has been discussed in the second Section (cf. Eq. (7)). 3. When the investigated volume is large enough, the apparent permeability tends to the geometric mean which is the effective permeability for a 2D lognormal permeability field [28]. 4. During the last period, the apparent permeability decreases sharply since the boundary conditions start to influence the well test. The corresponding time is typically dependent on a maximum radius of investigation [31] the value of which is equal to one third of the wellboundary distance: This value is in agreement with Oliver's approximation [4] (Indeed we have seen in the second Section that the kernel G (r, t) performs a permeability average within the range 0.12R app (t) , 2.34R app (t) ) Here, we pointed out the averaging process performed by the well test. Other effects, such as the influence of the permeability variability or of the correlation length can be found in Gautier [31].

Fast Evaluation of Well Tests
The fast evaluation of well tests is performed by first transforming the cartesian permeability grid into a homogeneous radial permeability field defined by a set of rings defined by (12). Using (7) we are then able to compute an estimation of the apparent permeability. In Figure 5, we show two examples. The permeability field description can be found in the Table 1 and the 50 ring radii used for the estimation are in geometrical progression with a minimum and a maximum radius, equal to ∆ x /2 and N x ∆ x /2 respectively. We can see that for both cases the difference between the apparent permeability and its estimation is small and fluctuates more for intermediate values of the investigated radius.
To check the global accuracy of the fast estimation, we computed the average of the estimation errors for well tests performed on 70 realizations of two permeability fields with σ 2 ln = 1 and l c = 50, 100 m ( Fig. 6 and 7). The bold line represents the errors (in %) and the thin line, the averaged radius of investigation. Globally, the errors remain less that 10%, and are smaller for short and long times. Both estimation errors have a small peak at the beginning. This peak appearing for an averaged radius of investigation close to half the size of a gridblock is due to simulator numerical errors and is independent of the permeability field properties. Of more importance is the second peak which is located at a distance from the well close to the correlation length. This is in accordance with the homogenization process performed by the well test. For short times (short radius of investigation), we have seen that both the apparent permeability and the estimation are close to the wellbore area permeability. For long times, if the reservoir is large enough, the apparent permeability tends to the effective permeability. The estimation tends to the permeability of the last ring which is  basically equal to the geometric mean if the last ring contains a sufficient amount of gridblocks. Therefore, the peak at the correlation length means that the apparent permeability is no longer sensitive to the well permeability and that it is not yet sensitive to the homogenization process. The fluctuations are bigger during this period and the approximation is slightly worse.

Geostatistical Parameter Estimation
We will now describe the GPs estimation using well test data. First we have to make some assumptions. As pointed out in the former section, the fluid flow simulator is not precise for short times (typically for R app (t) < ∆ x /2). For long times, the well test is sensitive to the boundary conditions which are independent of the permeability field structure. The last time step used in the optimization process will have a corresponding radius of investigation close to R max defined in (34). Therefore the set of time steps (t i ) i=1...N t will be selected such that: We will consider that the observation error matrix is the identity matrix multiplied by a proportionality factor σ 2 w . We will consider hereinafter that σ 2 w = (10%) 2 . 10% represent the tolerance of the discrepancies between the observation data and the estimations. Therefore R in (27) is defined by: Finally, we also include a parameter σ 2 c which adds an information (obtained by other measurements): our knowledge of the well permeability (3) . We will also consider that this parameter is fixed and that σ 2 c = (10%) 2 . Let us consider an example of a well test performed on a 799 × 799 permeability field (the other data are in Table 1). The variance and the correlation length of the variogram are σ 2 ln = 1 and l c = 50 m respectively. Figure 8 represents the apparent permeability with respect to the radius of investigation.
For this particular example, we have made a complete exploration of the GPs space and we have computed for each value of the GPs the negative log-likelihood function L σ 2 ln , l c . The L contour plot is represented in Figure 9. We can see that the minimum is for values of the GPs close to σ 2 ln = 1, l c = 40. The minimum of L is more precisely determined in Figure 10 where we used the optimization algorithm with different first guesses. We see that, in all the cases, the algorithm converges to a same point σ 2 ln = 1, l c = 42 in 5-10 iterations. The next two figures (11(a) and 11(b)) represent the GPs estimation for well tests performed over two other realizations of the same permeability field and show us that the estimation depends on the realization. A statistical analysis of the estimation is therefore necessary. Figure 12 represents the histograms of the correlation length and of the variance estimation based on well tests performed on 70 realizations of the permeability field that we have already considered. Statistical data are presented  Negative log-likelihood function L σ 2 ln , l c for the well test plotted on Figure 8.
in the fourth row of Tables 2 and 3. For this particular example, we found that the averages of the variance and of the correlation length estimations (resp 0.85 and 73 m) gives Minimization of the negative log-likelihood function L σ 2 ln , l c for the well test plotted on Figure 8 and for different first guesses.
a good order of magnitude of the theoretical values (resp. 1 and 50 m).
We performed the same study for other permeability fields, with σ 2 ln = 1, 2, 4 and l c = 30, 50, 100 m. In order to make comparisons between the results, we have plotted in Figures 13 and 14 the histograms of l c /l c and σ 2 ln /σ 2 ln ( l c and σ 2 ln represent the estimations). The statistics are presented in Tables 2 and 3. The main characteristics of these statistics are as follows: -The mean of the correlation length estimation overestimates the theoretical value. The difference lies between 20 − 30 m. The standard deviation is quite high (more that half of the mean value). -The mean of the variance estimation is close to its theoretical value. Nevertheless, the uncertainty remains high and is close to the mean value.
It is not clear if this overestimation is due to the estimation method or if it is simply due to an artifact which depends on the numerical scheme used to solve the well test Equation. In fact, the harmonic mean has been used to compute the transmissivity between two adjacent gridblocks. The well test performs a moving average of the transmissivities which are a function of the permeabilities. Therefore, to a certain extent, the apparent permeability reflects the statistical properties of the transmissivities [33]. A theoretical study would be necessary to relate the statistical properties of the transmissivities and of the permeabilities. We have also studied the influence of the control parameters (σ 2 w , σ 2 c , R min , R max ), and the evolution of the estimation statistics with respect to the number of tests [31]. For this latter effect, it has been shown, that, typically, 5-10 estimations of the GPs are needed to obtain a mean value close to the stable one.

CONCLUSIONS
Using a general Bayesian formalism, we built a method to estimate the dominant geostatistical parameters, l c and σ 2 ln , using well test data. The proposed method is rigorous, and general. To get practical results, we had to couple up-scaling techniques and inversion procedures.
Due to the rather low dependency of the well test on these GPs, the obtained estimations have important fluctuations. Accurate well test data are needed, and the results are improved when a lot of independent well test data are available.
The proposed technique was implemented in the lognormal case, and it would need to be generalized to discontinuous permeability distributions that occur often in practice. It could also be adapted to account for other fluid-flow data, such as interference tests, and multiphase flow which are likely to contain more information. Histograms of l c /l c from well tests performed on 70 realizations of 9 permeability fields with σ 2 ln = 1, 2, 4 and l c = 30, 50, 100.   Figure 14 Histograms of σ 2 ln /σ 2 ln from well tests performed on 70 realizations of 9 permeability fields with σ 2 ln = 1, 2, 4 and l c = 30, 50, 100.  Here, we will describe the different steps of the optimization algorithm. An important part of these results is based on the research of Kitanidis [17]. The first step is equivalent to minimizing −2 I (y) defined in (28) and appears in many papers. One can refer to Tarantola [34] for a thorough review. We use an iterative Gauss-Newton method: The (N r × N r ) matrix, H y n , is the Hessian of I (y) computed at the point y n : Exact computation of the Hessian requires the second order derivatives of I (y). But in the Gauss-Newton method, it is not necessary to have a very accurate evaluation of the Hessian, which is computed by differentiating twice Equation (28), ignoring second-order terms: where J y n is the (N t × Nr) Jacobian matrix of G (y) evaluated at y n (see its expression in Appendix B). By deriving Equation (28): ∂I (y) ∂y y n = −J y n R −1 (d − G (y n )) + C −1 (y n − y) (40) Therefore, including Equations (39) and (40) in Equation (37), we obtain: y n+1 = y n + J t y n R −1 J y n + C −1 −1 × J t y n R −1 (d − G (y n )) + C −1 (y n − y) Or using two classical relations (4) and a some matrix algebra: y n+1 = y + CJ t y n R + J y n C J t y n −1 Therefore the Hessian matrix Ξ n can be asymptotically estimated by its mean value given by the information Fisher matrix. After calculation: and finally to update Θ n , we use the following equation: Expressions of the Jacobian matrix of G (y), and of the Fisher matrix can be found in Appendix B.

APPENDIX B
Derivation of the Jacobian and the Fisher Matrices using Equations (7) and (26), and performing the differentiations: which is the expression of the Jacobian of G (y).

Expression of the Fisher Matrix
To compute the Fisher matrix (52) we have to find the derivatives of Σ and d with respect to Θ. Using (45) To compute the covariance matrix derivative with respect to the GPs, we decompose the sums in (25) into two different contributions: Let's assume that the correlation function is exponential with a variance, σ 2 ln and a correlation length l c , ie ρ (r) = σ 2 ln e −r/l c . Therefore Θ = σ 2 ln l c t and: and with D ij = (k∈Ωi),(l∈Ω j) ,k =l ω l ω k d kl ρ (d kl ) The derivative components of y defined by (24) with respect to the GPs, are: Since we assume in all our applications that the mean of the permeability field is constant: Copyright © 2004, Institut français du pétrole Permission to make digital or hard copies of part or all of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. Copyrights for components of this work owned by others than IFP must be honored. Abstracting with credit is permitted. To copy otherwise, to republish, to post on servers, or to redistribute to lists, requires prior specific permission and/or a fee. Request permission from Documentation, Institut français du pétrole, fax. +33 1 47 52 70 78, or revueogst@ifp.fr.