Introduction

The ingredients for a plurigaussian simulation are:

• The number of lithotypes to be simulated

• Their proportions - constant in the domain for a stationary case or variable in space for a non stationary case

• Their relationships given by the lithotype rule which involves one or two underlying Gaussian Random Function (GRF)

• The variogram models of the underlying GRF(s)

When all these parameters are defined, the non-conditional simulations can be performed using the classical parameters such as the seed, the number of bands (when the GRF simulation is performed using the Turning Bands method) and the number of simulations. Finally all these simulations can match the (gaussian) data when available: they are then called conditional simulations. This feature is not addressed in this case study.

Define the simulation grid

The non-conditional simulations will be stored in a Db (called outgrid1) organized as a 2-D grid with 200 nodes along X and 100 nodes along Y (square mesh of 1m side) and which is created in the next command:

outgrid1<-db.create(flag.grid=T,x0=c(0,0),nx=c(200,100),dx=c(1,1))

Number of lithotypes to be simulated

For this illustration we have chosen four lithotypes. Their colors are given in the following colormap:

pal4fac <- c("red","orange","green","blue")

Lithotype Rule

The relationships between lithotypes are summarized in a lithotype rule. It consists of a 2-D diagram (represented as a square) where the horizontal axis gives the range of values for the first GRF and the vertical axis for the second GRF. Both GRF vary between[). This square is subdivided into as many rectangles as there are lithotypes. The procedure rule.input allows the user to dene the lithotype rule interactively. The first question concerns the simulation method:

• The Plurigaussian simulation: this corresponds to the current illustration, where the two underlying GRFs are not correlated

• The Shifted simulation

• The Shadow simulation

The rule construction is performed sequentially, describing the rectangles starting from the lower left corner. The rule is read from left to right and from bottom to top. Each rectangle is defined by the ordered thresholds. Therefore S the user must choose which threshold is encountered: (if the next limit is T defined along the first GRF) or (if the next limit is defined along the second F GRF). When the rectangle is defined, the user chooses to dene the facies number.

rule4fac = rule.create(c("S","F1","T","F2","S","F3","F4"))

The visualization of the lithotype rule is obtained using the following command:

plot (rule4fac,col=pal4fac)

Variogram Models

In the scope of the Plurigaussian model, the lithotypes are obtained by thresholding the underlying GRFs. Consequently, the variograms of the lithotype indicators are related to the models of these GRFs, taking the proportions and the lithotype rule into account. The fitting procedure is a trial and error method: we postulate the model of the underlying GRFs and check that the quality of the t for the variograms of the corresponding lithotype indicators. Here, the model is assumed to be known and is defined interactively.

The model of each underlying GRF is defined using the command model.input. The model for the first GRF is an anisotropic Cubic structure and an isotropic Exponential structure for the second GRF. The models for both underlying GRFs are stored in the file modY1 and modY2.

modY1 = model.create(vartype="Cubic",sill=1,aniso.angles=-60,range=c(35,20))
modY2 = model.create(vartype="Exponential",sill=1,range=15)

We check the contents of these models:

modY1

## 
## Model characteristics
## =====================
## Space dimension              = 2
## Number of variable(s)        = 1
## Number of basic structure(s) = 1
## Number of drift function(s)  = 1
## Number of drift equation(s)  = 1
## 
## Covariance Part
## ---------------
## - Cubic
##   Sill        =      1.000
##   Anisotropy  =      1.000     0.571
##   Ranges      =     35.000    20.000
##   Rotation    =
##                 [,1]      [,2]
##       [1,]     0.500     0.866
##       [2,]    -0.866     0.500
##   Angle       =   -60.000 degrees
## Total Sill    =      1.000
## 
## Drift Part
## ----------
## Universality Condition

modY2

## 
## Model characteristics
## =====================
## Space dimension              = 2
## Number of variable(s)        = 1
## Number of basic structure(s) = 1
## Number of drift function(s)  = 1
## Number of drift equation(s)  = 1
## 
## Covariance Part
## ---------------
## - Exponential
##   Range       =     15.000
##   Theo. Range =      5.007
##   Sill        =      1.000
## Total Sill    =      1.000
## 
## Drift Part
## ----------
## Universality Condition

Stationary PluriGaussian Simulation

The simulation is performed and stored in the output grid. In the following example, a single non-conditional simulation is performed with:

• the lithotype rule (rule4fac)

• the models of the underlying GRFs (modY1 and modY2)

• the constant proportions: here, they are considered as all equal to 0.25 (as they must sum up to 1)

• 600 turning bands and a seed equal to 126543.

outgrid1 <- simpgs(dbout=outgrid1,rule=rule4fac,
                   model1=modY1,model2=modY2,props=c(0.25,0.25,0.25,0.25),
                   nbsimu=1,nbtuba=600,seed=126543)

The display of the simulated outcome is obtained with the following command:

plot(outgrid1,col=pal4fac,flag.scale=TRUE,
     title="Stationary PluriGaussian Simulation")

Non-stationary PluriGaussian Simulation

For illustrating the non-stationary case, a second output grid (called outgrid2) is created.

outgrid2<-db.create(flag.grid=TRUE,x0=c(0,0),nx=c(200,100),dx=c(1,1))

The proportions of each lithotype vary in the space: they are stored from the grid. In our example, four variables are defined (called Prop1, Prop2, Prop3 and Prop4) which contain the proportions. For instance:

outgrid2 <- db.add(outgrid2,Prop1=x1/200) 
outgrid2 <- db.add(outgrid2,Prop2=0.1*(1-Prop1)) 
outgrid2 <- db.add(outgrid2,Prop3=0.5*(1-Prop1-Prop2)) 
outgrid2 <- db.add(outgrid2,Prop4=1-Prop1-Prop2-Prop3)

The type (locators) of these variables (located in fields 4 to 7) must be defined as «p» :

outgrid2 <- db.locate(outgrid2,seq(4,7),loctype="p")
outgrid2

## 
## Data Base Characteristics
## =========================
## 
## Data Base Summary
## -----------------
## File is organized as a regular grid
## Space dimension              = 2
## Number of fields             = 7
## Maximum Number of attributes = 7
## Total number of samples      = 20000
## Grid characteristics:
## Origin :      0.000     0.000
## Mesh   :      1.000     1.000
## Number :        200       100
## Angles :      0.000     0.000
## 
## Variables
## ---------
## Field =   1 - Name = rank - Locator = NA
## Field =   2 - Name = x1 - Locator = x1
## Field =   3 - Name = x2 - Locator = x2
## Field =   4 - Name = Prop1 - Locator = p1
## Field =   5 - Name = Prop2 - Locator = p2
## Field =   6 - Name = Prop3 - Locator = p3
## Field =   7 - Name = Prop4 - Locator = p4

The different proportions can be displayed using:

plot(outgrid2,name="Prop1",pos.legend=1,
     title="Proportion of First Lithotype")

plot(outgrid2,name="Prop4",pos.legend=1,
     title="Proportion of fourth Lithotype")

The rule and the models for the two underlying GRF are kept the same as in the stationary case. The command for the non stationary simulation is:

outgrid2 <- simpgs(dbout=outgrid2, rule=rule4fac, dbprop=outgrid2,
                   model1=modY1,model2=modY2, nbsimu=1,
                   nbtuba=600,seed=126543,verbose=TRUE)

The corresponding result is displayed in the next figure:

plot(outgrid2,col=pal4fac,flag.scale=TRUE,
          title="Non Stationary PluriGaussian Simulation")

PluriGaussian Simulations

D. Renard