Sgems Manual

Geostatistical Earth Modeling Software:
Users Manual
Nicolas Remy
May 2004
Contents
1 General Overview 4
1.1 First Steps with GEMS . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.1.1 A quick tour to the graphical user interface . . . . . . . . . . . . . 4
1.1.2 A simple tutorial . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.1.3 Automating tasks in GEMS . . . . . . . . . . . . . . . . . . . . . 11
1.2 Dening a 3-D ellipsoid . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.3 Variogram Model Specication . . . . . . . . . . . . . . . . . . . . . . . . 15
1.4 File Formats . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.4.1 Objects Files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.4.2 Parameter Files . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
1.5 Data Set . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2 Estimation Algorithms 24
2.1 Common Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.2 Kriging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.3 Indicator Kriging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.4 CoKriging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3 Simulation Algorithms 35
3.1 Common Simulation Parameters . . . . . . . . . . . . . . . . . . . . . . . 36
3.2 SGSIM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.3 SISIM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.4 COSGSIM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
2
CONTENTS 3
3.5 COSISIM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
3.6 SNESIM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
Section 1
General Overview
GEMS, the Geostatistical Earth Modeling Software, is an example of software built from
scratch using the G
s
TL. The source code of GEMS serves as an example of how to use
G
s
TL facilities.
GEMS was designed with two aims in mind. The rst one, geared toward the end-
user, is to provide a user-friendly software which offers a large range of geostatistics tools:
the most common geostatistics algorithms are implemented, in addition to more recent
developments such as multiple-point statistics simulation. The user-friendliness of GEMS
mainly comes from its non-obtrusive graphical user interface, and the possibility to directly
visualize data sets and results in a full 3-D interactive environment. The second objective
was to design a software whose functionalities could conveniently be augmented. New fea-
tures can be added into GEMS through a system of plug-ins, i.e. pieces of software which
can not be run by themselves but complement a main software. In GEMS, plug-ins can
be used to add new (geostatistics) tools, add new grid data structures (faulted stratigraphic
grids for example) or dene new import/export lters.
1.1 First Steps with GEMS
1.1.1 A quick tour to the graphical user interface
The graphical user interface (GUI) of GEMS is divided into 3 main parts, see Fig. 1.1:
4
1.1. FIRST STEPS WITH GEMS 5
Figure 1.1: GEMSs graphical interface. The three main panels are highlighted in red. The
top-right panel is the AlgorithmPanel, the top-left is the Visualization Panel and the bottom
is the Command Panel
The Algorithm Panel The user selects in this panel which geostatistics tool to use and
inputs the required parameters (see Fig. 1.2). The top part of that panel shows a
list of available algorithms, e.g. kriging, sequential Gaussian simulation. When
an algorithm from that list is selected, a form containing the corresponding input
parameters appears below the tools list.
The Visualization Panel One or multiple objects can be displayed in this panel, e.g. a
Cartesian grid and a set of points, in an interactive 3-D environment. Visualization
options such as color-maps are also set in the Visualization Panel. The Visualization
Panel is shown in more details in Fig. 1.3.
6 SECTION 1. GENERAL OVERVIEW
Figure 1.2: The 3 parts of the Algorithm Panel highlighted in red. The top part displays the
list of available tools. The middle part is where the input parameters for the selected tool
are entered.
The Command Panel This panel is not shown by default when GEMS is started. It gives
the possibility to control the software from a command line rather that from the GUI.
It displays a history of all commands executed so far and provides an input eld
where new commands can be typed (see Fig. 1.4). See tutorial 1.1.2 for more details
about the Command Panel.
1.1.2 A simple tutorial
This short tutorial describes a GEMS session in which ordinary kriging is performed on
a 100 130 30 Cartesian grid. The data consist of a set of 400 points in 3-D space,
with a rock porosity value associated to each data point. This point-set object is called
porosity data.
Figure 1.3: The Visualization Panel. The left-hand side (highlighted in red) controls which
objects (e.g. grids) are visible. It is also used to set display options, such as which color-
map to use.
The steps involved are the following:
1. Load the data set
2. Create a Cartesian grid
3. Select the kriging tool and enter the necessary parameters
4. Display and save the result
Figure 1.4: The Command Panel
Loading the data set
The rst step is to load the data set into the objects database of GEMS. Click Objects | Load Object
and select the le containing the data. Refer to section 1.4.1 for a description of the avail-
able data le formats. When the object is loaded a new entry called porosity data
appears in the Objects section of the Visualization Panel, as shown in Fig. 1.5.
Figure 1.5: Object list after the data set is loaded
Click in the square before the point-set name to display it. Displayed objects have a little
eye painted inside the rectangle before their name. The plus sign before the square indicates
that the object contains properties. Click on the plus sign to display the list of properties.
Click in the square before the property name to paint the object with the corresponding
property (see Fig. 1.6).
Creating a grid
The next step is to create the grid on which kriging will be performed. The grid we will
create is a 3-D Cartesian grid with 100 130 30 cells.
Click Objects | New Cartesian Grid to open the grid creation dialog. Enter the di-
mensions of the grid, the coordinates of the origin of the grid, i.e. the lower left corner
The eye indicates that
the object is displayed
The eye indicates that
the property is displayed
Figure 1.6: Showing/hiding an object or a property
(i.e. the point of the grid with smallest x, y, z coordinates), and the dimensions of each grid
cell. Provide a name for the new grid, working_grid for example. Click Create Grid
to create the grid. A new entry called working_grid appears in the Objects panel of
the Visualization Panel, as shown in Fig. 1.7.
Figure 1.7: Object list after the Cartesian grid is created
The object data-base now contains two objects: a point-set with the rock porosity prop-
erty, and a Cartesian grid with not yet any property attached. We can proceed to the kriging
run.
Running the kriging algorithm
Select the kriging tool from the list in the Algorithm Panel. A form prompting for the
kriging parameters appears below the algorithms list. You can either type in the parameters
or load them from a le. Fig. 1.8 shows an example of parameter le for kriging (refer
to section 1.4.2 for a description of the parameter le format and to section 2.2 for details
on kriging parameters). Using the parameters of Fig. 1.8, ordinary kriging is performed
with an isotropic search ellipsoid of radius 50 (section 1.2 describes how to specify a 3-D
ellipsoid in GEMS) and an isotropic spherical variogram of range 30, sill 1, and a nugget
effect of 0.1 (section 1.3 explains how to specify a variogram).
<parameters> <algorithm name="kriging" />
<Grid_Name value="working_grid" />
<Property_Name value="krig_porosity" />
<Hard_Data grid="porosity data" property="porosity" />
<Kriging_Type type="Ordinary Kriging (OK)" >
<parameters />
</Kriging_Type>
<Max_Conditioning_Data value="12" />
<Search_Ellipsoid value="50 50 50
0 0 0 " />
<Variogram nugget="0.1" structures_count="1" >
<structure_1 contribution="0.9" type="Spherical" >
<ranges max="30" medium="30" min="30" />
<angles x="0" y="0" z="0" />
</structure_1>
</Variogram>
</parameters>
Figure 1.8: Kriging parameter le
Once all parameters have been entered, click the Run Algorithm button. If some pa-
rameters are not correctly set, they are highlighted in red and a description of the error will
appear if the mouse is left a few seconds on the offending parameter.
If kriging was run with the parameters shown on Fig. 1.8, the grid named working_grid
now contains a new property called krig_porosity.
Displaying and saving the result
The algorithmKriging created a newproperty krig_porosityin the grid working_grid.
Click on the plus sign before the working_grid entry in the objects list to show the list
of properties attached to the grid, and click in the square before the newly created property
to display it. To save the results, click Object | Save Object to open the Save Object dia-
log. Provide a le name, the name of the object to save, e.g. working_grid and the le
format to use (see section 1.4.1).
It is also possible to save all the objects at once by saving the project (File | Save Project).
1.1.3 Automating tasks in GEMS
Tutorial 1.1.2 showed how to perform a single run of the kriging algorithm. Next, one
would like to study the sensitivity of the algorithmto parameter Max Conditioning Data,
the maximum number of conditioning data retained for each kriging. The user would like
to vary that number from 1 to 50 in increments of 1. It would be very impractical to perform
such a study in successive sequences as explained in Tutorial 1.1.2.
GEMS provides a solution to this problem through its command line interface. Many
actions in GEMS can either be performed with mouse clicks or by typing a command in
the Command Panel. For example, loading the data set in step 1 of Tutorial 1.1.2 could
have been achieved by typing the following command:
LoadObjectFromFile /home/user/data_file.dat:All
Each command has the following format:
the name of the command, e.g. LoadObjectFromFile
a list of parameters, separated by a colon :. In the previous example two parameters
were supplied: the name of the le to load, and the le format All (meaning that
every available le format should be considered).
Every command performed in GEMS, either typed or resulting from mouse clicks, is
recorded to both the Commands History section of the Command Panel and to a le
called gstlappli history.log. Hence if one does not know a command name, one
can use the GUI to perform the corresponding action and check the command name and
syntax in the command history.
It is possible to combine several commands into a single le and have GEMS execute
them all sequentially. For the sensitivity study example, one would have to write a script
le containing 50 kriging commands, each time changing the Max Conditioning Data
parameter. Note that there are no control structures, e.g. for loops, in GEMS script les.
However there are many programming tools available (Perl, Awk, Shell,. . . ) that can be
used to generate such script le.
1.2 Dening a 3-D ellipsoid
Most of the algorithms in GEMS require the user to specify a 3-D ellipsoid, to represent
a search volume for example, or three anisotropy directions and afnity coefcients. In
GEMS a 3-D ellipsoid is represented by six parameters: the three dimensions of the ellip-
soid major axes r
max
, r
med
, r
min
, and three angles, , , positioning the ellipsoid in space
(see Figs. 1.10, 1.11 and 1.12).
Let (X, Y, Z) be the three orthogonal axes of a Cartesian coordinate system. The po-
sition of the ellipsoid is obtained by three successive rotations. All angles are measured
counter-clockwise. Initially, before any rotation, the major axis of the ellipsoid is aligned
with X, the medium axis with Y , and the minor axis with Z, see Fig. 1.9.
X
Z
Y
Figure 1.9: Starting point
1.2. DEFINING A 3-D ELLIPSOID 13
rst rotation about Z: The ellipsoid is rst rotated about axis Z, by angle , typically
called azimuth. Let (X
, Y
, Z
) be the rotated coordinate system, with Z
= Z, see
Fig. 1.10.
X
Z
Y
X'
Y'
Z'
Figure 1.10: First rotation about Z
second rotation about Y

: The ellipsoid is then rotated about axis Y

, by angle , called
the dip. Let (X
, Y
, Z
) be the rotated coordinate system, with Y

= Y
see
Fig. 1.11.
third rotation about X
: Last, the ellipsoid is rotated about axis X
by angle , called
the plunge or rake angle, leading to the nal coordinate system (X
, Y
, Z
), with
X
= X
, see Fig. 1.12.

The nal transformation matrix, product of three rotations, is given by:
T =
1 0 0
0 cos sin
0 sin cos
cos 0 sin
0 1 0
sin 0 cos
cos sin 0
sin cos 0
0 0 1

X
Z
Y
Y''=Y'
X'
Z''
Z'
X''
Figure 1.11: Second rotation about Y

Note that the order in which the rotations are performed is important. In the GEMS
parameter les, the six parameters dening an ellipsoid must be given in the following
order: r
max
, r
med
, r
min
, , , , where r
max
is the dimension of the major axis, r
med
the
dimension of the medium axis, and r
min
is the dimension of the minor axis.
Relation to GSLIBs conventions
GSLIB [Deutsch and Journel, 1992] uses a different set of three angles to position the
ellipsoid in space. They are related to GEMS angles as follows:
1
= 90
o
2
=
3
=
where
1
,
2
,
3
are the angles specied in the GSLIB parameter les, in that order.
1.3. VARIOGRAM MODEL SPECIFICATION 15
X
Z
Y
Y''
Z''
X''
Z'''
Y'''
Figure 1.12: Third rotation about X
1.3 Variogram Model Specication

Variogram models are called by most of GEMS algorithms. The current version of GEMS
allows four basic models and any positive linear combinations of them, that is the linear
model of regionalization, [Goovaerts, 1997a]. The four basic isotropic models are: the
nugget effect model, the spherical model, the exponential model and the Gaussian model.
nugget effect model
(h) =
0 if h = 0
1 otherwise
(1.1)
A pure nugget effect model for a variable Z(u) expresses a lack of (linear) depen-
dence between variables Z(u) and Z(u +h)
spherical model with actual range a
(h) =
3
2
h
a

1
2
(
h
a
)
3
if h a
1 otherwise
(1.2)
exponential model with practical range a
(h) = 1 exp
3h
a
(1.3)
Gaussian model with practical range a
(h) = 1 exp
3h
2
a
2
(1.4)
All these models are permissible in 3-D and have a covariance counterpart:
(h) = C(0) C(h)
For this reason, we will interchangeably refer to either the covariance or the variogram in
the remainder of the text.
In GEMS a variogram model: (h) = c
0
(0)
(h) +
L
l=1
c
l
(l)
(h) is characterized by
the following parameters:
a nugget effect c
0
(0)
with nugget constant c
0
0
the number L of nested structures. Each structure
(l)
(h) is then dened by:
a variance contribution c
l
0
the type of the variogram: spherical, exponential or Gaussian
six parameters: the three practical ranges and the three corresponding rotation
angles dening the anisotropy ellipsoid, see Section 1.2. Note that each nested
structure can have a different anisotropy.
Example:
Consider a variogram (h) = 0.3 + 0.4
(1)
(h) + 0.3
(2)
(h), with:
1.3. VARIOGRAM MODEL SPECIFICATION 17
a nugget constant 0.3

(1)
(h) an anisotropic spherical variogram with major range 40, medium range 20
and minor range 5, and angles = 45
o
, = 0, = 0

(2)
(h) an isotropic Gaussian variogram of range 200
Fig. 1.13 shows how this variogram model is input using the graphical interface.
Figure 1.13: Inputing using the GUI
The corresponding parameter le snippet would be:
<[Parameter_name] nugget="0.3" structures_count="2" >
<angles x="45" y="0" z="0" />
</structure_1>
<structure_2 contribution="0.3" type="Exponential" >
<angles x="0" y="0" z="0" />
</structure_2>
</[Parameter_name]>
Where [Parameter_name] is the name of the parameter. See Section 1.4.2 for more
details on the parameter le format.
1.4 File Formats
1.4.1 Objects Files
GEMS supports two le formats by default to describe grids and sets of points: the GSLIB
format and the GEMS binary format.
The GSLIB le format
It is a simple ASCII format used by the GSLIB software [Deutsch and Journel, 1992]. It is
organized by lines:
the rst line gives the title.
the second line is a single number n indicating the number of properties in the object.
the n following lines contain the names of each property (one property name per line)
each remaining line contains the values of each property (n values per line) separated
by spaces or tabulations. The order of the property values along each line is given by
the order in which the property names were entered.
1.4. FILE FORMATS 19
Note that the same format is used for describing both point sets and Cartesian grids.
When a GSLIB le is loaded into GEMS the user has to supply all the information that
are not provided in the le itself, e.g. the name of the object, the number of cells in each
direction if it is a Cartesian grid.
The GEMS binary le format
GEMS can use an uncompressed binary le format to store objects. Binary formats have
two main advantages over ASCII les: they occupy less disk space and they can be loaded
and saved faster. The drawback is a lack of portability between platforms. GEMS binary
format is self-contained, hence the user need nnot provide any additional information when
loading such a le.
1.4.2 Parameter Files
When an algorithm is selected from the Algorithm Panel (see step 3 of Tutorial 1.1.2),
several parameters are called for. It is possible to save those parameters to a le, and later
retrieve them.
The format of a parameter le in GEMS is based on the eXtended Markup Language
(XML), a standard formating language of the World Wild Web Consortium (www.w3.org).
Fig. 1.8 shows an example of such parameter le.
In a parameter le, each parameter is represented by an XML element. An element
consists of an opening and a closing tag, e.g. <tag> and </tag>, and one or several
attributes. Following is an example of an element called algorithm which contains a
single attribute name:
<algorithm name="kriging"> </algorithm>
Elements can themselves contain other elements:
<ranges max="30" medium="30" min="30"> </ranges>
<angles x="0" y="0" z="0"> </angles>
</structure_1>
</Variogram>
Here the element Variogram contains an element structure_1, which itself con-
tains two elements: ranges and angles. Each of these elements have attributes. Note
that if an element only contains attributes the closing tag can be abbreviated: in the previous
example, element ranges only contains attributes and could have been written:
The /> sequence indicates the end of the element
A GEMS parameter le always has the two following elements:
element parameters. It is the root element: it contains all other elements
element algorithm. It has a single attribute name which gives the name of the
algorithm for which parameters are specied.
All other elements are algorithm-dependent and are described in sections 2 and 3.
Such XML formated parameter le has several advantages:
Elements can be entered in any order
comments can be inserted anywhere in the le. A comment block starts with . They can span multiple lines. For example:


<Grid_Name value="working_grid" />
</parameters>
1.5. DATA SET 21
1.5 Data Set
Stanford V [Mao and Journel, 1999] is a synthetic uvial channel clastic reservoir data set.
Data on lithofacies, rock porosity, density and sonic velocity are available in the whole
reservoir. Thirteen vertical and horizontal wells, drilled from 5 plateforms, are shown on
Fig. 1.14.
The Stanford V reservoir is divided into three layers, each divided into 100x130x10
grid cells. In each layer, the channels have a different thickness and orientation. Fig. 1.15
shows some horizontal cross-sections of each layer. Fig. 1.16 shows the corresponding
porosity values, as well as the porosity histogram of each layer.
Figure 1.14: Wells and a horizontal cross-section of the Stanford V reservoir
(a) Layer 1
(b) Layer 2
(c) Layer 3
Figure 1.15: Cross-sections of each of the 3 layers showing the lithofacies
1.5. DATA SET 23
(a) Layer 1 (b) Layer 1
(c) Layer 2 (d) Layer 2
(e) Layer 3 (f) Layer 3
Figure 1.16: Porosity cross-sections and histograms of each of the 3 layers
Section 2
Estimation Algorithms
GEMS provides several tools for estimating a spatially distributed variable Z(u) from
a limited number of samples. All algorithms rely on the linear, least-square estimation
paradigm called kriging. The kriging estimator can be extended to deal with a variable
whose mean is not stationary, and to the case of multiple covariates (cokriging).
2.1 Common Parameters
All estimation algorithms in GEMS work on a 3-D object, loosely called the estimation
grid: in the current version, that object can either be a Cartesian grid or an unstructured
set of points. When an estimation algorithm is run on an object, a new property containing
the result of the estimation is attached to that object. All estimation algorithms require the
following two parameters:
Parameters description
Grid Name: The grid (or more generally, the object) on which the estimation is to be
performed
Property Name: The name of the new property resulting from the estimation
24
2.2. KRIGING 25
2.2 Kriging
Kriging is a 3-D estimation program for variables dened on a constant volume support.
Estimation can either be performed by simple kriging (SK), ordinary kriging (OK), kriging
with a polynomial trend (KT) or simple kriging with a localy varying mean (LVM). LVM is
a variant of kriging, in which it is assumed that the mean m(u) is not constant but is known
for every u. The kriging problem then becomes to nd the weights {
}, = 1, . . . , n
such that:
V ar
=1
[Z(u
) m(u
)] [Z(u) m(u)]
is minimum
Example:
Simple kriging is used to estimate porosity in the second layer of Stanford V. Porosity
is modeled by a stationary random function with an anisotropic spherical variogram. Its
expected value, inferred from the well data, is 0.17. The input parameters for algorithm
Kriging are reproduced in Fig. 2.1. Fig. 2.2 shows the estimated porosity eld.
Hard Data: The grid and the name of the property containing the hard data
Kriging Type: Possible types of kriging include: Simple Kriging (SK), Ordinary Krig-
ing (OK), Kriging with Trend (KT) and Simple Kriging with Localy Varying Mean
(LVM)
Trend: The trend components. Possible components are, with u = (x, y, z): t
1
(u) = x,
t
2
(u) = y, t
3
(u) = z, t
4
(u) = x.y, t
5
(u) = x.z, t
6
(u) = y.z, t
7
(u) = x
2
, t
8
(u) =
y
2
, t
9
(u) = z
2
.
The trend is coded by a string of 9 ags, for example the string 0 1 0 0 0 1 0 1 0
would correspond to a trend with components t
2
, t
6
and t
8
, i.e. T(u) = x + yz + y
2
.
Each ag is separated by a space.
26 SECTION 2. ESTIMATION ALGORITHMS
<Grid_Name value="layer 2" />
<Property_Name value="krig" />
<Hard_Data grid="layer 2 well data" property="porosity" />
<Kriging_Type type="Simple Kriging (SK)" >
<parameters mean="0.17" />
</Kriging_Type>
0 0 0 " />
<angles x="30" y="0" z="0" />
</structure_1>
</Variogram>
</parameters>
Figure 2.1: Simple kriging parameters
Local Mean Property: The property of the simulation grid containing the non-stationary
mean. A mean value must be available at each location to be estimated.
Max Conditioning Data: The maximum number of conditioning data to be used for
kriging at any location
Search Ellipsoid: The ranges and angles dening the search ellipsoid
Variogram: The variogram of the variable to be estimated by kriging
2.3 Indicator Kriging
Let Z(u) be a continuous random variable, and I(u, z
k
) the binary indicator function de-
ned at cutoff z
k
:
I(u, z
k
) =
1 if Z(u) z
k
0 otherwise
2.3. INDICATOR KRIGING 27
(a) 3D view
(b) Cross-section 1 (c) Cross-section 2
Figure 2.2: Porosity estimated by simple kriging
The aim of indicator kriging is to estimate the conditional cumulative distribution function
(ccdf) at any cutoff z
k
, conditional to data (n):
I
(u, z
k
) = E
I(u, z
k
) | (n)
= Prob
Z(u) < z
k
| (n)
(2.1)
I
(u, z
k
) is estimated by the simple kriging estimator:
I
(u, z
k
) E{I(u, z
k
)} =
n
=1
I(u
, z
k
) E{I(u
, z
k
)}
Estimating I
(u, z
k
) for different cutoffs z
k
, k =1,. . . ,K, yields a discrete estimate of
the conditional cumulative distribution function at the threshold values z
1
, . . . , z
K
.
The algorithm Indicator Kriging assumes that the marginal probabilities E{I(u
, z
k
)}
are constant and known for all u
and z
k
(Simple Indicator Kriging [Goovaerts, 1997a]).
Median IK
Estimating a conditional cumulative distribution function at a given location u requires
the knowledge of the variogram of each of the indicator variables I(., z
k
), k =1, . . . , K.
Inference of these K variograms can be a daunting task. Moreover a kriging system has
to be solved for each indicator variable. The inference and computational burden can be
alleviated if two conditions are met:
the K indicator variables are intrinsically correlated [Goovaerts, 1997a]:
Z
(h) =
I
(h, z
k
) =
I
(h, z
k
, z
k
) z
k
, z
k
where
Z
(h) is the variogramof variable Z(u) and
I
(h, z
k
, z
k
) is the cross-variogram
between indicator variables I(., z
k
) and I(., z
k
). All these variograms are standard-
ized to a unit sill
All vectors of hard indicator data are complete: there are no missing values as could
result from inequality constraints on the Z value.
When those two conditions are met, it is only necessary to infer one variogram
Z
(h) and
only one kriging system has to be solved for all thresholds z
k
(indeed, the data locations
and the variogram are the same for all thresholds). This simplication is called median
indicator kriging.
2.3. INDICATOR KRIGING 29
Ensuring the validity of the estimated cdf
Since kriging is a non-convex estimator, the cdf values estimated by kriging Prob
Z(u) < z
k
,
k = 1, . . . , K estimated by indicator kriging may not satisfy the properties of a cdf F, that
is:
k {1, ..., K} F(z
k
) = Prob(Z(u) < z
k
) [0, 1] (2.2)
z
k
z
k
F(z
k
) F(z
k
) (2.3)
If properties 2.2 and 2.3 are not veried, the program Indicator Kriging modies the
kriging values so as to ensure the validity of the ccdf. The corrections are performed in 3
steps.
1. Upward correction:
Loop through all the cutoffs z
1
,. . . , z
K
, starting with the lowest cutoff z
1
.
At cutoff z
k
, if the estimated probability I
(Z(u), z
k
) is not in [I
(Z(u), z
k1
), 1],
reset it to the closest bound
2. Downward correction:
Loop through all the cutoffs z
1
,. . . , z
K
, starting with the highest cutoff z
K
.
At cutoff z
k
, if the estimated probability I
(Z(u), z
k
) is not in [0, I
(Z(u), z
k+1
)],
reset it to the closest bound
3. The nal corrected probability values are the average of the upward and downward
corrections
Categorical Variables
Indicator Kriging can also be applied to categorical variables, i.e. variables that take a
discrete, nite, number of values (also called classes or categories): z(u) {0, ..., K1}.
The indicator variable for class k is then dened as:
I(u, k) =
1 if Z(u) = k
0 otherwise
and the probability I
(u, k) of Z(u) belonging to class k is estimated by simple kriging:

I
(u, k) E{I(u, k)} =

n
=1
I(u
, k) E{I(u
, k)}
In the case of categorical variables, the estimated probabilities must al be in [0, 1] and
verify:
K
k=1
P
Z(u = k) | (n)
= 1 (2.4)
If not, they are corrected as follows:
1. If I
(Z(u), z
k
) / [0, 1] reset it to the closest bound. If all the probability values are
lesser or equal to 0, no correction is made and a warning is issued.
2. Standardize the values so that they sum-up to 1:
I
corrected
(Z(u), z
k
) =
I
(Z(u), z
k
)
K
i=1
I
(Z(u), z
i
)
Hard Data Grid: The grid containing the hard data
Hard Data Property: The name of the properties containing the hard data. If there
are K indicators, K properties must be specied. Each property contains the indica-
tor values I(u, k), k = 1,. . . ,K
Categorical Variable Flag: A ag indicating whether the variable to be kriged
is categorical or continuous. If the ag is on (equal to 1), the variable is assumed
categorical
Marginal Probabilities: The marginal probabilities of each indicator. The prob-
ability values are separated by one or more space. The rst probability corresponds
to the probability of the rst indicator and so on.
Max Conditioning Data: The maximum number of conditioning data to be used for
each kriging
2.4. COKRIGING 31
Median Ik Flag: A ag indicating whether median IK should be performed
Full Ik Flag: Aag indicating whether full IKshould be performed. If Median Ik Flag
is on/off, Full Ik Flag is off/on
Variogram Median Ik: The variogram of the indicator variables used for median IK.
If Median Ik Flag is off this parameter is ignored
Variogram Full Ik: The variograms of each indicator variable. If Median Ik Flag
is on, this parameter is ignored
2.4 CoKriging
The aimof cokriging is to estimate a variable Z
1
(u) accounting for data on related variables
Z
2
, ..., Z
J+1
. The cokriging estimator is given by:
Z
1
(u) m
1
(u) =
Z
1
(u
) m
1
(u
+
J+1
j=2
N
j
Z
j
(u
j
) m
j
(u
j
)
(2.5)
Algorithm CoKriging allows to distinguish between the two cases of simple and ordi-
nary cokriging:
the means m
1
(u), ..., m
J+1
(u) are known and constant: simple cokriging
the means are locally constant but unknown. The weights
and
j
must then
satisfy:
= 1 (2.6)
N
j
j
= 0 j {1, ..., J} (2.7)
Solving the cokriging system calls for the inference and modeling of multiple vari-
ograms: the variogram of each variable and all cross-variograms between any two vari-
ables. In order to alleviate the burden of modeling all these variograms, two models
have been proposed: the Markov Model 1 and the Markov Model 2 [Almeida and Jour-
nel, 1994; Goovaerts, 1997b]
For text clarity, we will only consider the case of a single secondary variable (J = 1).
Markov Model 1
The Markov Model 1 (MM1) considers the following Markov-type screening hypothesis:
E
Z
2
(u) | Z
1
(u), Z
1
(u +h)
= E
Z
2
(u) | Z
1
(u)
(2.8)
i.e. the dependence of the secondary variable on the primary is limited to the co-located
primary variable.
The cross-variogram (or cross-covariance) is then proportional to the auto-variogram
of the primary variable [Almeida and Journel, 1994; Goovaerts, 1997b]:
C
12
(h) =
C
12
(0)
C
11
(0)
C
11
(h) (2.9)
Where C
12
is the covariance between Z
1
and Z
2
and C
11
is the covariance of Z
1
. Solv-
ing the cokriging system under the MM1 model only calls for knowledge of C
11
, hence the
inference and modeling effort is the same as for univariate kriging.
Although very congenial the MM1 model should not be used when the support of the
secondary variable Z
2
is larger than the one of Z
1
, lest the variance of Z
1
would be under-
estimated. It is better to use the Markov Model 2 in that case.
Markov Model 2
The Markov Model 2 (MM2) was developed for the case where the volume support of the
secondary variable is larger than that of the primary variable [Journel, 1999]. This is often
the case with remote sensing and seismic-related data. The Markov-type hypothesis is now:
E
Z
1
(u) | Z
2
(u), Z
2
(u +h)
= E
Z
1
(u) | Z
2
(u)
(2.10)
2.4. COKRIGING 33
i.e. the dependence of the primary variable on the secondary is limited to the co-located
secondary variable.
The cross-variogram is now proportional to the variogram of the secondary variable:
C
12
(h) =
C
12
(0)
C
11
(0)
C
22
(h) (2.11)
In order for all three covariances C
11
, C
12
and C
22
to be consistent, Journel (1999) pro-
posed to model C
11
as a linear combination of C
22
and any permissible residual covariance
C
R
. Expressed in term of correlograms (C
11
(h) = C
11
(0)
11
(h)), this is written:
11
(h) =
2
12

22
(h) + (1
2
12
)
R
(h) (2.12)
Primary Harddata Grid: The grid containing the primary hard data
Primary Variable: The name of the property containing the primary hard data
Assign Hard Data: A ag specifying whether the hard data should be re-located on
the simulation grid.
Secondary Harddata Grid: The grid containing the secondary hard data
Secondary Variable: The name of the property containing the secondary hard data
Kriging Type: Possible types of cokriging include: Simple Kriging (SK) and Ordi-
nary Kriging (OK). Note that selecting OK while retaining only a single secondary
datum (e.g. when doing co-located cokriging) amounts to completely ignoring the
secondary information : from Eq. (2.7), the weight associated to the single secondary
datum will be 0.
SK Means: The mean of the primary and secondary variables. If the selected kriging
type is Ordinary Kriging, this parameter is ignored
Cokriging Type: The cokriging option. Available options are: Full Cokriging, Co-
located Cokriging with Markov Model 1 (MM1), and Co-located Cokriging with
Markov Model 2 (MM2)
Max Conditioning Data: The maximum number of primary conditioning data to be
used for each cokriging
Search Ellipsoid 1: The ranges and angles dening the search ellipsoid used to
nd the primary conditioning data
Variogram C11: The variogram model for the primary variable
Max Conditioning Data 2: The maximum number of secondary conditioning data
used for each cokriging. This parameter is only required if Cokriging Type is
set to Full Cokriging
nd the secondary conditioning data. This parameter is only required if Cokriging Type
is set to Full Cokriging
Variogram C12: The cross-variogram between the primary and secondary variable.
This parameter is only required if Cokriging Type is set to Full Cokriging
Variogram C22: The variogram of the secondary variable. This parameter is only re-
quired if Cokriging Type is set to Full Cokriging
Correl Z1Z2: The coefcient of correlation between the primary and secondary vari-
able. This parameter is only required if Cokriging Type is set to MM1
MM2 Correl Z1Z2: The coefcient of correlation between the primary and secondary
variable. This parameter is only required if Cokriging Type is set to MM2
MM2 Variogram C22: The variogramof the secondary variable. This parameter is only
required if Cokriging Type is set to MM2. Note that the variogram of the sec-
ondary variable and the variogram of the primary variable must verify the condi-
tion (2.12)
Section 3
Simulation Algorithms
Stochastic simulations provides a depiction and a measure of uncertainty of the spatial vari-
ability of a phenomenon. This is done by generating multiple realizations of the stochastic
process modeling the spatial distribution under study. Several simulation techniques have
been developed, which can be classied into four categories: sequential simulation, p-eld
simulation, object simulation and optimization-based techniques.
All the simulation tools implemented in GEMS belong to the sequential simulation
category. Sequential simulation is a broad class of algorithms which can be used to solve a
very large spectrum of problems. It relies on the property that the joint spatial conditional
distribution of a variable Z(u) at several (N) locations in space can be decomposed into
the following product:
F
u
1
, . . . , u
N
; z
1
, . . . , z
N
| (n)
= F(u
N
; z
N
| (n + N1))
F(u
N1
; z
N1
| (n + N2))
.
.
.
F(u
2
; z
2
| (n + 1)) F(u
1
; z
1
| (n))
(3.1)
where F(u
j
; z
j
| (n + j)) is the probability that any single variable Z(u
j
) z
j
given n + j
conditioning data (n + j), j = 1, . . . , N. Hence to sample from the joint distribution
F
u
1
, . . . , u
N
; z
1
, . . . , z
N
| (n)
one could sequentially sample from N single-point con-

ditional distributions F(u
j
; z
j
| (n + j)), j = 1, . . . , N, each sampled value becoming a
conditioning datum for the next ccdf.
35
36 SECTION 3. SIMULATION ALGORITHMS
Notice that in Eq. 3.1 each ccdf F(u
j
; z
j
| (n + j)) involves all (n +j) data. Condition-
ing to all (n + j) data can be daunting and it is a common assumption that F(u
j
; z
j
| (n +
j)) = F(u
j
; z
j
| (n
+ j
)), where (n
+ j
) is a subset of the (n + j) data found in a neigh-

borhood of u
j
. This Markov-type assumption is often justied by the screening effect of
remote data by the closest (n
+ j
) data retained.
Another point worth noticing is that in order to actually sample fromF
u
1
, . . . , u
N
; z
1
, . . . , z
N
| (n)
all the conditional cumulative distribution functions in decomposition 3.1 must be exactly
known. That is generally not the case.
3.1 Common Simulation Parameters
The following parameters are common to all simulation algorithms:
Grid Name: The grid on which the simulation algorithm is to be performed
Property Name Sim: Each realization is stored as a new property attached to the ob-
ject specied by Grid Name. The properties are called Property Name Sim
followed by a unique index
Nb Realizations: The number of realizations to be generated
Seed: The seed of the random number generator used for the stochastic simulation
3.2 SGSIM
Let Z(u) be a multivariate Gaussian random function with 0 mean, unit variance, and a
given covariance model. Realizations of Z, conditioned to data (n) can be generated by
Algorithm 1.
It can be shown [Goovaerts, 1997a] that the resulting realizations will reproduce both
the rst and second moments (covariance) of Z. Note that the theory requires the variance
of each ccdf be estimated by the simple kriging variance. Algorithm SGSIM can account
3.2. SGSIM 37
Algorithm 1 Sequential Gaussian Simulation
1. Dene a random path visiting each node of the grid once
2. At each node u, the local conditional cumulative distribution function is Gaussian,
its mean is estimated by simple kriging and its variance is equal to the simple kriging
variance. The conditioning data consist of both neighboring original data (n) and
previously simulated values
3. Draw a value from that Gaussian ccdf and add the simulated value to the data set
4. Proceed to the next node in the path and repeat the two previous steps until all nodes
have been visited
for non-stationary behaviors by using other types of kriging (see 2.2) to estimate the mean
and variance of the Gaussian cdf. In that case however, variogram reproduction is not
theoretically ensured.
Histogram Transformation
The sequential Gaussian simulation algorithm as described above assumes that the variable
is Gaussian. If that is not the case, it is possible to transform the marginal distribution
of the variable into a Gaussian distribution and work on the transformed variable. The
transformation of variable Z(u) with cdf F
Z
into a standard normal variable Y (u) with
cumulative distribution function G is written:
Y (u) = G
1
Z(u)
(3.2)
This transformation does not ensure that Y is multivariate Gaussian, only its histogram
is. One should check that the multivariate (or at least bi-variate) Gaussian hypothesis holds
for Y before performing Gaussian simulation. If the hypothesis is not appropriate, other
algorithms that do not require Gaussianity, e.g. sequential indicator simulation (SISIM)
should be considered.
The sequential Gaussian simulation algorithm proceeds as follows:
1. Transform Z into a Gaussian variable according to Eq. (3.2)
2. Simulate Y as in Algorithm 1
3. Back-transform the simulated values y
1
, ..., y
N
into z
1
, ..., z
N
:
z
i
= F
1
G(y
i
)
i = 1,. . . ,N (3.3)
In SGSIM the distribution function F is infered from a set of values
1
. . .
L
which can either be the hard data or values read from a le. F is built such that: F(
1
) =
1
L
and F(
L
) = 1
1
L
. The back-transformation step of SGSIM requires in addition the
modeling of the tails of the target distribution F(z), for values outside the data interval:
z <
1
or z >
L
.
Two options are available in SGSIM:
Z is bounded: F(Z < z
min
) = F(Z > z
max
) = 0. The lower tail of F is then
modeled by:
F(
1
) F(z)
F(
1
)
=

1
z
1
z
min
3
z ]z
min
,
1
[ (3.4)
and the upper tail:
F(z) F(
L
)
1 F(
L
)
=
z
L
z
max

1
1/3
z ]
L
, z
max
[ (3.5)
Z is not bounded: the lower tail is then given by:
F(z) = F(
1
) exp
(z
1
)
2
z <
1
(3.6)
and the uppert tail:
1 F(z)
1 F(
L
)
=

L
z
z >
L
(3.7)
Example:
Algorithm SGSIM is used to simulate the shale porosity in the second layer of Stanford V.
The histogram of the shale porosity data is shown in Fig 3.2. It is used as the reference his-
togram in the normal score transformation of the data (see the discussion about histogram
transformation). Since the data are too sparse to allow a reliable inference of the variogram,
3.2. SGSIM 39
the variogram is arbitrarily modeled by an isotropic spherical variogram, of range 40. Dif-
ferent views of a realization are show in Fig. 3.3. The parameter le corresponding to the
run is reproduced on Fig. 3.2.
Fig. 3.2 and 3.2 show that both the target distribution and variogram are correctly re-
produced.
Figure 3.1: Histogram of shale porosity
<parameters> <algorithm name="sgsim" />
<Property_Name value="sgsim_shale_porosity" />
<Nb_Realizations value="5" />
<Seed value="14071789" />
<Hard_Data grid="layer 2 shale data" property="porosity" />
<Assign_Hard_Data value="1" />
<Use_Target_Histogram value="1" />
<Target_Hist_From_Harddata value="1" />
<Target_Hist_From_File value="0" />
<Use_Min_Max value="1" />
<Target_Hist_Min value="0.01" />
<Target_Hist_Max value="0.27" />
<Kriging_Type value="Simple Kriging (SK)" />
0 0 0 " />
<Variogram nugget="0" structures_count="1" >
<structure_1 contribution="1" type="Spherical" >
<angles x="0" y="0" z="0" />
</structure_1>
</Variogram>
</parameters>
Figure 3.2: Sequential Gaussian simulation parameters
3.2. SGSIM 41
(a) 3D view
Figure 3.3: shale porosity simulated with SGSIM
Figure 3.4: Reproduction of the target histogram
Figure 3.5: Reproduction of the input variogram. The continuous black line is the input
variogram, the red crosses the variogram of the SGSIM realization
3.2. SGSIM 43
Hard Data: The grid and the name of the property containing the hard data
Use Target Histogram: A ag indicating whether to perform steps 1 and 3 of Al-
gorithm 1. If the ag is set, distribution F (use in relations (3.2) and (3.3)) must be
specied. The distribution can either be infered from the hard data, if any, or dened
explicitely (see parameters Target Hist From Harddataand Target Hist From File).
Setting this ag is necessary if the hard data do not follow a standard normal distri-
bution, as dicussed in Section 3.2
Target Hist From Harddata: If this ag is on (equal to 1), the target histogram
corresponding to cdf F will be infered from the hard data
Target Hist From File: If this ag is on, the histogram is read from a le. The le
name is specied by Target Hist File. Note that Target Hist From Harddata
and Target Hist From File are mutually exclusive : both can not be set to 1
at the same time.
Target Hist File: The name of the le containing the histogram description. This
is an ASCII le containing values drawn from the target distribution. The values
can be separated by spaces, tabulations of return carriage. The target distribution is
computed from that list of values. Any number of values can be specied in the le.
This parameter is ignored if Target Hist From File is set to 0
Use Min Max: Aag indicating whether variable Z is bounded: z
min
Z(u) z
max
u.
Target Hist Min: The minimum of the target distribution. This parameter is ignored
if Use Min Max is set to 0
Target Hist Max: The maximum of the target distribution. This parameter is ignored
if Use Min Max is set to 0
Kriging Type: The possible types of kriging are: Simple Kriging (SK), Ordinary Krig-
ing (OK), Kriging with Trend (KT) and Simple Kriging with Localy Varying Mean
(LVM). Because of Gaussian simulation theory, simple kriging should be the pref-
ered option unless deemed inappropriate.
Trend: The trend components. The possible components are, with u = (x, y, z) being
the coordinates vector:
t
1
(u) = x, t
2
(u) = y, t
3
(u) = z, t
4
(u) = x.y, t
5
(u) = x.z, t
6
(u) = y.z, t
7
(u) = x
2
,
t
8
(u) = y
2
, t
9
(u) = z
2
.
The trend is coded by a string of 9 ags, for example string 0 1 0 0 0 1 0 1 0
corresponds to a trend with only components t
2
, t
6
and t
8
, i.e.: T(u) = x + yz + y
2
.
Each ag is separated by a space.
Local Mean Property: The property of the simulation grid containing the non-stationary
mean. A mean value must be available at each location to be simulated.
Max Conditioning Data: The maximum number of conditioning data used to infer
the local conditional distribution at each simulation node
Variogram: The variogram model for the variable to be simulated
3.3 SISIM
The algorithm SISIM relies on indicator kriging (see Section 2.3) to infer the conditional
cumulative distribution function s F
Z
(u
k+1
| (n + k)), k = 1, . . . , N 1 of Eq. (3.1).
Let Z(u) be a continuous variable. SISIM implements Algorithm 2.
Step 2(c) of Algorithm 2 calls for the solution of K kriging systems. Computations
can actually be alleviated by using the same indicator variogram for all thresholds (see
discussion on median indicator kriging in section 2.3).
3.3. SISIM 45
Algorithm 2 SISIM - continuous variable
1. Choose a discretization of the range D of Z(u): z
1
, . . . , z
K
2. Dene a path visiting all locations to be simulated
3. For each location u along the path:
(a) Retrieve the neighboring conditioning data of location u: z(u
), = 1, . . . , N
(b) Turn each conditioning datum z(u
) into a vector of indicator values: v(u
) =
z(u
), z
1
, . . . , i
z(u
), z
K
. The indicator variable is dened as:

i(u, z
k
) =
1 if z(u) z
k
0 otherwise
(c) Estimate the indicator random variable I(u, z
k
) for each of the K cutoffs by
solving a kriging system. The simple kriging estimator is given by:
I
(u, z
k
) E{I(u, z
k
)} =
N
=1
I(u
, z
k
) E{I(u
, z
k
)}
(d) The estimated values i
(u, z
k
) = Prob
(Z(u) z
k
), after correction of order-
relation problems (see Section 2.3), dene an estimate of the cumulative distri-
bution function F
Z(u)
of the variable Z(u). Draw a value from that ccdf and
assign it as a datum at location u.
(e) Loop until all locations are visited
4. Repeat the previous steps to generate another simulated realization
Sampling the estimated distribution function
At each location to be simulated, the ccdf F
Z(u)
is estimated at values z
1
, . . . , z
K
. However
sampling from F
Z(u)
as described in step 3(e) of Algorithm 2 requires the knowledge of
F
Z(u)
(z) for all z D (see Algorithm 2). In SISIM F
Z(u)
is interpolated as follows:
F
Z(u)
(z) =
Z(u)
(z
1
) e
(z
1
z)
2
if z z
1
F
Z(u)
(z
1
) +
zz
k
z
k+1
zk
Z(u)
(z
2
) F
Z(u)
(z
1
)
if z
k
z zk + 1
1
1 F
Z(u)
(z
K
)
z
K
z
if z
k+1
z, > 0
(3.8)
where F
Z(u)
(z
k
) = i
(u, z
k
) is estimated by indicator kriging.
Categorical Variables
If Z(u) is a categorical variable which only takes the integer values {1, . . . , K}, Algo-
rithm 2 is modied as described in Algorithm 3.
Implementation limitations
Algorithm SISIM does not allow for inequality constraints or soft data. If such data must
be accounted for, use algorithm COSISIM instead, see section 4.
Example:
Porosity in the second layer of Stanford V is bi-modal (see Fig 1.16): it is higher in the
sand facies than in shale. Its spatial distribution is also variable: the high porosity region
has a 45
o
azimuth anisotropy, while the low porosity is isotropic. Full indicator simulation
is used to simulate such a complex eld. Five thresholds are selected such that each of the
two modes of the porosity distribution is described by the same number of thresholds. The
cut-offs z
1
, . . . , z
5
retained and the corresponding probabilities P(Z z
1
), . . . , P(Z z
5
)
are:
threshold z
i
0.06 0.15 0.23 0.27 0.3
P(Z z
i
) 0.22 0.54 0.63 0.69 0.9
The variograms corresponding to the three rst thresholds are spherical isotropic mod-
els, of range 40, while the last two variograms have a 45
o
azimuth anisotropy and an
anisotropy ratio of 3. The parameter le for the simulation is reproduction in Fig 3.6.
3.3. SISIM 47
Algorithm 3 SISIM - categorical variable
1. Dene a path visiting all locations to be simulated
(a) Retrieve the neighboring categorical conditioning data: z(u
), = 1, . . . , N
(b) Turn each datum z(u
) into a vector of indicator data values: v(u
) =
z(u
), z
1
, . . . , i
z(u
), z
K
. The indicator variable is dened by:

i(u, k) =
1 if Z(u) = k
0 otherwise
(c) Estimate the indicator random variable I(u, k) for each of the K categories by
solving a kriging system. The simple kriging estimator is given by:
I
(u, k) E{I(u, k)} =

N
=1
I(u
, k) E{I(u
, k)}
(d) The estimated values i
(u, k) = prob
(Z(u) = k), after correction of order-

relation problems (see 2.3), dene an estimate of the discrete probability density
function (pdf) f
Z(u)
of the categorical variable Z(u). Draw a realization from
f and assign it as a datum at location u.
(e) Loop until all locations are visited
3. Repeat the previous steps to generate another realization
Different views of a sample realization are shown on Fig 3.7.
Example:
As noted in the previous example, the complexity of the porosity eld of Stanford V layer
2, comes from its non-stationarity: porosity is high in the sand regions and low in the
shale regions. Full indicator simulation provides the necessary exibility to simulate such
a complex eld. Another possible approach is to rst separate the eld into two stationary
regions, channel sand and shale, and then simulate porosity in each region independently.
In this example, algorithm SISIM is used to simulate the two sand and shale facies. Full
indicator simulation is chosen to simulate sand bodies preferentially oriented along a 45
o
direction and isotropic shale bodies. A long range is used for the sand variogram to mimic
the high connectivity of the actual channels of layer 2 (see Fig. 1.15). The input parameters
are reproduced in Fig 3.8. Different views of a sample realization are shown in Fig 3.9.
It can be noticed that while the connectivity of the simulated sand channels is high, the
sinuosity of the reference channels is not reproduced.
3.3. SISIM 49
<parameters> <algorithm name="sisim" />
<Property_Name value="sisim" />
<Seed value="14071789" />
<Hard_Data_Grid value="layer 2 well data" />
<Hard_Data_Property value="porosity" />
<Categorical_Variable_Flag value="0" />
<Nb_Indicators value="5" />
<Thresholds value="0.06 0.15 0.23 0.27 0.3" />
<Marginal_Probabilities value="0.22 0.54 0.63 0.69 0.9" />
<Cdf_Min value="0" />
<Cdf_Max value="0.4" />
0 0 0 " />
<Full_Ik_Flag value="1" />
<Variogram_Full_Ik nugget="0" structures_count="1" >
<angles x="0" y="0" z="0" />
</structure_1>
</Variogram_Full_Ik>
<Variogram_Full_Ik_2 nugget="0" structures_count="1" >
<angles x="0" y="0" z="0" />
</structure_1>
</Variogram_Full_Ik_2>
<angles x="0" y="0" z="0" />
</structure_1>
<angles x="30" y="0" z="0" />
</structure_1>
<angles x="45" y="0" z="0" />
</structure_1>
</parameters>
Figure 3.6: Sequential indicator simulation parameters
(a) 3D view
Figure 3.7: Porosity simulated with full indicator simulation
3.3. SISIM 51
<parameters> <algorithm name="sisim" />
<Property_Name value="sisim" />
<Seed value="14071789" />
<Hard_Data_Grid value="well data 2 facies" />
<Hard_Data_Property value="facies " />
<Marginal_Probabilities value="0.63 0.37" />
<Cdf_Min value="0" />
<Cdf_Max value="0.4" />
0 0 0 " />
<Median_Ik_Flag value="0" />
<Full_Ik_Flag value="1" />
<Variogram_Full_Ik nugget="0" structures_count="1" >
<angles x="0" y="0" z="0" />
</structure_1>
</Variogram_Full_Ik>
<angles x="45" y="0" z="0" />
</structure_1>
</parameters>
Figure 3.8: Sequential indicator simulation parameters - Facies simulation
(a) 3D view
Figure 3.9: Sand (in blue) and shale (in white) facies simulated with full indicator simula-
tion
3.3. SISIM 53
Hard Data Grid: The grid containing the hard data property
Hard Data Property: The name of the property containing the hard data
Categorical Variable Flag: A ag indicating whether the variable to be simu-
lated is categorical or continuous. If the ag is on (i.e. equal to 1), the variable is
categorical
Thresholds: The threshold values z
1
, . . . , z
K
dening each indicator variable. This is
only needed for continuous variables. The threshold values are separated by spaces.
This parameter is ignored if Categorical Variable Flag is set to 1
ability values are separated by spaces. In the case of a continuous variable (ag
Categorical Variable Flag set to 0) the marginal probability values p
1
, . . . , p
K
must dene a valid cumulative distribution function, that is:
p
k
[0, 1] k
p
k
p
k
, k k
If the variable is categorical (ag Categorical Variable Flag set to 1),

p
1
, . . . , p
K
must dene a valid pdf, that is:
p
k
[0, 1] k
K
k=1
p
k
= 1
Max Conditioning Data: The maximum number of conditioning data used for each
kriging
Full Ik Flag: Aag indicating whether full IKshould be performed. Note that Median Ik Flag
and Full Ik Flag are mutually exclusive
Variogram Median Ik: The single indicator variogram model used for median IK. If
Median Ik Flag is off this parameter is ignored
Variogram Full Ik: The indicator variogram models (one for each indicator). If
Median Ik Flag is on this parameter is ignored
3.4 COSGSIM
COSGSIM allows to simulate a Gaussian variable accounting for secondary information.
Let Z
1
(u) and Z
2
(u) be two correlated random variables. Z
1
(u) is called the primary vari-
able, and Z
2
(u) the secondary variable. The implementation of the COSGSIM simulation
of the primary variable Z
1
conditioned to both primary and secondary data is described in
Algorithm 4.
Accounting for non stationarity
Reproduction of the covariance of Z
1
is ensured in theory only if simple cokriging is used
in step 2(b) of Algorithm 4. Simple cokriging requires that the expected values of Z
1
and
Z
2
are known. If the means are not known, COSGSIM allows the use of ordinary cokriging
in step 2(b) of Algorithm 4. However, the reproduction of the covariance of Z
1
is not
guaranteed anymore.
The Markov Models for coregionalization
As discussed in section 2.4, it is possible to simplify the inference and modeling of the three
covariances Cov
Z
1
(u), Z
1
(u+h)
, Cov
Z
1
(u), Z
2
(u+h)
, Cov
Z
2
(u), Z
2
(u +h)
by
only retaining the co-located secondary datum at each location u and assuming one of
two Markov models. The algorithm COSGSIM allows to choose between solving a full
cokriging system, or using either Markov Model 1 (MM1) or Markov Model 2 (MM2).
3.4. COSGSIM 55
Algorithm 4 COSGSIM
1. If necessary, transform Z
1
and Z
2
into Gaussian variables Y
1
and Y
2
, according to
Eq. (3.2). Make sure that the transformed variables Y
1
and Y
2
are at least bigaussian.
If they are not, another simulation algorithm should be considered, COSISIM for
example (see 3.5).
2. Simulate variable Y
1
:
(a) Dene a path visiting each node of the grid
(b) At each node u, the local conditional cumulative distribution function is Gaus-
sian, its mean is estimated by the simple cokriging estimate and its variance
by the simple cokriging variance. The conditioning data consist of the original
primary and secondary data and the previously simulated values
(c) Draw a value from the previous ccdf and add it to the data set
(d) Proceed to the next node in the path and repeat the two previous steps until all
nodes are visited
3. Back-transform the simulated values y
1,1
, ..., y
1,N
into z
1,1
, ..., z
1,N
:
z
1,i
= F
1
1
G
1
(y
1,i
)
i = 1,. . . ,N (3.9)
where F
1
is the target distribution function of Z
1
and G
1
is the standard normal cdf.
Example:
Seismic impedance in Stanford V is correlated to porosity with a coefcient of correla-
tion of 0.48, see Fig. 3.10. Hence it is interesting to use seismic impedance as a secondary
data to simulate porosity. Seismic information is available at every grid node, and in order
to simplify the covariances modeling only the co-located seismic information is used to
simulate a given node. The covariance inference problem is further simplied by assuming
that the MM2 screening property (Eq. 2.10) applies. The variogram of seismic impedance
is rst modeled by an anisotropic exponential variogram of ranges 40/15/5, the longest
range is oriented along the X direction. The variogram of porosity is then modeled accord-
ing to Eq. (2.12), with
R
a spherical model of ranges 40/40/5. The complete list of input
parameters to COSGSIM is reproduced in Fig. 3.11. Only the second layer of Stanford V
is simulated.
Fig. 3.12 shows different views of the simulated eld, along with the corresponding
seismic impedance and true porosity elds.
Figure 3.10: Porosity vs. seismic impedance
3.4. COSGSIM 57
<parameters> <algorithm name="cosgsim" />
<Property_Name value="cosgsim_porosity" />
<Primary_Harddata_Grid value="layer 2 well data" />
<Primary_Variable value="porosity" />
<Assign_Hard_Data value="1" />
<Secondary_Harddata_Grid value="layer 2" />
<Secondary_Variable value="seismic impedance" />
<Seed value="14071789" />
<Kriging_Type value="Simple Kriging (SK)" />
<Cokriging_Type value="Markov Model 2" />
<Max_Conditioning_Data_1 value="12" />
<Search_Ellipsoid_1 value="70 70 6
0 0 0 " />
<Variogram_C11 nugget="0" structures_count="2" >
<structure_1 contribution="0.25" type="Exponential" >
<angles x="0" y="0" z="0" />
</structure_1>
<angles x="0" y="0" z="0" />
</structure_2>
</Variogram_C11>
<MM2_Correl_Z1Z2 value="-0.5" />
<MM2_Variogram_C22 nugget="0" structures_count="1" >
<structure_1 contribution="1" type="Exponential" >
<angles x="0" y="0" z="0" />
</structure_1>
</MM2_Variogram_C22>
<Transform_Primary_Variable value="1" />
<Prim_Histogram_From_Harddata value="1" />
<Primary_Histogram_From_File value="0" />
<Target_Hist_Min value="0" />
<Target_Hist_Max value="0.4" />
<Transform_Secondary_Variable value="1" />
<Sec_Histogram_From_Harddata value="1" />
<Secondary_Histogram_From_File value="0" />
</parameters>
Figure 3.11: Sequential Gaussian co-simulation parameters
(a) 3D view
(b) Cross-section 1 - Seis-
mic impedance
(c) Cross-section 1 - COS-
GSIM realization
(d) Cross-section 1 - True
porosity eld
(e) Cross-section 2 - Seis-
mic impedance
(f) Cross-section 2 - COS-
GSIM realization
(g) Cross-section 2 - True
porosity eld
Figure 3.12: Gaussian co-simulation of porosity, conditioned to seismic impedance
3.4. COSGSIM 59
Primary Variable: The name of the property containing the primary hard data
Secondary Harddata Grid: The grid containing the secondary hard data
Secondary Variable: The name of the property containing the secondary hard data
Transform Primary Variable: A ag indicating whether steps 1 and 3 of Algo-
rithm 4 should be performed on the primary variable. Setting this ag is necessary
if the primary hard data do not follow a standard normal distribution
Prim Histogram From Harddata: If this ag is on, the distribution function of the
primary variable is inferred from the primary hard data. This parameter is ignored if
Transform Primary Variable is set to 0
Primary Histogram From File: If this ag is on, the distribution of the primary
variable is inferred from values read from a le. The le name is specied by
Primary Histogram File. See parameter Primary Histogram File for
a description of the le format and its expected content. Note that Prim Histogram From Harddata
and Primary Histogram From File are mutually exclusive
Primary Histogram File: The name of the le containing the histogram descrip-
tion of the primary variable. This is an ASCII le containing values drawn from
the distribution of the primary variable. The values can be separated by spaces, tab-
ulations of return carriage. The distribution (cdf) of the primary variable is computed
fromthat list of values. This parameter is ignored if Primary Histogram From File
is set to 0
Use Min Max: A ag indicating whether we want to specify the min/max of the distribu-
tion of the primary variable. If not set, the min/max of the distribution is the min/max
of the hard data if Target Hist From Harddata is set, or the min/max of the
values in le Target Hist File
Target Hist Min: The minimum of the distribution of the primary variable. This pa-
rameter is ignored if Use Min Max is set to 0
Target Hist Max: The maximum of the distribution of the primary variable. This
parameter is ignored if Use Min Max is set to 0
Transform Secondary Variable: Aag indicating whether step 1 of Algorithm4
should be performed on the secondary variable. Setting this ag is necessary if the
secondary hard data do not follow a standard normal distribution.
Sec Histogram From Harddata: If this ag is on, the distribution function of the
secondary variable will be inferred from the secondary hard data. This parameter is
ignored if Transform Secondary Variable is set to 0
Secondary Histogram From File: If this ag is on, the distribution of the sec-
ondary variable will be inferred from values read from a le. The le name is speci-
ed by Secondary Histogram File. See parameter Secondary Histogram File
for a description of the le format and its expected content. Note that Sec Histogram From Harddata
and Secondary Histogram From File are mutually exclusive
Secondary Histogram File: The name of the le containing the histogram de-
scription. This is an ASCII le containing values drawn from the distribution. The
values can be separated by spaces, tabulations of return carriage. The target distribu-
tion is computed fromthat list of values. This parameter is ignored if Secondary Histogram From File
is set to 0
Kriging Type: The possible types of kriging include: Simple Kriging (SK) and Ordi-
nary Kriging (OK)
Cokriging Type: The cokriging option. Available options are: Full Cokriging, Co-
located Cokriging with Markov Model 1 (MM1), and Co-located Cokriging with
Markov Model 2 (MM2)
3.5. COSISIM 61
Max Conditioning Data 1: The maximum number of primary conditioning data
used for each cokriging
nd the primary conditioning data
Variogram C11: The variogram of the primary variable
Max Conditioning Data 2: The maximum number of secondary conditioning data
to be used for each cokriging. This parameter is only required if Cokriging Type
nd the secondary conditioning data. This parameter is only required if Cokriging Type
Variogram C12: The cross-variogram between the primary and secondary variables.
This parameter is only required if Cokriging Type is set to Full Cokriging
Variogram C22: The variogram of the secondary variable. This parameter is only re-
quired if Cokriging Type is set to Full Cokriging
Correl Z1Z2: The coefcient of correlation between the primary and secondary vari-
ables. This parameter is only required if Cokriging Type is set to MM1
MM2 Correl Z1Z2: The coefcient of correlation between the primary and secondary
variables. This parameter is only required if Cokriging Type is set to MM2
MM2 Variogram C22: The variogramof the secondary variable. This parameter is only
required if Cokriging Type is set to MM2. Note that the variogram of the sec-
ondary variable and the variogram of the primary variable must verify Eq. (2.12)
3.5 COSISIM
Algorithm SISIM simulates a variable Z conditioned only to z-values. The indicator for-
malismactually allows the conditioning to very diverse types of data: besides equality data,
inequality constraints or probabilistic constraints can be honored. Algorithm COSISIM al-
lows to account for such types of information.
Coding information as indicator values
Consider a continuous variable Z(u) and a discretization of its range by the K threshold
values z
1
, . . . , z
K
.
Different types of information can be coded into a vector of K indicator values I =
[i
1
, . . . , i
K
].
Equality data
The value of Z at a given location u is known, equal to z(u). The corresponding K indica-
tor values are given by:
i
k
=
1 if z(u) z
k
0 otherwise
k = 1, . . . , K
Example:
If K = 5, and z
1
= 1, z
2
= 2, z
3
= 3, z
4
= 4, z
5
= 5, then z = 3.7 would be coded as:
I =
0
0
0
1
1
Inequality data
There are cases when the value of z(u) is not known precisely; all that is known is an
interval that contains the value, e.g z(u) [a, b], or z(u) ]a, +[. The indicator vector
3.5. COSISIM 63
I corresponding to z(u) [a, b[ is given by:
i
k
=
0 if z
k
< a
not dened if z
k
[a, b[
1 if z
k
b
k = 1, . . . , K
Example:
If K = 5, and z
1
= 1, z
2
= 2, z
3
= 3, z
4
= 4, z
5
= 5, then z ]2.1, 4.9] would be coded
as:
I =
0
0
?
?
1
where ? denotes a value that is not dened (missing data value)

Soft data
Information that are only indirectly related to variable Z can also be coded into indicator
values, or more precisely, into prior probability values.
Consider a variable Y (u) correlated with Z. Y could be the permeability of a rock and
Z its porosity. If the conditional prior probability function F(Z(u) z
k
| Y(u)) is known,
it can be used as an indicator coding of the inuence of Y on Z:
I
k
(u) = F
Z(u) z
k
| Y(u)
k = 1, . . . , K
I
k
(u) is not anymore equal to 0 or 1 only but is a probability value, valued in [0, 1], reect-
ing the uncertainty on Z(u) given the information Y (u).
The prior probability function F could be estimated in different ways. One possibility
might be to use the experimental conditional distribution: If Y is a continuous variable and
y
1
, . . . , y
K
a discretization of its range, then
F
Z z
k
| Y [y
l
, y
l+1
]
=
1
uU
i(u
, [y
l
, y
l+1
])
uU
i(u
, z
k
) i(u
, [y
l
, y
l+1
])
U is the set of locations where both Z and Y are known and:
i(u
, z
k
) =
1 if Z(u
) z
k
0 otherwise
i(u
, [y
l
, y
l+1
]) =
1 if Y (u
) [y
l
, y
l+1
]
0 otherwise
If variable Y is categorical, taking the categorical values {1, . . . , K
}, then the prior

distribution function is given by:
F
Z z
k
| Y = k
=
1
uU
i(u
, k
uU
i(u
, z
k
) i(u
, k
) (3.10)
with i(u
, k
) dened as:
i(u
, k
) =
1 if Y (u
) = k
0 otherwise
Accounting for soft data: the Markov-Bayes algorithm
Indicator coding provides a exible means of accounting for different types of information.
The indicator coded soft data can be interpreted as a secondary variable Y and accounted
for by solving a cokriging system. The (ordinary) cokriging estimator is given by:
I
(u, z
k
) =
I(u
), z
k
) +
, z
k
with Y
, z
k
= F
Z(u) z
k
| Y
. In all rigor, the cokriging estimator should include

indicator and soft data at thresholds other than z
k
. However Goovaerts (1994) showed that
in practice those terms could be neglected.
Solving this cokriging system calls for the knowledge of two covariance functions and
one cross-covariance function for each threshold z
k
. The inference and modeling work can
be alleviated if the following Markov property is assumed:
Prob
Z(u
) z
k
| i(u, z
k
), y(u, z
k
)
= Prob
Z(u
) z
k
| i(u, z
k
)
u, u, z
k
(3.11)
3.5. COSISIM 65
With this screening (Markov) hypothesis, it can be shown [Zhu and Journel, 1993]
that:
C
Y
(h, z
k
) =
|B(z
k
)| C
I
(h, z
k
) if h = 0
B(z
k
)
2
C
I
(h, z
k
) if h > 0
(3.12)
C
IY
(h, z
k
) = B(z
k
) C
I
(h, z
k
) (3.13)
where B(z
k
) is the difference between two conditional expectations:
B(z
k
) = m
1
(z
k
) m
0
(z
k
) [1, 1]
dened as:
m
1
(z
k
) = E
Y (u, z
k
)|I(u, z
k
) = 1
m
0
(z
k
) = E
Y (u, z
k
)|I(u, z
k
) = 0
Updating of the prior probabilities Y (u, z

k
) under that Markov hypothesis is referred to
as the Markov-Bayes algorithm.
Example:
Section 3.3 showed an example of indicator simulation of sand and shale facies in the
layer 2 of Stanford V. In that example, the simulation was only constrained to facies well
data, and the sand bodies, while having the long continuity of the true channels, failed
to reproduce their sinuosity. The result of the indicator simulation can be improved by
integrating additional information, such as seismic impedance.
As mentioned in the discussion on information coding into indicator values, seismic
impedance has to be transformed into a prior probability P(I(u) = 0 | Y(u)) of observing
the shale facies given a seismic impedance datum. The prior probabilities are computed
according to Eq. (3.10). Fig. 3.14 shows the obtained probability maps for two cross-
sections.
In order to simplify the variogram inference, the Markov property (3.11) is assumed.
Parameters B(0) and B(1) (Eq. (3.12)), inferred from the well data, are both equal to 0.5.
The complete parameter le for the run is reproduced in Fig. 3.13.
Fig. 3.14 shows two cross-sections of a sample realization. Due to the extremely high
quality of the seismic information (the prior probability maps clearly identify the channels),
the channels are very well reproduced.
<parameters> <algorithm name="cosisim" />
<Property_Name value="cosisim_facies" />
<Seed value="14071789" />
<Marginal_Probabilities value="0.57 0.43" />
<Kriging_Type value="Simple Kriging" />
<Median_Ik_Flag value="1" />
<Primary_Harddata_Grid value="well data" />
<Primary_Indicators count="2" value="shale indicator; sand indicator" />
<Max_Conditioning_Data_Primary value="12" />
<Search_Ellipsoid_Primary value="60 60 7
0 0 0 " />
<Secondary_Harddata_Grid value="layer 2" />
<Secondary_Indicators count="2" value="P(I=0|Y); P(I=1|Y)" />
<Max_Conditioning_Data_Secondary value="12" />
<Search_Ellipsoid_Secondary value="20 20 4
0 0 0 " />
<Bz_Values value="0.5 0.5" />
<Variogram_Median_Ik nugget="0" structures_count="1" >
<angles x="45" y="0" z="0" />
</structure_1>
</Variogram_Median_Ik>
</parameters>
Figure 3.13: Sequential indicator co-simulation parameters
3.5. COSISIM 67
(a) Seismic impedance (b) P(I = 0 | Y ) (c) COSISIM realization
(d) Seismic impedance (e) P(I = 0 | Y ) (f) COSISIM realization
Figure 3.14: Indicator co-simulation of facies, conditioned to seismic impedance. Each line
of plots corresponds to a different cross-section. Each line shows the seismic impedance,
the corresponding prior probability P(I = 0 | Y ) (I = 0 corresponds to the shale facies)
and the COSISIM realization. The shale facies is in white, sand in blue.
Primary Indicators: The name of the properties containing the indicator values for
the primary variable. Exactly K properties must be specied.
Secondary Harddata Grid: The grid containing the soft data. This parameter is
optional
Secondary Indicators: The name of the properties containing the indicator values
for the soft data. Exactly K properties must be specied. This parameter is optional
Categorical Variable Flag: A ag indicating whether the variable to be simu-
lated is categorical or continuous. If the ag is on (equal to 1), the variable is cate-
gorical
Nb Indicators: The number K of indicator variables or the number K of categories
Thresholds: The threshold values z
1
, . . . , z
K
dening each indicator variable. This is
only needed for continuous variables. The threshold values are separated by spaces.
This parameter is ignored if Categorical Variable Flag is set to 1
ability values are separated by one or more space
Kriging Type: The possible types of kriging are: Simple Kriging (SK) and Ordinary
Kriging (OK)
Full Ik Flag: Aag indicating whether full IKshould be performed. Note that Median Ik Flag
and Full Ik Flag are mutually exclusive
Max Conditioning Data Primary: The maximum number of primary condition-
ing data to be used for each cokriging
3.6. SNESIM 69
Search Ellipsoid Primary: The ranges and angles dening the search ellipsoid
used to nd the primary conditioning data
Max Conditioning Data Secondary: The maximum number of soft conditioning
data to be used for each cokriging
Search Ellipsoid Secondary: The ranges and angles dening the search ellip-
soid used to nd the soft conditioning data
Bz Values: The B values for each indicator variogram (see the Markov-Bayes algo-
rithm). The values are separated by spaces. This parameter is not required if no soft
data are used
Variogram Median Ik: The single indicator variogram model used for median IK. If
Median Ik Flag is off this parameter is ignored
Variogram Full Ik: The indicator variogram models (one for each indicator) of each
indicator variable. If Median Ik Flag is on this parameter is ignored
3.6 SNESIM
Algorithms based on two-points statistics, e.g. all algorithms based on kriging, fail to model
complex curvilinear heterogeneities. It is possible to produce images with very different
structures (e.g. some with channels stretching from one end of the image to the other, and
some featuring non-overlapping ellipsis of small dimensions) which all share the same
two-point statistics (i.e. covariance) [Strebelle, 2000]: two point-statistics are not enough
to fully characterize complex (low-entropy) structures, and one must rely on higher order
statistics to model those structures.
AlgorithmSNESIM implements a sequential simulation algorithm, called single normal
equation simulation [Strebelle, 2000], that simulates categorical variables using multiple-
point statistics. In SNESIM, the multiple-point statistics are not provided by an analytical
function but are directly borrowed from a training image. The training image depicts the
type of structures the simulated realizations should exhibit, and it needs not be constrained
to any local data such as sampled values or well data. Training images can conveniently
be generated by unconditional object-based simulations [Tjelmeland, 1996; Holden et al.,
1998; Chil` es and Delner, 1999].
Algorithm 5 describes the simplest version of the snesim algorithm.
Step 2(b) of Algorithm 5 is critical: it requires the search of all locations u
of the
training image where D
j
(u
) = D
j
(u), 1 j J. Performing an exhaustive search for
each location to be simulated is not practical. Instead SNESIM uses a tree data structure
(called search-tree) to efciently perform the search. The search-tree stores all the differ-
ent data events D
j
(u
), j = 1, . . . , J, found at all locations u
of the training image, and

allows to retrieve the conditional probability distribution of step 2(c) in O(log(J)). This
speed comes at the cost of a possibly large RAM memory demand.
Let N
TI
be the total number of locations in the training image. Since for a given data
event size j, there can not be more than N
TI
different data events in the training image, a
(crude) upper-bound of the memory demand of the search tree is:
Memory Demand
J
j=1
min(K
j
, N
TI
)
Multiple grid simulation
Reproducing large scale structures would call for using a large template. However because
of memory constraints, it is not practical to work with too large templates. A work-around
is provided by the multiple-grid approach [Tran, 1994].
Consider a 3-D Cartesian grid G. Dene G
(l)
as the sub-set of G such that: G
(0)
= G
and G
(l)
is obtained by down-sampling G
(l1)
by a factor of 2 in all 3 coordinate directions:
G
(l)
is the sub-set of G
(l1)
obtained by retaining every other node of G
(l1)
. G
(l)
is often
called the l
th
level grid.
Algorithm 6 describes how multiple grids are used in program SNESIM
Marginal distribution reproduction
It is sometimes desirable that the histogram of the simulated variable be close to a given
target distribution, e.g. the sample histogram. There is however no constraint in SNESIM
3.6. SNESIM 71
Algorithm 5 Simple Snesim
Consider a categorical variable Z(u) valued in {1, . . . , K}.
1. Dene a geometrical template (or window) T
J
= {h
1
, . . . , h
J
}, h
j
is a vector in
R
3
, and a path visiting all locations to be simulated
(a) Retrieve the conditioning data event D
J
(u) dened by template T
J
:
D
J
(u) = {z(u +h
1
), . . . , z(u +h
J
)}
(b) Find all locations u
, = 1, . . . , n, of the training image such that D

J
(u
) =
D
J
(u), i.e. such that:
z(u +h
j
) = z(u
+h
j
) j {1, . . . , J}
If the number n of such locations is lesser than a xed threshold n
min
, retrieve
the smaller data event D
J1
dened by T
J1
and repeat the search. This step is
repeated until n n
min
. Call J
the largest template size for which n n

min
.
(c) The conditional probability distribution P(Z(u) = k | D
J
(u)) =
E(I(u, k) | D
J
(u)), k = 1, . . . , K is modeled by the experimental conditional
expectation E
(I(u, k) | D
J
(u)) :
P
Z(u) = k | D
J
(u)
Z(u) = k | D
J
(u)
(I(u, k) | D
J
(u))
=
1
n
n
=1
i
z(u
), k
where i
z(u
), k
is the indicator function:

i
z(u
), k
1 if z(u
) = k
0 otherwise
Draw a simulated value from that conditional distribution and add it to the data
set
3. Repeat step 2 until all locations in the path are visited
Algorithm 6 SNESIM with multiple grids
1. Choose the number L of multiple grids to consider.
2. Start with the coarsest grid G
(L)
.
3. Build a new geometrical template T
(L)
J
by re-scaling template T
J
:
T
(L)
J
= {2
L1
h
1
, . . . , 2
L1
h
J
}
Template T
(L)
J
has the same number of vectors as T
J
but has a much greater extent
than T
J
, hence will be able to capture large-scale variations more efciently without
increasing the search tree size.
4. Build the search tree using the training image and template T
(L)
J
5. Simulate all nodes of G
(L)
as in Algorithm 5
6. All the nodes of G
(L)
are now part of the data set. Repeat steps 3 to 5 on the next
ner grid G
(L1)
7. Loop until all nested multiple grids have been simulated
as described in Algorithm 5 or Algorithm 6 to ensure that any target distribution will be
reproduced. SNESIM uses an updating of the conditional distribution function computed
at each location (see step 2(c) of Algorithm 5) to control the histogram of the simulated
values.
Let p
c
k
, k = 1, . . . , K, denote the proportions of values in class k simulated so far,
and p
t
k
, k = 1, . . . , K, be the target proportions. Step 2(c) of Algorithm 5 is modied as
follows:
1. Compute the conditional probability distribution as originally described in step 2(c)
of Algorithm 5.
2. Update the computed probabilities P
Z(u) = k | D
J
(u)
into probabilities P
Z(u) =
k | D
J
(u)
with:
P
Z(u) = k | D
J
(u)
= P
Z(u) = k | D
J
(u)
+

1
(p
t
k
p
c
k
)
3.6. SNESIM 73
where [0, 1[ is a control parameter. If = 0, no updating is performed. Con-
versely, if 1, reproducing the target distribution entirely controls the simulation
process.
If P
Z(u) = k | D
J
(u)
/ [0, 1] it is reset to the closest bound. All updated

probability values are then rescaled to sum up to 1:
P
Z(u) = k | D
J
(u)
=
1
K
k=1
P
Z(u) = k | D
J
(u)
P
Z(u) = k | D
J
(u)
Conditioning to soft data

SNESIM can also account for a secondary variable Y (u), for example a seismic amplitude
pattern centered on u. The rst step is to turn Y (u) into a vector of conditional probabilities
(see the paragraph on indicator coding of soft data in section 3.5):
Z(u) = 1 | Y(u)
.
.
.
P
Z(u) = K | Y(u)
The conditional distribution function of step 2(c) of Algorithm5, P
Z(u) = k | D
J
(u)
,
is then updated into P
Z(u) = k | D
J
(u), Y(u)
using the formula proposed by Journel

(2002). Dene the four quantities:
x =
1 P
Z(u) = k | D
J
(u), Y(u)
Z(u) = k | D
J
(u), Y(u)
a =
1 P
Z(u) = k
Z(u) = k
b =
1 P
Z(u) = k | D
J
(u)
Z(u) = k | D
J
(u)
c =
1 P
Z(u) = k | Y(u)
Z(u) = k | Y(u)
Under the hypothesis of conditional independence:

P
D
J
(u), Y(u) | Z(u)) = P
D
J
(u) | Z(u)) P
Y(u) | Z(u)) (3.14)

the updated probabilities are given by:
x =
b.c
a
(3.15)
Step 2(c) of Algorithm 5 is then modied as follows:
1. Estimate the probabilities P
Z(u) = k | D
J
(u)
as described in Algorithm 5
2. Compute the updated probabilities P
Z(u) = k | D
J
(u), Y(u)
using Eq. (3.15)

3. Draw a realization from the updated distribution function
Accounting for local non-stationarity
The training image should be reasonably stationary so that meaningful statistics can be in-
ferred from it. It is however possible to introduce some non-stationarity in the simulation
by accounting for local rotation information. Call R(u) the rotation angle about the (verti-
cal) Z-axis at location u. In SNESIM the rotation axis is always the Z-axis. The angle R(u)
is restricted to take only a nite number of different values r
1
, . . . , r
M
.
Step 2(a) of Algorithm 5 could be modied as follows to account for R:
1. Dene a new geometrical template RT
J
(u) by rotating T
J
by angle r(u):
RT
J
= {(u)h
1
, . . . , (u)h
1
}
where (u) is the rotation matrix:
(u) =
cos(r(u)) sin(r(u)) 0
sin(r(u)) cos(r(u)) 0
0 0 1
2. Retrieve the conditioning data event D

J
(u) dened by template RT
J
(u)
This solution, while elegant, yields disappointing results because of aliasing issues
caused by the rotation of template T
J
.
Those aliasing issues are alleviated if instead of rotating template T
J
the entire training
image is rotated and scanned with T
J
.
3.6. SNESIM 75
Algorithm 7 Snesim with locally varying azimuth
Consider a categorical variable Z(u) valued in {1, . . . , K}.
1. Dene a geometrical template T
J
= {h
1
, . . . , h
J
} and a path visiting all the locations
to be simulated
2. Compute the rotated training images RTI
m
, m=1,. . . ,M corresponding to each angle
r
1
, . . . , r
M
. Build the corresponding M search trees, using template T
J
3. For each location u in the path:
(a) Retrieve the conditioning data event D
J
(u) dened by template T
J
:
D
J
(u) = {z(u +h
1
), . . . , z(u +h
J
)}
(b) Find all locations u
, = 1, . . . , n, of the training image rotated by angle r(u),

RTI
r(u)
, such that D
J
(u
) = D
J
(u), i.e. such that:
z(u +h
j
) = z(u
+h
j
) j {1, . . . , J}, u
RTI
r(u)
If the number n of such locations is lesser than a xed threshold n
min
, retrieve
the smaller data event D
J1
dened by T
J1
and repeat the search. This step is
repeated until n n
min
. Call J
the largest template size for which n n

min
.
(c) The conditional probability distribution P(Z(u) = k | D
J
(u)), k = 1, . . . , K is
modeled as follows:
P
Z(u) = k | D
J
(u)
= P
Z(u) = k | D
J
(u)
=
1
n
n
=1
i
z(u
), k
where i
z(u
), k
is the indicator function:

i
z(u
), k
1 if z(u
) = k
0 otherwise
Draw a realization from that conditional distribution, which is then added to the
data set
4. Repeat step 2 until all locations in the path have been visited
The new snesim algorithm is described in Algorithm 7
This solution which is the one implemented in algorithm SNESIM can be very memory
demanding as one search tree has to be built for each rotation angle r
1
, . . . , r
M
. However,
practice has shown that it is possible to generate fairly complex models while using a very
limited number of angles: about 5 angle values should be sufcient in most cases
Example:
The sand/shale facies simulation of section 3.3 is repeated using SNESIM. Cross-sections
of the 3-D, 100x130x8 training image are shown on Fig. 3.16. Four multiple grids are
used, and the template contains 100 nodes. The complete set of parameters is reproduced in
Fig. 3.15. Fig. 3.17 shows different views of a sample realization. Cross-sections (e) and (h)
of Fig. 3.17 correspond to the bottom of the channels of cross-sections (f) and (i), which
explains the noticeable lack on continuity of the channels in those sections. The channel
proportion in the realization of Fig. 3.17 is 35%, compared to the 30% target proportion.
Notice the sharp difference with the SISIM realization (Fig. 3.9): the sand bodies are
crisp are reproduce the continuity and sinuosity of the reference channels (Fig. 1.15).
Inuence of the template size
One of the key SNESIM parameter is the size of the template: the bigger the template,
the more structure it can capture. Fig. 3.18 shows the results obtained by increasing the
number of template nodes. In all cases, no multiple grids are used, and the other parameters
are those reproduced on Fig. 3.15.
Even with 200 nodes in the template, the channel structures are very poorly reproduced
if no multiple grids are used. 200 nodes indeed constitute a template of very limited extent,
e.g. 7x7x4 grid cells, making it impossible to capture structures that are continuous across
100 grid cells.
Inuence of the number of multiple grids
The template size is limited by the amount of RAM available, and it is necessary to use
the multiple-grid simulation technique to keep the template size low, yet capture large scale
structures. Fig. 3.19 shows the inuence of the number of multiple-grids on the quality of
3.6. SNESIM 77
<parameters> <algorithm name="snesim" />
<GridSelector_Sim value="layer 2" />
<Property_Name_Sim value="snesim_facies" />
<Seed value="211175" />
<PropertySelector_Training grid="TI" property="facies (2)" />
<Nb_Facies value="2" />
<Marginal_Cdf value="0.7 0.3" />
<Max_Cond value="100" />
0 0 0 " />
<Hard_Data grid="well data" property="facies (2)" />
<Use_ProbField value="0" />
<Use_Rotation value="0" />
<Use_Affinity value="0" />
<Cmin value="1" />
<Constraint_Marginal_ADVANCED value="0.5" />
<Nb_Multigrids_ADVANCED value="4" />
<Subgrid_choice value="1" />
<Previously_simulated value="4" />
</parameters>
Figure 3.15: SNESIM parameters
the realizations. The parameter used are those shown on Fig. 3.15, except for the number
of multiple grids and the template size, which was set to 30.
The use of multiple grids drastically improves the reproduction of large scale structures,
even with small, 30-nodes template.
The same sensitivity analysis is repeated with a 100-nodes template in Fig. 3.20, and a
200-nodes template in Fig. 3.21. In this case study, increasing the template size has little
inuence on the quality of the reproduction of large scale structures.
(a) XY cross-section (b) XY cross-section (c) XY cross-section
(d) XZ cross-section (e) XZ cross-section
Figure 3.16: Cross-sections of the training image
3.6. SNESIM 79
(a) 3D view (b) XZ section (c) XZ section
(d) section 1 (e) section 2 (f) section 3
(g) section 5 (h) section 8 (i) section 9
Figure 3.17: One sample SNESIM realization
(a) section 1 (b) section 1 (c) section 1 (d) section 1
(e) section 2 (f) section 2 (g) section 2 (h) section 2
(i) XZ section (j) XZ section (k) XZ section (l) XZ section
Figure 3.18: Inuence of the template size. Column 1: 30 nodes. Column 2: 60 nodes.
Column 3: 100 nodes. Column 4: 200 nodes.
3.6. SNESIM 81
(a) section 1 (b) section 1 (c) section 1 (d) section 1 (e) section 1
(f) section 2 (g) section 2 (h) section 2 (i) section 2 (j) section 2
(k) XZ sec-
tion
(l) XZ
section
(m) XZ sec-
tion
(n) XZ sec-
tion
(o) XZ sec-
tion
Figure 3.19: Inuence of the number of multiple grids, with a 30-nodes template. Column
1: 1 multiple grid. Column 2: 2 multiple grids. Column 3: 3 multiple grids. Column 4: 4
multiple grids. Column 5: 5 multiple grids.
(k) XZ sec-
tion
(l) XZ
section
(m) XZ sec-
tion
(n) XZ sec-
tion
(o) XZ sec-
tion
3.6. SNESIM 83
(k) XZ sec-
tion
(l) XZ
section
(m) XZ sec-
tion
(n) XZ sec-
tion
(o) XZ sec-
tion
PropertySelector Training: The grid and property containing the training im-
age
Nb Facies: The number K of different values that the simulated categorical variable Z
can take
Marginal Cdf: The marginal pdf of the variable. The probability values must be given
in the following order: P(Z = 0), . . . , P(Z = K 1)
Max Cond: The maximum number J of conditioning data to be retained for each simu-
lation
Search Ellipsoid: The ranges and angles dening the ellipsoid used to search neigh-
boring conditioning data. Template T
J
will be automatically built from the search
ellipsoid.
Hard Data: The grid and property containing the hard data
Use ProbField: This ag indicates whether the simulation should be conditioned to
prior local probability information (e.g. the indicator coding of soft data). Set the
ag to 1 to condition to prior local probability information
ProbField properties: The grid and properties containing the prior probability
values. One property must be specied for each possible value of Z(u). The k
th
property corresponds to P
Z(u) = k | Y(u)
. This parameter is ignored if Use ProbField

is set to 0
Use Rotation: This ag indicates whether to use additional angle information.
Use Global Rotation: This ag indicates whether to use a global rotation angle. If
set to 1, a single angle must be specied: R(u) = R is constant for all locations u
Use Local Rotation: This ag indicates whether to use a local rotation angle: if set
to 1, a rotation angle must be specied for each node of the simulation grid. Note that
Use Global Rotation and Use Local Rotation are mutually exclusive
3.6. SNESIM 85
Global Angle: The global rotation angle. The training image will be rotated by that
angle before the simulation begins. This parameter is ignored if Use Global Rotation
is set to 0
Rotation property: The property of the simulation grid containing the local angle
category. Angle categories range from 0 to M 1 where M is the total number of
angle categories. The angles corresponding to each angle category are specied by
Rotation categories
Rotation categories: The angles, expressed in degrees, corresponding to each cat-
egory. The angles are separated by spaces. This parameter is ignored if Use Global Rotation
is set to 0
Cmin: The minimal number of training replicates of a given data event that must be found
in order for that data event to be retained. This is referred to as n
min
in step 2(b) of
Algorithm 5.
Constraint Marginal ADVANCED: Anumber in [0,1] indicating the extent to which
the marginal pdf of the variable should be matched. See the constant in the discus-
sion on marginal distribution reproduction
Nb Multigrids ADVANCED: The number of multiple grids
Subgrid choice: A ag indicating whether sub-grids should be used. Sub-grids
is a structured simulation path which possibly allows to use a smaller geometrical
template D
J
Previously simulated: The number of nodes belonging to other sub-grids that
should be used to simulate a given sub-grid. Typically 4 nodes are used. This param-
eter is ignored if Subgrid choice is set to 0
Bibliography
Almeida, A. and Journel, A.: 1994, Joint simulation of multiple variables with a markov-
type coregionalization model, Mathematical Geology 26(5), 565588.
Chil` es, J. and Delner, P.: 1999, Geostatistics: Modeling spatial uncertainty, John Wiley
& Sons, New York.
Deutsch, C. and Journel, A.: 1992, GSLIB: Geostatistical Software Library and Users
Guide, Oxford University Press, New York.
Goovaerts, P.: 1994, Comparative performance of indicator algorithms for modeling con-
ditional probability distribution functions, Mathematical Geology 26, 389411.
Goovaerts, P.: 1997a, Geostatistics for natural resources evaluation, Oxford University
Press, New York.
Goovaerts, P.: 1997b, Geostatistics for natural resources evaluation, Oxford University
Press, New York.
Holden, L., Hauge, R., Skare, O. and Skorstad, A.: 1998, Modeling of uvial reservoirs
with object models, Mathematical Geology 30(5), 473496.
Journel, A.: 1999, Markov models for cross covariances, Mathematical Geology 31.
Journel, A.: 2002, Combining knowledge from diverse sources: An alternative to tradi-
tional data independence hypotheses, Mathematical Geology 34(5), 573596.
86
BIBLIOGRAPHY 87
Mao, S. and Journel, A.: 1999, Generation of a reference petrophysical/seismic data set: the
stanford v reservoir, Report 12, Stanford Center for Reservoir Forecasting, Stanford,
CA.
Strebelle, S.: 2000, Sequential simulation drawing structures from training images, PhD
thesis, Stanford University, Stanford, CA.
Tjelmeland, H.: 1996, Stochastic models in reservoir characterization and Markov ran-
dom elds for compact objects, PhD thesis, Norwegian University of Science and
Technology, Trondheim, Norway.
Tran, T.: 1994, Improving variogram reproduction on dense simulation grids, Computers
& Geosciences 20(7), 11611168.
Zhu, H. and Journel, A.: 1993, Formatting and interpreting soft data: Stochastic imaging
via the markov-bayes algorithm, 1, 112.

Sgems Manual

Uploaded by

Copyright:

Available Formats

Sgems Manual

Uploaded by

Document Information

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Sgems Manual

Uploaded by

Copyright:

Available Formats

Geostatistical Earth Modeling Software:

) be the rotated coordinate system, with Z

) be the rotated coordinate system, with Y

: Last, the ellipsoid is rotated about axis X

, see Fig. 1.12.

14 SECTION 1. GENERAL OVERVIEW

Figure 1.11: Second rotation about Y

1.3 Variogram Model Specication

(u, k) of Z(u) belonging to class k is estimated by simple kriging:

(u, k) E{I(u, k)} =

one could sequentially sample from N single-point con-

) is a subset of the (n + j) data found in a neigh-

) into a vector of indicator values: v(u

. The indicator variable is dened as:

(d) The estimated values i

) into a vector of indicator data values: v(u

. The indicator variable is dened by:

(u, k) E{I(u, k)} =

(d) The estimated values i

(Z(u) = k), after correction of order-

If the variable is categorical (ag Categorical Variable Flag set to 1),

where ? denotes a value that is not dened (missing data value)

}, then the prior

. In all rigor, the cokriging estimator should include

Updating of the prior probabilities Y (u, z

), j = 1, . . . , J, found at all locations u

of the training image, and

, = 1, . . . , n, of the training image such that D

the largest template size for which n n

is the indicator function:

/ [0, 1] it is reset to the closest bound. All updated

Conditioning to soft data

The conditional distribution function of step 2(c) of Algorithm5, P

using the formula proposed by Journel

Under the hypothesis of conditional independence:

Y(u) | Z(u)) (3.14)

using Eq. (3.15)

2. Retrieve the conditioning data event D

, = 1, . . . , n, of the training image rotated by angle r(u),

the largest template size for which n n

is the indicator function:

. This parameter is ignored if Use ProbField

You might also like