An implementation of conditional Markov mesh simulation
Contents
1. zi fi 25 Jat fe s 1 In the following two sections we explain how the f functions have been chosen and which input ranges that are to be determined by the user 4 7 1 Two point interactions For two point interactions we consider the strength between the variables x and all the variables in the two point interaction neighborhood For each facies k we i 2point assign an indicator function f x 2poin ft apene 1 if I tZk which yields e 0 otherwise where as ane the cell Figures 3 a c display the variables included in the two point interaction is a variable in the two point interaction neighborhood j labeling neighborhood The figures show one plane at the time where l ly and l is the An implementation of conditional Markov mesh simulation with parameter estimation 19 L y Figure 3 Illustration of the two point interaction neighborhoods for all three planes a x y plane b x z plane and c y z plane extension of the neighborhood in the corresponding direction Note that we only consider variables from the planes orthogonal to eachother For the total num ber of N cells in the two point interaction neighborhood of cell i we get the following set of explanatory variables 2poin 2poin 2poin 2poin on TE a a In the model file the user must give the parameters l ly and I for instance TW2. 1983 An implementation of conditional Markov mesh simulation with parameter estimation 11 For each cell i this expression is used to find the predictors Z xi xj lt i xw and Zi ges of 6 2 5 Seismic conditioning If seismic data is present this should also contribute to Eq 5 We assume that the seismic data is available in the form of a full facies probability cube i e that the seismic probabilities p x m with m being the seismic raw data are known for all cells i in the grid Equation 5 asks for the likelihood and we therefore calculate e PY o mos 0 where p x m is taken from the seismic facies probability cube and p x is the marginal facies probability found from the training image A weight factor 5 for the seismic potential is manually included which allows for some user s adjust ment of the importance of the seismic This is common procedure when taking into account seismic data With seismic data being present in addition to well data the adjusting factor Y of 6 thus is replaced by V zilzjai tw gt La valia 2w 10 with 3 0 1 If no well data is present V x 2 1y U x 2 6 Local update The modified Markov mesh model 7 is well suited for iterations A main mo tivation for establishing an iterative method is to render possible local update in an already existing grid configuration typically as a result of new well data or modified seismic potentials Iterations coul
3. 22 69 76 60 nr nr no Title Author Date Publication number Abstract An implementation of conditional Markov mesh simulation with parameter estimation Marita Stien Heidi Kj nsberg Odd Kolbj rnsen Petter Abrahamsen 19th June 2008 SAND 08 04 The main purpose of these notes is to provide documentation of version 0 1 of the Markov mesh implementation done for the Multipoint project The report is mainly meant for internal use in the project serving as a user s manual for running the program and as a guide for further improvements and modifications of the code Keywords Target group Availability Project Project number Research field Number of pages O Copyright Markov Mesh models sequential simulation parameter estimation data conditioning All employees MPS project partners Open Multipoint 808002 37 Norwegian Computing Center Contents 1 Introduction 0 0 222208 7 2 Main methods 2 2048 8 2 1 Markov mesh models 2244 8 2 2 000000 s s s g 225 46 8 4 B60 e ee ee oe Bo 9 2 2 1 Parameter reduction 4 10 23 Simulation algorithm 2 s oe wa amp BS w g 10 2 4 Well conditioning 2248 10 2 5 Seismic conditioning 2 222 a 12 2 6 Localupdate 4 2 repr ae GE A 12 3 Program structure LL 22222 ra aars 14 4 Input parameters 222 va varar a 16 4 1 Main f
4. 27 6 Results This section shows a few results that have been obtained by using the program The main objective is to illustrate possible usage not to discuss in detail the va lidity of the results All illustrations were created from the data of the program s output files using MatLab Unconditional simulation Figure 6 displays two different 2D training images one characterized by chan nels the other by more irregular sand objects These training images have been used to estimate model parameters as described in these notes and hence de fine two different Markov mesh models These Markov models have in turn been used for simulation be it unilateral unconditional or conditional or local update Figure 6 Two training images Figure 7 shows one random sample for each model Visual comparison to the training images shows that the implemented method for parametrization and parameter estimation in these cases is able to reproduce the characteristics of the training image quite well Figure 7 Results from unconditional simulation of estimated Markov mesh model 2 ME An implementation of conditional Markov mesh simulation with parameter estimation Unilateral conditioned simulation Figure 8 shows two random samples from conditional simulations using the mod ified unilateral Markov mesh model of Equation 7 Conditioning is done with respect to a well in the middle of each figure Comparing these two samples with the sa
5. all of size 12 4 4 Training image The training image must be either a 2D rectangle or a 3D box Under the com mand TI the first input to be specified is the name and location of the training image file The file format is a text file with no header containing only the facies values The values are separated either by space or line break and the order of cell listings are fastest in x direction and slowest in the z direction Next the size of the training image must be given see example TI C TrainingImages channelTI txt IDirectory and filename of TI 40 40 50 IX Y and Z dimension 3 4 5 Direction Path The program supports three main simulation paths determined by the command DIRECTION The first main path is the usual left to right top to bottom simulation path The second path is from left to right and right to left for every other row for each vertical layer Each of these two paths is suitable for sequential simulation The third choice of path is a random path which is useful for the kriging only version of the program see Section 4 1 In addition there is a fourth path choice This choice is useful for local update The subpath traversing the cells to be updated then varies randomly between a left to right top to bottom simulation path NorthWest and a right to left bottom to top SouthEast simulation path in each vertical layer The vertical layers are traversed from top to bottom TOP DIRECTION 10 One w
6. neighborhood A snapshot of a simulation is displayed and the grey cells have not yet been sampled We write the conditional probability for the facies at cell i as m xi tj lt i m zilzg 1 where zx is the set of facies values for the cells in the sequential neighborhood Markov mesh models are fully specified through the conditional probabilities in 1 ie the joint probability is N 21 Lo EN z 7Ceilzs 2 i 1 and in order to define the conditional probabilities we turn to the class of gener alized linear models The concept in generalized linear models is that a linear combination of inde pendent variables z is linked to the expected value y of the dependent variable 1 P McCullagh and J A Nelder Generalized Linear Models Chapman amp Hall 1989 gs HE An implementation of conditional Markov mesh simulation with parameter estimation y through a non linear function 270 n 1 We use the link function n m log 4 3 1 p We introduce the dependent variable y 1 when z k and O otherwise and let z be a vector consisting of functions of the sequential neighbourhood of node i and a constant term 2 f 25 fp as 17 The functions are discussed in more detail in Section 4 7 The expected value j can be written e 1 yt 1 2 00 s a for each facies k 0 1 K 1 The vectors 0 0 are of size P 1 x 1 and represent coefficients f
7. number of parameters is reduced by performing a matrix reduction of Z A matrix Z is the reduced matrix of size N x P where P lt P and is the product Z y where V isa P 1 x P matrix of singular vectors which correspond to the P largest singular values of the matrix Z The columns of Z are then linear combinations of the columns in Z The set 0 0 7 is the corresponding reduced set of coefficients We write ZO ZVO ZO Thus the coefficients 0 0 are estimated using Z and then the true coefficients 0 0 7 are computed from 0 VO 2 3 Simulation algorithm Simulation from the Markov mesh model is performed by following the path i 1 2 N throughout the grid For each cell the facies value is drawn accord ing to the conditional probability 7 x s It is crucial to the algorithm that when updating cell i only information about previously visited cells is taken into ac count no attention is payed to any cell along the future path The grid is scanned once and the resulting grid configuration follows the joint probability distribu tion N Til EN reilzs 4 i 1 2 4 Well conditioning In the case of unconditional simulation the algorithm of Section 2 3 produces grid configurations that are consistent with the Markov mesh probability distribution 4 This is not necessarily true if the simulation is to be conditioned on well data The reason is that updat
8. the program If Output prefix ends with backslash V it refers only to the output folder and all output files will be given their default names see Section 5 The file folder must be manually created before simulation REALIZATIONS 10 INumber of realizations 1 IHeader 1 true 0 false C MarkovMeshSimulations realization Output prefix 3 4 3 Simbox definition The simbox definitions command SIMBOXDEFINITION are for defining the area for simulation Only two dimensional rectangles or three dimensional boxes are supported The user must set the corner point x0 yO and z0 which is the upper left corner at the top of the cube Cell sizes are given by the parameters dx dy and dz and the number of cells by the parameters nx ny and nz Padding to avoid the edge effects is included in the program with default parameters Note that for two dimensional grids nz is set to 1 see example 3 H Soleng A R Syversveen and O Kolbj rnsen Comparing Facies Realizations Defining Met rices on Realization Space Proceedings of the 10th European Conference in the Mathematics of Oil Recovery 2006 An implementation of conditional Markov mesh simulation with parameter estimation 17 SIMBOXDEFINITION 000 1x0 yO zO 111 ldx dy dz 100 100 1 nx ny nz 2 The example above describes a simulation rectangle in two dimensions where the upper left corner is located at the origin There are 100 cells in each direction and they are
9. IGHBOURHOOD 8 X 8 Y 1 IZ 1 XY Diagnoals 1 XZ Diagonals 1 YZ Diagnoals 1 XYZ Diagonals TWOPOINT 6 6 XY plane 1x ly 1 1 XZ plane 1x lz 1 1 YZ plane ly lz 2 CONDITIONING 1 variogram type 0 spherical 2D 1 empiric 2D and 3D 30 Ivariogram range Rx empiric variogram range spherical 30 Ivariogram range Ry empiric variogram subrange spherical 10 variogram range Rz empiric variogram angle spherical 3 Imaximum number of kriging neighbours 11 lImaximum number of kriging neighbours forward used during iterations 11 maximum number of kriging neighbours backward used during iterations WELLDATA welldata txt Name and location of the welldata file 3 Example 2 WORKFLOW 0 estimate O only_kriging_ 1 luse_wells_ O use_seismic_ 1 do_iterate_ REALIZATIONS 1 INumber of realizations 0 IHeader 1 true 0 false C Project Output directory and prefix 34 TE An implementation of conditional Markov mesh simulation with parameter estimation SIMBOXDEFINITION 000 1x0 yO zO 1 1 1 dx dy dz 300 300 1 nx ny nz 2 TI C TrainingImages channelTI dat 250 250 1 Size x y z direction DIRECTION 10 IOne way simulation 11 IBack and forth simulation top to bottom 12 Random 3 Main path is one way subpaths vary between NWTOP and SETOP FACIES 2 Inumber of facies CONDITIONING 1 Ivariogram type 0 spherical 2D 1 empiric 2D and 3D 30 variogram rang
10. OPOINT 2 2 XY plane 1x ly 1 1 XZ plane 1x lz 1 1 YZ plane ly 1z 2 2 ME An implementation of conditional Markov mesh simulation with parameter estimation L Figure 4 Illustration of the strips of cells where higher point interactions are considered Arrows indicates the directions and in which order the number of interaction terms in creases x 3 point interaction x 4 point interaction x 5 point interaction x L point interaction Figure 5 Example of a strip of cells and the interactions that are included in the explana tory variables 4 7 2 Higher point interactions In addition to two point interactions three to L point interactions are also im plemented The size of L must be given by the user for each direction x y and z Not all possible higher point interactions within a certain neighborhood range are considered Certain strips of cells are selected see Figure 4 where each strip in one plane is indicated by arrows All three planes are similar Four additional strips are also included going from the center cell i diagonally out one step in each direction Figure 5 shows for one of the strips which interactions that are taken into account We let x be the set of cells included in a l point interaction for strip j with the indicator function lpoint E k 0 otherwise p alpont 1 if all cellsinx For the total number of R strips root
11. Section 4 1 An implementation of conditional Markov mesh simulation with parameter estimation 23 only kriging 1 filling in a coarse grid In each shell of the shell search the program identifies cells that are either future wells or have been filled in before The single parameter N under CONDITIONING defines the upper limit for how many cells to include only_kriging 0 use wells 1 do_iterate_ 0 sequential simula tion with hard data conditioning Before shell search is started the program accepts as kriging data future wells that are inside the box defined by the variogram range parameters If the number of these future wells exceeds N a subset of the wells is randomly chosen If the number of wells is lower than N the shell search algorithm starts and previsouly simulated cells are included until the limit N is reached or shell search is beyond variogram range only_kriging 0 use wells 1 do_iterate 1 local update of the initial configuration Two parameters are used for the iterative phase of the simulation N1 and N2 are maximal limits for the number of conditioning cells along the forward and backward path respectively Shell search is used for each direction independently 4 8 2 Well input file Well positions and well facies must be specified in an input file The name of the file is specified as in this example WELLDATA welldata txt Name of the welldata file 2 The file forma
12. ay simulation 11 Back and forth simulation top to bottom 12 Random 3 Subpaths vary between NW_TOP and SE_TOP 18 ME An implementation of conditional Markov mesh simulation with parameter estimation Be aware that regardless of simulation path only one set of parameters is esti mated This parameter set refers to the one way simulation path The use of any other path specification under DIRECTION implicitly assumes a symmetric model 4 6 Facies The number of facies in the training image and realizations must be given as an input If the only kriging version of the program is being used also the volume fractions of the facies must be provided Volume fractions are in all other cases automatically calculated from the training image FACIES 2 number of facies 0 723 Ivolume fraction facies 0 used iff only_kriging 1 0 277 Ivolume fraction facies 1 used iff only_kriging 1 3 4 7 Neighborhood If the estimate parameter is set to true neighborhood parameters must be pro vided If estimate is set to false and the parameters are needed for simulation only_kriging is set to false they are automatically read from the coefficient file see Section 5 We have divided the neighborhood into two categories two point interac tions and higher point interactions under the command names TWOPOINT and NEIGHBORHOOD respectively The explanatory variables in the conditional prob abilites are functions of neighborhood variables i e
13. cribes the data in terms of output files that are returned by the program A few results are for illustrative purposes included in Section 6 and a short summary and closing remarks are given i Sec tion 7 Appendix A provides examples of typical model files The user can on request to the authors be provided with examples of input files such that some simple models can be run right away An implementation of conditional Markov mesh simulation with parameter estimation 7 2 Main methods This section gives a description of the main methods used in the program It will be helpful for those interested to know the details of the various algorithms The main parts of the section involve model specification estimation data condition ing and local updating 2 1 Markov mesh models Consider a finite regular grid in two or three dimensions and let the one dimensional index i label the cells of the grid The set of all cells is 9 1 2 N where cell value x can take K different facies values i e x 0 1 K 1 Markov mesh models follow a sequential path Probabilities are constructed such that they depend on a subset of cell values from earlier in the path and this subset is called the sequential neighborhood Figure 1 gives an illustration of a sequential neighborhood on a two dimensional grid Sequential neighborhood mr Figure 1 Illustration of sequential
14. d also be used to clear unwanted kriging effects during the initial establishment of a reservoir configuration to en sure that it is consistent with the statistics of the prior model We propose to use a Metropolis Hastings algorithm where proposal configurations are established via block update based on the modified Markov model 7 and the accept prob ability ensures that sampling is done from the prior model conditioned on data In the following we describe one step in the Markov chain of the Monte Carlo simulation Let y be a label for the existing configuration of the grid i e the grid configu ration is x The proposal configuration will be denoted 1 Pick according to some rule a connected set of cells B C 9 where 9 is the full grid If the set B includes data cells or cells that for some other reason are supposed to have fixed facies throughout the iterations let the set of these cells be denoted Bo Let Bj BA Bo and define 08 as the set of cells in 9 B that according to the prior model may be affected by a change in the set B That is j OB iff dk 0 k By Then 2 Me An implementation of conditional Markov mesh simulation with parameter estimation B9QUOBi N j j gt i is the union of the set of future cells from a fixed location i that are within B and that we want to condition on and the set of future cells from position i that according to the prior model are affected by a change in the set B T
15. d in iterations 1 N2 the max number of kriging neighbours backward used in iterations 3 2 NE An implementation of conditional Markov mesh simulation with parameter estimation In the following we describe the details of the algorithms controlled by the pa rameters and thereby explain what each parameter means Correlation structure Two different variogram types are supported by our implementation spherical variogram and empiric variogram The choice of variogram type is controlled by the first parameter listed under CONDITIONING If empiric variogram is chosen the program automatically calculates the variogram in the x y and z direction from the input training image see Section 4 4 The spherical variogram is all synthetic with 2D range subrange and angle relative to the x axis specified as input parameters The variograms are subsequently used for the correlations that enter the kriging matrix and kriging vector with elements given by X j Cov j1 j2 and c Cov i j respectively The meaning of the three next parameters under CONDITIONING depends on the choice of variogram type If empiric variogram is chosen the parameters de note range R in x direction range R in y direction and range R in z direction For any fixed cell i these parameters define a box centered at cell i with side lengths 2R 1 2R 1 and 2R 1 This box defines the cells that are candi dates for being used as conditioning cells when doin
16. del Far away from any wells we expect V 1 and hence the new model based on equations 6 and 7 gives statistics similar to the prior Markov mesh model If cell i has a non zero correlation with a future well the kriging function V af fects the overall probability P A positive correlation implies that Z x j lt i w increases the probability for x to be updated to the same facies as the well a neg ative correlation decreases the well s contribution to this probability The overall effect of the well is modified by the past cells influence on the predictors Simple kriging is suitable when the mean volume fraction of each facies is given Let u be the mean volume fraction of facies k and let D be the set of cells on which to condition in the kriging algorithm Our choice for how to iden tify D is explained in Section 4 8 1 Let X be a matrix whose elements are the covariances between the different data points j j2 D Sjj Cov ji j2 and c is a vector of the covariance functions between the fixed update location i and the data points c Cov i j j D Define in addition the n x 1 vectors y vi Y YR and e 1 1 1 where n is the number of cells in D and Y k 1 when z k and 0 otherwise The predictor for facies k is then found via the linear expression Z xi klagen p eE gt ute 8 2 AG Journel Nonparametric Estimation of Spatial Distributions Math Geol 1503
17. e Rx empiric variogram range spherical 30 Ivariogram range Ry empiric variogram subrange spherical 10 Ivariogram range Rz empiric variogram angle spherical 3 Imaximum number of kriging neighbours 11 Imaximum number of kriging neighbours forward used during iterations 11 maximum number of kriging neighbours backward used during iterations WELLDATA welldata txt Name and location of the welldata file 2 LOCALUPDATE maskslocalupdate txt IName and location of the local update masks file ITERATIONS 1000 the number of iterations if do iterate true 25 Ilxmin 25 lymin 0 1zmin 50 1xmax An implementation of conditional Markov mesh simulation with parameter estimation 35 50 lymax 50 llzmax variable box to be used during local update size 2 1x 1 2 ly 1 2 1z 1 lonly cells listed in input file set in LOCALUPDATE will actually be updated 2 36 NE An implementation of conditional Markov mesh simulation with parameter estimation An implementation of conditional Markov mesh simulation with parameter estimation 37
18. e input parameters in the model file We use a model file structure which is divided into commands Each command starts with the command name written in capital letters and ends with a semicolon The input parameters to be specified by the user are written between the com mand name and the semicolon The exclamation mark is used in front of com ments i e everything written on the line behind an exclamation mark is not read by the program In some cases input files in addition to the model file are needed by the pro gram These must all be located on the root directory of the program 4 1 Main flow parameters The first parameters to be set determine which workflow to run using command WORKFLOW We refer to Figure 2 which shows various workflow actions and how one follows the other There are five boolean parameters to be set where 1 means true and 0 false 1 estimate Set to true if the coefficients 0 should be estimated or false if the coefficients should be read from an existing file Further parameter input is described in sections 4 4 4 5 4 6 and 4 7 2 only_kriging Set to true if coarse grid values from file should be filled into a finer grid The coarse grid values must be available via file See Section 4 11 for further details If the parameter is set to false sequential simulation is run either with or without conditioning See sections 4 2 4 3 4 4 and 4 5 for further parameter settings 3 use_wells Set to t
19. e of any cell i does not condition on data located along the future path from cell i Hence it may happen that when the simulation hits a well data point w the fixed value x may be inconsistent with the probability Fl 25 The probability that should be used for update of x is p i xj lt i tw Here W is the set of well data along the future path from cell i W being the set of all w Me An implementation of conditional Markov mesh simulation with parameter estimation well data cells We rewrite this as P TilZj lt i Tw p Tilzjei tw Playa 5 a plzil j lt i K Equation 5 can be considered a factorization of the posterior probability into well likelihood leftmost factor and prior probability rightmost factor The un conditional Markov mesh probability is the prior information whereas we have chosen to approximate the likelihood via indicator kriging That is for p x j lt i xw we use the approximation Z XilZj lt is Z1 Z TilTj lt i P axilrj lt i 2w m 2 5 V ailzjen tw m zijzs 6 where Z z xj lt is the predictor for x found by indicator kriging conditioned on cells in the past of i Z xi xj lt i xw is the predictor conditioned also on fu ture data points and 7 is the unconditional Markov mesh probability The joint probability J Peile 2w 7 is then an approximation to the true Markov mesh model and defines what we will refer to as the modified Markov mesh mo
20. ed at cell i we get the following sets of ex An implementation of conditional Markov mesh simulation with parameter estimation 21 planatory variables ENN a ae nd a EE A onl OO EN In the model file the user must give the parameter L for the various directions both along the three main axes and the diagonal lengths for instance NEIGHBOURHOOD 2 X direction 2 Y direction 1 Z direction 1 XY 2D diagnoals 1 XZ 2D diagonals 1 YZ 2D diagnoals 1 XYZ 3D diagonals 3 4 8 Well conditioning There are two aspects of the kriging algorithm for hard data conditioning wells that are controlled by parameters the correlation structure used to set up X and c see Eq 8 and the number of cells from which the predictor is calculated These parameters are set under the model file command CONDITIONING In addition the locations and facies of wells must be provided via an input file that is specified under the command WELLDATA 4 8 1 Correlation structure and shell search The parameters that determine the kriging details command CONDITIONING are specified in the model file as CONDITIONING 1 variogram type 0 spherical 1 empiric 30 variogram range Rx empiric variogram range spheric 30 variogram range Ry empiric variogram subrange spheric 10 variogram range Rz empiric variogram angle spheric 1 N the max number of kriging neighbours 1 N1 the max number of kriging neighbours forward use
21. es LOCALUPDATE provided the cell is not a well 2 Pick randomly the size parameters l ly and l such that each parameter is in the interval I Imin 1max 3 Define a box with center in the picked cell and sides of lengths 21 1 21 1 and 2l 1 The set of cells within the box constitute B 4 The set of cells that are updated during this iteration step is given by Bj UC 4 11 Fill in from coarse grid If the main flow parameter only kriging is set to true coarse grid values from file should be filled into a finer grid The coarse grid values must be available in a file called coarsegrid txt The format of this file should be identical to the header less output file from sequential simulation see sections 4 2 and 5 and the number of values should in 3D be 1 8 of the size of the present simbox see Section 4 2 The program will krige unsampled cells in the finer grid 2 NE An implementation of conditional Markov mesh simulation with parameter estimation 5 Program output Depending on the main flow parameters the program produces various output information in the form of files The names of these files all start with the prefix specified in the command REALIZATIONS All relevant file directories must be set up before simulation Apart from the prefix the output file names are If estimate true coeff txt Estimated model parameters The file header contains the se quential neighborhood parameters If only_
22. g kriging for cell position i Cells inside the box are taken into account provided they fulfil various crite ria explained below cells outside the box are not taken into account Hence the ranges R Ry and R should have values that are roughly the same as the correla tion lengths of the TI The output file correlations txt see Section 5 contains the correlations of the training image and can be used to identify appropriate values for R Ry and R If spherical variogram is chosen the three parameters signify 2D range sub range and angle relative to the x axis From these parameters the program auto matically calculates a box that encompasses the variogram range This box de fines cells that are considered candidates for playing the role as kriging data Shell search To identify the cells D used to calculate the predictor Z x k z jeny in Equation 8 the basic algorithm is a shell search algorithm That is we search outwards from cell in cubic shells if 2D in quadratic shells and identify in each shell the cells that either are future wells or have already been simulated though with some minor requirements that are listed below We continue until the shell is entirely outside the box defined by the range parameters described above or until the number of found cells reaches a set limit whichever occurs first The details of the shell search algorithm differs slightly depending on the flow parameters of the program see
23. he new configuration y is established as follows scan through the part of the simulation path that is within B and for each cell i B draw its new facies value x for configuration u from the probability P xilEj lt i TW ULBUOBi NLD T Lal Lo 11 for all cells in 9 B let the facies be as in the state v This defines the proposal configuration u Acceptance or rejection of the configuration v is done with probability a min 22 i 12 qual where q is the probability for suggesting the new configuration y starting from V q is the probability for suggesting the old configuration y starting from y and r x and 7 x are the prior Markov mesh probabilities Since all suggested states conform with data this accept probability ensures that sampling is done from the prior model conditioned on data The accept probability can be rewritten min VU za v T zi u o 1 Ta A T x l 1 13 In the notation of equation 13 the conditional dependencies in the function ar guments are suppressed for readability An implementation of conditional Markov mesh simulation with parameter estimation 13 3 Program structure The main task supported by the present implementation are Based on an input training image estimate parameters for a Markov mesh model Create a number of independent realizations from the Markov mesh model by using sequential simulation Each realization can be conditio
24. kriging true 0 txt Resulting configuration of the fine grid If only_kriging false and use wells false true and use_seismic false true lt realization number gt txt Resulting grid configuration for the indicated realization number If use_wells true correlations txt Correlations of the training image If do_iterate true lt realization number gt txt Final grid configuration for the indicated real ization number lt realization number gt _ lt iteration number gt after txt Intermediate grid con figurations written every 10th iteration lt realization number gt _faciesCount_n lt iteration number gt _N lt total number of iterations gt txt Facies count for each facies for each cell Written every 1000th iteration and at the end of the iterations The listing is according to the rule first run over cell index then run over facies accept txt Statistics for the accept rate Updated each 1000th iteration and at the end of the iterations each line being on the format lt realization number gt lt number of tried iterations gt lt number of accepted iterations gt In addition for future reference all relevant input files are automatically copied to the chosen output directory The possible files are model txt coeff txt welldata txt seismicattributes txt initialgrid txt and maskslocalupdate txt An implementation of conditional Markov mesh simulation with parameter estimation
25. low parameters LL 2222 arr 222 16 4 2 Realizations ios a amp Bode de me r 17 4 3 Simbox definition 17 4 4 Trainingimage e a 4 eee eS Ga 18 45 Direction Path 2 gue RO a we eS 18 46 Faces 64 wee FENG ee Eee s EM ow OR we ow 19 4 7 Neighborhood aa aa 19 4 7 1 Two point interactions 19 4 7 2 Higher point interactions 21 48 Well conditioning 2 22 aar 22 4 8 1 Correlation structure and shell search 22 48 2 Well input file 24 4 9 Seismicconditioning 24 4 10 Localupdate 2 0 ee ew 0 ee er AG 25 4 11 Fill in from coarse grid S ss ke be Ke ee 26 5 Program output aa aaa a 27 6 Results hs s bw de E See e Ek O ene 28 7 SUMMA lt 4 coke ewe bed e dog SEA GE A 32 A Examples of model files 33 An implementation of conditional Markov mesh simulation with parameter estimation 5 1 Introduction These notes contain documentation of a prototype implementation of a Markov mesh model The model is developed as part of a larger project for development of models for multipoint statistics and is attractable due to its fast sequential simulation scheme The program implementation can perform the following tasks 1 estimate a statistical model from a training image 2 simulate realizations from a given model
26. mand specifies the number of iterations and the possible sizes of the update boxes used during the iterations The input file must be specified as in this example LOCALUPDATE maskslocalupdate txt IName of the local update mask file The file format is assuming n cells are allowed updated x y1 z1 x2 y2 z2 x3 y3 z3 xn yn zn The coordinates x y z must refer to the same unit system as is used to define the simulation box see Section 4 2 An implementation of conditional Markov mesh simulation with parameter estimation 25 Input parameters that control the iteration algorithm are specified as follows ITERATIONS 1000 Ithe number of iterations if do iterate true 25 1xmin 25 lymin 0 1zmin 50 1xmax 50 lymax 50 1zmax The first parameter describes the number of iteration steps that are to be car ried out The remaining parameters are understood as follows Each step of the Markov chain that defines the iterative process uses a connected set of cells B C 9 where is the full grid See Section 2 6 for a detailed description of the theory be hind the method The present implementation uses a rectangular box for B The six parameters 1xmin 1zmax defines this box according to the following algo rithm which is used once for each step of the Markov chain 1 Pick randomly one of the cells that are allowed to change facies during the it erative phase i e a cell from the set specified by the input fil
27. mples from the unconditional Markov mesh models see figure 7 there is no observable statistical difference between the samples of the prior models and the modified models The only difference is that the modified Markov mesh models locate sand objects correctly around the well position Figure 8 Results from conditional simulation using indicator kriging and estimated prior model Conditioned on object facies well in cell the center of the figure the well being indicated with a circle We can investigate the statistical properties in more detail by simulating many realizations from one model and then summarize the results In particular the statistics of the modified Markov mesh model can be compared to the statistics re sulting from rejection sampling of realizations created by the prior Markov mesh model The approximation done when using indicator kriging is to be consid ered satisfactory if rejection sampling and modified Markov model give the same statistical results Figure 9 displays the conditional marginal probabilities for the channel model conditioned on object well top row and background well bottom row In each case the well is located in the middle of the grid The results from rejection sam pling are shown in the left column The right column displays the results of con ditioned simulation using the modified Markov mesh model with indicator krig ing Comparison to the rejection sampling gives the impression that the
28. ned on wells and or seismic data In addition the program can Do local update on an existing realization Instead of estimation and subsequent sequential simulation the user can pro vide a grid realization on a coarse scale then ask the program to fill in the grid on a finer scale Figure 2 illustrates the logical structure of the program from a user s perspective Parameters that control the behavior of the program are described in the next section 14 Me An implementation of conditional Markov mesh simulation with parameter estimation read coarse grid from file fill in using default method lt facies fractions true coarse grid or a set of wells gt exit read estimated parameters from file do estimation lt neighbourhood parameters Tl gt read wells from file read data from file read existing realisation from file false1 true2 false3 false1 true2 true3 or or truet false2 false3 truet false2 true3 falset false2 false3 false1 false2 true3 or or truet true2 false3 truet true2 true3 do sequential unconditioned iterate do sequential conditioned iterate conditioned simulation unconditioned simulation lt path gt lt path gt lt path gt lt path gt exit exit exit exit Figure 2 Graphical illustration of the program structure An implementation of conditional Markov mesh simulation with parameter estimation 15 4 Input parameters In this section we explain how to set th
29. nn m Norsk Regnesent or NORWEGIAN COMPUTING CENTER An implementation of conditional Markov mesh simulation with parameter estimation Overview and user s manual for version 0 1 SAND 08 04 Marita Stien Heidi Kj nsberg Odd Kolbj rnsen Petter Abrahamsen 19th June 2008 En N Norsk Regnesentral e NORWEGIAN COMPUTING CENTER Norwegian Computing Center Norsk Regnesentral Norwegian Computing Center NR is a private indepen dent non profit foundation established in 1952 NR carries out contract research and development projects in the areas of information and communication tech nology and applied statistical modeling The clients are a broad range of indus trial commercial and public service organizations in the national as well as the international market Our scientific and technical capabilities are further devel oped in co operation with The Research Council of Norway and key customers The results of our projects may take the form of reports software prototypes and short courses A proof of the confidence and appreciation our clients have for us is given by the fact that most of our new contracts are signed with previous customers Norsk Regnesentral Bes ksadresse Telefon telephone Internett internet Norwegian Computing Center Office address 47 22 85 25 00 www nr no Postboks 114 Blindern Gaustadall en 23 Telefaks telefax E post e mail NO 0314 Oslo Norway NO 0373 Oslo Norway 47
30. ns The figure illustrates that along the edges of the update rectangle A the probabilities follow the conditioning provided by the cells out side the rectangle Also conditioning to the well is good The other areas of the local update rectangle provide evidence that the mixing of the iterative process is satisfactory This conclusion follows from the fact that the cell wise statistics tend to smear out and match the marginal probability p x 1 only broken by areas with a higher channel probability areas that to a large extent obviously follow from the edge and well conditioning Figure 11 Left Initial grid configuration before local update middle snapshot of the grid configuration during iterative process right cell wise probability p x 1 after 3000 iterations The well position is indicated but not conditioned to in the initial grid The iterative process conditions on the well which has facies x 1 The rectangle marks the area of local update An implementation of conditional Markov mesh simulation with parameter estimation 31 7 Summary These notes provide documentation for a prototype implementation of a Markov mesh model The implementation can estimate a statistical model from a training image simulate realizations from a given model either unconditioned or condi tioned to well s and or seismic data update a given realization locally or down scale a realization The theoretical formulation of the implemen
31. or each facies corresponding to the functions in z Combining the above result with 3 yields the following expression for the con ditional probabilities Tok rly lzi 8 ae aa J ka Zo exp 2 0 2 2 Estimation The model coefficients are estimated by the maximum likelihood estimator A training image contains the data basis for the estimation The likelihood function is given by L 0 0 Z J I explzt0 ut LN where Z is an N x P 1 matrix of the full set of observations z and Y is an N x K matrix of the full set of observations y for each facies Simple derivations yield the following system of equations to be solved Z y ZT p 0 a ot for k 0 1 K 1 The vector u 0 0 is the N x1 vector of u 09 057 while y yf y3 y The above expressions are solved using the standard approach of iterated weighted least squares For each iteration the following ex pressions are computed OF 0k ZWEZ IZ yt 0 0 for k 0 1 K 1 Here W diag 1 18 0m it Om 1 p Om uk Om is an N x N matrix of weights An implementation of conditional Markov mesh simulation with parameter estimation 9 2 2 1 Parameter reduction The number of coefficients to estimate is P for each facies It is often favorable to reduce the number of parameters We do this by removing the parameters that have least influence on the model More specifically the
32. overall results of the modified Markov mesh model are satisfactory Figure 10 displays results for conditioning on a channel edge i e a transition from an x 0 well to an x 1 well Comparison to rejection sampling left shows that the main features are well reproduced by the modified model right In particular the results are good close to the wells An implementation of conditional Markov mesh simulation with parameter estimation 29 1 80 80 os 90 90 os 400 100 o7 110 110 os 120 120 os 130 130 o4 140 140 os 450 150 02 160 160 01 170 170 80 90 100 90 120 130 140 160 160 170 80 90 100 110 120 130 140 150 160 170 gt 1 80 20 os 90 so os 100 100 o7 110 110 os 120 120 os 130 130 o4 140 140 o3 450 150 02 160 160 o1 170 170 so 99 100 110 120 130 140 150 160 170 80 90 100 110 120 130 140 160 160 170 e Figure 9 Conditional marginal statistics displayed as the cell wise probability p x 1 Left column rejection sampling from prior right column indicator kriging with respect to the well and two cells in the past path Top row well has facies x 1 bottom row well has facies x 0 19 60 70 80 90 100 110 120 130 140 15 ka 60 70 amp 90 100 110 120 130 140 150 Figure 10 Conditioning on two neighboring wells one with x 0 in position 101 100 the other with x 1 in position 101 101 Cell wise probability p x 1 Left rejection sampling Right Modified unilate
33. ral Markov mesh model Local Update The iterative model can be used to perform local update of an existing realization We illustrate this in Figure 11 The left hand side of the figure shows an example of an initial grid configuration created by unconditional simulation This realiza tion is to be updated with respect to a new well in grid position 126 126 The well position is indicated with a circle at the center of the figure The value of the well is supposed to be x 1 which is not in agreement with the initial grid configuration Local update is to be carried out in an area A defined by the rect angular box shown in red in the initial grid All cells outside A are fixed The iterative method using block update is now applied to this problem The region B used to generate a suggestion state for the Metropolis Hastings algorithm is for each element of the Monte Carlo chain a randomly drawn rectangle and the cells 3 HE An implementation of conditional Markov mesh simulation with parameter estimation that satisfy i BN A are updated according to the algorithm The middle grid of Figure 11 displays a snapshot of the grid during the itera tive process The channels have now been adjusted such that they are consistent with the well data In addition the channels nicely match the fixed part of the grid outside the marked rectangle At the right hand side of Figure 11 we display the cell wise statistics of the grid after 3000 iteratio
34. rue if wells should be taken into account and false if there is no well conditioning Further parameter input is described in Section 4 8 4 use_seismic Set to true if seismic should be taken into account and false if there is no seismic conditioning Seismic input information is described in Section 4 9 5 do_iterate Set to true if local update is to be performed If so the initial grid must be provided in a file named initialgrid txt Further parameter input is described in Section 4 10 If false then no local update is to be done 6 Me An implementation of conditional Markov mesh simulation with parameter estimation Example of workflow parameters WORKFLOW O lestimate_ 1 true 0 false O only_kriging_ 1 true 0 false 1 luse wells 1 true 0 false 1 use_seismic_ 1 true 0 false O do_iterate_ 1 true 0 false 2 4 2 Realizations The number of realizations the output file format for writing the grid configu rations to file and the location and filename prefix for the realizations must be given under the REALIZATIONS command in the model file There are two options for the file format either a header is included or not A header must be included if the realizations are to be run through the analysis program FaProp gt otherwise only facies values from 0 to K 1 are printed The program adds realization number and a txt ending to the output prefix The prefix is also used for all other output files from
35. t is assuming n well positions with facies values ky k2 kn x1 y1 zi Kl x2 y2 z2 k2 x3 y3 z3 k3 xn yn zn fn The coordinates x y z must refer to the same unit system as is used to define the simulation box see Section 4 2 4 9 Seismic conditioning If seismic conditioning is to be done the name of a file specifying the seismic prob ability cube must be provided In addition the seismic weight factor must be set see Eq 10 The following example illustrates how to set the seismic parameters 24 ME An implementation of conditional Markov mesh simulation with parameter estimation SEISMICATTRIBUTES 0 5 10 lt geis factor lt 1 seismicattributes txt Name of the seismic attributes file 3 The file format is assuming K facies values p x1 0 p x2 0 p xN 0 p xi 1 p x2 1 p xN 1 p x1 K 1 p x2 K 1 p xN K 1 Here x1 2 Ty are facies variables one for each cell of the grid N is the total number of grid cells listed in the file and this number must be consistent with the simbox definitions of Section 4 3 The indices 1 2 N label the cells according to the ordering first run over x then over y last over z 4 10 Local update Local update requires two commands to be specified in the model file LOCALUPDATE and ITERATIONS The first command simply specifies the name of the required in put file listing the cells where update is allowed to happen the second com
36. ted methods is provided as well as a detailed overview of the program s input parameters and main work flows A few results were included for illustrative purposes 32 NE An implementation of conditional Markov mesh simulation with parameter estimation A Examples of model files We give here two examples of model files In the first example the program is supposed to estimate parameters based on a training image then make 100 re alizations conditioned on well data In the second example we show the model file that will take a grid realization as input then perform local update on this realization In reality a likely situation is that Example 2 is carried out some time after Example 1 and that one of the realizations created by Example 1 is desired updated due to new well data Example 1 WORKFLOW 1 lestimate_ O only_kriging_ 1 luse wells O luse_seismic_ O do_iterate_ 3 REALIZATIONS 100 Number of realizations 0 Header 1 true O false C Project Output directory and prefix SIMBOXDEFINITION 000 x0 yO zO 111 dx dy dz 300 300 1 nx ny nz 3 TI C TrainingImages channelTI dat 250 250 1 ISize x y z direction DIRECTION 0 One way simulation 11 IBack and forth simulation top to bottom 12 Random 13 Main path is one way subpaths vary between NWIOP and SETOP FACIES 2 number of facies An implementation of conditional Markov mesh simulation with parameter estimation 33 3 NE
37. unconditional conditioned to well s and or seismic data 3 update a given realization locally 4 down scale a realization Facies grids can be in two or three dimensions and contain multiple facies In principle there are no limitations to the methodology but in practice it is limited by the memory on the computer Local updating is a post process that involves resimulation in local areas of the grid by using an McMC algorithm and is typ ically used to modify existing realizations when new well data or new seismic potentials are available Down scaling is performed to create a model on a finer level than the present realization In order to run the program the user needs to specify a model file containing all necessary model parameters The parameters are explained in Section 4 with the basic theory of the applied methods being presented in Section 2 It is not crucial to read and understand the theory of Section 2 in order to properly set the input parameters as Section 4 is attempted being self explanatory The organization of the documentation is as follows Section 2 Main methods presents the theory related to the model construction estimation conditioning and the local updating This establishes the notation Section 3 Program struc ture illustrates the main tasks enumerated above showing how they are con nected Section 4 Input parameters explains the input parameters needed to run the program Section 5 Program output des
Download Pdf Manuals
Related Search