Generating Multiple

History-Matched Reservoir-Model
Realizations Using Wavelets
Isha Sahni and Roland N. Horne, SPE, Stanford U.

Summary information, and seismic surveys may provide 3D impedance dis-

This paper focuses on an automated way to generate multiple tributions that can be inverted and used as “soft-conditioning data”
history-matched reservoir models with the inclusion of both geo- at the corresponding locations. These different sources of data can
logical uncertainty and varying levels of trust in the production be combined together with different approaches (e.g., Bayesian
data, using wavelet methods. As opposed to previously developed probability techniques) to give a single set of probabilities.
automated history-matching algorithms, this methodology not only Production data (and well-testing data), however, are of a fun-
ensures geological consistency in the final models but also in- damentally different nature and can be looked at as reservoir re-
cludes uncertainty in the production data. sponse data. If the fluid-flow model is thought of as a function/
A data distribution, such as a permeability field, can be (re- operator and the reservoir features are parameters, then production
versibly) transformed into wavelet space in which it is fully de- data would be the result of applying this function to a given input
scribed by a set of wavelet coefficients. It was found that different or stimulus (for example, a well test or oil production). The func-
subsets of the collection of wavelet coefficients can be constrained tion linking the production data and reservoir properties is based
separately to (a) the production history and (b) the geological on flow equations and simulation, which renders the task of con-
constraints. This means that the history match need be performed ditioning reservoir models to production data much harder than
only once, after which multiple realizations can be generated by direct conditioning to “hard data.” Automated history-matching
adjusting only the second subset of coefficients. algorithms usually require iterative optimization techniques to
The ability to include both geological and production-data un- match or honor production data (Carter et al. 1974; Chen et al.
certainty into the reservoir model automatically is of great conse- 1974; Jacquard and Jain 1965). A collection of production data is
quence to reservoir modeling and, hence, to reservoir management, not easy to transform into a probability distribution in the Bayesian
risk analysis, and making key economic decisions. A more com- framework. This is the reason why the integration of these data
plete and realistic reservoir model will lead to better reservoir types is difficult to achieve simultaneously (de Marsily et al. 1984;
production and development decisions. Datta-Gupta et al. 1997; Bissel 1996).
In general practice, a reservoir model is first built using all the
Introduction other information available, and the production data are then su-
Reservoir modeling is an important step in forecasting the perfor- perimposed on the existing model by way of history matching. It
mance of a reservoir, forming the basis for reservoir management, has been shown (Sahni and Horne 2005; Sahni 2003) that some
risk analysis, and making key economic decisions. A history methods of history matching might destroy (or remove) previously
match, however, is not a sufficient condition for a reservoir to integrated geologic information and/or produce artifacts that are
make better predictions for future production. The model should at nongeologic in nature (Fig. 1). The resulting reservoir models will
least conform to all the available data and the geologist’s prior then, as a result, match production data but may no longer be
conception of the reservoir. Thus, the purpose of reservoir mod- consistent with the geologic data that were integrated previously.
eling is to use all available sources of information to develop such As mentioned before, the purpose of reservoir modeling and his-
a reservoir model. This model then can be used to forecast future tory matching is not limited to building a model that is consistent
performance and optimize reservoir-management decisions. with the production data currently available, but one that gives
It is essential to integrate all the different sources of data to good predictions of its future behavior. Reservoir models that are
provide the most complete reservoir model or models (Landa and inconsistent with the geology are not likely to give good forecasts.
Horne 1997; Landa 1997; Wang 2001). Our model certainty is Therefore, it is essential to develop reservoir models that conserve
always limited by the data available to us. As such, it is never geologic information while being consistent with production-
possible to infer or develop a reservoir model with full certainty. history data at the same time.
However, the optimal use of all consistent data available will yield The methodology presented in this paper uses wavelets and
reservoir models that are less and less uncertain. Herein lies the flow simulation to interpret the production data as influencing a
significance of methodologies that can realistically and efficiently spatial distribution of wavelet coefficients, a form that is easier to
integrate different sources of reservoir information. integrate with other sources of reservoir data such as seismic or
Reservoir data are, generally speaking, divided into two cat- well logs. As a result of the above transformation, if we study the
egories: production data (e.g., pressure and water-cut histories prior permeability model in terms of its wavelet coefficients, we
from wells) and all other sources of data (e.g., core samples, seis- find that production data constrain a particular subset of the wave-
mic, and well logs). This second category of data depends on let coefficients of the given reservoir model. This is analogous to
reservoir properties like porosity and permeability in a relatively the way that hard data from cores and soft data from seismic
direct way. Core samples can be used to provide porosity and surveys constrain regions of the reservoir at different scales.
permeability measurements at specific locations (well locations); The particular advantage of using wavelets is the ability to
semivariograms (Deutsch and Journel 1998; Isaaks and Srivastava constrain multiple scales of data simultaneously. This is useful
1989) obtained from outcrops, for example, act as spatial statistics because different sources of data provide information about the
reservoir at different scales. In terms of data integration using
wavelets, this implies that different types of data will potentially
constrain sets of wavelet coefficients that describe the reservoir at
Copyright © 2006 Society of Petroleum Engineers
different resolutions (Mallat 1987). This technique is superior to
This paper (SPE 89950) was first presented at the 2004 SPE Annual Technical Conference purely pixel-based techniques (Sahni and Horne 2005; Sahni 2003)
and Exhibition, Houston, 26–29 September, and revised for publication. Original manuscript
received for review 7 June 2004. Revised manuscript received 6 March 2006. Paper peer
because, besides providing the power to change the model at the
approved 16 March 2006. highest resolution (pixel level), it also provides a more realistic

June 2006 SPE Reservoir Evaluation & Engineering 217

 variable V (in our case, the logarithms of permeabilities) is a
Gaussian random variable with mean⳱m. The variance is
␴2 = E兵V 2 其 − m2, . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (1)
and the semivariogram is ␥ (h).
Now, consider linear combinations of these random variables,
such as W, where ␻i are weights corresponding to each V.
n

W= 兺␻ V . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (2)
i=1
i i

Note that a linear combination of Gaussian random variables

yields a Gaussian random variable.

再兺 冎 兺
Fig. 1—(a) Reference permeability field; (b) history-matched
n n
model using streamline algorithm shows streamline artifacts
(from Wang 2001). E兵w其 = E ␻iVi = ␻iE兵Vi其. . . . . . . . . . . . . . . . . . . . . . (3)
i=1 i=1

higher-level support by enabling the data-integration algorithm to

constrain large-scale averages and contrasts in permeabilities di-
rectly. In other words, the technique provides more degrees of
Var兵W其 = Var 再兺 冎 i=1
n

␻iVi =
n

i=1 j=1
n

兺兺␻ ␻ Cov兵V ,V 其. . . . . . . . . (4)

i j i j

freedom that can be modified independently for the purpose of Consider

constraining to geostatistical and production data (Lu and Horne Cov兵V1,V2其 = E兵共V1 − E兵V1其兲共V2 − E兵V2其兲其. . . . . . . . . . . . . . . . . (5)
2000; Lu 2001).
In a previous paper (Sahni and Horne 2005), we outlined a Because the random function is stationary, we have
methodology using wavelets to generate multiple history-matched E兵V1其 = E兵V2其 = m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (6)
reservoir models. The algorithm dealt with the task of integrating
multiple sources of data at different scales for the purpose of Moreover, we consider an isotropic field. For this case, for a
reservoir modeling. This process was stochastic, producing a num- distance h as separation between the two random variables
ber of different, and equally likely, results. The current paper not V1⳱V(x1) and V2⳱V(x2 (that is, for |x1−x2|⳱h), we get
only provides the theoretical justification of the methodology, it Cov共h兲 = Cov兵V共x兲, V共x + h兲其 . . . . . . . . . . . . . . . . . . . . . . . . . . . (7)
also proposes a new method that is many times more efficient than
the previous one and addresses a number of issues that were not = Cov兵0其 − ␥共h兲 = ␴2 − ␥共h兲. . . . . . . . . . . . . . . . . . . . . (8)
dealt with in the previous work.
To better understand the algorithm proposed in this paper, a Now, from Eqs. 4 and 8, we can say that
summary of the formulation of the problem, as well as some key
statistical concepts, are described in the following section. Var兵W其 = Var 再兺 冎 i=1
n

␻iVi =
n

i=1 j=1
n

兺兺␻ ␻ Cov兵V ,V 其 i j i j

Problem Formulation/Theory n n n n

Logarithm Permeability Model. Consider an n×n reservoir per-

meability distribution that has already been matched to history.
= ␴2 兺兺
i=1 j=1
␻i␻j − 兺兺␻ ␻ ␥ 共h 兲. . . . . . . . . . . . . . . (9)
i=1 j=1
i j i,j

Permeability exhibits numerical behavior akin to that of a Jeffrey’s Further, for two different sets of points denoted by V1 and V 2,
parameter (Tarantola 2003) because it can take values between we can calculate the covariance between their linear combina-
zero and infinity. The proper way of evaluating the contrasts and tions using

再兺 冎
averages of such parameters is to work with their logarithm. Tak-
n m n m
ing the logarithm of this set of Jeffrey’s parameter yields Cartesian
parameters that may range anywhere from –⬁ to ⬁. This formu- Cov
i=1
␻iV1i , 兺␻ V
j=1
j
2
j = 兺兺␻ ␻ Cov兵V ,V 其. . . . . . . . (10)
i=1 j=1
i j
1
i
2
j
lation, besides being appropriate for computations involving the
permeability parameter, also lends itself very well to the Haar These equations suggest that there exists a way of calculating
wavelet transformation. This is because if we use permeability the statistics of a linear combination of random variables given the
values directly for wavelet analysis, we find that it is possible that statistics of the random variables themselves. That is, we see that
some combinations of evaluated wavelet coefficients, upon inver- the mean of a linear combination of random variables is a linear
sion, may yield negative values. It is hard to condition the sets of combination of the means of the random variables; the variance
wavelet coefficients so that they would yield only positive values and covariance of the linear combination depend on the corre-
upon inversion. Using logarithms of permeabilities for wavelet sponding variance and covariance (or semivariogram) of the ran-
analysis ensures that the wavelet inversion yields values that are all dom variable. Also note that each set of wavelet coefficients is a
valid (being within the range of –⬁ to ⬁). linear combination of the reservoir parameters (log permeabilities)
Haar wavelets are mathematical functions that allow us to rep- that are random variables. The weights (␻i ) associated with each
resent the parameter distributions at different resolutions in terms linear combination are given by the wavelet function used (the
of linear combinations of the parameter itself (for a more detailed Haar wavelet function, in our case). Thus, if we are given the
description of wavelet functions, refer to Lu and Horne 2000; Lu mean, variance, and variogram of the reservoir parameter (log
2001; Sahni and Horne 2005; Sahni 2003; and Appendix A). Given permeability) we can, under the assumption of Gaussianity, com-
the fact that the Haar wavelet coefficients are linear combinations pute the mean, variance, and variogram for each of the sets of
of the log permeabilities, we can compute the statistical parameters wavelet coefficients separately.
describing these coefficient sets with the formulation shown in the With the problem formulation and these mathematical devel-
following section. The log-permeability distribution is represented opments, we can proceed to describe the data-integration/history-
by a random function V(x) such that each location xi of the distri- matching approach proposed in this paper.
bution is composed of a random variable V(xi ).
Proposed Approach
Gaussian Variables and the Haar Wavelet Transform. Suppose As can be seen from the formulation of the Haar wavelet coeffi-
that we have a random function, V(x) composed of a (stationary) cients (Sahni and Horne 2005; Sahni 2003) (and the formulas
spatial distribution of random variables. Assume that the random listed in Appendix A), different sets of wavelet coefficients are

218 June 2006 SPE Reservoir Evaluation & Engineering

essentially different linear combinations of the original parameters
(log permeability, in our case). Thus, using the development
shown in the previous section, given the a priori statistics of the
original parameter, it is possible to evaluate the a priori statistics of
the sets of wavelet coefficients directly.
Consider a prior reservoir model (of permeabilities, for in-
stance) consisting solely of hard data and geostatistical informa-
tion in the form of a Gaussian histogram (parameterized by mean
and variance) and a variogram of the parameter. Consider also a
history-matched reservoir model that is not constrained to the sta-
tistical parameters of geology (histogram and variance). This
study shows how the latter model can be constrained to geo-
statistical information without losing the history match. In other
words, the algorithm shows that by holding a subset of the wave-
let coefficients of the model fixed and modifying the rest to honor
the geostatistical data, we can obtain, stochastically, reservoir
models that honor both geological and production-history data.
This can be done as many times as we choose without redoing the
history match.
To start with, the Haar wavelet transformation is performed on
the history-matched model, and the problem is mapped to the
wavelet domain. The subsequent step is to determine the subset of
wavelet coefficients that are most sensitive (Lu and Horne 2000;
Lu 2001) to production data. This is done by computing the gra-
dient of the production data (e.g., for pressure and water-cut pro-
files) with respect to each wavelet coefficient of the reservoir
model (Anterion et al. 1989; Chu et al. 1995]. This is the most
CPU-consuming (slowest) step of the algorithm. The higher the
sensitivity coefficient for a particular wavelet coefficient, the more
the production history depends on its magnitude. In other words,
changing a wavelet coefficient with high sensitivity will lead to a
greater deviation in the simulated production data as compared to
Fig. 2—Permeability distributions with oriented artifacts caused
changing one with a lower sensitivity. We showed previously by modifying sets of wavelet coefficients constraining only the
(Sahni and Horne 2005; Sahni 2003) that a threshold can be set on corresponding orientations.
the minimum number of wavelet coefficients required to be held
constant to provide a satisfactory history match. This set forms the
history-matching wavelet coefficients that correspond to produc- entations and the scales of wavelet coefficients in a random order
tion-data information. After the key aspects of the reservoir per- based on a random seed. The resultant reservoir models obtained
meability field (and corresponding wavelet coefficients) that in- with random traversal had relatively uniform statistics in all direc-
fluence the production match have been identified, the remaining tions (Fig. 3).
wavelet coefficients corresponding to the geological model may be As mentioned previously, the method of simulated annealing in
modified without significantly affecting the history match. The wavelet space is based on honoring statistical constraints (mean,
orthonormality of the wavelet basis function used is a key property variance, and variogram) in permeability space. Therefore, this
that makes this decoupling possible. To incorporate geostatistical method requires frequent conversions from the permeability grid
information into the model, what remains to be done is to find a to the corresponding wavelet-coefficient grid and back. Given that
way of modifying these “free” wavelet coefficients. the wavelet transform is a linear operation, the cost of frequent
The method used in our previous paper (Sahni and Horne 2005)
was based on an undirected iterative optimization technique
(simulated annealing, Gill et al. 1981) that sought to match the
variogram of the model to the prescribed model variogram. This
simulated-annealing algorithm visits each of the free wavelet-
coefficient nodes and perturbs the magnitude of the coefficient.
Because the geostatistical constraints are on the permeability val-
ues themselves, the wavelet coefficients need to be inverted back
to permeability values to compute the objective function. Based on
the change in objective function, the perturbation is either retained
or ignored. A second node is then picked and perturbed in a similar
fashion. This process, as presented in our previous paper, had an
issue with sometimes generating artifacts based on the method of
traversal of the wavelet coefficients (Sahni and Horne 2005; Sahni
2003). This problem was because it is only a smaller subset of the
free wavelet coefficients that are required to incorporate the geo-
statistical constraints into the model. Moreover, different sets of
wavelet coefficients constrain the reservoir-model statistics in dif-
ferent orientations. The a priori statistics of the reservoir model
could be honored even if the wavelet coefficients belonging to
only a particular orientation got perturbed, based on the traversal
path within the simulated-annealing module. This potentially could
lead to the production of oriented artifacts (Fig. 2) in the resultant
reservoir model obtained, whereas the reference model is thought
to be isotropic. This problem subsequently was solved by using a Fig. 3—Reservoir-model results obtained using random tra-
method of random traversal that visited the nodes across all ori- versal to avoid oriented artifacts.

June 2006 SPE Reservoir Evaluation & Engineering 219

 Fig. 4—Data integration: new methodology.

inversions is not very high. However, the wavelet coefficients are Fig. 5—Binary wavelet mask. Probability of perturbation is 0 for
perturbed randomly, and at every iteration, the objective function the blue wavelet coefficients and 1 for the red wavelet coefficients.
(i.e., the variogram) needs to be computed to check if that particu-
lar perturbation was successful. This procedure was found to be Using different random seeds for the sequential Gaussian simu-
wasteful because more than half the perturbations turned out to be lation of wavelet coefficients will yield different permeability
unsuccessful. As such, a better way of constraining to geostatistical fields, all of which will be constrained to all sources of data for the
parameters was sought. Given the fact that the Haar wavelet co- reservoir. Interestingly, because the simulations are independent
efficients of a parameter distribution are mere linear combinations from one set of wavelet coefficients to another, a combination of
of the parameters themselves, a more efficient, noniterative tech- wavelet coefficients from across the cases with different random
nique for geostatistical data integration is proposed in this paper. seeds will yield yet more permeability-field models.
Using the fact that the sets of wavelet coefficients are linear Note also that sgsim visits each node using a random path based
combinations of Gaussian parameters (log permeability), we can on a random seed. This eliminates the chance of traversal-based
compute the statistical properties of each set (mean, variance, and artifacts appearing in the results as was a possibility in the previous
variogram). The analysis of the sensitivity coefficients yields the algorithm (Sahni and Horne 2005; Sahni 2003).
set of wavelet coefficients that need to be “fixed” in order to “fix”
the production history of the model. For each set of wavelet co- Probabilistic History Matching
efficients at each scale, we can then perform sequential Gaussian In the procedures of data integration described in the previous
simulation (sgsim) using the fixed wavelet coefficients as hard data section, the production data were the first to be integrated and
and constraining to the computed statistics for that particular scale. fixed deterministically before moving on to the integration of other
Thus, we perform sgsim in wavelet space to re-evaluate the free sources of data (e.g., geostatistical data). In other words, based on
coefficients with statistics that guarantee that wavelet inversion what is thought to be good history match, a threshold of sensitivi-
will yield a distribution that is constrained by the reference geo- ties was determined, and the corresponding wavelet coefficients
statistical data. The fact that we fix the wavelet coefficients cor- were fixed. This step ensures that for all the subsequent models
responding to production data as hard data in the simulations ensures generated by integrating geology, the history match is always al-
that the history of the field is also preserved upon wavelet inversion. most exactly as good as that fixed at the thresholding stage. Im-
The result of performing sgsim in wavelet space yields, upon wavelet plicitly, we are making an assumption of perfect history informa-
inversion, the log of a permeability distribution that is constrained tion and then integrating as much geology into the model as is
to both geological and production data. Fig. 4 describes the meth- consistent with that assumption.
odology in the form of a flow chart; the emphasized box shows the In realistic situations, we are never perfectly sure of the pro-
difference between the current approach and the algorithm used duction-history data. The data are susceptible to many types of
previously (Sahni and Horne 2005; Sahni 2003). errors from different sources, from equipment malfunction and
measurement uncertainty to random noise. That being the case,
deterministically fixing a particular set of wavelet coefficients cor-
responding to the fixed production data would be incorrect because
it also means that we are limiting the integrated geological infor-
mation to be consistent with this fixed set. This is equivalent to
giving too much weight to the available production data and con-
straining the model very strongly to match history.
For a more realistic picture of the uncertainty, it is important to
consider some degree of freedom in matching the production data
as well. This uncertainty in production data can be integrated
easily into the approach by modifying the sensitivity mask that sepa-
rates the sets of wavelet coefficients constraining production data
and geology. Fixing the wavelet mask constraining the history-
matching wavelet coefficients is equivalent to assigning all the
wavelet coefficients a probability of either 0 or 1 to be perturbed
to make the model geologically consistent (Fig. 5). This is the
approach we used previously (Sahni and Horne 2005; Sahni 2003).
Uncertainty in production data can be included in the model by
replacing this black and white probability mask with a grayscale
Fig. 6—Grayscale wavelet mask. Probability of keeping a wave- mask so that each wavelet coefficient has a probability between 0
let coefficient fixed for history match may lie between 0 and 1. and 1 to be perturbed in order to match geology (Fig. 6). As can

220 June 2006 SPE Reservoir Evaluation & Engineering

 Fig. 7—Reference log-permeability field. Fig. 8—Histogram of permeabilities in reference distribution.

be seen from Fig. 6, most of the coefficients that earlier had a shows the random traversal of the simulated-annealing algorithm
probability of zero of being modified (Fig. 5) still have a very low and the number of times each node is visited. The actual reservoir-
probability, and most of the coefficients that had a probability of model results are depicted in Fig. 11 (log permeabilities) and
one of being modified still have a very high probability. In the Fig. 12. Fig. 13 shows the decrease of the objective function cor-
grayscale sensitivity mask, however, there exist wavelet coeffi- responding to the variogram match (Fig. 14). We see that the
cients that have intermediate probabilities of being perturbed to simulated production is still close to the original production-
match geology. Thus, using this methodology there is a chance of history data, and the geologic constraints are met. Being less con-
constraining different wavelet coefficients to history as well as strained than the previous models, these results are expected to
geology for each resulting reservoir model generated. show more variability in their prediction of future production. This
The grayscale approach to constraining wavelet coefficients is more realistic than assuming that the history data are perfect,
was applied to a log-permeability distribution, as shown in Fig. 7 which would give an unrealistically low value of uncertainty.
(Example 1 from Sahni and Horne 2005; Sahni 2003). Fig. 8
shows the a priori histogram for the permeabilities in the reservoir Results
model. Starting from the reference reservoir permeability model The noniterative algorithm proposed in this study was applied in
shown in Fig. 7, we calculate sensitivity coefficients in wavelet Example 1 from Sahni and Horne (2005) and Sahni (2003) for the
space and evaluate and fix those coefficients that the production purpose of comparison with the previous (iterative) algorithm (see
history is most dependent on (that is, probability of perturbation is Fig. 7 for the corresponding reference log-permeability distribu-
zero for these coefficients). Keeping the values of these coeffi- tion and Fig. 8 for the corresponding histogram). This example
cients fixed and setting the rest of the coefficients values to zero, consists of a reference 2D Gaussian permeability distribution that
we perform an inverse wavelet transform to get a permeability matches the production data of four wells. Starting with a history-
model, as depicted in Fig. 9. This permeability distribution, de- matched model that is not constrained to geological parameters
rived from the history-matched model, shows the parameters to (Fig. 9) and using an iterative algorithm (simulated annealing), it
which the production history is most sensitive—that is, permeabil- was shown in Sahni and Horne (2005) and Sahni (2003) that a
ity details are retained in the regions close to wells, whereas further number of history-matched, variogram-constrained algorithms can
away from the wells, we see that large block averages of the be generated.
parameter are constrained by production history. The nonzero We applied the new noniterative algorithm to the same example
wavelet coefficients corresponding to this reservoir model are the for comparison. Similar to the previous method, at the first stage,
input to our geostatistics integration algorithm, which modifies history-constraining wavelet coefficients were determined using a
these zero-valued coefficients to get a better representation of ge- sensitivity coefficient mask. The mean, variance, and variogram of
ology in our final reservoir models. The results from using a gray- the wavelet-coefficient sets were then computed using the property
scale sensitivity map are shown in Figs. 10 through 14. Fig. 10 that they are linear combinations of permeabilities at appropriate

Fig. 10—Random traversal showing number of visits to a par-

Fig. 9—Starting history-matched log-permeability field. ticular wavelet-coefficient node.

June 2006 SPE Reservoir Evaluation & Engineering 221

Fig. 11—Reservoir-model result (log permeabilities) using gray-
scale sensitivity coefficients.
Fig. 12—Reservoir-model result (as in Fig. 11) shown as perme-
ability in darcies. Compare with Fig. 3.

Fig. 13—Objective function vs. iteration number for matching

geostatistical parameters.
Fig. 14—Variograms for reference, history-matched model, and
the improvement of variogram match as the algorithm proceeds
from iteration to iteration.

Fig. 15—Log-permeability field: Result 1.

Fig. 16—Reservoir-model result using new approach (as in Fig.

scales (or volume of support). In the new algorithm, sgsim was 15), shown as permeabilities in darcies. Compare with Figs. 3
performed in wavelet space using these statistics, constraining to and 12.
the hard data of the fixed history-matching parameters. Wavelet
inversion of these wavelet coefficients gave a permeability reser-
voir model (Figs. 15 and 16) that matched both production data production data constant while modifying the rest to match geo-
(Figs. 17 and 18) and geostatistical constraints (mean, variance, statistical constraints. The two distributions (reference and result)
hard data, and semivariogram; see Figs. 19 and 20). Fig. 21 shows are distinct from each other in the areas in which there are no
how the algorithm keeps wavelet coefficients corresponding to the production data to constrain the reservoir models (see Fig. 22). In

222 June 2006 SPE Reservoir Evaluation & Engineering

 Fig. 17—Water-cut history: reference and Result 1.

the regions around the wells, the new result is similar to the his- Conclusions
tory-matched model, whereas far away from the wells, these real-
izations are quite different. This is a more accurate model of our A methodology for generating multiple history-matched, geologi-
uncertainty about the reservoir in the regions in which we have cally constrained reservoir models was developed. This method-
limited information. ology takes advantage of the fact that wavelet coefficients are
By changing the random seed used in sgsim, many different linear combinations of the actual reservoir parameter (log perme-
reservoir models can be obtained by inversion. All of these match ability). As a result, the statistics (mean, variance, and variogram)
both production data and geological constraints. Only a few sec- of the wavelet coefficients can be computed from the correspond-
onds of CPU time are required to generate each new model. ing statistics of the reference. Given the statistics for each wavelet

Fig. 18—Pressure history: reference and Result 1 (Curves P1 through P4).

June 2006 SPE Reservoir Evaluation & Engineering 223

 Fig. 19—Variogram match: reference and Result 1.

set and using history constraints as hard data, we can use sgsim in Fig. 20—Histogram of permeabilities in Result 1. Compare with
wavelet space to generate equiprobable reservoir models. This is a the distribution in Fig. 8.
more intuitive approach for reproducing the statistics of a distri-
bution than the iterative procedure of simulated annealing that we
used in an earlier paper. h ⳱ distance vector
Applying the proposed new approach, we generated a number h(n) ⳱ Haar wavelet scaling weights
of reservoir models using only approximately 5% of the CPU time j ⳱ scale for wavelet analysis
required for the simulated-annealing approach. This is because the k ⳱ displacement in time domain (integer value)
realizations are generated by drawing from the theoretically com- m ⳱ mean of random variable
puted statistics directly. The simulated-annealing approach was n ⳱ size of grid for wavelet decomposition
based on random perturbation, followed by wavelet inversion and p ⳱ spatial coordinates in permeability grid
checking the objective function to match the variogram. Using
u ⳱ location vector for geostatistical grid
grayscale sensitivity fields aids in better capturing the uncertainties
of production data and provides a more complete set of possible V ⳱ random variable
alternative reservoir models given the data. V(x) ⳱ random function
Var{V} ⳱ variance of random variable V
Nomenclature wi(t) ⳱ 1D Haar wavelet function
aj(p) ⳱ 2D scaling coefficient at scale j and spatial W ⳱ linear combination of random variables
location p x ⳱ location in space
C(h) ⳱ covariance for vectors separated by distance h ␥(h) ⳱ variogram as a function of distance
C(0) ⳱ variance ␴2 ⳱ variance
Cov{V1,V2} ⳱ covariance between random variables V1 and V2 ␸(t) ⳱ scaling function for 1D Haar wavelets
dj(p) ⳱ two dimensional wavelet coefficient at scale j ␺(t) ⳱ 1D Haar wavelet basis function
and spatial location p ␻i ⳱ weight corresponding to ith random variable Vi
E{V} ⳱ expected value of random variable V
f (t) ⳱ continuous function of time Subscripts
fo(t) ⳱ piecewise continuous sampling of function f (t)
i,j ⳱ indices for matrix or vector
g(n) ⳱ Haar wavelet weights

Fig. 21—Difference between wavelet coefficients of reference perme- Fig. 22—Difference between history-matched permeability dis-
ability distribution and Result 1, also showing the wavelet mask. tribution and Result 1.

224 June 2006 SPE Reservoir Evaluation & Engineering

Acknowledgments Wang, Y. 2001. Streamline Approaches for Integrating Production History
We would like to thank the Stanford Graduate Fellowship (SGF) With Geologic Information in Reservoir Models. PhD dissertation,
and the Stanford U. Petroleum Engineering Inst.—Innovations in Stanford U., Stanford, California.
Well Testing (SUPRI-D) for making this research possible. Wickerhauser, M.V. 1994. Adapted Wavelet Analysis From Theory to Soft-
ware. Wellesley, Massachusetts: A.K. Peters.
References
Appendix
Anterion, F., Eymard, R., and Karcher, B. 1989. Use of Parameter Gradi-
ents for Reservoir History Matching. Paper SPE 18433 presented at the
Haar Wavelets. The Haar wavelet transform used in this study is
SPE Symposium on Reservoir Simulation, Houston, 6–8 February.
a linear transform (Boggess and Narcowich 2001; Daubechies
1988, 1992; Jaffard and Meyer 2001; Wickerhauser 1994). We
Bissel, R. 1996. History Matching a Reservoir Model by the Positioning of
start by describing the wavelet transform in one dimension and
Geological Objects. Paper presented at the 5th European Conference on
then move on to 2D descriptions. Given the piecewise approxima-
the Mathematics of Oil Recovery, Leoben, Austria, 3–6 September.
tion fo(t) of a function f (t) taken at k (∈ integers) discrete points,
Boggess, A. and Narcowich, F.J. 2001. A First Course in Wavelets With
we define the function fj (t) as represented by blocks of width 1/2 j
Fourier Analysis. Upper Saddle River, New Jersey: Prentice-Hall.
( j being an integer scaling index) as
Carter, R.D., Kemp, L.F. Jr., Pierce, A.C., and Williams, D.L. 1974. Per-
formance Matching With Constraints. SPEJ 14 (2): 187–196; Trans.,
AIME, 257. SPE-4260-PA.
fj共t兲 = 兺a ␸共2 t − k兲, . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (A-1)
k⑀Z
k
j

Chen, W.H., Gavalas, G.R., Seinfeld, J.H., and Wasserman, M.L. 1974.
A New Algorithm for Automatic History Matching. SPEJ 14 (6): 593– where the Haar scaling function ␸(t) is given by
608; Trans., AIME, 257. SPE-4545-PA.
Chu, L., Reynolds, A.C., and Oliver, D.S. 1995. Computation of Sensitivity
Coefficients for Conditioning the Permeability Field to Well-Test Pres-
sure Data. In Situ 19 (2): 179–223.
␸共t兲 = 再 1
0
if
otherwise
0ⱕt⬍1
. . . . . . . . . . . . . . . . . . . . . . . . . . . (A-2)

Datta-Gupta, A., Vasco, D.W., and Long, J.C.S. 1997. On the Sensitivity The Haar wavelet basis function ␺(t) is written in terms of the
and Spatial Resolution of Transient Pressure and Tracer Data For Het- scaling function as
erogeneity Characterization. SPEFE 12 (2): 137–144. SPE-30589-PA. ␺共t兲 = ␸共2t兲 − ␸共2t − 1兲. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (A-3)
Daubechies, I. 1988. Orthonormal Bases of Compactly Supported Wave-
lets. Comm. On Pure and Appl. Math. 41 (7): 909–996. For j ⳱ 1, f1(t) is the piecewise continuous approximation of f (t)
Daubechies, I. 1992. Ten Lectures on Wavelets. CBMS-NSF Regional defined by blocks of thickness 1⁄2. We would like to express
Conference Series in Applied Mathematics, Soc. for Industrial and
f1共t兲 = f0共t兲 + w0共t兲, . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (A-4)
Applied Mathematics (SIAM) 61: 115–137.
de Marsily, G., Lavedan, G., Boucher, M., and Fasanino, G. 1984. Inter- where

兺a ␺共2 t − k兲 . . . a
pretation of Interference Tests in a Well Field Using Geostatistical
Techniques to Fit the Permeability Distribution in a Reservoir Model. wj共t兲 = k
j
k ∈ℜ. . . . . . . . . . . . . . . . . . . . . (A-5)
k⑀Z
In Geostatistics for Natural Resources Characterization, Part 2: 831–849.
Deutsch, C.V. and Journel, A.G. 1998. GSLIB: Geostatistical Software Li- In this way, the original sampled function f (t) can be successively
brary and Users Guide, second edition. New York City: Oxford U. Press. decomposed in terms of the unique sum
Gill, P.E., Murray W., and Wright, M.H. 1981. Practical Optimization.
New York City: Academic Press. fj共t兲 = wj−1共t兲 + wj−2共t兲 + . . . + w0共t兲 + f0共t兲. . . . . . . . . . . . . . . (A-6)
Isaaks, E. and Srivastava, M. 1989. An Introduction to Applied Geostatis- The process of expressing a function in this form is called wavelet
tics. New York City: Oxford U. Press. decomposition.
Jacquard, P. and Jain, C. 1965. Permeability Distribution From Field Pres- For a 2D matrix a, the Haar wavelet decomposition is carried
sure Data. SPEJ 5 (4): 281–294; Trans., AIME, 234. SPE-1307-PA. out similarly, as follows. The scaling coefficient at location p is
Jaffard, S. and Meyer, Y. 2001. Wavelets: Tools for Science and Technol- defined denoted by aj+1(p), and the wavelet coefficient is given by
ogy. Philadelphia, Pennsylvania: Soc. for Industrial and Applied Math- dj+1(p), where j denotes scale.
ematics (SIAM).
⬁

兺h共n − 2p兲a 共n兲 = h共0兲a 共2p兲 + h共1兲a 共2p + 1兲,

Landa, J.L. 1997. Reservoir Parameter Estimation Constrained to Pressure
Transients, Performance History and Distributed Saturation Data. PhD aj+1共p兲 = j j j (A-7)
n=−⬁
dissertation, Stanford U., Stanford, California.
Landa, J.L. and Horne, R.N. 1997. A Procedure to Integrate Well Test where
Data, Reservoir Performance History and 4D Seismic Data Into a Res-

冦
1
ervoir Description. Paper SPE 38653 presented at the SPE Annual if n = 0,1
Technical Conference and Exhibition, San Antonio, Texas, 5–8 October. h共n兲 = 公2 . . . . . . . . . . . . . . . . . . . . . . . . . . (A-8)
Lu, P. 2001. Reservoir Parameter Estimation Using Wavelet Analysis. PhD 0 otherwise
dissertation, Stanford U., Stanford, California.
Lu, P. and Horne, R.N. 2000. A Multiresolution Approach to Reservoir Substituting these coefficients, we obtain
Parameter Estimation Using Wavelet Analysis. Paper SPE 62985 pre-
1
sented at the SPE Annual Technical Conference and Exhibition, Dallas, aj+1共p兲 = 关aj共2 p兲 + aj共2 p + 1兲兴, . . . . . . . . . . . . . . . . . . (A-9)
1–4 October. 公2
Mallat, S. 1987. An Efficient Image Representation for Multiscale Analy-
sis. Paper presented at the Machine Vision Conference, Lake Tahoe, or in general,

冉公 冊 冋兺 册
February. 2 j−1
1 j
Sahni, I. 2003. Multiresolution Wavelet Analysis for Improved Reservoir aj共p兲 = a0共2 j p + i兲 . . . . . . . . . . . . . . . . . . . (A-10)
Description. MS report, Stanford U., Stanford, California. 2 i=0
Sahni, I. and Horne, R.N. 2005. Multiresolution Wavelet Analysis for
Improved Reservoir Description. SPEREE 8 (1): 53–69. SPE-87820- The corresponding wavelet coefficients are given by the expression
PA. ⬁
Tarantola, A. 2003. Inverse Problem Theory and Methods for Model Pa-
rameter Estimation. Philadelphia, Pennsylvania: Soc. for Industrial and
dj+1共p兲 = 兺g共n − 2p兲a 共n兲 = g共0兲a 共2p兲 + g共1兲a 共2p + 1兲,
n=−⬁
j j j

Applied Mathematics (SIAM). . . . . . . . . . . . . . . . . . . . . . . . (A-11)

June 2006 SPE Reservoir Evaluation & Engineering 225

 where

冦
1
− if n=0
公2
g共n兲 = 共−1兲1−n h共1 − n兲 = 1 . . . . . . . (A-12)
if n=1
公2
0 otherwise
Substituting for g(n) in Eq. A-11, we obtain:
1
dj+1共p兲 = 关−aj共2p兲 + aj共2p + 1兲兴. . . . . . . . . . . . . . . . . . (A-13)
公2
Wavelet reconstruction is the inverse process of wavelet decom-
position. The Haar wavelet reconstruction equations are derived by
solving the decomposition equations A-9 and A-13 for aj and aj+1:
1
aj共2p兲 = h共0兲aj+1共p兲 + g共0兲dj+1共p兲 = 关aj+1共p兲 − dj+1共p兲兴 .
公2
. . . . . . . . . . . . . . . . . . . . . . . (A-14)
1
aj+1共2p + 1兲 = h共1兲aj+1共p兲 + g共1兲dj+1共p兲 = 关aj+1共p兲 + dj+1共p兲兴 .
公2
. . . . . . . . . . . . . . . . . . . . . . . (A-15)
Because the wavelet transform is linear, any wavelet coefficient
aj (p) is a linear combination of the original signal a0.
Isha Sahni is a PhD student in the Dept. of Petroleum Engineer-
ing at Stanford U. She completed her MS degree in petroleum
engineering, also from Stanford U., in 2003. Sahni currently
holds the Stanford Graduate Fellowship (SGF). She won the SPE
International Student Paper Contest at the 2003 SPE Annual
Technical Conference and Exhibition and the 2001 Miller
Award for academic performance in the Petroleum Engineer-
ing Dept. at Stanford. She holds a B.Tech. degree in chemical
engineering (computer science minor) from the Indian Inst. of
Technology, New Delhi, India. Roland N. Horne is Professor and
Chairman of Petroleum Engineering at Stanford U. He holds BE,
PhD, and DSc degrees from the U. of Auckland, New Zealand,
all in engineering science. Horne has been an SPE Distin-
guished Lecturer and has been awarded the SPE Distinguished
Achievement Award for Petroleum Engineering Faculty, the
Lester C. Uren Award, and the John Franklin Carll Award. Horne
is a member of the U.S. Natl. Academy of Engineering and is
also an SPE Distinguished Member.

226 June 2006 SPE Reservoir Evaluation & Engineering