A Multidimensional Markov Chain Model for Simulating Stochastic Permeability Conditioned by Pressure Measures

In this paper, we are interested in simulating a stochastic permeability distribution constrained by some pressure measures coming from a steady flow (Poisson problem) over a two-dimensional domain. The permeability is discretized over a regular rectangular gird and considered to be constant by cell but it can take randomly a finite number of values. When such permeability is modeled using a multidimensional Markov chain, it can be constrained by some permeability measures. The purpose of this work is to propose an algorithm that simulates stochastic permeability constrained not only by some permeability measures but also by pressure measures at some points of the domain. The simulation algorithm couples the MCMC sampling technique with the multidimensional Markov chain model in a Bayesian framework.


INTRODUCTION
Reservoir characterization aims to model the heterogeneous elements in a reservoir using subsurface and outcrop data.In the recent years, many techniques have been developed to describe and quantify the heterogeneity in a reservoir with great precision.These methods are divided into two categories: the kriging methods and the Markov chain methods.Both of these two categories are interpolation methods that seek to estimate the values of a stochastic permeability field over the reservoir domain by interpolating some permeability measures at some points of the domain.
Kriging methods are linear least squares methods.
For more details on the kriging method see ( [8]).The drawback of these methods is that they assume some linear properties of the permeability and a linear system is solved at each point where the permeability has to be estimated.Markov chain interpolation methods have the advantage of being less expensive in computations than the kriging methods because the permeability is simulated directly from the transition probabilities (these methods are described later in details).Another advantage of these methods is that the transition probabilities are easier to interpret than the variograms used for the kriging methods and they give more flexibility to be extended to more general interpolation problems.
In this paper, we are interested in simulating a stochastic permeability distribution over a domain given an underground water flow.More specifically, we consider a bounded domain Ω and the following Poisson problem: (1) where p is the pressure in the domain, g is a source term, h is the Diriclet boundary condition and K(x) is the permeability field over the domain Ω.To solve numerically this problem, a finite volume procedure is used and the domain is discretized by a regular rectangular grid of nodes.Let K ~be the discretized permeability over the grid.K ~is a two-dimensional array of values.The discretized permeability is considered to be stochastic but it can take a finite number of values and it is constant inside each cell.Let S l , l = 1 … n be the set of possible values of each element of K ~and there are n possible values.
We are interested in simulating such permeability field, K ~, constrained by some hard data (permeability measures κ) and soft data (pressure measures π) at a set of points ξ of Ω.The inverse problem of estimating such permeability is also considered.It is known that there is no guarantee on the uniqueness of this inverse problem, in the least squares sens, given that the unknown permeability function K(x) is not a smooth function neither its derivatives.Therefore, there is a non-singleton set K of the solutions of the least squares inverse problem.K is the set of the potential permeabilities that correspond to the pressure measures.The purpose of this work is to propose an algorithm that simulates stochastic permeability K ~that belong to K i.e. constrained by permeability and pressure measures (κ, π).The multidimensional Markov chain permeability model can be used to simulate K ~constrained only by permeability measures (κ).Our algorithm couples this model with the Markov Chain Monte Carlo (MCMC) sampling method in order to simulate K ~constrained by permeability and pressure measures (K ~from K).
The efficiency of this coupling algorithm is studied by computing the uncertainty in the permeability given the permeability and the pressure measures.The uncertainty in the permeability at cell (i, j) of the mesh grid is represented by the probability Pr [K ~ij = S l |κ, π], which is computed using our algorithm.
The paper starts by discussing in section (2) the general concept of the one dimensional Markov Chain for modeling one-dimensional permeability and how it has been extended to the multidimensional case to simulate stochastic permeability constrained by some permeability measures.Next, in section (3) the Bayesian framework is introduced followed by the coupling algorithm to simulate stochastic permeability constrained by both permeability and pressure measures.Finally, some numerical experiments are studied to show the performance of the coupling algorithm.

THE ONE-DIMENSIONAL CASE
We start by presenting the one-dimensional Markov chain model that is used to describe a stochastic permeability distribution over a one-dimensional domain (for example a vertical borehole).The domain is discretized using evenly spaced cells.Consider a sequence of random variables indexed by the space (Z 0 , Z 1 ,..., Z N ) where all these random variables can take a finite number of values.The random variable (Z k ) represents the permeability inside the cell k of the discretized domain (see Fig. l).Let S be the set of all the possible values of the permeability along the domain and n = card (S).(Z k ) is called a Markov chain if it verifies the following property: (2) In other words, if the permeability values Z k are simulated iteratively from k = 0…N the value of the permeability at the current cell k does not depend on all the previously simulated values (Z i , i = 0…k − 1) but only on its nearest simulated neighbor, Z k − 1 .
Therefore, to be able to simulate a one-dimensional Markov chain one has to define the transition probabilities, P ij = Pr(Z k + 1 = S j |Z k = S i ), the probability that the permeability jumps from the value S i to the value S j between two consecutive cells.These transition probabilities can be written together in a matrix form called transition matrix: (3) One must verify the constraint on the lines of this matrix, They can be estimated experimentally by counting the number of transitions occurred from value S i to value S j , T ij , along some boreholes and defining (4) These transition probabilities are called single step transitions because they relate the permeability values of two consecutive cells.One can consider also N-steps transitions and compute P N ij = Pr(Z k +N = S j |Z k = S i ), the probability of having the value S j at cell k + N if its nearest simulated neighbor value is S i at cell k.One can easily demonstrate that this probability is the (ij) th element of the matrix P (N ) .The N -step transition probabilities can be obtained by multiplying the single-step transition probability matrix by itself N times.

THE TWO-DIMENSIONAL CASE
In the previous section the one-dimensional Markov chain model was presented to describe the correlation between the value of the random permeability and its nearest known neighbor  value.This section is a brief introduction to the two-dimensional Markov chain model used to simulate a random permeability over a regular rectangular grid of cells.It is a generalization of the preceding model for the case where the value of the permeability is correlated by many known neighbors in different directions.The simulation algorithm is based on the conditional independence of the known nearest neighbors at a given location.
This assumption states that to simulate the value Z of the permeability at some cell of the grid given that it has m known neighbors in the four cardinal directions Z i ,i = 1, m it is assumed that each Z i has only an effect on Z and not on the other neighbors of Z: (5) This assumption does not have a mathematical basis in the theory of Markov fields (see [3]), but given its practicality and that its satisfactory results, it is used in many geostatistical applications.
The probability of Z given its known nearest neighbors Z 1 ,…, Z m is derived as follows: with .Using the conditional Independence as stated above, we have: and finally: Let h i be the number of cells that separates Z from its neighbor Z i .Each of the linear sequences Z to Z i constitute a one-dimensional Markov chain with a specific transition matrix (see [10]).Given that each cell can have maximum four neighbors in the four cardinal directions, four transition matrices must be defined P m , m = 1,4.The probability Pr(Z |Z 1 , …, Z m ) in equation ( 6) can be written in terms of N-steps transition probabilities as follows: ) .
,( ) ,( Pr( , , ) Pr( , , ) where refers to the element (k, l m ) of the matrix P m powered to h m .Equation ( 7) is normalized to eliminate the constant C: For more details on the mathematical basis of this model see ([4] and [2]) Now, it is possible to simulate a two-dimensional Markov chain over a regular rectangular grid of points that represents a two-dimensional stochastic permeability distribution.To illustrate the algorithm of simulating such a Markov chain, consider the domain described in Fig. 2. We suppose that the permeability is known in the dark brown cells from some experimental measurements.It is called the hard data.The cells of the grid are scanned sequentially form up to down and alternatively from right to left and from left to right.In the case of this particular domain, each cell Z to be simulated has three known neighbors: Z 1 and Z 2 which were simulated before (light brown) and Z 3 which is a hard data.Using Equation ( 8) with m = 3, one can have the probability law of Z given its neighbors and then simulate it (see [5]).The three dimensional case can be easily derived from the preceding equations.

BAYESIAN INFERENCE
In this section, the mathematical expression of the probability density function (pdf) of Pr(K ~|κ, π) is derived.This pdf represents the constraints on the permeability distribution by taking into account the permeability and the pressure measures.If a particular realization of the random permeability distribution K ~has a high probability density then its corresponding simulated pressure values are close to the measured ones.
The error ε in the pressure measures π is considered to be an independent white noise with variance s: where p ~is the numerical solution of the flow problem (1) obtained using the finite volume method.From Bayes's formula we have: (10) The joint event (K ~, κ) is identical to the single event K ~because K ~is an interpolation of κ using the tow-dimensional Markov chain so all the values of κ are in K ~: (11) Equation ( 11) is the basis of the Bayesian inference (see [7]).It involves three different probability densities: • pdf (K ~), the prior probability density function, includes all the information on the permeability distribution K ~without considering the pressure measures.In our case, it is equal to one if K ~is a possible realization of a Markov chain interpolating κ and zero elsewhere.• pdf (π |K ~), the likelihood probability density function, represents the mathematical model of the underground water flow which maps the permeability K ~into experimental pressure measures.The density of this probability can be deduced from the one of the noise in the pressure measures: Hence, the density function of this probability is deduced from the one of ε: • Pr (K ~|π ), the posterior probability density function, is the probability density function of our interest.It represents our knowledge about the permeability K ~after making the observations on the pressure.From above, we have: The expression of this probability density shows the link it has with the least squares inverse problem discussed above: the K ~that maximizes pdf (K ~|κ, π ) minimizes ||π − p ~(K ~, ξ)|| 2 so it is a solution of the least squares problem.Note that each evaluation of pdf(K ~|κ, π ) require solving numerically the Poisson problem in (1).

THE UNCERTAINTY IN THE PERMEABILITY
We recall that K(x) is the permeability over the domain Ω.K ~is the discretized permeability over the regular rectangular grid of points discretising Ω, so K ~is a two-dimensional array.Each element of K ~is considered to be stochastic but it can take n different values.S is the set of these possible values.First, the inverse problem of estimating such permeability is discussed.The data of the inverse problem are: κ, a vector of measured values of K ~at some points ξ i of the domain and π a vector of pressure measures at the same points ξ i .K ~is the stochastic permeability to be simulated conditioned by (κ, π).The inverse problem can be defined as finding the distribution of each value S l of the permeability over the domain.
A classical formulation of this inverse problem would be the least squares formulation: (13) where p ~(K ~, ξ) is the numerical solution of problem (1) at points ξ i .In order to guarantee the uniqueness of the solution of this optimization problem, one must impose some regularity conditions on K ~.In [1], for example, K ~is sought in BV, some bounded value space defined by: ( for a fixed lower bound (K ~min ) and upper bound (K ~max ) of K ~.Such formulation gives a smoothed function K ~for the permeability which can not include jumps in the values of K ~.
Thus, the least squares formulation does not admits a unique solution for the case where K can only take a finite number of values.Let K be the set of solutions of the problem (13) and (15) be the set of all the permeabilities K ~that are solution of problem (13) and have K ~ij = S l in common.Given that K is not a singleton, this explains our motivation to simulate K ~∈ K and not estimate it from (κ, π).
The probability represents the relative occurrence of the value S l at in cell (i, j) between all the potential solutions in K.This probability can be written as follows: (16 where 1 is the indicator function and E is the expectation operator.A numerical approximation of this probability, using a Monte Carlo method, would be: (18) with (K ~m, m = 1 … M ) a sequence of M samples of K ~following the probability Pr(K ~|κ, π) (see [9]). and This probability represents the uncertainty in the permeability.If , for some value S l , this means that all the potential solutions have K ij = S l in common.One can deduce that in this case there is no uncertainty in the value of the permeability at the cell (i, j), the value at this cell is S l .If for some value S l , this means that κ and π do not influence the value of K ~ij and it is not possible to deduce any information about the value of K ~ij .The computation of this probability is requires sampling K ~from the probability law Pr(K ~|κ, π).We propose to do this by coupling the MCMC sampling algorithm with the multidimensional Markov chain model as it is described in the next section.

THE COUPLING ALGORITHM
The Markov Chain Monte Carlo (MCMC) algorithm is a general method for generating random variables that follow a defined probability density (see [6]).It is used to generate the sequence of random realizations of K ~following pdf (K ~|κ, π), the probability density that conditions K ~by the permeability and the pressure measures (see section 3).A general description of MCMC is summarized in the following steps: 1.
At iteration i, generate a candidate C using q(C |K ~i).

2.
Generate a random number u following a uniform distribution over [0,1].If Under some conditions, the distribution of the sequence K ~i converges to pdf(K ~|κ, π), and then we have K ~i ∈ K. Using this sequence the probability Pr(K ~ij = S l |κ, π) representing the uncertainty in the permeability is computed.
One can see that in step one of MCMC, a candidate C has to be generated from the permeability K ~i using a proposal distribution q(.|K ~i ).This proposal distribution defines how the candidate C is generated.It has to insure that the candidate C does not differ a lot from K ~i.A few number of cells in C are allowed to take different values from their corresponding cells in K ~i.The difficulty of implementing the MCMC algorithm is the definition of a suitable proposal distribution q(C|K ~i).The candidate C generated by q must be a Markov chain having the same transition matrices as K ~and interpolating the same values κ.To generate such a candidate, we propose to couple the two-dimensional permeabilities C with K ~i by considering that both of them constitute a unique threedimensional Markov chain (see Fig. 3).The values of C and K ~i are coupled such that C ij is not only correlated with its neighbors in the four cardinal directions but also with a fifth neighbor which is K ~ij .A transition matrix in the third dimension of this multidimensional Markov chain must be defined.It can be in the form: One can see that if α = 1 the two permeability fields C and K ~i are identical, and if the two permeability fields are completely decoupled.For some value of , C and K ~i differ by a fixed number of pixel, in average.The choice of α is tuned manually by checking the efficiency of the MCMC sampling algorithm.Some methods to study the efficiency of the sampling algorithm are proposed in ( [6]).

NUMERICAL EXPERIMENT
We recall that the Poisson problem defined in (1) is considered.The domain Ω is a square discretized by a grid of 50 × 50 nodes.The source term, g, is taken equal to zero and the boundary condition, h, is equal to 10 on the left and lower sides of Ω and zero on the other sides.The discretized permeability K ~takes its values in {1,10} and its modeled using a twodimensional Markov chain having its four transition matrices int the four cardinal directions equal to: (20) with β = 0.80.The set of measurements points, ξ, constitute a grid of n s × n s points placed at the center of Ω.Finally, the transition matrix in the third direction used in the coupling algorithm to generate the candidates with the MCMC is: (21) with α = 0.9955.
In the following numerical experiments, the multidimensional Markov chain model is compared with the the multidimensional Markov chain model coupled with MCMC.The comparison is based on computing Pr(K ~ij = S 2 |κ) with the first one (using only permeability measures) and Pr(K ~ij = S 2 |κ,π) with the second one (using permeability and pressure measures) for all (i, j).As explained in section [4], these probabilities represent the uncertainty at cell (i, j) of having the permeability equals to S 2 when it is conditioned by the pressure and permeability measures.One would expect to have the second probability higher then the first one because the additional information in the pressure measures must constrain more the permeability.The uncertainty at cell (i, j) of having the permeability equals to S 1 can be deduced from Pr(K ~ij = S 1 |κ,π) = 1 − Pr(K ~ij = S 2 |κ,π).
Both algorithms are run with the same measure grid points (n s = 5, 7 and 9). Figure ( 4) shows the permeability field used to generate the "experimental" permeability and pressure measures (κ,π) and it represents the exact solution of the inverse problem. Figure (5) shows The permeability distribution to be estimated.
the unnormalized posterior probability density (Pr(K ~i| π), i = 1…M) of the samples generated by MCMC.One can notice that these density values are not centered around one fixed mean value.This shows the existence of multiple local maximum in the posterior probability density and thus multiple solutions of the inverse problem.Figures (10, 8 and 6) correspond to the values of Pr(K ~ij = S 2 |κ) using the permeability measures (the multidimensional Markov chain) and figures (11, 9 and 7) correspond to the values of Pr(K ~ij = S 2 |κ,π) using the permeability and pressure measures (the coupling algorithm).We recall that if then there is no information about the value of the permeability at the cell (i, j ) (the uncertainty is maximal) and if Pr(K ~ij = S 2 |κ) = 1 then it is sure that the cell (i, j) has the value S 2 (no uncertainty).The first remark is that, with both algorithms, these probabilities take higher values when the number of measures increases.This means that the uncertainty in the estimated permeability is directly affected by the number of measures.
The second remark is that, for the same number of measure points, the probability takes higher values with the coupling algorithm using permeability and pressure measures (figures in the right column).This shows that there is less uncertainty in the permeability when it is conditioned by the pressure measures.In figures ( 6) and ( 7) for example, it is shown how the pressure measures improve our knowledge on the permeability distribution.The regions of probability less than 0.8 in the first figure become of probability 1 in the second figure.These       regions are in good agreement with Figure (4), the exact permeability distribution.Thus, the coupling algorithm is able to generate efficiently stochastic permeabilities by taking into account information from permeability and pressure measure.One can also deduce from figure (7) that the region of probability one is the distribution of the value S 2 over the domain.
As a concluding remark, the coupling algorithm generates a sequence of stochastic permeabilities having less variability in its random values when pressure measures are used.These generated stochastic permeabilities have more cells with same value in common.This shows that the coupling algorithm is able to give efficiently the inference of the pressure measure on the permeability allowing to have a better idea on the distribution of the values of the permeability over the domain.

CONCLUSION
In this paper, an algorithm is presented for simulating stochastic permeability conditioned by some permeability and pressure measures.The permeability is considered to be stochastic but can take a finite number of values.The simulation algorithm is based on coupling the multidimensional Markov chain model with the MCMC sampling technique.Some numerical experiments were carried out with binary valued permeability.They showed a good accuracy in the results for different numbers of permeability and pressure measures.The proposed coupling algorithm generates stochastic permeabilities with less uncertainty then the classical multi-dimensional Markov Chain one meaning that it is able to efficiently constrain the permeability with the pressure measures.
This work is funded by the German federal ministry of education and research.

Figure 1
Figure 1 One dimensional Markov chain for simulating permeability distribution over a borehole.

Figure 2
Figure2Simulating the permeability of a cell (Z ) constrained by its nearest known simulated neighbors (Z 1 , Z 2 ) and and nearest had data (Z 3 ).

nFigure 3
Figure 3 Coupling the two-dimensional permeability fields C and K ~i by considering that both of them constitute a unique three-dimensional Markov chain.

Figure 5
Figure5The sequence of the MCMC samples and their corresponding posterior probability in the case of 9 × 9 measures.

Figure 6
Figure 6 Probability distribution over the domain of having the permeability equals the second value, without conditioning by the pressure, in the case of 9 × 9 measures.

Figure 7
Figure 7 Probability distribution over the domain of having the permeability equals the second value, conditioned by the pressure, in the case of 9 × 9 measures.

Figure 8
Figure 8 Probability distribution over the domain of having the permeability equals the second value, without conditioning by the pressure, in the case of 7 × 7 measures.

Figure 9
Figure 9 Probability distribution over the domain of having the permeability equals the second value, conditioned by the pressure, in the case of 7 × 7 measures.

Figure 10
Figure 10 Probability distribution over the domain of having the permeability equals the second value, without conditioning by the pressure, in the case of 5 × 5 measures.

Figure 11
Figure 11 Probability distribution over the domain of having the permeability equals the second value, conditioned by the pressure, in the case of 5 × 5 measures.