Suscripción a Biblioteca: Guest
Portal Digitalde Biblioteca Digital eLibros Revistas Referencias y Libros de Ponencias Colecciones
International Journal for Uncertainty Quantification
Factor de Impacto: 4.911 Factor de Impacto de 5 años: 3.179 SJR: 1.008 SNIP: 0.983 CiteScore™: 5.2

ISSN Imprimir: 2152-5080
ISSN En Línea: 2152-5099

Acceso abierto

International Journal for Uncertainty Quantification

DOI: 10.1615/Int.J.UncertaintyQuantification.2012003480
pages 323-339


Richard Booth

Schlumberger BRGC, Rio de Janeiro, Brazil

Kirsty Morton

Schlumberger BRGC, Rio de Janeiro, Brazil

Mustafa Onur

Technical University of Istanbul, Istanbul, Turkey

Fikri Kuchuk

Schlumberger Ribould Product Centre, Clamart, France


In any subsurface exploration and development, indirect information and measurements, such as detailed geological description, outcrop studies, and direct measurements (such as seismic, cores, logs, and fluid samples), provide useful data and information for static reservoir characterization, simulation, and forecasting. However, core and log data delineate rock properties only in the vicinity of the wellbore, while geological and seismic data are usually not directly related to formation permeability. Pressure transient test (PTT) data provide dynamic information about the reservoir and can be used to estimate rock properties, fluid samples for well productivity, and dynamic reservoir description. Therefore, PTT data are essential in the industry for the general purposes of production and reservoir engineering as well as commonly used for exploration environments. With the need for improved spatial resolution of the reservoir parameters, grid-based techniques have been developed in which the reservoir properties are discretized over a fine grid and characterization of the probable state of the reservoir is sought using the Bayesian framework. Unfortunately, for the exploration of hydrocarbon-bearing formations, the available prior information is often limited: in particular, unexpected geological features, such as fracture and faults, may be present. There are two groups of recent methods for dynamic characterization of the reservoir: (i) data assimilation techniques, e.g., ensemble Kalman filter (EnKF) and (ii) maximum-likelihood techniques, such as gradient-based methods. The EnKF is designed to produce a set of realizations of the reservoir properties that fit the PTT data; however, the method often fails to honor the data when unexpected features are not captured by the prior model. The alternative gradient-based methods do provide a good fit to the PTT data. They can also be made efficient for high-dimensional problems by using an adjoint scheme for determining the gradient of the log-likelihood function. However, as a maximum likelihood technique, this method only yields a single realization of the reservoir. It is important to maintain a model of the uncertainty of the reservoir characterization after PTT data assimilation, so that the risk associated with future decisions is understood.We therefore present and investigate a stochastic, gradient-based method that allows for proper sampling of realizations of the reservoir parameters that preserve the fit with the PTT data. The results indicate that our proposed method is quite encouraging for efficiently generating realizations of rock property distributions conditioned to PTT data sets and a given prior geostatistical model.

KEYWORDS: Bayesian inference, geological applications, transient well test interpretation, reservoir characterization


Pressure transient testing is a long established procedure [1–3] for determining the productivity of a well and the properties of the formation (reservoir) from downhole and/or surface pressure and flow-rate measurements. In a pressure transient test (PTT), a well is produced (or injected) for a known period of time (drawdown test) and subsequently may be shut in for a buildup test, or for a sequence of production and buildup periods. During the PTT, the flow rate at the production well (surface or downhole) and transient pressure in the wellbore and/or in a observation well are measured. Using conventional interpretation methods [Horner or Semilog (MDH) and/or type-curve matching with the aid of derivative plots], reservoir pressure, an effective ‘average’ permeability of the reservoir, the skin factor, the wellbore storage, etc., are obtained from the PTT data.

The permeability obtained from a single drawdown or buildup test is an average permeability over a radius of investigation [4]. In other words, the mild permeability heterogeneity around the wellbore within a radius of investigation is averaged by the transient pressure diffusion during a transient well test. However, if in addition to the well test data acquired at the production well, the transient pressure data acquired at other spatial locations are available (interference tests performed in several wells [5]), it is possible to establish a more complete picture of the spatial distribution of the reservoir parameters. For instance, vertical interference tests conducted with multiprobe wireline formation testers can be used to determine the spatial distribution of the permeability vertically around the wellbore [6]. However, whether using a single transient test or multiwell (multiple spatial location) interference tests, it is well known that the sensitivity of PTT to small-scale spatial variations of the permeability or other formation parameters is fairly low. On the other hand, significant geological features, such as a region of low permeability (sealing faults) or high permeability (conductive fractures and faults), might still be identifiable from the pressure transient data.

At the exploration phase during the early years of the well, when PTTs are used for reservoir characterization and determination of well productivity, the prior reservoir knowledge is usually minimal at the well test scale. Larger scale reservoir properties might be available from seismic and geology and/or outcrop analogies, and at a small scale, the prior formation knowledge might be available from logs or cores. Therefore, for exploration and development pressure transient testing, the existence of and locations of many subseismic features, such as strong spatial permeability variations, faults, fractures, pinch outs, etc., are not known a priori. Furthermore, we are unlikely to have a detailed geostatistical description of the reservoir.

The method of Backus and Gilbert as also described in [4, 7], provides one method for determining the spatial variation of reservoir properties. However, the assumption of small variations in the properties makes the method of Backus and Gilbert inappropriate for geological feature identification. Moreover, its lack of robustness when it is used with inaccurate data suggests the need for some sort of prior modeling.

We therefore need to follow a Bayesian framework to use PTT data for spatial resolution of reservoir parameters, with a probabilistic prior description of the reservoir updated by the measurements to give a (probabilistic) posterior description of the reservoir. To capture the spatial variability of the reservoir parameters, we turn to “pixel” methods, with the reservoir parameters discretized over a grid. In previous works, these grids have typically been regular [8, 9], principally because the prior geostatistical model required a regular grid. Irregular grids with a finer resolution close to the wellbore are useful for the numerical simulation of the well test, and moreover, we expect that the resolution of the inversion should be highest closest to the wellbore. For this reason, we use the same irregular grid that is also used for numerical simulation of the PTT data to discretize the reservoir.

In this paper we first describe a prior model that can efficiently describe the geostatistics of the reservoir properties on an irregular grid.We then obtain an expression for the posterior probability distribution of the reservoir parameters. The probability density function is expensive to evaluate, and therefore, a more simple description is required. One of the most useful pieces of information that can be extracted from the posterior probability density is the maximum likelihood parameters. However, this information alone does not tell us how likely these parameters are and, thus, some measure of the uncertainty is also usually desirable, typically either random samples from the distribution or a Gaussian approximation to the posterior probability distribution.

One method that has received much recent attention [10, 11] is the ensemble Kalman filter (EnKF), in which the probability distribution is represented by an ensemble of samples that are updated as each new measurement becomes available. The ensemble Kalman filter assumes that at each step the probability distribution is Gaussian, and failure of either the Gaussian model to capture the true probability distribution or the ensemble to accurately capture the Gaussian model can lead to failure of EnKF. This is potentially more of a problem for inversion of PTT data than for other applications of EnKF, because in well testing, the early time measurements give important information about the near-wellbore permeability distribution, which cannot easily be recovered from measurements at later times. In [12], we found that, although the EnKF method can give excellent results, its success is dependent on how accurately the prior described the true model and, therefore, is not always appropriate for identifying and characterizing unexpected reservoir features. The EnKF method does not always accurately match the PTT data.

In contrast, gradient-based methods for determining the maximum a posteriori estimation of parameters are able to accurately match with PTT data and are applicable even when the probability distribution is not adequately approximated by a Gaussian model. In this paper, we review how to efficiently calculate gradients using adjoint variables and then show how the priors that we have introduced are useful in improving the progress of gradient-based descent methods.

Because these methods only determine the most likely parameters, we then investigate methods to produce samples from the posterior distribution. Existing methods once again rely on the assumption of a Gaussian model for the posterior probability distribution, but in this paper, we detail a method for proper sampling from the posterior distribution by means of the Langevin equation.


The prior model is a probability density functional π0[u], which describes the likelihood of a particular set of parameters (reservoir or formation properties), u (e.g., the spatially variable porosity and components of permeability), representing the true state of the reservoir. The prior model can always be written in the form

where for commonly used Gaussian priors the function R[u] will be quadratic. The function R[u] will typically itself be described by parameters that characterize our prior knowledge of the reservoir, such as a typical correlation length or the local variance of the permeability. These parameters that describe the prior are known as “hyperparameters”, to distinguish them from the parameters of the model under analysis. Given a particular set of parameters, we need to model the likelihood of a particular set of pressure measurements ℳ being made, or the “posterior update”. The pressure measurements will never be exact due to various sources of noise, particularly when oil is being produced from the well. When the errors of each measurement are considered as independent and Gaussian, the likelihood of a particular set of measurements ℳ is

where Nw, Nt are respectively the number of wells and pressure measurements; pwm(u; tr) and wm(tr) are respectively the simulated pressures in the well (of the model described by parameters u) and the measured pressures in the well at times tr; and σ2wm(tr) is the error variance of the measurement. Using Bayes′ theorem, we are able to deduce that the posterior distribution of the parameters as

For this application, it is important that the prior is independent of the grid and that it provides a sufficient description of the small-scale spatial correlation of the parameters. At the same time, only very limited geostatistical information may be available and, thus, ideally the prior should be described by a small number of hyperparameters. Local Gaussian random fields can satisfy all of these requirements. For a parameter u(x), where x is a spatial vector, a local Gaussian random field prior is described by

where 𝓛 is a linear differential operator, such as

for arbitrary constants h0, h1, h2. This random field has a useful variogram, as shown in [13], and is best described by rewriting the linear operator as


where ρ is a correlation length, S (which can vary between −1 and 1) controls the shape of the variogram, and σ*2 is related to the local covariance σ2 by

in two dimensions, and for three-dimensional random fields by

The key advantage of employing this type of random field is that, once the problem is discretized, the inverse covariance matrix is readily obtained—it is simply the finite difference matrix for the differential operator described in Eq. (1)—and it will be sparse.

Each component of the permeability is approximately log-normally distributed, and thus,

can be modeled by a local Gaussian prior as above. The remaining components of the permeability distribution may be treated similarly, defining uy and uz. It is possible to treat the porosity in the same way,

However, such a model allows the possibility of φ > 1 and therefore can only be applied when either it is clear beforehand that the porosity is small (i.e., φ << 1) or that the variation of the porosity is small.

The prior model should not only describe the spatial correlation of the parameters, but also the local correlation between distinct parameters (e.g., permeability and porosity). Then, the complete prior model can be written as


where u = (ux, uy, uz, uφ)T and A is a 4 × 4 matrix, the inverse of the covariance matrix between the parameters. The operators 𝓛x, etc., are of the form (1), although the h0 term can be included within A. This model is described by at most 22 hyperparameters: four average values; four correlation lengths; four structure hyperparameters; and ten entries in A. However, it is easy to significantly reduce this number of hyperparameters by, for example, assuming that the correlation lengths and structure hyperparameters are all equal, and/or by assuming that the permeability is either fully isotropic or instead that there are only horizontal and vertical components of permeability, with the permeability identical in the x and y directions. If the permeability is known to be isotropic with u = log klog k, then the prior model is simply



3.1 Forward Model

In this paper, we consider pressure diffusion for a single-phase, constant-viscosity, and slightly compressible (the rate of change of density with pressure is linear) fluid in a three-dimensional porous media. In the absence of any sources or sinks, the pressure field in the reservoir is described by the diffusivity equation as


where ct is the total compressibility, K is the diagonal permeability tensor, and µ is the fluid viscosity. Note that gravity is neglected throughout, because for a weakly compressible fluid, the only effect of gravity is to introduce a hydrostatic term to the pressure, which is normally termed as potential.

We do not intend for our numerical grid to capture the geometry of the wellbore. Therefore, we apply Peaceman′s model [14] and treat the well as a line source within the grid blocks through which it penetrates. The wellbore pressure in the mth well then evolves according to


where Cwm is a property of the well and the fluid, qm is the flow rate of fluid produced from the well, m is the set of grid blocks that well m penetrates, and in the jth grid block WIj is the Peaceman well index and pj is the grid block pressure. The discretized model for the evolution of the pressure in grid block i is


where Vi is the volume of the ith grid block, Aijpj is the (symmetric) finite-volume discretization of the right-hand side of (4), and δij is the Kronecker δ function. The initial pressure in both wells and reservoir is known (although it would also be possible to include this initial pressure as an additional unknown parameter to be determined from the pressure data). The reservoir permeability enters the forward model [Eqs. (5) and (6)] via Aij and WIj .

To determine the sensitivity of the pressure transient data to changes in the reservoir parameters, we must determine the gradient of H[u]. The gradient of R[u] can be easily determined yielding a linear function of u; however, determining the gradient of I[u] is more difficult. The most direct method of determining this gradient is simply to use finite differences with small variations made to the parameter values. This approach is impractical because we would need to vary the parameters in each grid cell, leading to an extremely large number of model evaluations. Fortunately, we are able to obtain the gradient much more effectively by employing an adjoint method, which will require only a single evaluation of the forward model [Eqs. (5) and (6)] and a single evaluation of an adjoint model.

3.2 Adjoint Model

The difficulty in finding the gradient of I with respect to u arises because we are unable to find a closed form for pwm[u]; however, by introducing adjoint variables to the pressure in the grid blocks and the wells, λi, λwm, we need never explicitly calculate the gradient of pwm with respect to u. We introduce

and do not strictly enforce that pwm and pi satisfy (5) and (6), but instead note that when these equations are satisfied, I[u] = J[u, pi, pwm, λi, λwm], ∂J/∂λwm = 0, and ∂J/∂λi = 0. We can also seek the necessary conditions which λi and λwm must satisfy if we are to ensure that ∂J/∂pwm = 0 and ∂J/∂pi = 0, and find that we require




with the final conditions λi(T) = 0 and λwm(T) = 0. The adjoint problem [Eqs. (7) and (8)] is almost identical to the forward model for the pressure; however, unlike the forward problem, it must be solved backward in time, starting at a time t = T, after the last of the pressure measurements. Also in solving the adjoint problem, the rate of fluid production from the well must be replaced by the mismatch between the measured values of pressure in the wells and the values obtained by the forward model.

When the pressure and its adjoint solve (5)–(8), we have (writing Uj for the value of a generic parameter in grid block j),

and moreover, with (5) and (6) satisfied, I[u] = J[u, pi, pwm, λi, λwm] for all u, and so the derivative of I[u] is

For example, the derivative of I with respect to the porosity in grid block i is


It should be noted that this is not the same as the value in grid block i of the derivative of I with respect to the parameter, which is instead


This second representation of the gradient is more useful because it is not affected by the size of the grid block. The derivative with respect to ux, uy, uz may be obtained in a similar way, although the matrix Aij and the well indexes must be differentiated.

The approach that we have outlined thus far is applicable to determining a fully anisotropic permeability distribution and a porosity distribution. Henceforth, to aid the exposition of the problem, we restrict ourselves to a reservoir with an isotropic permeability distribution, with log-permeability denoted by u, and with known porosity. The more general methodology follows along similar lines to the simplified approach that we now describe.


Our aim in studying this inverse problem is primarily to find the true reservoir parameters; however, we must recognize that the information that we provide will almost never be sufficient to conclusively answer this question. A complete understanding of the posterior probability distribution would be the most complete answer that we could hope to give, but even when the reservoir is described by a modest number of parameters, a complete description of the posterior probability distribution is often numerically intractable.

A starting point for describing the posterior probability distribution is to find the most likely set of parameters. The simplest method for finding the most likely set of parameters is to use the steepest descent method. Crucially, the steepest descent method only uses the value of H[u] and its gradient, which, as we have described in the previous section, may be evaluated inexpensively. The steepest descent method is usually given by

where αn is a parameter to be determined and un is the estimate of the minimum after n iterations of the method. For sufficiently small values of αn, it is guaranteed that H[un+1] will be lower than H[un]; however, to make good progress with as few iterations, and therefore model simulations, as possible we need to choose αn as large as possible. A line search scheme is therefore used to determine an optimal value of αn. As noted in the previous section, when working with an irregular grid it is appropriate to use the discretization of the gradient as given by (10) rather than (9), as this expression is not affected by the size of the grid blocks.

One interpretation of the steepest descent scheme is as the explicit Euler method discretization of


A minimizer of H[u] is a stable, steady-state solution to the above equation, and thus for large values of τ, we expect this equation to find a local minimizer of H. The explicit Euler discretization of Eq. (11) is only conditionally stable. In practice, the line search ensures that the τ-step is sufficiently small to ensure stability. The maximum permissible τ-step for which the explicit Euler discretization is stable will be small whenever the problem is stiff. Because from (3)

the local Gaussian regularization that we are proposing will yield a stiff problem at this stage, becoming stiffer as the grid is refined. This stiff behavior is likely to be a generic property of all scale-independent prior models. If we are forced to take small steps, then we will have slow convergence and be required to perform many simulations of the model and its adjoint.

To avoid the restrictions of the stability criterion, we should consider implicit schemes for the discretization of (11), e.g.,

where θ lies between 0 and 1 (this scheme is unconditionally stable provided θ ≥ 1/2). However, because ∂H/∂u is both nonlinear and nonlocal, it is not easy to solve for un+1, unless θ = 0. We shall therefore treat the term in H arriving from the prior model differently from the term arriving from the update and consider the scheme


The linear part of the scheme shown in Eq. (12) is still unconditionally stable for θ ≥ 1/2. Thus, this scheme avoids the stiffness problem caused by the prior model while still allowing us to easily solve for un+1, because the update is treated explicitly. When the nonlinear and nonlocal term provides the dominant contribution toward stiffness, it may be appropriate to use a predictor-corrector scheme; however, the fast convergence of such schemes is still dependent on the use of a reasonably small step size.

The semi-implicit scheme (12) permits larger steps to be taken, although the scheme is still not guaranteed to be unconditionally stable. To determine an optimal value of αn, a line search is still required. We ensure that the line search satisfies the Wolfe conditions [15] given as


where 0 < c1 < c2 < 1 are chosen to control the accuracy of the line search. The second Wolfe condition (14) requires the derivative of the likelihood function with respect to the line-search parameter. Because we are using the scheme (12), rather than the standard steepest descent scheme, this is no longer simply given by the scalar product of the current gradient and the descent direction. Nevertheless, it is easily evaluated and given by


found by differentiating (12) with respect to α.

The steepest descent method will converge toward the most likely reservoir parameters. The adaptation to the steepest descent method (12) allows us to work with appropriately stiff priors on fine grids.

4.1 Application to an example

We seek to identify a zone of low-permeability across a single-layered reservoir, as shown in Fig. 1, with PTT data taken from an active production well and two observation wells. The true permeability distribution contains a zone of low permeability, separating one of the observation wells from the producer. The permeability distribution used as a starting point for steepest descent is homogeneous, with a permeability identical to the permeability away from the zone of low permeability. This permeability field is also used as the mode in our prior model. A drawdown test, in which fluid is produced from a well, is followed by a buildup test, where the well is shut in and no fluid is produced. Observation wells are present on either side of the band of low permeability as shown in Fig. 1. All the wells are vertical and fully completed within the reservoir. The observed pressure measurements are shown in Fig. 2, along with the pressure response from both the initial homogeneous model of the reservoir and the final model suggested from the gradient descent scheme.

FIG. 1: Plots of the permeability distribution. The first image represents the true permeability distribution, from which synthetic pressure measurements were produced. The second image is an estimate of the most likely permeability distribution consistent with the synthetic pressure measurements, as generated by the gradient descent scheme. The location of the producer well is marked with PROD, and the observation well locations are marked with OBS1/OBS2.

FIG. 2: Comparison of the observed pressure measurements (points) with the pressure response generated by the reservoir model of Fig. 1 (solid line). The pressure response of the initial homogeneous permeability distribution is shown as a dotted line.


The steepest descent scheme can only provide information about the most likely realization of the reservoir parameters; however, the information provided from well testing will never be sufficient to completely determine the reservoir properties and, thus, we should expect our result to be subject to a great deal of uncertainty. There are two broad approaches that can be applied to dealing with this uncertainty. One possibility is to describe the likelihood of realizations that are close to the most likely realization and thereby attempt to quantify the level of uncertainty in the realization.

This could be done by approximating the second derivative of the log-likelihood H[u] for the most likely realization, which can be done efficiently by making use of the gradients used as part of the minimization procedure. An alternative approach is to seek actual realizations from the posterior, and even with this approach the gradient can be used to help match the pressure transient data. We present, in this section, a method for each of these two approaches.

5.1 Quasi-Newton Methods

Quasi-Newton methods attempt to build an increasingly accurate description of the second derivative of H[u] as part of the minimization of H. This is achieved by using the many gradients of H that are calculated as part of the minimization scheme. The availability of the second derivative also improves the convergence of the scheme, beyond that of steepest descent. One of the most effective quasi-Newton schemes is L-BFGS scheme [15], an adaptation of the Broyden-Fletcher-Goldfarb-Shanno (BFGS) scheme. Both the BFGS and L-BFGS schemes use the values of the parameters and their gradients at previous iterations to approximate the Hessian matrix; however, the L-BFGS scheme only retains a limited number of the previous iterates. This allows the Hessian matrix to be represented by a low-rank update to an initial approximation.

If the updates to the parameters and the gradients are written as

then the Hessian matrix of the BFGS scheme can be represented [15] in the following compact form:

where B0 is a simple initial approximation,

and Ln and En are defined by

To obtain the L-BFGS scheme, one simply must discard the earliest entries from Sn and Yn to ensure that no more than a limited number m of updates are stored. Once an approximation of the Hessian matrix is obtained, the new reservoir parameters are found by solving

where αn should be determined by a line search. As the scheme begins to converge and the approximation to the Hessian improves, we expect that αn should approach 1. Because the L-BFGS scheme will only make a low-rank update to B0, it is sensible to include the prior covariance matrix as a term in B0, e.g.,

where Bprior is the discretized covariance matrix for the prior and χ is a parameter that can be scaled to accommodate the new information obtained from the calculation of gradients.

With a local approximation to the posterior covariance, one can try to obtain samples close to the parameters of maximum a posteriori likelihood; however, the samples produced by this approach may not be reliable for nonlinear problems, failing to accurately reproduce the data [16]; which an alternative approach, the randomized maximum likelihood method (RML) is suggested to produce more reliable samples of the posterior distribution. This method replaces the prior mode with a random sample drawn from the prior, and adds noise to the true data measurements. A sample is then drawn by finding the parameters of maximum a posteriori likelihood for this new distribution. This method seems to be a reasonable approach and will clearly produce samples that appear to honor both the prior model and the data. However, in addition to being computationally intensive to produce samples, it is only guaranteed to draw samples correctly from the posterior if it is Gaussian. We instead propose a method based on the Langevin equation in which noise is added as part of the minimization process, rather than a priori, and which, in the long term, is guaranteed to draw true samples from the posterior, regardless of its distribution.

5.2 Langevin Method

In the context of history matching of geological models, the Langevin method has previously been suggested in [13] and applied in [17, 18].We now extend the use of adjoint-derived gradients and semi-implicit schemes to the Langevin method. As noted in [13], the Langevin equation given by


where η(x, τ) represents white noise, allows one to generate samples from the posterior distribution. The distribution of u as given by (15) is governed by the Focker-Planck equation [19]. The Focker-Planck equation shows that the equilibrium distribution, achieved for large values of τ, is equal to the posterior distribution. Therefore, for large values of τ, realizations of (15) are all samples of the posterior distribution and it is possible to produce many samples from the posterior (although these will only become independent when separated by large values of τ).

One can note the similarity between the Langevin equation and the continuous form of steepest descent (11), and thus, we approximate (15) using a similar scheme to (12), which can be written as


where η(x) represents white noise only in the variable x and the √αn factor is a consequence of integration of white noise. Using a partially implicit scheme for simulating the Langevin equation has been considered before [20], although applied to a much simpler application than considered here.

The introduction of random noise in (16) prevents a line search being used to optimally choose αn. Indeed, we must allow the objective function to increase sometimes. However, it is still important that a check is made because the scheme (16) is only partially implicit.We cannot expect to generate samples if the value of αn has been chosen so large that the underlying deterministic scheme is numerically unstable.We can apply a check based on the Metropolis- Hastings algorithm, in place of a line search, thereby guaranteeing that we correctly draw samples from the posterior distribution.

The Metropolis-Hastings algorithm [21] requires that, to determine the next sample un+1, we have a method of picking a proposed sample un (the “proposal”) from a probability distribution (the “proposal distribution”) with density function Q[un, un], which may depend on the current sample un. Once a proposal has been drawn, we calculate

i.e. , the ratio between relative likelihoods of un and un in the posterior distribution and the Langevin-derived proposal distribution. If paccept > 1, then the proposal un should always be accepted as the new value for un+1; otherwise, it should only be accepted with probability paccept, and on rejection un+1 should take the value un. In the long term, the Metropolis-Hastings algorithm will produce samples from the target posterior distribution provided that the proposal distribution allows all possible values of u to be reached. However, in practice, convergence can be very slow, either because many of the proposals are rejected or because there is slow mixing (i.e., un is strongly correlated to un). In the ideal case, where the proposal distribution is identical to the posterior distribution, we always find that paccept = 1 and, thus, acceptance is guaranteed and there is perfect mixing because the proposal distribution is independent of un.

The effectiveness of the Metropolis-Hastings scheme lies in choosing a good proposal distribution for which rejections are rare and the level of mixing is high (i.e., proposals are not strongly correlated to the current s). Simple proposal distributions, such as setting un equal to un plus some random noise, often fail to simultaneously meet these two requirements.

To use the Metropolis-Hastings algorithm as a check on our results from the Langevin equation, we use (16) to define the proposal distribution as


It is reasonable to expect that, for the same rate of rejection, the use of the gradient should allow for greater mixing than is possible when the proposal only updates the previous sample with random noise. The parameter αn controls the rate of mixing, with greater mixing for large values of αn. However, because the deterministic scheme fails for large αn we should expect a high probability of rejection if αn is chosen to be too large. If the proposal is accepted, then the Metropolis-Hastings check does not require any additional expensive evaluations of the forward model or adjoint model.

We can investigate the rate of acceptance by first calculating paccept when there is no update to the posterior, i.e. , when I = 0, for which


Equation (17) shows that when we are sampling from a Gaussian distribution, the proposal generated by the Langevin method is guaranteed to be accepted provided we take θ = 1/2. This is, therefore, the optimal choice of θ in (15) and shows that this method can quickly sample from Gaussian probability distributions (although other methods will often be even more efficient). For a more abstract nonlinear Brownian bridge sampling problem considered in [20], the choice of θ = 1/2 proved crucial to ensure that the optimal step size did not tend to zero as the number of parameters to be sampled increased. As noted in [20] it is likely that such behavior holds more generically when sampling from a distribution that is a nonlinear update of a Gaussian distribution, and our experience supports using θ = 1/2 for the realization of samples conditioned to pressure transient testing data.

With θ = 1/2, we find that the acceptance probability satisfies

For small values of αn, we have

and the conditions under which each of these two terms will be positive are analogous to the first and second Wolfe conditions (13) and (14). Moreover, for a descent direction, each of these terms would be positive for sufficiently small αn in the absence of noise, and so acceptance is likely for small αn.

The step size, αn, serves as a mixing parameter with small values of αn giving slow mixing, but ensuring acceptance, and large values of αn giving enhanced mixing, but risking rejection. The efficiency of the scheme depends carefully on an optimal choice of αn to obtain a balance, avoiding both slow-mixing and frequent rejection. To date we have only employed ad hoc methods to determine αn based on the rejection likelihoods obtained.

5.2.1 Application to the Example

We have applied the Langevin method to the problem considered in Section 4.1. The initial permeability distribution used was identical to the homogeneous permeability distribution used to initialize the steepest descent method. This allows us to demonstrate the ability of the Langevin method to invert for the pressure response, although in practice it might be more efficient to initialize the Langevin method with the final state of a deterministic scheme, such as steepest descent.

The step size required to achieve a reasonable acceptance rate must be small for the first few steps at the start of the inversion procedure, but we are soon able to increase it to αn = 100 while ensuring an acceptance rate which remained consistently high. A further increase of the step size was feasible for many iterates, but the consistent application of a higher step size eventually caused the algorithm to become stuck at a certain set of parameters with all generated proposals having an extremely low acceptance probability. The importance of taking θ = 0:5, and therefore avoiding any bias in simulating the prior, was also demonstrated: significantly smaller step sizes were required to maintain the acceptance rate with θ ≠ 0:5

As can be seen in the first plot of Fig. 3, the log-likelihood function approaches its equilibrium fairly rapidly, despite the lack of an efficient method of choosing the search step. However, the rapid convergence of the objective function does not necessarily confirm stationarity of the distribution, and it can be seen from the second plot of Fig. 3 that it takes longer before the initial state is forgotten. The progress toward providing a history match and towards sampling from the stationary posterior distribution is shown in the permeability samples of Fig. 4, whereas samples of the permeability distribution from later iterates are shown in Fig. 5. The pressure responses corresponding to these samples are shown in Fig. 6.

FIG. 3: Plots of the progress of the objective function (log likelihood) and the norm of the vector of parameters against the increase in the Langevin τ parameter.

FIG. 4: Samples of the permeability distributions obtained from the initial stages of the Langevin method. In this period an increasingly good match with the measured pressure response is obtained, but it is not clear that independence from the initial state has been achieved.

FIG. 5: Samples of the permeability distributions obtained from the later stages of the Langevin method. The match with the measured pressure response is maintained, samples are more likely to be representative of the stationary distribution, although mixing is still slow with strong correlation between successive samples.

FIG. 6: Comparison of the observed pressure measurements (points) with the pressure responses generated by the reservoir models of Fig. 5 (solid lines). The pressure response of the initial homogeneous permeability distribution is shown as a dotted line.

5.2.2 Improvement of Performance of the Langevin Method

It is critical that the mixing rate is improved and/or that the work required for each proposal is reduced if the Langevin method is to be applicable for real-world problems. We have found that using a semi-implicit method allows one to increase the mixing rate while avoiding frequent rejections, and clearly, the adjoint method significantly reduces the numerical effort per proposal, but these steps alone are unlikely to be sufficient to make the method feasible in practice. In [17] the number of parameters to be determined was reduced by the use of the Karhunen-Loeve expansion to represent the prior model of the permeability distribution. In [17] the gradient was also first calculated using a coarse grid to quickly establish the likely acceptance probability of the true permeability distribution; although the benefit of this strategy is somewhat diminished here by our use of adjoint methods for determining the gradient.

Because the use of a semi-implicit method for the prior helped significantly, we might try to apply a similar idea for the nonlinear nonlocal likelihood update (the I[u] term). Although it is possible to use an explicit predictor-implicit corrector scheme to allow this term to be treated implicitly, even when we choose θ = 1/2, the proposal will not always be accepted because the objective function is not generally quadratic.

Improved methods [22] of optimally choosing the value of the search distance would also improve the mixing rate.


In this paper we have developed methodologies for estimating the likely spatial variation of reservoir parameters based on pressure transient testing data from multiple well locations. The resolution of the grid over which reservoir parameters are discretized and, in particular, the resolution of the grid near the wellbore, need not be an obstacle to the matching of these reservoir parameters to pressure transient testing data. Provided that the adjoint method is employed, there is no increase in the numerical effort required to obtain the gradient as the number of parameters increases. Some prior information is always required, and as the resolution of the grid increases, the local correlation implied by the prior model becomes important. Prior models that have a strong local correlation will typically lead to stiff problems when gradient descent schemes or quasi-Newton schemes are used to search for likely posterior parameters. The stiffness can be overcome by using a semi-implicit scheme to integrate the prior gradient.

The information provided by the pressure transient test data is usually not sufficient to eliminate the uncertainty in the estimation of the reservoir parameters. We have shown that the Langevin approach may be used to directly produce realizations of the posterior distribution. The use of a semi-implicit scheme proved to be crucial for the Langevin method to work efficiently with a large number of parameters. However, more work is needed to further improve the rate of mixing of the Langevin method.


We are grateful to Schlumberger and the Technical University of Istanbul for giving permission to publish this paper.


1. Earlougher, R., Advances in Well Test Analysis, 1st ed., Monograph Series 5, Society of Petroleum Engineers of AIME, Dallas, 1977.

2. Horner, D., Pressure buildup in wells, Proc. of 3rd Petroleum Congress, Brill, E. J., ed., Leiden, pp. 503–521, 1951.

3. Matthews, C. and Russell, D., Pressure Buildup and Flow Tests in Wells, 1st ed., Monograph Series 1, Society of Petroleum Engineers of AIME, Dallas, 1967.

4. Oliver, D. S., Estimation of radial permeability distribution from well test data, SPE Form. Eval., 7:290–296, 1992.

5. McKinley, P. M., Vela, S., and Carlton, L. A., A field application of pulse-testing for detailed reservoir description, J. Pet. Technol., 20(3):313–321, 1968.

6. Kuchuk, F. J. and Onur, M., Estimating permeability distribution from 3D interval pressure transient tests, J. Pet. Sci. Eng., 39:5–27, 2003.

7. Backus, G. E. and Gilbert, J. F., Numerical applications of a formalism for geophysical inverse problems, Geophys. J. R. Astron. Soc., 13:247–276, 1967.

8. Oliver, D. S., Multiple realizations of the permeability from well test dat, SPE J., 1(2):145–154, 1996.

9. Oliver, D. S., Reynolds, A. C., Bi, X., and Abacioglu, Y., Integration of production data into reservoir models, Pet. Geosci., 7:65–73, 2001.

10. Liu, N. and Oliver, D., Critical evaluation of the ensemble Kalman filter on history matching of geological facies, SPE J., 8(6):470–477, 2005.

11. Li, G., Han, M., Banerjee, R., and Reynolds, A., Integration of well-test pressure data into heterogeneous geological reservoir models, SPE J., 13(3):496–508, 2010.

12. Morton, K., Booth, R., Onur, M., and Kuchuk, F., Grid-based inversion methods for spatial feature identification and parameter estimation from pressure transient tests, Vienna, SPE EUROPEC, Paper No. 142996, 2011.

13. Farmer, C. L. Bayesian field theory applied to scattered data interpolation and inverse problems, Algorithms for Approximation, Iske, A. and Levesley, J., eds., Springer: Verlag, Berlin, pp. 147–166, 2007.

14. Peaceman, D., Interpretation of well-block pressures in numerical reservoir simulation with nonsquare grid blocks and anisotropic permeability, SPE J., 23(3):531–543, 1983.

15. Nocedal, J. and Wright, S. J., Numerical Optimization, 2nd ed., Springer, New York, 2006.

16. Oliver, D. S., Reynolds, A. C., and Liu, N., Inverse Theory For Petroleum Reservoir Characterization and History Matching, Cambridge University Press, Cambridge, UK, 2008.

17. Dostert, P., Efendiev, Y., Hou, T. Y., and Luo, W., Coarse-gradient Langevin algorithms for dynamic data integration and uncertainty quantification, J. Comput. Phys., 217:123–142, 2006.

18. Ma, X., Al-Harbi, M., Datta-Gupta, A., and Efendiev, Y., An efficient two-stage sampling method for uncertainty quantification in history matching geological models, SPE J., 13(1):77–87, 2008.

19. Gardiner, C. W., Handbook of Stochastic Methods, 3rd ed., Springer, New York, 2004.

20. Beskos, A., Roberts, G., Stuart, A., and Voss, J., MCMC methods for diffusion bridges, Stochas. Dyn., 3:319–350, 2008.

21. Hastings, W., Monte Carlo sampling methods using Markov chains and their applications, Biometrika, 57(1):97–109, 1970.

22. Roberts, G. O. and Rosenthal, J. S., Optimal scaling for various metropolis-hastings algorithms, Stat. Sci., 16(4):351–367, 2001.