IPTC 11119

Using Production Data to Mitigate Reservoir Connectivity Uncertainty

Hong Tang, SPE, Louisiana State University

Copyright 2007, International Petroleum Technology Conference

This paper was prepared for presentation at the International Petroleum Technology
Conference held in Dubai, U.A.E., 46 December 2007.
This paper was selected for presentation by an IPTC Programme Committee following review
of information contained in an abstract submitted by the author(s). Contents of the paper, as
presented, have not been reviewed by the International Petroleum Technology Conference
and are subject to correction by the author(s). The material, as presented, does not
necessarily reflect any position of the International Petroleum Technology Conference, its
officers, or members. Papers presented at IPTC are subject to publication review by Sponsor
Society Committees of IPTC. Electronic reproduction, distribution, or storage of any part of this
paper for commercial purposes without the written consent of the International Petroleum
Technology Conference is prohibited. Permission to reproduce in print is restricted to an
abstract of not more than 300 words; illustrations may not be copied. The abstract must
contain conspicuous acknowledgment of where and by whom the paper was presented. Write
Librarian, IPTC, P.O. Box 833836, Richardson, TX 75083-3836, U.S.A., fax 01-972-952-9435.

Abstract
Even though production data have direct responses of reservoir
heterogeneity and connectivity, they are rarely incorporated into
reservoir modeling workflow among the geological community. In
this paper, a designed simulation method is proposed to mitigate
reservoir connectivity uncertainty. This proposed method is more
accurate and efficient by integrating production data with reasonable
computational cost. This method, applied on a North Africa shallow
marine reservoir, includes following four steps.
1. A Folded Plackett-Burman design (FPBD) is used to generate
typical reservoir facies models using multipoint geostatistics method.
Then reservoirs properties (i.e. porosity, permeability and saturation)
are populated among different facies. The different scenarios of
reservoir properties represent geological uncertainties.
2. The recovery factor of each model is computed by using a
commercial black oil simulator. A polynomial response surface of
recovery factor and modeling parameters is generated as a proxy of
flow responses. A Monte Carlo Simulation Bayes Method (MCSBM)
estimates the uncertainty of recovery factors and derives the posterior
probabilities for modeling parameters.
3. The facies models are changed by using a probability
perturbation method, which reduces the error between observed
pressure data and simulation prediction. Different from traditional
methods, the objective of this step is to mitigate connectivity
uncertainty by integrating production data instead of forcing history
match.
4. Similar to step 2, another polynomial response surface between
updated recovery factor and modeling parameters is generated by
using perturbed reservoir models. MCSBM derived the conditional
probability of modeling parameters and probability distribution of
recovery factors. Both derived probabilities are compared with initial
models.
Sensitivity analysis shows that recovery factor heavy hitters of
updated models are significantly different from initial models.
MCSBM shows that the production updated modeling parameters
have smaller uncertainty ranges than the prior uncertainty ranges,
which means the production update reduces parameters uncertainty.
The updated recovery factors have wider uncertainty ranges (18
percent) compared with initial recovery factors (15 percent). This is
understandable because the updated geological models honoring
production data may be more heterogeneity and have tortuous

connectivity. The updated models help capturing the full uncertainty

ranges for recovery factor.
The updated recovery uncertainty as well as reduced modeling
parameters uncertainties help correct business decisions. Furthermore,
the updated response surfaces show confounded higher order
nonlinear effects. Therefore, further study based on more detailed
design is desired. The method in this paper honors available data
without over-tuning geological parameters for history matching. It
also provides guidance for data acquisition to mitigate production
prediction risks.

1. Introduction
The goals of this work are 1) reducing modeling parameters
uncertainty by integrating dynamic data; 2) changing the reservoir
connectivity to capture correct uncertainty ranges; 3) developing a
flexible, integrated method to model sedimentary facies connectivity.
Stochastic reservoir modeling technology provides geologists and
engineers an efficient way to integrate various types of data for
predicting inter-well reservoir properties. At the early stage of
reservoir exploration, well data are usually limited, which pose a big
challenge of inter-well prediction. Those reservoir properties are then
used to compute oil in place (OIP) on which business decisions are
usually based. In the study area, existing geological data are limited
with poor quality. One major object is then to integrate all available
data to accurately estimate the reservoir properties using accurate
modeling parameters.
Another objective of this study is to capture the connectivity
uncertainty ranges by integrating production data. The original oil in
place (OIP) is a good proxy for recoverable oil. However, this proxy
might not be true in reservoirs with low net-to-gross ratio or high
heterogeneity such as carbonate reservoirs. Reservoir connectivity
analysis provides a better alternative to estimate for recoverable oil.
The reservoir connectivity is also important for infill well planning,
enhance oil recovery, and current production optimization. However,
the production data such as well transient test data and production
history are less practically used in connectivity prediction. There are
two major directions for dynamic data integration. Direct analytical
estimation of reservoir properties (such as KH, well distance to faults
etc.) from well transient test or PLT logs is informative. Since these
data can be tested in short time, they are mostly used to estimate
average reservoir properties. However, they are seldom used to
update reservoir geological models. History match/ assisted history
match (HM/AHM) are the champion of inversely integrating dynamic
data for model updating2,4-6,9,13,15. The major limitation of most
current history match methods is that they failed to change the
reservoir model in a way consistent with geological concepts. Thus
the matched model may not be able to correctly predict the
production. Furthermore, because of non-uniqueness and
nonlinearality of solutions, the convergence of objective function
strongly depends on the initial models; wrong geological models may
either converge to local minimum or fail to converge9. Some
algorithms (such as metropolitan algorithm, Simulated Annealing, the
Tunneling Method, and Genetic Algorithms 2,9,13) help find the global
minimum of objective function. But these methods usually require

www.petroman.ir

IPTC 11119

extensive simulations runs, which are prohibitively time consuming

for complex reservoirs or immature reservoirs such as study area that
has limited well data and short production history. A practical
reservoir modeling workflow should be able to mitigate uncertainties
of reservoir parameters by integrating all available data instead of
pushing for history match.
In this paper, a recently developed Probability Perturbation Method
(PPM) is used in integrating production data by adjusting reservoir
models honoring geological concepts 1,5-7. An experimental designed
Monte Carlo Simulation Bayes Method (MCSBM) is used to detect
the factor sensitivity, quantify reservoir parameters uncertainties.
Two sets of models are compared in this paper. One is based on a
FPB Design and direct flow responses (initial model); another model
is based on the same initial model. But each model is updated by
production data using PPM method, which is called updated model.
This paper starts with an introduction of methods and proposed
workflow. It is followed by an application on a North Africa shallow
marine reservoir. Finally, the application results are discussed.

2. Methods
Probability Perturbation Method (PPM)
Probability Perturbation method (PPM) is developed to integrate
production data without creating any artifacts in geological models1,57
. PPM treats production data as a probability cube which is a soft
constraint for facies modeling and can be manipulated in a similar
manner as seismic derived probability cube. However, unlike seismic
data, production data are difficult to be converted into 3D probability
cubes because production data are only available at several well
penetrated reservoir cells. Furthermore, the dynamic data are mostly
controlled by nonlinear partial differential equations of fluid flow in
porous media, which needs to be solved globally. Based on Bayes
theorem, Caers1 introduced a simple formulation to convert
production data into conditional probability P(A|D):
P(A|D) = (1-rD) i(0)(u)+ rD P(A)
(1)
Where P(A|D) is the probability of alternative facies providing
production data D, which is a function of facies prior probability
P(A) and initial realization i(0)(u). i(.)(u) is an indicator function at
reservoir location u. The dimensionless perturbation factor rD[0,1]
controls how the initial facies are changed. i(0)(u) is adjusted to
match the production data. To understand how the formulation works,
two extremities of rD are discussed here. If rD =1, the P(A|D)=P(A),
which means the alternative facies is not informative and the prior
facies probability is honored. If rD =0, the P(A|D)= i(0)(u). The initial
alternative facies is retained in its entirety. By gradually adjusting rD,
facies realization is adjusted between two end members (initial facies
and prior probabilities) until production data are matched or
maximum iteration number is exceeded. A gradient based method
(Brent Algorithm) is used to minimize the error between observation
and simulation prediction.
The PPM algorithm has two loops for facies updating. The intra-loop
represents small adjustments of facies realization among two end
members. If satisfying match cannot be achieved by these two end
members, a more dramatic facies adjustment (another equi-probably
facies model) is added by changing random seed of facies simulation
algorithm. Based on previous study, the iteration upper limit for intraloop is 6 and outer-loop is 20. The assumption is that if the
production data cannot be approximated within limited runs, the
starting geological model is questionable and not likely to be matched.
The model that has lowest mismatch error is recorded for further
experimental design analysis.
In summary, the PPM uses a simple formulation to convert dynamic
data into reservoir properties. To reduce prediction error, two
iteration loops are used to gradually change the facies realizations.

An inner loop locally adjusts facies realizations. An outer loop

continuously adds new information till history matched.

Folded Placket-Burman Design (FPBD) and Response

Surface Method
Placket-Burman designs (PBD) are very useful economically for
detecting large main effects, assuming all interactions are negligible
when compared with the few important main effects2,3,8,12. However,
when such assumption is not satisfied, the lack of ability to handle
interactions and curvatures may give rise to significant prediction
error12. A FPBD or mirror-image fold-over design is used to increase
the resolution of PBD. It is obtained by reversing the signs of all the
columns of the original design matrix. The original design runs are
combined with the mirror-image fold-over design runs. This
combination can then be used to estimate all main effects clear of any
two-factor interaction. A FPBD is a method to build a design with
resolution for interactions from a design with resolution for main
factors only.
Selection of design is a subjective decision. It depends on the purpose
of design, number of factors, cost of design runs, and nonlinearity of
responses, For example, for 2 to 4 factors, randomized block design
is used to compare different models; full or fractional factorial
Central composite design or Box-Behnken design is used to
screening models; and Central composite design (CCD) and BoxBehnken design (BBD) is used to construct response surfaces12.
For this study, a design is used to construct response surface to
allow us to estimate interaction effects, and therefore give us an idea
of the (local) shape of the flow response surface. Four design factors
are investigated. Based on above rule of thumb, a CCD or BBD
should be used. However, considering a large number of simulation
runs may be involved, the more efficient FPBD is selected for this
study.
The complexity of recovery factors is expected to require a second
order of polynomial regression model. Due to properties of FPBD,
the curvature and higher order effects are confounded. The second
order model is expressed as

y = 0 + i xi + i j xi x j +
i

(2)

Where . is the regression coefficient. It is computed using least

square method by minimizing the prediction and observation error.

x.

is the designed factors.

y is

predicted value or response. All

curvature and higher order effects are confounded into the error
term .
The performance of Regression model is evaluated by R2, which
is a fraction of model interpreted variance vs. total variance. R2
ranges from 0 to 1. For multivariate regression, the model interpreted
variances may increase by adding variables or factors, although some
of this increase in R2 would be simply due to chance variation in that
particular sample. To remove the effects of number of variables
included, the adjusted R2 is used, which is adjusted based on the
number of variables increased or number of freedom decreased.

Monte Carlo Simulation and Bayes Theorem

Monte Carlo simulation has been widely used to propagate
uncertainty of input parameters through a performance model (i.e.
through a polynomial regression responses surface model into
uncertainty of recovery factor)14,17-19. Monte Carlo simulation
requires large number of sampling (104-105) for moderate complex
problems. Many other sampling methods such as qausi-Monte Carlo
Harmmardy sequence and Latin Hypercubes have been discussed in
literatures. These methods reduced required sampling size without
losing Monte Carlo desired properties. However, these alternative
methods need the specific sampling strategy, which is a subjective

www.petroman.ir

IPTC 11119

decision and may cause extra errors in uncertainty estimation. In this

study, because the computation cost for response surface method is
trivial, the uniform sampled MCS is used. Another concern of MCS
is that the input variables should be independent, i.e. orthogonal
assumption. The orthogonal assumption could be safely justified
because four input geological modeling parameters represent
different geological properties and there is no obvious correlation
between each other. The MCS derived flow responses can be used to
compute posterior probability distribution using Bayes theorem.
Bayes theorem is widely used due to its abilities of integrating
different data, including prior probability and inferring event
probabilities. One of the major objectives of this study is to use
production data to infer reservoir factors and subsequently mitigate
the uncertainty range of reservoir factors. For instance, in this paper,
the Bayes theorem is used to update the distribution of reservoir
modeling parameters by integrating production data17-19. The Bayes
theorem is expressed as:
P(F|R)=P(R|F)*P(F)/P(R)
(3)
The P(F|R) is called posterior probability distribution of modeling
factors conditioned to production data. P(F) is the prior probability of
modeling factors. P(R|F) is the conditional probability of responses
given factors. P(R) is the boundary probability of responses. The
prior distribution P(F) can be estimated from the sample data or
existing geological knowledge. Alternatively, if there is no or limited
information of prior distribution, a uniformly distributed prior
probability distribution could be used.
The P(R|F) is computed from MCS results. The factor values are
drawn randomly from their prior distributions. The response surface
is used to compute recovery factor given factor combination. The
number of times this joint response is observed is then tabulated into
discrete bins for each factor, which is proportional to the product of
P(R|F) and P(R)18,19. After a large number of simulations, the P(R|F)
can be computed.
Using Equation 2, the P(F|R) is computed in simple EXCEL
spreadsheets.

Proposed Workflow
The proposed workflow updates geological models by using
production data and generates dynamic response surface for the
updated models. These geological models and response surfaces of
updated models are compared with those of initial models. The
proposed workflow includes 4 relating steps.
1. Designed Geological Modeling: Based on geological
knowledge, uncertainty factors are selected. A Folded
Plackett-Burman design (FPBD) generates different
combination of modeling factors. Using these factors
combination, a Multiple-Point Statistic (MPS) method
integrate geological and geophysics data into generating
multiple initial reservoir facies models16. The training
images are based on neighboring reservoir models.
2. Flow responses of the initial geological models are
computed by using a black oil simulator. The Monte Carlo
Bayes Method (MCBM) uses a polynomial response
surface to estimate the recovery uncertainties and infer the
conditional probability for factors providing production
data.
3.
Each initial reservoir model is perturbed to matching the
production data by using the Probability Perturbation
Method (PPM). This step decreases the uncertainty of
sedimentary connectivity in a way consistent with
geological models.
4. Similar to step 2, the new response surfaces oil recovery
are recomputed based on the perturbed models. Both
responses and conditional probability distribution are
compared.

This workflow is different from traditional HM/AHM methods in

several aspects:
1) It iteratively integrates dynamic data for all designed runs instead
of limited initial models (i.e. three selected P10-50-90 (OOIP)
models). Several methods ensure such workflow is computationally
reasonable. First, an efficient FPBD is selected to represent
geological uncertainty. The FPBD also has good resolution for
interaction effects, which can describe a moderate complex dynamic
response surface. Second, because each initial model is independent
from each other, they are perfect candidates for parallel computation,
which dramatically decrease the total simulation time. Finally, the
small intra and outer loop maximum iteration number 7 and 20
guarantee the efficiency without forcing history match.
2) This workflow is able to capture the full uncertainty range of
responses. FPBD regularly samples the geological parameters and
generates diversified geological models. The designed initial
geological models ensure the perturbed final models will not
converge to local minima. As a result, the perturbed response surface
could capture the full uncertainty range for correct business decision
making.
3) The probability perturbation method adjusts facies connectivity
through the perturbation number (rD); by gradually add alternative
facies realization, all reservoir models are equally probable as the
initial models. PPM honors the geological conceptual model and
avoids traditional HM artifacts such as locally changing the reservoir
properties
4) Because the workflow only changes geological parameters, once
the initial simulation template deck is set up, a simple EXCEL
VBA script helps geologist quickly QC and update reservoir models.

3. Application on a shallow marine reservoir North

Africa
Geological Setting
Based on regional geological study and log analysis, the sediments
include a poorly sorted conglomeratic section caused by movement
along the boundary faults. It was deposited during continued marine
transgression and has a pronounced fining upward grain-size profile.
The sedimentary system is a combination of both continental and
marine systems: (1) swamp-alluvial-fan system with three lithofacies swamp facies, paleosoil facies, alluvial-fan facies; (2)
transgressive shallow marine facieswith two litho-facies:
beach/barrier shoreline facies and lagoonal/tidal flat facies.
Petrophysical analysis indicates best reservoirs are beach/barrier
facies followed by alluvial fan facies and lagoon tidal facies (Fig. 1)
(Table 1).

Flow Simulation Model Construction

Available data include 17 well log data, 2D structure maps and
reservoir isochors, which are used to build a faulted structure maps.
A general property workflow is followed to populate reservoir
properties and manage the geological uncertainties (Fig. 2). To
reproduce the geological conceptual model, which is important
reservoir with spares or poor quality data, a prior probability cube of
reservoir quality are created by using well log facies proportion and
horizontal geological maps. Multiple Point Statistic Method (MPS) is
used to populate reservoir facies models. A simplified binary facies
(reservoir and non-reservoir) are modeled. However, the extension to
multiple facies is straightforward. Moreover, uncertainty
management plan determines the four major uncertain parameters are
oil water contact, fault sealing, net to gross ratio, and tar mat. FPBD
generates multiple facies realizations. Then the petrophysics
properties of each zone are modeled separately. Empirical
relationships derived from core data are used for permeability and
saturation transformation. Production bottom-hole pressure (BHP)
data from four wells are available in each fault block. Probability
perturbation method is used to integrate them into geological model.

www.petroman.ir

IPTC 11119

The final geological model is converted into a corner point grid for
numerical simulation. The gridded model has 983,280 cells, of which
472,715 are active. The 3D grid preserves the geological details while
still keeps the model small enough to be manageable.

Flow Simulation
Flow was simulated with a commercial simulator. Four wells
naturally deplete the reservoir with reservoir pressure above the
bubble point. Four wells penetrate reservoirs with 50 feet above oil
water contact. Based on previous simulation report, the PVT, relative
permeability data are assigned and the Kv/Kh is set to 0.2. The PPM
and an Excel script drive the simulator to automatic approach the
production pressure profile. Simulation results are automatically
retrieved for response analysis.
Flow simulation fence plots of the fifth run illustrate (Fig. 3) the
change of water saturation during the depletion process. The warm
color represents high oil saturation and cold color indicates low oil
saturation. Fig. 3 (a) is the initial saturation distribution; (b) is
saturation distribution at the water breakthrough time. (c) is for
saturation after the 1 year production; (d) is after one pore volume
depletion. The depletion process shows the water finger through the
high permeability formation 2 (c). The vertical flow is impeded by
the sedimentary layering which is modeled by vertically proportional
gridding and associated Kv/Kh ratio.

Integrating Production History

The production data integration of well 4 is plotted in Fig. 4. Green
curves represent initial models; and red curves represent perturbed
facies models. This result may not be true for other reservoirs. Ideally,
facie perturbation may visually separate reservoir models into two
groups: 1) easy or possible to converged models and 2) difficult to
converged models. Since all these models are equally probable, they
are all included in uncertainty analysis in this paper. It is a subjective
choice to differentiate these two models from each other. However, it
will be reasonable to choose the models closer to the production
history for final history match. Those models could be assumed to
better represent real geology of subsurface reservoirs.
Another limitation in this study is that visually quantification of
reservoir connectivity or continuity is intuitive, further study on
quantitatively estimate the reservoir connectivity will provide more
insight on the reservoir static-dynamic relationships.
The determination of maximum outer loop iteration number is a
subjective decision, which is a trade off between total computational
cost and prediction accuracy. Too much iterations waste
computational time without improving the history match. Less
iteration runs may miss the correct geological models before reaching
the history match. The maximum iteration number of inner loop is
seven which is based on previous study6. The outer loop iteration
number is 20 for this study. So the maximum total iterations for one
model are 140 iterations. This iteration number is acceptable because
in this study reservoir is reasonably small. The average simulation
time for one run is 5 to 10 minutes.

Change of Facies Realizations and Recovery Uncertainties

The perturbation changes the reservoir facies realizations (Fig 5). The
reservoir facies is in yellow color; the blue color represents nonreservoir facies. To match the production data, though the modeling
parameters are the same, the perturbation dramatically changed the
facies realizations. The prior facies probability cube tends to
distribute reservoir facies around north paleo-high (Fig. 5 (g)). The
production data makes some of the updated facies realizations less
continuous, particularly in the southwest block where four producers
locate. The change of reservoir continuity might increase the
reservoir uncertainty ranges, which may finally increase the reservoir
recovery factor range. The predicted recovery for all designed

simulation indicates the uncertainty range increase from 15 percent to

18 percent after perturbation.

Experimental Design and Response Surface Model

In this study, four designed factors are Net-Gross-Value (NTG), Fault
Sealing (Fault), Oil Water Contact (OWC), and Tar Mat (Tar). The
recovery factor is the only response. Three dynamic responses are
initial recovery factor (RFini); updated recovery factor (RFupd) and
predicted recovery factor at 1PV depletion (RpD) (Table 2). A full
three level four factor design will have (81 runs). FPB has reduced
total runs from 81 in full factorial to 17 runs. Tar Mat acts as a flow
barrier right above the oil water contact. Both Tar mat and Fault
sealing are modeled as transmissibility multiplier. The all factors are
scaled linearly with xi = ( f i f i min ) /( f i max f i min ) where
the scaled factor,

fi

xi

is

is the unscaled factor. The flow response is the

recovery factor. Recoveries from both initial model and updated

model are included. Based on the updated model, the forecasting
recovery factor after 1 pore volume injection is also recorded. Factor
combinations are specified by a 17 Folded Plackett-Burman design
(Table 3).
Analysis of Variance
Response surface model are constructed by using polynomial
regression method11,17-19. As mentioned in equation 3, FPBD cannot
differentiate the quadratic term (or curvature). A polynomial
regression with major factors and interactions are modeled
(APPENDIX A). The regression factors are illustrated in Table.
The integration of production data changed the flow response Table
4 and Fig. 6. Several results are observed. The fact that adjusted
regression R2 adjust (0.96) means 96 percent of total variance are
interpreted by response surface model. Three out of four linear
factors are significant! This indicates the recovery factor are highly
sensitive to the change of three modeling factors in a straightforward
manner (liner relationship). However, after the integration of
dynamic data, the response surface model can only interprets 72
percent of total variance with R2 adjust (0.72). The production data
integration indicates the reservoir dynamic responses are more
complex than our initial guess. The simple response surfaces model is
not adequate to interpret the simulation responses. A higher order
nonlinear response surface is needed to represent the reservoir
behavior.
Sensitivity Analysis
The regression coefficient can be used to estimate the sensitivity of
the factors. The interpretation of the coefficient is the increment in
response when the factor increase one unit by assuming the rest
factors constant. Usually, a tornado chart is used to visually check the
uncertainties of factors (Fig. 6). In the initial model, three factors are
significant; however, only OWC is significant in the updated model.
The interaction between OWC and NTG is significant in initial model;
but not significant in updated model. Another intuitive way to
visualize the interaction terms is using 2D response surface plot (Fig
7). The response surface of NTG and OWC interaction illustrate that:
a) The production updated recovery factor is not sensitive to
the NTG and OWC interaction. However, the initial model
has strong NTG and OWC interaction.
b) In the updated model, recovery factor is only sensitive to
OWC. Raising the OWC increases the final recovery. This
could be explained by rising the OWC will have more high
permeability formations connect to aquifer that may
provide better pressure support and higher final recovery
for the reservoir.
Similar response surfaces can be plotted for other factors.
Because other interaction terms are not significant in both updated
and initial models, it will be not interesting to exam other interaction
terms.
In summary, unlike the initial model, updated model
indicates current response surface model is not adequate to explain

www.petroman.ir

IPTC 11119

the variation of dynamic responses. A more complex model and more

runs are required to fully explore the recovery responses.

Application of Monte Carlo Simulation and Bayes method

Uncertainty of Recovery Factors

The response uncertainty includes measure errors, upscaling errors,
oversimplified reservoir models, and biased sampling. Failure in
capture right recovery uncertainty may cause mistakes in business
decision. However, the high nonlinearality in flow responses makes
the uncertainty estimation difficult. Monte Carlos Simulation helps to
transfer the factor errors into responses uncertainty estimation. The
uncertain factors are drawn from specific distributions assigned by
geologist. Besides the MCS computed uncertainties, a measurement
error of recovery factors are added. The MCS computed recovery
uncertainty for both models indicates that dynamic updated model
has larger uncertainty ranges than the initial model. A student t test
assuming unequal variances indicates the two recovery distribution
are significantly different (with t value= 65 with confidence interval
>95%). Practically the uncertainty range for initial model is also
different (Fig. 8). This means the production data update is
informative and should be considered in decision analysis.
Posterior analysis
Besides analyzing the uncertainty of recovery factor and, by
combining Monte Carlo Simulation and Bayes theorem (MCSBM),
the MCSBM is used to infer the posterior probability of reservoir
factors condition to responses. For example, in this paper, production
data is used to derive the posterior probability distribution of
geological factors such as NTG and Fault. The Bayes formula shows
that the probability density function (PDF) of factors posterior
probability distribution is a function of prior probability P(F) and
condition probability of responses give factors P(R|F).
First, the initial prior probability is assumed uniformly distributed. If
the curved changes to be sharply distributed (less variation), the
MCSBM is informative. Updated NTG has a skewed distribution
based on RSM from updated model Fig 9. (a). The initial NTG and
Fault distributions are more sharply distributed Fig 9. (b). This is
because recovery is more sensitive to these factors in the initial
response surface than the updated response surface. By adding a
response error, the factor posterior probability distributions have the
similar trend as the previous models, but the distribution is less
pronouncedly distributed. In other words, the MCSBM becomes less
informative Fig 9 (c) (d). This is not the defect of this method. The
predictability reflects the data quality.
By adding more confined prior probability (triangular distributions)
to all factors, and using the same response surfaces, the factors
posterior is strongly affected by prior probability. To generation of
triangular prior distribution, a triangular probability density function
is converted to cumulative probability function (CDF) by integration.
Then the CDF converts a randomly drew probability into
corresponding factor values. Similar to PDFS using uniform priori,
MCSBM has more significant effects on the initial models Fig. 10 (b)
(d) than updated model fig. 10 (a) (c). The deviation from prior
probability distribution is more pronounced for models without
response error added fig. 10 (a) (b) than those with response error (c)
(d).
In summary, the MCSBM provides uncertainty estimation for
responses as well as modeling factors.

Acknowledgement
The Author would like to extend thanks to Dr. X. -H Wen for
providing useful guidance for this work. PPM code is provided by Dr.
Hoffman Todd. Dr. N. Liu provides insights and constructive
suggestions. Dr. B. Kurnianwan helps proof read this paper.

Nomenclature
Roman Symbols
fi,fimin,fimaxunscaledfactor
pprobability density
Pprobability
  Dimensionless perturbation factor
N recovery factor at breakthrough
x scaled factor
Greek Symbols
regressioncoefficients
uspatial coordinate
error
References
1.

2.

3.

Summary

A designed simulation method (Folded Plackett-Burman

design) and probability perturbation method (PPM) are
proposed to mitigate reservoir connectivity uncertainty
with reasonable computational cost.
A Monte Carlo Simulation and Bayes Method (MCSBM)
exam the response surfaces of recovery factors for both
initial and updated models.

A case study on a North Africa reservoir indicates the

production data updated reservoir model has different
heavy hitters for recovery factor from initial models.
MCSBM and posterior study show the modeling
parameters conditioning to production data have smaller
uncertainty ranges for both initial and updated models,
which means the production updating reduces parameters
uncertainty. The updated recovery factors have wider
uncertainty ranges (18 percent) compared with initial
recovery factors (15 percent). This is understandable
because the updated geological models honoring
production data may have wider heterogeneity ranges,
which capture the full uncertainty ranges for recovery
factor.
The updated recovery uncertainty as well as reduced
modeling parameters uncertainties help to make correct
business decision. Furthermore, the updated response
surfaces shows confounded higher order nonlinear effects.
A further study based on more detailed design is needed.
The above method honors available data without overtuning geological parameters for history matching. It also
provides more accurate prediction for risk analysis with
acceptable computational cost.

4.

www.petroman.ir

Caers, J.and Hoffman, T. 2006, The Probability

Perturbation Method: A New Look at Bayesian Inverse
Modeling, Mathematical Geology, Vol. 38, No. 1
Emanuel, A.S. and Milliken, W.J.: History Matching
Finite Difference Models With 3D streamlines, paper SPE
49000 prepared for presentation at the 1998 SPE Annual
Technical Conference and Exhibition, New Orleans, 2730
September.
Friedmann, F., Chawathe, A., Larue, D. K.: Assessing
Uncertainty in Channelized Reservoirs Using Experimental
Designs, SPE Reservoir Evaluation and Engineering,
August, 2003.
Gomez, S., Gosselin O., and Barker, J.W.: GradientBased History-Matching With a Global Optimization
Method, paper SPE 56756presented at the 1999 SPE
Annual Technical Conference and Exhibition, Houston, 3
6 October.

IPTC 11119

5.

6.

7.

8.

9.

10.

11.

12.
13.

14.

15.

16.

17.

18.

19.

20.

21.

Hoffman, B. T., and Caers, J., 2003, Geostatistical history

matching using the regional probability perturbation
method, in Proceedings of the SPE Annual Conference and
Technical Exhibition, Soc. Petrol. Eng., Paper no. 84409,
Denver, CO, October 58, 16 p.
Hoffman, B. T., Geologically Consistent history matching
while perturbing facies Ph.D. Thesis at Stanford
University, 2005.
King, P.R., Harding, T.J., and Radcliffe, N.J.:
Optimization of Production Strategies using Stochastic
Search Methods, Proc., 5th European Conference on the
Mathematics of Oil Recovery, Leoben, Austria,36
September 1996.
Journel, A., 2002, combining knowledge from diverse
sources: An alternative to traditional data independence
hypotheses: Mathematical Geology, v. 34, no. 5, p. 573
596.
Kabir, C. S., Chawathe, A., Jenkins, S. D., Olayomi, A. J.,
Aigbe, C., Faparusi, D. B.: Developing New Fields Using
Probabilistic Reservoir Forecasting, paper SPE 77564,
presented at the SPE Annual Technical Conference and
Exhibition, San Antonio, Texas, USA, 29 September 2
October, 2002.
Milliken, W.J., Emanuel, A.S., and Chakravarty, A.:
Applications of 3D Streamline Simulation to Assist
History Matching, paper SPE 63155 presented at the 2000
SPE Annual Technical Conference and Exhibition, Dallas,
14 October.
Narayanan, K., White, C. D., Lake, L. W., Willis, B. J.:
Response Surface Methods for Upscaling Heterogeneous
Geologic Models,
NIST/SEMATECH e-Handbook of Statistical Methods,
http://www.itl.nist.gov/div898/handbook/, 2004.
Liu, N. and D. S. Oliver, Evaluation of Monte Carlo
methods for assessing uncertainty, SPE 84936, SPE Journal,
8(2), 188195, 2003b.
Liu, N. and D. S. Oliver, Automatic history matching of
geologic facies, SPE 84594, Proceedings of the 2003 SPE
Annual Technical Conference and Exhibition, pp. 115,
2003a.
Ouenes, A. et al.: A New Algorithm for Automatic
History Matching: Application of Simulated Annealing
Method to Reservoir Inverse Modeling, paper SPE 26297
available from SPE, Richardson, Texas (1983).
Strebelle, S., and A. Journel, 2001, Reservoir modeling
using multiplepoint statistics: Presented at the Society of
Petroleum Engineers Annual Technical Conference and
Exhibition, SPE Paper 71324, 10 p.
Wang, F., White, C. D.: Designed Simulation for a
Detailed 3D Turbidite Reservoir Model, paper SPE 75515,
presented at the SPE Gas Technology Symposium held in
Calgary, Alberta, Canada, 30 April-2 May, 2002.
White, C.D. et al.: Identifying Controls on Reservoir
Behavior Using Designed Simulations, paper SPE 62971
presented at the 2000 SPE Annual Technical Conference
and Exhibition, Dallas, 14 October.
White, C. D., Willis, B. J., Narayanan, K., Dutton, S. P.:
Identifying and Estimating Significant Geologic
Parameters with Experimental Design, SPE Journal
(September, 2001).
Willis, B.J. and White, C.D.: Use of Quantitative Outcrop
Description for Reservoir Modeling, J. Sedimentary
Research (2000) 70, No. 4, 788.
Mantica, S. Cominelli, A. and Mantica, G.: Combining
Global and Local Optimization Techniques for Automatic

History Matching Production and Seismic Data,

Journal, June 2002

SPE

Appendix
Response Surfaces for initial and updated models
RFupdated = 0.0627 + NTG * 0.0029 + Fault * -0.00332 + OWC * 0.0042 + Tar * -0.00023 + NTG * Fault * 0.002278 + NTG * OWC *
0.00065 + NTG * Tar * -0.000829 + Fault * OWC * 0.002153 +
Fault * Tar * -0.00208 + OWC * Tar * -0.0017
RFinit = 0.0595+ NTG * 0.0021 + Fault * -0.0027 + OWC * 0.00099 + Tar * 0.00072 + NTG * Fault * -0.0006 + NTG * OWC * 0.0021 + NTG * Tar * -0.000028 + Fault * OWC * 0.00152 + Fault *
Tar * -0.000153 + OWC * Tar * 0.000347

www.petroman.ir

IPTC 11119

Table1. Sedimentary Facies Description by Zones

Formation

F1

F2
F3

Litho-Facies Desciption
Swamp alluvial fan system with swap
facies, paleosoil facies, and alluvial
facies
Transgressive shallow marine system
including beach/barrier facies
Shallow marine system including
lagoonal/tidal flat facies

Porosity (%)
Avg. Stdev.

Sw (%)
Permeability (md)
Avg.
Stdev.
Avg.
Stdev.

14.0

3.0

48.0

5.0

124.0

41.0

24.0

4.0

33.0

5.0

159.0

6.0

16.0

5.0

16.0

3.3

35.8

2.8

Note: the F1 is the youngest zone, follows by F2 and F3. F2 has highest permeability and may cause uneven oil depletion.

Table2. Experimental Designed Factors

F ac tor R anges and S c aling
F ac tor M appings
F ac tor N am e
S y m bols U nits
-1
0
1
Index
1
N et to G ros s R atio
N TG
%
0.2
0.5
0.8
2
F ault S ealing *
F ault
0
0.5
1
3
O il W ater C ontac t
OW C
ft
-10000
-10250
-10500
4
Tar m at *
Tar
0
0.5
1
* F ault s ealing and Tar m at are defined as different trans m is s ibility m ultiplier
0 m eans c om pletely s eal, 1 m eans no barrier, 0.5 m eans partially s eal

Table3. Folded Placket-Burman Design

Factors
Fault

CPR

Lower Fold

Folded PB

Upper Fold

NTG
1
1
1
-1
1
-1
-1
-1
0
1
1
1
-1
1
-1
-1
-1

OW C
-1
1
1
1
-1
1
-1
-1
0
1
1
-1
1
-1
-1
-1
1

Tar
-1
-1
1
1
1
-1
1
-1
0
1
-1
1
-1
-1
-1
1
1

1
-1
1
1
1
1
-1
-1
0
1
1
-1
-1
-1
-1
1
-1

Response
RF_ini
RF_upd
RF_1pv
0.070
0.074
0.176
0.060
0.069
0.149
0.055
0.051
0.153
0.059
0.052
0.155
0.062
0.063
0.288
0.053
0.055
0.105
0.058
0.057
0.256
0.059
0.067
0.180
0.060
0.074
0.188
0.060
0.061
0.206
0.060
0.065
0.144
0.059
0.061
0.279
0.052
0.054
0.123
0.069
0.066
0.164
0.060
0.071
0.191
0.060
0.055
0.231
0.057
0.053
0.167

Note: RF_ini=Recovery Factor of initial model; RF_upd= RF of updated model;

RF_1pv=RF after 1 pore volumn injection for updated model; CPR=Center Point
Run

Interaction Terms Linear Terms

Table 4. Sensitivity Analysis of Response Surface

Recovery Factor
Simulation Model
Initial
Updated
R2 Adjusted
0.96
0.72
10
10
Number of Regressors
Parameter Estimated
Constant
0.0595
0.0627
NTG
0.0029
0.0021
Fault
-0.0033
-0.0027
OWC
-0.0010
-0.0042
Tar
0.0007
-0.0002
NTG*Fault -0.0006
0.0023
NTG*OWC -0.0021
0.0007
NTG*Tar
0.0000
-0.0008
Fault*OWC 0.0015
0.0022
Fault*Tar
-0.0002
-0.0021
Note:

OWC*Tar
0.0003
-0.0017
Significant coefficients are shown
in Bold type, and insignificant
coefficients are italicized.

www.petroman.ir

IPTC 11119

Sedimentary System Co nceptual M o del

MARINE SEDIMENT

ACTIVE
FAULT
MARGIN
TRIASSIC

PALAEOZOIC
PASSIVE MARGIN

BASEMENT

(From Sghair and A lami)

Figure 1.The sedimentary system includes both continental and marine systems.
The petrophysical study shows the best reservoirs are shore face /barrier facies
followed by alluvial fan facies and lagoon tidal facies.

Probability Cube+MPS
Vertical (from well) and Area (geological prior)

Experimental Designed Modeling

to determine the range of reservoir values

Porosity Permeability Transformation

Using empirical correlation

Probability Perturbation Method to

Integrate Production data

Figure2. Reservoir Property Modeling Workflow

(a)

(b)

Vertical Exaggeration is 2

(c)

(d)

Figure 3. Reservoir Depletion Process Fence Plot

Warm color represents high oil saturation; cold color represents low oil saturation.
(a) initial saturation (b) saturation distribution at water breakthrough (c) Saturation
distribution at 1 year depletion (d) saturation depletion at 1 pore volume depletion.

www.petroman.ir

IPTC 11119

5000

4500

Initial Presure

History Data

4000

WBHP pressure (psi)

3500

3000

2500

2000

Updated Presure

1500

1000

500

0
0

200

400

600

800

1000

1200

1400

1600

time (days)

Figure 4. PPM integration of Production history from well 4.

Red curves are production profile from initial geological models. Green curves are PPM updated
geological models. Production history is plotted in a dotted solid line The Updated model
separate into two groups; one converges or possible to converge to history data. These models
represent the correct connectivity and are possible for further History match. Another group drifts
away from production history after perturbation, which may not correctly represent the geological
connectivity and is less possible be matched later.

Initial Facies Distribution

Updated Facies Distribution

Run5 high

80%

(a)

(b)
Run6 low

20%

(c)

(d)

Run9-Base

50%

(e)

(f)

Prior probability

(g)
Figure 5. Comparison of facies realizations before and after production update and facies prior probability
Yellow represents reservoir, blue represents non-reservoir. Left column is for initial facies distribution with different NTG;
right column represents updated facies distribution with different NTG. Bottom row is the initial facies prior probability cube.
(a) (b) are high initial/updated facies with 80 % NTG. (c) (d) are low initial/updated facies with 80 % NTG. (e) (f) are middle
initial/updated facies with 80 % NTG. The production update changes the original reservoir facies distribution and
connectivity, which may subsequently change the final recovery.

www.petroman.ir

10

IPTC 11119

Fig. 6 Plato-Chart for sensitivity analysis

The sensitivity analysis indicates updated models even using the same modeling parameters, the major
hitter is OWC (a). Rest of factors is less significant. However the major hitters for initial response surfaces
are Fault, NTG, and interaction of NTG and OWC (b). The significant factors are marked with a *.

0.068

0.068

0.066

0.066

0.064

0.064

0.062

0.062

0.06

0.06

0.058

0.058
0.5
0.44
0.38
0.32 NTG

NTG

-10250

0.2
-10300

OWC

-10350

0.29
-10500

-10250

-10350

-10300

-10400

0.38

0.052
0.05

0.2

OWC

0.47

0.054

0.26
-10450

-10500

0.052
0.05

0.056

-10450

0.054

-10400

0.056

Fig. 7 Response Surfaces

The interaction term in the updated model response surface model is not significant (a); however,
it is significant in the origin model.

Probability Density Function

0.7
0.6

Original Recovery
Uncertainty

Updated Recovery Ucertainty

0.5
0.4
0.3
0.2
0.1
0
0.05

0.055

0.06
0.065
Recovery Factor (%)

0.07

0.075

Fig. 8 Comparison of recovery factors uncertainty

Updated uncertainty of recovery factor is larger than initial recovery uncertainties. This may be
because of the change of connectivity increase the uncertainty range of connectivity.

www.petroman.ir

IPTC 11119

11

0.16

0.2

0.14

0.18

(a)

0.12

0.14

Probability Density

Probability Density

0.16

NTG

0.12

Fault

0.1

OW C

0.08
0.06

NTG

(b)
Prior

0.1

OW C

0.08

Fault

0.06
0.04

0.04

Tar
0.02

0.02

0
-1

-0.5

0.5

Tar

-1

-0.5

Factor Value

Initial factor pdf conditional to produc tion data

0.2

0.12

0.18

(c)

0.16

0.1

NTG

0.14

Probability Density

Probability Density

0.5

Factor Value

Updated factor pdf conditional to production data

NTG

0.12

Fault

0.1

OW C

0.08
0.06

0.08

Prior

OW C

0.06

Fault

0.04

0.04

Tar

0.02

(d)

Tar

0.02
0

-1

-0.5

0.5

-1

-0.5

Factor Value

0.5

Factor Value

Updated fac tor pdf including response error

Initial factor pdf inc luding response error

Figure 9.The probability density function of factors posterior probability distribution.

If the curved is sharply distributed, the MCSBM is informative. (a) Updated NTG has a skewed distribution. b)
initial NTG and Fault distributions are more sharply distributed. This is because recovery is more sensitive to these
factors in the initial response surface. (c) (d) by adding the response error, the factor posterior probability
distributions have the similar trend as the previous models, but the distributions less pronouncedly distributed.
0.3
0.3

(a)

Tar

Fault

0.15

(b)

0.25

OWC

0.2

prior
0.1

Probability Density

Probability Density

0.25

0.2

NTG
Fault
0.15

Tar

0.1

NTG
0.05

0.05

OWC
0

-1

-1

-0.5

-0.5

0.5

0.5

Factor Value

Factor Value

Updatedl factor pdf without response error

Initial factor pdf without response error

0.25
0.3

Fault

Fault
0.25

OWC

(c)

Probability Density

Probability Density

0.2

0.15

0.1

Tar
NTG

Tar

(d)

0.2

NTG
0.15

OWC
0.1

0.05
0.05

-1

-0.5

0.5

-1

Factor Value

-0.8

-0.6

-0.4

-0.2

0.2

0.4

0.6

0.8

Factor Value

Updatedl factor pdf including response error

Initial factor pdf including response error

Figure 10 the probability density function of factors posterior probability distribution.

By adding more confined prior probability (triangular distributions), the factors posterior is strongly affected by
prior probability. MCSBM skewed the sensitive factor poster from symmetric triangular distribution (the prior
distributions). Updated PDFs are skewed but less uncertain (less sharply distributed) (a). Initial PDF of Fault and
NTG are skewed with sharp distribution;

www.petroman.ir