Multi-point Reservoir Modelling

Field of the invention

[0001] Embodiments of the present invention relate generally to methods for building statistic models of geological reservoirs, to an improved method for storing and building such models. Related Art

[0002] Reservoir simulation is a technique used by major oil and gas companies to predict oil and gas flow through the porous media of the reservoir. Such simulations typically use 2D or 3D models of the reservoir that includes a grid of a large number, often in excess of a million, of individual cells. These models typically represent the values of one physical property, often the underlying geological facies expressed as a discrete- valued category variable. The goal of modelling is to obtain an image of the underground reservoir that reflect the knowledge of the general geological structures presents on the site, restrained by conditioning values measured at certain points, for example corresponding to exploration and production wells, and by other auxiliary information, like seismic data. [0003] Oil reservoirs present several different morphologies, possess a large amount of heterogeneity, and typically span a large extension, both in surface and in volume. Conditioning data derived from assay wells, on the other hand are almost always quite sparse and other data, like seismic data, have, in general, a relatively coarse resolution. In order to make the most of available data, reservoir modelling must necessary use probabilistic, or geostatistical methods. [0004] Multiple-point statistics is a known technique allowing the simulation of heterogeneous geological facies images. The simulated images (either in 2D or in 3D) respect the structures and morphology of a given reference training image and adhere to the conditioning data. This technique has been described, among others in the following documents, which are incorporated hereby by reference: "Sequential Simulation Drawing Structures from Training Images"

(2000) by Strebelle S. and Journel A. in Kleingeld and Krige (eds.), Geostatistics 2000 Cape Town, Vol. 1;

"Multivariate geostatistics: beyond Bivariate moments"(1993) by Guardiano F. and Srivastava M. in Soares A. (ed.) Geostatistics Troia'92 vol. 1 Kluwer Acad. Publ., Dordrecht. "Conditional simulation of complex Geological structures using multiple-point statistics", Strebelle S, Mathematical Geology, 35,

No 8, pp. 915-925.

[0005] Following the algorithm proposed by Strebelle, the training image is scanned with a search template and the multiple-point statistics are stored in a dynamically allocated tree structure. Further on, the statistics so gathered are used to fill a simulated model.

[0006] Document WO2006/023602 teaches a method for creating a reservoir facies model based on multiple-point statistics. FR2905181 also describes a similar method that further includes steps to take non- stationarity into account. [0007] A common drawback of the known techniques is that the amount of memory required for storing the tree structure holding the multi-point statistics can quickly become prohibitive as the size of the model increases. Particularly for 3D images, the need to keep the size of the tree within manageable limits imposes the use of a small search template, with the result that complex structures are much worse represented in 3D than in 2D. Brief summary of the invention

[0008] It is an aim of the present invention to propose a method to create a reservoir model that is free from the aforementioned limitation and, in particular, to propose a method that can be scaled to models of larger size

[0009] Another aim of the invention is to propose a method to create a 3D reservoir model that allows a better representation of complex structures. [0010] According to the invention such objects are obtained by the method that is the object of the main independent claim and, in particular, Short description of the drawings

[0011] The invention will be better understood with the aid of the description of an embodiment given by way of example and illustrated by the figures, in which:

Figure 1 and 2 show an example of a training image and of a search template.

Figure 3 is a graphical representation of a tree structure issued from the training image and the template of figures 1 and 2. Figure 4 is a graphical comparison of the memory usage in the method of the invention and in a known method.

Figure 5 illustrates an aspect of the method of the invention.

Figure 6 shows possible weight distribution for selection and penalization in presence of an auxiliary variable representing non- stationarity

Detailed Description of possible embodiments of the Invention [0012] The basic principles of multipoint statistics for simulating a facies at each cell of a grid are now recalled summarily. We assume that one disposes of a Training Image (Tl) representing a model geologic structure as a 2D or 3D vector. The present invention, however, is not limited to two- dimensional or three-dimensional problems and could be applied also to systems having a different dimensionality. It is not essential for the invention that the model space is defined in term of spatial coordinates either. The model space could in fact also include a temporal coordinate or any kind of other variable. The method of the invention could also be advantageously applied to the modelling of non-geologic systems. [0013] A search template is defined as a set of relative nodes locations h _v ...h _N . For a given reference node u, the search template at u selects the sets of nodes τ(u) = {u + h _u ...u + h _N } and, if s(v) denotes the facies at node v , the vector d(u) = {s(u + h,),...s(u + h _N )} defines the data event at u. Note that one can consider also a data event with undefined components (for a node that has not yet been simulated). [0014] To attribute a facies at a node u in the simulated image, we consider the nodes v in the Training Image where the data event c/(v ) has the same components as those of d{u) . Then, the number of occurrences of all the facies at the nodes v is counted. This provides a probability distribution function (pdf) that is used to draw a facies at the node u

randomly. More precisely, if the positions of the known components in d(u) are i, ≤ ... ≤ i _m (with O ≤ m < N) ₁ then the probability to draw the facies k at the node u is forO < k < M , where M is the number of facies. Note that if the number of matching data events (the denominator of the formula above) is too small, (for example less than a given minimal number of replicates), we consider the last defined d{u) component as undefined, and this operation is repeated until the number is acceptable. [0015] Preferably, the multigrid approach, as proposed in the above- identified art is used to capture structures within the training image that are defined at different scales. Let us introduce some terminology. A 2D or

3D grid is a box-shaped set of pixels

G = {(n _x ,n _y (,n _z )) : 0 ≤ n _x ≤ N _x - 1,0 ≤ n _y ≤ N _y - %0 ≤ n _z ≤ N _z - D}, (2) [0016] where N _x , λ/ _y , (λ/ _z ) are the dimensions in the x-axis, y-axis (and z-axis) directions respectively. In the main grid G, the subgrid with a step of 2' between neighbouring nodes is defined as

SG ₁ = {(n _x ,n _y (,n _z )) : (n _x ,n _y (,n _z )) e G and n _x ,n _y (,n _z ) = 0 mod 2'}, (3) where / is a positive or zero integer. [0017] For a given situation, m levels of multigrid G ₀ ,...,G _m _ _λ are defined as

G _m _, = SG _m _i and G, = SG ₁ \ SG ₁₊ , for 0 ≤ / ≤ m - 2. (4)

They form a partition of the main grid G.

[0018] The simulation proceeds successively with the simulation of all the nodes in the multigrid level G _m _ _λ . The lag vectors in the search template are then scaled such that their value is 2 ^m~1 -/i, instead of /i, . The process continues with the nodes in the multigrid level G _m _ ₂ and repeats for all the multigrid levels (in decreasing order) similarly. The simulation starts with the coarsest multigrid level and finishes with the finest one. [0019] According to one aspect of the present invention, the multiple- point statistics inferred from the training image and provided by a search template are stored in a dynamically allocated list. An element of the list is

a pair of vector (d,c) where d = (S ₁ , ...,s _w ) defines a data event and c = (c ₀ ,...,C^ ₁ ) is a list of counters summing up the number of occurrences for each facies: c, is the number of data events d{v) equal to d found in the training image with facies / at the reference node v. Note that the total number of elements in the list is the number of different data events in the training image.

[0020] The list is sorted such that the identities of the data events (vector d) are in increasing order according to the following relation: c/ ⁽¹⁾ = (sj ¹⁾ ,...,s< _/ ¹⁾ ) < c/ ^<2) = (sj ²⁾ ,...,s<, ²⁾ ) if and only if there exists y e {1,...,λ/} such that s, ^<1) = s, ^<2) for 1 ≤ / < y andsj ¹⁾ < sj ²⁾ .

[0021] To build the list, the training image is scanned moving the search template such that:

1) either the search template is entirely included in the training image (scan in a window of the Tl) ₁ 2) or the reference node of the search template covers all the nodes of the training image (whole scan of the Tl); in this case, a new value for "lacking" facies should be introduced. [0022] To simulate a facies in a node u, for each facies k (0 ≤ k < M) we compute the sum Q of the occurrence counter c _k of each compatible element in the list with the data event d(u) . Then, the conditional probability distribution function knowing the data event d(u) , used to draw the facies in the node u, is given by

ψ(s(u) = k \ d(u)) = ^, 0 ≤ k < M, (5)

where C = ^ ¹ Q . [0023] In practice, a given value C _mm represents the minimal number of replicates that must be found in the training image for validating the corresponding pdf, i.e. the pdf above is validated only if C ≥ C _mm . [0024] Assume that the positions of the known components (corresponding to the nodes that has been already simulated) in d(u) are /., < ... < /„ (with O ≤ n < N). In this case, we compute the sum C _k ^J) of the occurrence counters c _k of each compatible element in the list with the data event d ^{i) (u) whose defined components are s(u + h _h ),..., s(u + h _: ) (only the first j simulated nodes in the search template are taking into account), for

O ≤ / ≤ n. Then, we choose the greatest index/ such that C ^ω = ∑y _o C _k ^{J) ≥ ^c _m , _n and the corresponding pdf

P(s(u) = k \ d(u)) = 0 ≤ k < M ₁ to draw the facies in the node u. [0025] Note that for the simulation of each multigrid, the Training Image is scanned once for building the corresponding list. Then this list is used for simulating the nodes of the considered multigrid. [0026] In order to illustrate the method of the invention, a full example will be carried out with reference to the Figures 1 and 2. This example is extremely simplified, and it must be understood that the invention is applicable to much larger 2D or 3D reference images and search templates, and to variables having more than two discrete values. Yet it describes the essential steps of the inventive method. [0027] Figure 1 shows a training image, in this case a 2D 6x 6 grid containing values of a binary variable, assuming the values 0 or 1, represented as white and black squares. A possible corresponding search template is shown in figure 2, containing four cells in the neighbourhood of the considered cell u. [0028] In conventional methods, a tree structure is used instead of a list for storing the multiple-points statistics inferred from the training image. The tree has a depth of N+1 (where N is the size of the search template, N=4 in the presented example). It is made up of cells divided in M (number of facies, M=2 in our simple example) subcells. The levels of the tree are numbered from 1 to N+1 and the subcells in a cell from 0 to M-1. Each subcell holds a counter and can have a child cell (in the next level of the tree): the tree is an M-ary tree. The counter in a subcell is defined as follows. Let /(1), /(2), ., i{k + 1) a path in the tree where i(j ) is the identification number of a subcell in a cell of level j. The counter in the last subcell of the path is the number of data events found in the training image with facies i(j ) at node v + h _j of the search template τ (v) for

1 < j ≤ k , and with facies i(k+l) at the reference node v. In each cell, there is at least one subcell with a strictly positive counter, i.e. empty cells are not stored.

[0029] If the training image is scanned such that the search template is always entirely included in the training image for building the tree, the obtained tree for the training image of figure 1 and the search template of figure 2 is presented in figure 3. [0030] According to the present invention, instead of using the complete search tree, the multiple-point statistics inferred from the training image is stored in a list structure as explained above. Only the data events corresponding to the last level in the search tree of figure 3 are stored in the list of table 1. The statistics stored in the list and the search tree are identical. The knowledge of the list allows reconstructing exactly the search tree.

Elements of the list Corresponding nodes

Li = (0,0,1,1), (0,1) (0,{(2,2)}}

L ₂ = (0,1,0,1), (0,2) {0,{(3,1),(3,4)}}

L ₃ = (0,1,1,0), (0,2) (0,{(1,4),(4,2)}}

L ₄ = (0,1,1,1), (1,0) {{(3,2)},0}

L ₅ = (1,0,0,0), (1,0) {{(3,3)},0}

L ₆ = (1,0,0,1), (0,2) (0,{(4,1),(4,4)}}

L ₇ = (1,0,1,0), (1,1) {{(4,3)},{(1,3)}}

L ₈ = (1,0,1,1), (1,0) {{(2,3)},0}

L ₉ = (1,1,0,1), (0,2) {0,{(1,2),(2,4)}}

Lio = (1,1,1,0), (1,1) {{(1,1)},{(2,1)}}

Table 1: list for training image and template of figures 1 and 2

[0031] The size of the search tree and the list depends on the size λ/ of the search template, the number M of facies and the entropy of the training image. Moreover, the size of the search tree can also depend on the order of the nodes numbering in the search template.

[0032] The memory requirement of the method of the invention and of the known tree method can be compared either by combinatorial considerations or by numerical simulation. The latter method has shown that, in typical cases of realistic size, the memory occupation of the list structure of the invention is lower of a factor 10-30 than the size of the tree used in known methods. Figure 4 shows the relative memory gain simulated for a 3D training image having 600'0OO nodes and 4 facies (curve 125), and for a 3D training image having 14'128'385 nodes and 6 facies (curve 120).

[0033] One possible method of simulating a model with the treatment of conditioning data according to the present invention is now described summarily. It must be kept in mind that this particular algorithm for the treatment of the conditioning data is provided simply as an example to better understand the invention, but is not in itself an essential feature, and the treatment of conditioning data could be, in the frame of the present invention also performed otherwise.

[0034] Assuming, for example, that the simulation is carried out on a 2- dimensional or 3-dimensional space, a node location (n _{χ l} n _y (,n _z )) in the simulation grid is the pixel [n _{χ l} n _x + 1[x[n _y ,n _y + 1[(x[n _z ,n _z + 1[) (the notation [...[ indicates a inferiorly closed and superiorly open real interval) and corresponds to the region

[o _x + n _x - d _{χ l} o _x + (n _x + \)d _x [ x [o _y + /yc/ _y ,o _y + (n _y + 1)c/ _y [ (6)

(x [o ₂ + n _z -d _z ,o _z + (n _z + 1)c/ _z [ ) of the whole area [o _x ,o _x + N _x -d _x [x[o _y , o _y + λ/ _y -d _y [ (x[o _z ,o _z + N _z -d _z [), (7) where (o _{χ l} o _y (,o _z )) is the origin

, d _{χ l} d _y (,d _z ) are the node spacing in each direction, and N _x x /V _y (x/V _z ) are the dimensions of the simulation grid. A conditioning data is a point (x,y(,z)) with an attributed facies. The coordinates x,y(,z) are real and must be in the area defined by (7).

[0035] Assume that m levels of multigrid are used for the simulation, in order to account for multiple-scale structures. Before starting the simulation, the conditioning data ought to be spread into all subgrids SG ₀ ,...,SG _m _i, in order to take them in account from the beginning of the simulation. This is done as follows.

1) For each conditioning data, the attributed facies is assigned to the node in simulation grid whose corresponding region (6) contains the point of the data. (If more than one conditioning data lead to the same node location, the conditioning data whose point is the closest to the centre of the corresponding region is retained.)

2) Each conditioning node u _c in G obtained by the step above is spread in the subgrid SG^ as follows: a. We select all the closest nodes in the subgrid SG^ to u _c , i.e. the realize the minimum where d is the dimension (2 or 3) and u(j) (resp. u _c (j) are the integer coordinates in G of the node u (resp. u _c ). b. If all the selected nodes are not yet simulated, we choose one randomly and simulate a facies for the chosen node using the search template τ ₀ , corresponding to the scale of the subgrid SG ₀ , otherwise nothing is done (no need to spread this data).

3) The previous step is repeated for the sub-grids SG ₂ ,...,SG _m _^ successively. The conditioning nodes in G obtained in the first step are spread in the subgrid SG ₁ using the search template τ,^ corresponding to the scale of the sub-grid SG ₁ ^ . [0036] The simulation continues with the simulation of all the un- simulated nodes in the multigrid level G _m _ _v G _m _ ₂ ...,G ₀ successively. [0037] Figure 5 shows, by way of example, a 2D grid G defined by the origin (o _{χ l} o _y ) = (0.0,0.0) , node spacing d _x = d _y = 1 and dimensions N _x = 15 and λ/ _y = 10. Assume that 3 levels of multigrid are used and that there are four conditioning data whose location are O ₁ = (4.4,6.4), Q ₂ = (8.5,8.8) , O ₃ = (9.2,4.7) and O ₄ = (11.4,7.3) . These points are marked as black dots 130 in figure 5. The conditioning nodes obtained in step 1) above are u _{c 1} = (4,6), u _{c 2} = (8,8), u _{c 3} = (9,4) and u _{c 1} = (11,7) , and marked as white squares 140 in figure 5. The nodes selected in step 2)a) above are marked as white circles 141 in figure 5 for the subgrid SG^ ( (4,6) for u _{c 1} , (8,8) for u _{c 2} , (8,4) and (10,4) for u _{c 3} , and (10,6), (10,8), (12,6) and (12,8) for u _cA ) and marked as crosses 142 in figure 5 for the subgrid SG ₂ ((4,4) and (4,8) for u _{c l} , (8,8) for u _{c 2} , (8,4) for u _{c 3} , and (12,8) for u _{c λ} ).

[0038] In step 2)b) above, one of the two nodes (8,4), (10,4) and one of the four nodes (10,6), (10,8), (12,6), (12,8) will be simulated for the subgrid SG^ . For the subgrid SG ₂ , one of the two nodes (4,4), (4,8) will be

simulated, and the nodes (8,4) and (12,8) will be also simulated if they are not chosen during the spreading of the data in SG^ . [0039] For each multigrid level, a path covering all the un-simulated nodes should be chosen. A possible strategy is to build, in a given multigrid level, a path from the simulated nodes in this multigrid level (provided by conditioning data and the steps described above). For each of these nodes, we consider the first layer of nodes (in the considered multigrid level) that just encompasses it, then the second layer of nodes that just encompasses the first one, and so on. A layer of nodes is simply made up of nodes in the border (edges of faces) of a box. The proposed path visits all the nodes in the first layers, then all of them in the second layers, and so on. Thus, such a path moves away from the conditioning data.

[0040] The examples presented above are valid for stationary training images. It is possible, however, to extend the method of the invention, by known techniques, to account for auxiliary continuous variables describing non-stationarity. In the following, a preferred method for dealing with non-stationary training images is presented in some detail. Other methods are however possible and included in the scope of the invention. [0041] In this case, a training image 77 containing facies (primary variable s) and an auxiliary training image Tl _aux with a continuous value attributed in each node are used. These two training images have the same dimensions and the continuous variable t of the second one describes the non-stationarity of the first one. For simulating a grid G, we need also an auxiliary grid G _aux of same dimensions (as G) with the same continuous variable as in the auxiliary training image.

[0042] In presence of an auxiliary variable, a vector m is appended to each element of the list described above. An element of the list is then a triplet of vectors (d,c,m) where c/ = (s _1f ...,s _w ) defines a data event, c = (c _0l ...,c _M _ _λ ) is a list of occurrence counters for each facies and /T7 = (m _o ,...,/T7 _M _ ₁ ) is a list of mean values of the auxiliary variable in each case: c, is the number of data events d{v) equal to d found in the training image with facies / at the reference node v and m, is the mean of the auxiliary variable at these nodes v. [0043] Before starting the simulation and building the list, the auxiliary variable t is normalized in the interval [0,1], via the linear transformation

[a,b] -> [0,1], t \→ (t - a) l (b - a), where a and b are respectively the minimum and the maximum of the auxiliary variable t(v), v in 77 _aux u G _aux . [0044] To simulate a facies in a node u of the grid G ₁ knowing the data event d(u) and the auxiliary variable t(u) (provided from the auxiliary grid G _aux ), we proceed as follows:

[0045] A tolerance error δue [0,1] is fixed. For each facies k, we retain the set E of the elements in the list that are compatible with the data event d(u) and that satisfy

\ m _k - t(u) \≤ δu. (8) [0046] Then, we compute the sum C _k of the occurrence counter c _k of the elements of E and the resulting mean M _k for the auxiliary variable:

Q =£c< ^e) and M, =-L£, (e) ml

[0047] where the exponent (e) denotes the corresponding element in the list. The resulting means are used to penalize the corresponding counters: for each facies k, a weight y _k is defined as where ω is a positive real parameter (weight) controlling the penalization. [0048] Then, the conditional pdf knowing the data event d(u) and the auxiliary variable t{u), used to draw the facies in the node u, is given by P(s(u) = /r | d(i/) _f t(i/)) = ^-, 0 < /r < M _f (10) γ where J = ∑ ^M _κ j' _k .

[0049] In summary, the presence of an auxiliary variable adds a step of selection (8) and a step of penalization (9) for retrieving the conditional pdf. The weight function y _k are illustrated, for three values of ω, on Fig. 6. [0050] In practice, we check if the sum C = ∑ ^{™ 1} Q is greater or equal to a given value C _mιn . We reduce the data event d(u) as much as needed so that it is the case, as described previously. Then, the pdf is computed as described above. Note that if the minimal number of replicates C _mιn is not reached with the empty data event centred at u, we set the tolerance δu to 1 temporarily (i.e. no selection by the auxiliary variable) and we do not take in account the data event: we use the pdf given by

(marg)

F(s(u) = k \ t(u)) = íj^y, O ≤ k < M, (1 1)

where

M-1 yjmarg) ₌ N _ | _M <™φ> _ _t(u) | ^ω | _C < _m aφ) _{a nd} y(marg) ₌ £γ<™φ> _{(1 2} )

/t=0 with c| ^maλ9) the number of nodes in the training image with the facies k and M _k ^mar9) the mean of the auxiliary variable at these nodes. [0051] Although the above example deals with one single auxiliary variable, the invention can be extended to more complex cases requiring several auxiliary variables. For example, in presence of n auxiliary variables, the definition (9) of y _k is replaced by where t ^ω (u) denotes the y-th auxiliary variable at u, M _k ^J) the resulting mean for the y-th auxiliary variable, and δu ^{J) and ω ^<y) the tolerance and the weight respectively for the y-th auxiliary variable. Similarly, j _k ^mar9) (see (12)) is replaced by γ [ ^marg) = J! f 1 - - t ^{J) (u)[ ^{ω 1} )c _k ^{mar9) ,

where M _k ^J){marg) is the mean of the y-th auxiliary variable at the nodes in the training image with the facies k.

[0052] The multiple-point statistics provide a pdf used for attributing a facies at a given node u. To compute this pdf, the elements of the list matching with the data event d(u) are retained, as previously explained, for simulation, eventually also considering auxiliary variables expressing non-stationarity. The sum C of all the occurrence counters of these elements is the number of replicates found in the training image that are compatible with the node u to be simulated. For validating the pdf provided by the multiple-point statistics, this number of replicates should be large enough, i.e. C > C _mιn . If it is not the case, the data event is reduced (dropping its last simulated nodes) as much as needed for obtaining this condition, i.e. the complete data event d(u) is incompatible with the reference training image.

[0053] Preferably, the method of the invention comprises steps for eliminating or reducing such incompatibilities. In a preferred variant, a method called synchronized post-processing is adopted. It is a real-time process, for avoiding the propagation of inconsistencies during the simulation.

[0054] According to this variant of the invention, The synchronized postprocessing algorithm consists in the following steps at each node u of the simulation grid:

1) Compute the pdf provided by the multiple-point statistics for the node u as above, and check if this pdf is acceptable. The criterion of acceptance is based on the number of dropped nodes in the data event d(u) : the pdf is considered as acceptable if the ratio of the number of nodes in d{u) considered for retrieving the pdf over the total number of simulated nodes in d{u) is greater or equal to a given threshold (in [0,1 ]).

2) If the pdf is unacceptable the simulation of the facies at node u is postponed, and one of the simulated nodes in d(u) , N ^~ \u), is un-simulated (i.e. its attributed facies is removed) and simulated again (re-simulation). For each re-simulation, the criterion of acceptance is considered. For the simulation at the node λT ^D (u) : a. if the multiple-point pdf is unacceptable for N ^~D (u) , the simulation of the facies at the node N ^~D (u) u is postponed, and one of the simulated nodes in d(N ^~D (u)) , N ^~D~ \u), is un-simulated and re-simulated; b. if the multiple-point pdf is acceptable for N ^~D (u) , the facies at this node is simulated, and a new attempt at the node N ^~D+ \u) is done. [0055] Note that D is the number of currently un-simulated (removed) nodes from the initial node u , called depth. To ensure the convergence of this recursive algorithm, a maximum total number of un-simulations for each initial node is preferably fixed. A maximal depth can also be fixed. If one of these maximal values are reached, the simulation continues using the current ((un)acceptable) multiple-point pdf for the current step.

[0056] In this process, the node that will be un-simulated could be chosen randomly, but for saving CPU time the following strategy is preferably adopted. Assume that τ(u) = {u + h _v ... ,u + h _N } is the search template at u , and the simulated nodes of the data event d(u) are L/ + /i, ,...,u + /i, _n . Assume that the nodes u + h, , J ≤ j ≤ n were dropped for computed the multiple-point pdf, i.e. the condition C ≥ C _mιn is satisfied considering the data event with the first 7 - 1 simulated nodes (d ^υ~ ^(u)) but not satisfied if the next simulated node (u + h _: ) is taken in account in the data event. In this case, the node u + h _h is un-simulated and re- simulated.

[0057] Preferably the method of the invention contains steps to ensure that some particular nodes not be un-simulated. It is the case, for example, of the conditioning nodes or the nodes that are not included in the multigrid currently simulated. To this end, before the simulation of a multigrid, these nodes might be marked as "fixed nodes" ( for instance with a Boolean flag). Then in the case above, we re-simulate the node u + h _lk where k is the greatest index in {1,...,7} such that u + is not flagged as "fixed node" instead of the node u + h _ιr If no such index exists, the simulation might continue using the current (unacceptable) multiple- point pdf for the current step.

[0058] The synchronized post-processing can sometimes un-simulate and re-simulate the same nodes in a cyclic manner. Such a situation is a waste of time because the simulation is not improved. Moreover, these cycles generate biases: a facies could be simulated many times at a node with the same multiple-point pdf. Preferably the method of the present invention includes steps to avoid useless re-simulation cycles. For avoiding these cycles, for example, all the un-simulations (location, facies and depth) from the initial node are recorded in a log. Before un-simulating a node, the log is scanned to determine if the simulation grid is found in the same state as in the past and if the node that would be un-simulated is the same one in the both situations. If it is the case, no un-simulation is done and the simulation continues using the current (unacceptable) multiple-point pdf for the current step (or another simulated node in the data event could be un-simulated, for breaking the cycle).

[0059] The following embodiment of the invention presents a possible parallel implementation of the multiple-point simulation algorithm of the present invention that, advantageously, can reduce the simulation time. Other manner of parallelising the method of the invention are however possible.

[0060] According to this embodiment the simulation is carried out by a number p of parallel processes, numbered from 0 to p-7, that may, or may not, be executed by several physically distinct hardware processors. The list containing the multiple-point statistics is then divided in p parts of balanced size. Each part is called a partial list. Every process contains one different partial list.

[0061] More precisely, if K is the number of elements in the list and K = qφ + r is the Euclidean division of K by p (q and r are integers such that 0 ≤ r < p), then there will be /" partial lists with g + 1 elements and p - r partial lists with q elements.

[0062] Let R(O) ,..., L(K - 1)} be the entire list, sorted as previously explained. The partial lists stored on the processes are then given in table 2. [0063] For computing the pdf provided by the multiple-point statistics, each process retrieves the needed information (counters, auxiliary variable's mean) scanning their local list. Then, the retrieved information by each process are grouped together (using parallel computing communication), and the pdf is computed. So, the simulation of each node is parallelized. [0064] The building of one list remains serial: one process build a list and then distributes it to all the processes as above. If n lists must be built for the simulation, and p processes are used, and n = qφ + r is the Euclidean division of n by p, the process/ build q + 1 lists for 0 < j < r, and q lists for r ≤ j < p .

Proc. Elements of the partial list

0 L(O) Up) U(q-D-P) L(q-p)

1 UD Up + 1) U(g - D-P + 1) Uq-P + 1)

M L(r - 1) Up + r-1) ... Lϋq-λ)- _P + r-λ) Uq-P + r-1) r Ur) L(p + r) U(q-D-p + r) r+ ^* \ L{r + 1) Up + r + 1) ... U(q-D-p + r + 1)

p-1 L(p - D Up + p-D ... U(q-D-P + P-1)

Table 2: partial lists