METHOD FOR DETERMINING THE GEOMETRY OF AN AREA OF A

RESERVOIR

FIELD OF THE INVENTION

The invention relates to a method for determining the geometry of an investigated area of a reservoir, in particular an oil reservoir.

STATE OF THE ART

In the petroleum industry, numerous data acquisition activities are known for broadening the knowledge about an underground reservoir. In particular, it is known to drill one or several wells in the reservoir, and to perform a so called well test to gather data about the connected volume of the reservoir, and to assess the equivalent permeability of the reservoir or specific areas thereof.

Another important step that can be performed before beginning a production phase of the reservoir is a reservoir modelling, which consists of a first determination of the geometry of the reservoir, or of its investigated areas.

To this end, the amount of data that can be provided by the exploration wells is quite limited, and complementary techniques have been developed, including acquisition and processing of seismic data to identify and characterize the various layers of rocks in the soil and the potential geometry of the oil reservoir.

However these techniques are limited as well; for instance seismic data is insufficient to interpret whether a structural trap has been sealed.

A method is also known from document WO2013/144458 which allows precisely modelling a soil’s sedimentary facies and mineralogical composition, by simulating the flowing of particles through the soil and the resulting changes in the soil’s mineralogical composition. The modelling of the soil obtained by this method can be an additional tool to provide indications regarding the ability of the soil to trap hydrocarbons, and therefore to locate possible parts of the reservoir.

However, an inconvenience of this method is that it is difficult to determine when the current state of the model corresponds to the real state of the soil, and hence when to stop the simulation of particles flowing through the soil.

It is therefore necessary to provide an improved solution for determining the geometry of an investigated area of a reservoir. DESCRIPTION OF THE INVENTION

The present invention aims at providing an improved method for determining the geometry of an investigated area of a reservoir.

To this end, the invention proposes a computer implemented method for determining a geometry of an investigated area of a reservoir, wherein a volume of the investigated area is known, the method comprising:

receiving a geological model of the soil in which said area is located, the geological model comprising a plurality of adjacent cells each having respective spatial coordinates and respective permeability values, and two wells being simulated in said geological model,

determining a preferred flowing path between the two wells, the preferred flowing path being formed by a series of adjacent cells of the model, and determined based on the respective permeability values of the cells, evaluating a volume of the preferred flowing path and,

if the volume of the preferred flowing path is lower than the volume of the investigated area, updating the preferred flowing path until its volume is superior or equal to the volume of the investigated area.

According to an embodiment, the updating of the preferred flowing path comprises adding to the current flowing path at least one cell, each added cell being adjacent to at least one cell of the current flowing path and having a permeability value suitable for increasing the permeability of the current flowing path.

In an embodiment, the preferred flowing path is determined by implementing a fast- marching algorithm, which may be configured to maximize the permeability of the preferred flowing path.

In an embodiment, said method is performed after a well testing preliminary phase during which a volume of the investigated area and an equivalent permeability between two wells drilled in the reservoir are determined, and the two wells of the geological model correspond to said two wells of the well testing preliminary phase. The method may further comprise a step of defining a search zone of the geological model, comprising:

delimiting a three-dimensional zone in the geological model such that said zone has a volume corresponding to that of the investigated area and the two wells of the geological model are contained within said zone,

establishing two planes within the geological model, each plane extending vertically and orthogonally to a segment linking the two wells,

defining the search zone as the part of the three dimensional zone contained within the planes,

wherein the search zone defines the cells that may be included in the preferred flowing path.

In embodiments, the method may further comprise, if the volume of the preferred flowing path is superior or equal to the volume of the investigated area, steps of computing the permeability of the preferred flowing path and comparing said permeability to an equivalent permeability between the wells.

Each cell of the geological model may further comprise at least one local mineralogical composition parameter, and the method may further comprise, if the permeability of the preferred flowing path is different from the equivalent permeability between the two wells, steps of:

simulating a stochastic movement of a particle in the geological model, said particle having coordinates in said model and having an aggressiveness parameter,

modifying the local mineralogical composition parameter while taking into account at least :

o the coordinates of the particle in said model

o the aggressiveness of the particle and

o the local mineralogical composition parameter,

modifying the aggressiveness of the particle while taking into account at least the modification of the local mineralogical composition, and

computing an updated permeability of the preferred flowing path taking into account the modification of the local mineralogical composition, and comparing the updated permeability of the preferred flowing path to the equivalent permeability. The method may then terminate if the updated permeability reaches the equivalent permeability. Otherwise the method may comprise a step of detecting an error in the geological model if the updated permeability of the preferred flowing path becomes more remote from the equivalent permeability than the initial permeability of the preferred flowing path.

According to another embodiment, a computer program product is disclosed, comprising a set of instructions adapted to implement the method according to the previous description, when it is executed by a processor.

A non-transitory computer readable storage medium is also disclosed, having thereon a computer program comprising program instructions, the computer program being loadable into a processor and adapted to cause the processor to carry out, when the computer program is run by the processor, the method according to the previous description.

A device for determining the geometry of a reservoir is also disclosed, the device comprising a processor configured for implementing the method according to the previous description.

A reservoir characterization method is also disclosed, comprising:

performing a well test during which a volume of an investigated area of the reservoir and an equivalent permeability between two wells drilled in the reservoir are determined, and

performing the method according to the preceding description

In an embodiment, the well test comprises:

causing a pressure disturbance at a first one of the wells drilled in the reservoir,

measuring a radius of investigation, and

inferring from said radius of investigation a volume of the investigated area.

In an embodiment the well test may further comprise: measuring at the second of the wells drilled in the reservoir a pressure variation in response to the pressure disturbance, and

inferring from said pressure variation an equivalent permeability value between the two wells.

The present invention makes use of a geological model of the soil in which the reservoir is located, an estimate of the volume of an investigated area of the reservoir, and a value of equivalent permeability between two wells drilled in the reservoir, to infer the geometry of the investigated area. This geometry is obtained by computing a preferred flowing path between the two wells and updating said paths until the equivalent permeability between the wells is reached.

For a better precision of the method, the geological model of the soil may be updated by simulating the flowing of particles between the wells to infer a change in the local mineralogical composition of the soil, and hence a change in the permeability between the wells.

The method also allows relevantly simulating the mineralogical composition of the reservoir.

BRIEF DESCRIPTION OF THE DRAWINGS

Other characteristics and advantages of the invention will become apparent upon reading the following description. This description is purely illustrative and it should be read in conjunction with the attached Figures, which show the following:

Figures 1a to 1 d schematically show the main steps of a method according to an embodiment of the invention.

Figure 2 is a functional diagram according to one embodiment of the invention,

Figures 3a and 3b schematically show two successive states of a preferred flowing path between two wells.

Figure 4 shows an example of a device for evaluating the geometry of a reservoir according to an embodiment of the invention.

Figure 5 shows an exemplary embodiment of a geological model according to one embodiment of the invention Figures 6a to 6c show the phenomenon of changing the mineralogical composition of grid cells of a model according to an embodiment of the invention.

DETAILED DESCRIPTION OF AT LEAST ONE EMBODIMENT OF THE INVENTION

With reference to figure 2, a method for evaluating the geometry of an investigated area of a reservoir is preferably performed after a well test 100 has been operated, in order to evaluate the volume V of the investigated area and the equivalent permeability K _eq between two wells P1 and P2 drilled in the reservoir.

The well test comprises a first step 110 of causing a pressure disturbance in the reservoir at a first one P1 of the wells. The pressure disturbance can for example be the injection of a high-pressure fluid into the well.

Accordingly, the pressure disturbance propagates within the reservoir and causes measurable pressure changes in the reservoir.

The radius of investigation, or also known as transient drainage radius, well known by the skilled person, is defined as the maximum distance from the point where a pressure disturbance is caused (i.e. P1 ), at which a pressure change caused by the pressure disturbance can be measured. In other words, the radius of investigation is the furthest distance at which the effect of the pressure disturbance is detectable.

The radius of investigation r _inv is defined by the following equation:

948 fme

Where :

t is the time elapsed since the pressure disturbance,

k is the average permeability of the soil,

f is the porosity of the soil,

c is the compressibility of the soil, and

m is the viscosity of the fluid in which the pressure disturbance propagates.

In order to compute the radius of investigation of the investigated area, the pressure at the well P1 is monitored in order to detect the time at which the pressure value at the well P1 comes back to what it was before the pressure disturbance. The radius of investigation is then computed by applying the above formula to the detected time.

The volume of the investigated area is then estimated during a step 120 by approximating the area as a cylinder having a radius equal to the radius of investigation, and a height estimated from geological data.

Regarding the equivalent permeability K _eq between the wells P1 and P2, it can be determined in step 130 by monitoring the pressure change at the well P2 caused by the pressure disturbance, by using the following Darcy equation:

Where:

Ap is the diminution of pressure change between the point where the pressure disturbance is caused, for instance the well P1 , and the point where it is measured, for instance the well P2,

K _eq is the equivalent permeability between the two points,

Q is the flow rate of the fluid flowing between the two points, determined from the injected flow rate,

L is the distance between the two points, and

S is the average section of a cylinder of volume V of the investigated area and of length L.

Once the volume V of the investigated area and the equivalent permeability K _eq between the wells P1 and P2 have been determined, the method 200 of determining the geometry of the investigated area is performed.

With reference to figure 4, this method is preferably implemented by a device 500 for determining the geometry of an investigated area of a reservoir. In an embodiment, this device comprises a central processing unit 503 which comprises a processor 501 adapted to execute instructions stored in a memory 502 for implementing this method, the memory 502 also storing the geological model. The device 500 can further comprise a human-machine interface 504, comprising for instance a display for displaying information regarding the current state of the geological model and the preferred flowing path, and input means for an operator to enter instructions. As schematically shown on figure 5, the method 200 of determining the geometry of the investigated area is performed based on a geological model of the soil in which the reservoir is located, said model being loaded by the processing unit at the beginning of the method.

The geological model is a three-dimensional grid comprising a number of cells, each cell C being associated with respective geographical coordinates (x,y,z) and a respective permeability value K _xyz .

Preferably, each cell of the model is further associated with a mineralogical composition parameter (designated as CM _xyz ). The variables x,y and z indicate the spatial positions of the cells. The term“mineralogical composition parameter” (or simply“mineralogical composition”) is understood to mean a parameter representing the ratios of minerals (such as sulfur, galena, cassiterite, fluoride, calcite, colemanite, chalcanthite, legrandite, etc.) within a portion of the soil, but also the ratios of other non-mineral chemical compounds such as molecules (for example C02, 02, etc.) or ions (Ca<2+>, HC03<->, etc.).

Furthermore, in the geological model, the two wells P1 and P2 are simulated as vertical axes extending at edges of the model, at locations corresponding to their real respective locations. The method is initialized by simulating a three-dimensional zone, representing the investigated area, said zone being in the shape of a cylinder of volume V centered on the well P1.

The geometry of the investigated area is evaluated by determining a preferred flowing path between the wells P1 and P2. To this end, during a first step 210 a search zone is defined; said search zone defining the cells of the geological model that may be added to the preferred flowing path. Thus, a cell may be added to the preferred flowing path only if it belongs to the search zone.

With reference to figure 1 b, the definition of the search zone is performed by deforming the three-dimensional zone mentioned above so that it incorporates both wells P1 and P2 and remains of volume V.

With reference to figure 1c, two vertical planes 1 and 2 are defined, each plane being orthogonal to a horizontal segment linking the wells P1 and P2, and the axis representing each well P1 , P2 being included in the respective plane 1 and 2.

The search zone Z is defined as the part of the three-dimensional zone that is comprised between the planes P1 and P2. Back to figure 2, the method then comprises a step 220 of determining a preferred flowing path between the wells P1 and P2. The preferred flowing path is formed by a series of adjacent cells of the geological model that link together the wells P1 and P2, and which maximize the permeability between P1 and P2.

The cells that are part of the preferred flowing path are chosen according to their respective permeability values.

The preferred flowing path may be determined by computing a percolation path or applying a method based on connectivity threshold.

According to a preferred embodiment, the preferred flowing path is determined by implementing a fast-marching algorithm weighted by the permeability values of the cells. The fast-marching algorithm is a well-known method based on the solving of an Eikonal equation as follows:

The skilled person may refer to the publication by J.A.Sethian,“A Fast Marching Level Set Method for Monotonically Advancing Fronts”, Proc. Natl. Acad. Sci., 93, 4, pp. 1591-1595, 1996 for the detailed implementation of this algorithm.

In the present case, the fast-marching algorithm is applied to the resolution of an equation representing the evolution of a front propagating radially from the well P1. The position of this front is obtained by the computation of directional arrival times, which depend on the distance from the well P1 and the permeabilities that are encountered. The propagation of the front is governed by the following Eikonal propagation equation:

u(x, f) = T(x)— t

Where u(x,t) represents the propagation front, c(x) is the propagation speed of the front and T(x) represents the duration for arriving at location x at the time t. T(x) is obtained by integrating on each volume unit the stationary equation below.

According to a preferred implementation of the fast-marching algorithm, each cell added to the preferred flowing path is selected based on an estimated time for a fluid to reach said cell, starting from the last cell added to the path. This time can be computed according to the distance between the centers of the adjacent cells, the respective permeabilities of the adjacent cells, the viscosity of the fluid, and possibly other parameters. The time can in particular be computed by a numerical method such as finite elements method.

As shown on figures 1d and 3a, at the end of step 220 a preferred flowing path is obtained, linking P1 and P2.

Back to figure 2, during a step 230 the volume V _p of the preferred flowing path, as defined by the sum of the volumes of its cells, is compared to the volume V of the investigated area.

If the volume V _p of the preferred flowing path is less than the volume V, then during a step 235 at least one cell is added to the preferred flowing path.

Referring to figure 3b, each cell added to the preferred flowing path must be adjacent to at least one cell already belonging to the preferred flowing path. As shown with the crossed cell in figure 3b, a cell not neighboring any cell of the preferred flowing path may not be added to the latter.

Moreover, each cell is added according to its permeability value, in order to increase the permeability of the preferred flowing path. In one embodiment, the selected cell is the cell exhibiting the highest permeability value among the cells neighboring the preferred flowing path. According to another embodiment, the selected cell is determined according to both its permeability value and position in the path, as the one increasing most the equivalent permeability between P1 and P2.

In an embodiment, during step 235, a single cell is added to the flowing path.

The step 330 of comparing the volume V _p of the updated preferred flowing path to the volume V is then repeated until the volume V _p of the updated preferred flowing path becomes superior or equal to the volume V, which acts therefore as a threshold.

In an embodiment, when the volume V _p of the updated preferred flowing path becomes superior or equal to the volume V, the process terminates at step 240 and the geometry of the investigated area is determined as being that of the flowing path.

In another embodiment, the process then comprises a step 250 during which the permeability K _p of the preferred flowing path is computed, during a step 250 based on the permeability values of the cells of the flowing path. The computed permeability is then compared during a step 260 to the equivalent permeability K _eq between the wells P1 and P2. This comparison allows a tolerance margin. Thus, if the permeability is comprised in the range [Keq-m, Keq+m], where m is the tolerance margin, the method terminates at step 240 and the geometry of the investigated area corresponds to the preferred flowing path.

If the permeability is not comprised in the range indicated above, the method further comprises a step 270 of updating the geological model, by simulating changes in the petrophysical characteristics of soil (in particular through the mineralogical composition of the soil) which result from infiltration of water. These changes are called diagenesis and can be the consequences of phenomena such as dissolution, precipitation or change of lithology, which can be simulated for an improved precision of the reservoir characterization.

In this embodiment, the movement of particles representing water infiltrated in the soil is simulated during step 270 in the geological model.

According to this embodiment, the geological model can further simulate fractures in the soil, and to this end comprises edge parameter values, for example the diameter of conduits d213v (vertical edge between two cells C1 13 and C213 in Figure 5), d123h (horizontal edge between two nodes N113 and N123), etc. The probability of stochastic movement of a particle in the second environment is calculated while taking into account these conduit diameter values, so as to simulate flow through the fractures.

The particles may be introduced at a given cell, for example C123, or at several cells. The introduction of particles can be performed at predetermined periodic intervals.

In this embodiment, it is assumed that the particles are subjected to two types of movement:“advective movement” and“dispersive movement”.“Advective movement” (or, for the latter,“advective direction”) is the most likely movement of the particle (or, for the latter, the most likely direction of movement of the particle). For a given cell, the advective movement is likely to take place in one direction and orientation, determined by the hydraulic gradient corresponding to the modeled region based on whether this region is or is not saturated. In addition, limestone lithologies can be characterized by a reaction index I Rij (or dolomitization index IDij in the case of limestone) describing the capacity of their mineralogical compositions to be modified (or respectively the capacity of the limestone to be transformed into dolomite).

One can thus predict that the particles introduced into the model may cause the mineralogical composition parameters of the first model (matrix) to evolve during the course of the simulation. Precipitating elements in significant proportions may cause a change of lithology and thus of the facies of the rock (for example a dolomitization phenomenon), in particular based on the aggressiveness of the particles, the current (or initial) mineralogical composition CMxyz of the matrix, the current reaction index IRxyz, the chemical content of the particle, etc.

The diagenesis of the rock can be also taken into account, namely the type of action of the particles on the rock during their movement. For example, the type of action of the particles on the rock is selected from among the following:

dissolution, whereby the water dissolves the rock,

precipitation, whereby the ions dissolved in the water are deposited on the rock,

change in lithology related to a change (increase or decrease) in the porosity.

Using dissolution as an example, the parameters taken into account can be the following:

the minerals to be dissolved (chosen for example from a list including limestone, dolomite, limestone dolomite, dolomitic limestone, argillaceous sandstone, and argillaceous limestone, in predefined ratios), a maximum porosity value for the matrix and/or a maximum conduit diameter value which can be obtained with dissolution,

a reactivity index (between 0 and 1 ), which can be defined for each of the minerals,

a transformation of the facies, defining the values both before and after dissolution,

an inhibitor of the dissolution reaction,

the kinetics of the dissolution reaction (taking into account, for example, factors such as the temperature, activation energy, etc.). In the case of precipitation, similar parameters can be defined to take into account the particular characteristics of precipitation in the rock in question. For example, the minimum porosity and diameter values can be defined instead of the maximum values defined for dissolution.

In the case of a change in lithology, it is possible to determine for example: the mineral to be transformed,

the mineral created,

the transformation of the facies, defining the values both before and after the changes in lithology,

the minimum and maximum rates of change of the mineral,

inhibitors or kinetics as mentioned above.

These parameters can be defined, where applicable, for either one or for both environments.

In this case, the particle aggressiveness describes the capacity of the particle for transforming the lithology or the chemical composition of the affected rock.

The particles can have variable characteristics, for example depending on their coordinates in the model, the simulation time, the distance traveled by a particle, the type of fluid, etc. In particular, the particle may lose aggressiveness in the case of dissolution, while its aggressiveness may be increased in the case of precipitation.

An example of the stochastic movement simulation model for particles is described in patent application PCT/FR201 1/052099.

FIGS. 6a to 6c illustrate the phenomenon of modifying the mineralogical composition of cells of a simplified, two-dimensional geological model for the sake of clarity.

In these Figures, in order to simplify the description, it is assumed that the rocks represented in the model are limestone rocks, but other compositions such as shale and sandstone are also possible.

Limestone rocks are readily soluble in water and consist mostly of calcium carbonate (CaC03). By way of an example, several minerals/rocks which may be included in the composition of a rock are listed below:

Calcite: a mineral consisting of calcium carbonate (CaC03), with traces of Mn, Fe, Zn, Co, Ba, Sr, Pb, Mg, Cu, Al, Ni, V, Cr, Mo; Dolomite: a mineral consisting of calcium carbonate and magnesium having the chemical formula CaMg(C03)2 and with traces of Fe, Mn, Co, Pb, Zn;

Clay: a sedimentary rock containing at least 50% alumina silicate, to which are added other minerals (quartz, feldspar, calcite, iron oxides).

FIG. 6a shows a two-dimensional gridded geological model at simulation time t. By way of simplification, this model includes 9 cells. This model may for example have been obtained from simulation software by electronic means, via a connection interface (for example a USB interface, an Ethernet interface, or a bus connection to a hard drive).

Each of the cells in this model has a mineralogical composition that is provided and fixed prior to the simulation. The mineralogical composition of cell C21 is indicated as CM21 (t), the mineralogical composition of cell C22 is indicated as CM22(t), and the mineralogical composition of cell C23 is indicated as CM23(t). The model is therefore described as having a local mineralogical composition parameter which can depend on the cell in question (i.e. the coordinates of the cell concerned in the model).

At time t-i, a particle 103 is present in cell C21 and has an aggressiveness a(t1 ) (or aggressiveness parameter). The aggressiveness of the particle can be represented by a vector of numeric values, where each numeric value of this vector represents a capacity of the particle for modifying a given mineral. By way of an example, it is assumed that this particle 103 has an aggressiveness vector which is equal to [calcite=0.7; dolomite=0.3; clay=-0.5](or simply [0.7; 0.3; -0.5]). In the present case, this vector indicates that the particle 103 has a high capacity for dissolving calcite (i.e. 0.7), a low capacity for dissolving dolomite (i.e. 0.3), and an average capacity for precipitating clay (i.e. -0.5).

In addition, the mineralogical composition of the cell C21 includes, for example, a large amount of dolomite and a small amount of calcite and clay: it is thus said that the proportions of dolomite, calcite and clay are the components of the mineralogical composition of cell C21. These components can be written down in vector form as [0.8; 0.1 ; 0.1]=CM21 (t). By extension of this concept, the mineralogical composition of the cell can also include chemical compounds which are not explicitly considered to be minerals; for example, the mineralogical composition of a cell may also include molecules (such as C02, 02, etc.), or ions (Ca<2+>, HC03<->, etc.).

It is possible to simulate the stochastic movement of the particle between times t-i and t-i+1 , particularly in relation with the methods presented in patent application PCT/FR2011/052099. For simplicity in the examples in FIGS. 6a to 6c, the simulated movement of the particle is considered to be vertical, in direction -y.

Between times t-i and t-i+1 , it is also possible to simulate the fact that the presence of the particle in the matrix of cell C21 induces a modification to the mineralogical composition of the cell (or local mineralogical composition parameter). In fact, if the mineralogical composition of the cell includes dolomite and the particle also has the capacity to dissolve dolomite, the ratio of dolomite in the cell will consequently be decreased while the ratio of other minerals will be increased. The same will be true for the precipitation capacity of a given mineral, for which the ratio of the given mineral will be increased.

In the example in FIGS. 6a to 6c, we have a(t ₁ )=[0.7; 3; -0.5] and CM21 (t ₁ )=[0.8; 1 ; 0.1]. In one embodiment of the simulation, the new ratios for the mineralogical composition (after modification) of the cell can be determined as follows:

Step 1 : intermediate mineralogical composition values are determined by multiplying the old ratio by the corresponding aggressiveness component and by subtracting the result of the old ratio calculation;

Step 2: once these intermediate values have been determined for each component, each of these intermediate values is divided by the sum of the determined intermediate values.

In the present case, the following result will be obtained in step 1 : [0.8 ^* (1 -0.7); 0.1 ^* (1 -0.3); 0.1 ^* (1 +0.5)] which is [0.24; 0.07; 0.15]. In addition, the result in step 2 will be [— 0,46 ;— 0,46 ;— 0,46] (because 0.24+0.07+0.15=0.46) or about [0.52; 0.15, 0.33]

Other calculation methods are of course possible. In particular, each of the new components of the mineralogical composition may be a function of a plurality of old mineralogical composition components or of aggressiveness components so as to represent the underlying chemical reactions in more detail.

The aggressiveness of the particle may evolve in response to this change in the mineralogical composition. Indeed, if the particle induces strong calcite dissolution, it is likely that its capacity to dissolve such a mineral will actually be reduced. In an embodiment of the simulation, the new aggressiveness components of the particle (after modification) can be determined by multiplying the old aggressiveness component by the associated ones' complement of the old mineralogical component (i.e. corresponding to the same mineral).

In this case, the following result would be obtained: [0.7 ^* (1 -0.8); 0.3 ^* (1—0.1 ); — 0.5 ^* (1—0.1 )] which is [0.14; 0.27; -0.45]. Other calculation methods are of course possible. In particular, a new aggressiveness component can be a function of a plurality of old aggressiveness components or local mineralogical composition components so as to represent the underlying chemical reactions in more detail.

Indeed, if the aggressiveness of the particle depends on its “ionic” composition (or, in inaccurate terms, its mineralogical composition), for example the concentration of acids or of Ca<2+> or even of HC03<->), the chemical mechanisms of precipitation or dissolution can cause the mineralogical composition of the particle to evolve and thus change its aggressiveness.

For example, based on the premise that the mineralogical composition of the particle 103 includes a high concentration of Ca<2+> ions and that the mineralogical composition of cell C21 includes HC03<-> ions in suspension, precipitation of calcium carbonate may occur (which is a frequent situation in aquatic environments, in particular in marine environments). Once the ions have been precipitated, the concentration of Ca<2+> ions in the particle is decreased and the particle thus loses its“aggressiveness”. It should be noted that in the latter case, the aggressiveness of the particle does not depend solely on the particle but also on the environment in which the particle is found: namely cell M21.

Therefore, in FIG. 6b and at time t-i+1 , the particle 103 is located in cell C22. The aggressiveness of the particle a(t-i+1 ) at this moment has been modified between times t-i and t-i+1 as a function of the modifications introduced into the mineralogical composition of cell C21. The mineralogical composition of other cells (i.e. C22 and C23) was not modified, as no particles were located inside these cells between simulation times t-i and t-i+1.

Also, at time t-i+2 in FIG. 6c, a repetition of the steps presented above can be carried out. In this manner, between times t+1 and t+2, the particle 103 moves from cell C22 to cell C23. The aggressiveness a(t-i+2) has been modified between time t-i+1 and time t-i+2 as a function of the modifications made to the mineralogical composition CM22 of cell C22. The mineralogical composition of other cells (i.e. C21 and C23) is unchanged, as no particles were located inside these cells between simulation times t-i+1 and t-i+2.

Obviously, although FIGS. 6a through 6c describe the movement of one particle, similar processing can be performed for a plurality of particles, for example 100,000 particles.

Once the movement of one or several particles have been simulated, and the geological model has been updated accordingly, back to figure 2 the method comprises a step 280 of computing the permeability of the preferred flowing path, taking into account the changes in the mineral composition of the cells forming the preferred flowing path, and hence of their respective permeability values.

Then the permeability value is compared to the equivalent permeability between the wells P1 and P2 considered as a threshold, during a step 290.

This means that, if the initial permeability value of the preferred flowing path was below - resp. above - the equivalent permeability, the method terminates at step 240 if the updated permeability value becomes superior or equal to - resp. inferior or equal to - the equivalent permeability.

Otherwise, if the updated permeability value has become closer from the equivalent permeability, while remaining below - resp. above - the latter, then the method comprises another step 270 of updating the preferred flowing path by further simulating movement of particles through the reservoir to induce changes in the mineral composition thereof. The method then repeats until the updated permeability of the flowing path reaches the threshold of the equivalent permeability between the wells.

However if the updated permeability value has become more remote from the equivalent permeability, then an error on the geological model used for the investigated area is detected at step 299.

The flowing path obtained according to the above method has the volume and equivalent permeability that have been determined during a well test. It also takes into account diagenesis simulation for an improved precision in its mineralogical composition. Furthermore, the method also provides an indicator for stopping diagenesis simulation, said indicator being the equivalent permeability of the preferred flowing path.

Last, it also provides an indicator of the reliability of the geological model that is used for performing diagenesis.