Stochastic Propagation of Discrete Fracture Networks

in: Proc. International Petroleum Technology Conference

Abstract

Fractures are ubiquitous structures occuring in a wide variety of rock types and tectonic settings over a broad range of scales. The permeability of these structural heterogeneities may, on average, be a few orders of magnitude higher or lower than those of the surrounding matrix rocks. Consequently, fractures and faults in rocks are essential elements of hydrocarbon migration and entrapment (Aydin, 2000). For that reason, the knowledge and modeling of the geometry and properties of individual fractures and of the fracture pattern are critical to better assess oil recovery and to optimize production of naturally fractured reservoirs. Two approaches can be selected to represent fractures : (1) using an homogenization method of the properties of the intact rocks, plus the individual fractures and its pattern, into a stochastic continuum, (2) or by a discrete fracture network (DFN) representation that models the characteristics of each individual fractures. We focus here on the DFN approach. The technique was created in the early 1980s (e.g. Priest and Hudson 1976, Dershowitz 1984, Chiles 1988) and has been continuously developed afterwards with many applications in civil, environmental and reservoir engineering and other geoscience fields. It emerges as a peerless technique to simulate explicitly 3D fracture patterns. However, the DFN approach based on marked-point process suffers from some restrictions and consequently does not "look" realistic. One reason is that it relies on a simplistic description of the geometry of fractures as a planar surface of predefined shape. Indeed, planar fractures or faults are rather the exception than the rule. Mean fracture orientation can change throughout the reservoir depending on the local driving stresses. Mechanical interactions between adjacent fractures, material properties of the rock and local heterogeneities may also influence the trajectory of fractures and, in some complex manner, it is the interplay of all these factors that is responsible of the resultant fracture pattern. Another limitation is the independence between point positions, between points and marks, and among marks trough the simulation process. Variations of the general Poisson point processes can take into account some spatial correlation among points, for instance clustering effect with parent-daughter model or minimal spacing between points with hard-core process (Stoyan et al., 1995, Chiles, 2004), but none of these can alone reproduce the whole complexity of fracture systems as observed on outcrops. This is a significant problem for dynamic reservoir studies, especially because the connectivity of the DFN realizations is not consistent. Jointly to the development of stochastic fracture models at reservoirs scales, the mechanical processes controlling the initiation, propagation and interactions of fractures have been better and better understood. This leads to the development of process-based algorithms (Renshaw and Pollard, 1994) which provide more realistic synthetic images of fracture patterns. They are, however, prohibitively time consuming for reservoirs scale applications and have been so far limited to simple 2D cases. Hybrid approaches can therefore be very useful. The present work motivation is to integrate in a stochastic framework insights from fracture mechanics theory. The proposed approach borrows an original idea of Srivastava et al. (2004) which is to limit the contribution of the marked-point process to the simulation of initial fracture seeds. Each one is then propagated to simulate the growth of fractures. We extend this method to full 3D propagation of fractures. We also change the way to determine the orientation of fractures at each propagation step. The procedure begins with the seeding of points, each corresponding to the center of a fracture. As in the DFN approach the point's positions follows a Poisson distribution law that expresses the probability of a number of points occuring in a volume, given the expected number of fractures per unit volume. In a second step, an initial fracture plane is generated for each central point previously simulated. Their geometrical properties are sampled in the probability density functions of fracture parameters (orientation and sizes) formulated according to available data. Once the distribution of initial fracture has been determined, the geometry of all fractures is simulated by sequential growth. For each propagation step, the fracture are visited by decreasing sizes. Considering fracture mechanics, this is consistent because longer fracture have the greatest impact on the final fracture pattern. By propagating first the longer fractures, we reproduce the effect of differential velocity propagation. Consequently, longer fractures having their process zones larger affect much more the propagation of smaller fractures and not the inverse. During propagation, the positions of new fracture extremities are determined. The first step is to examine if there is mechanical interaction between the growing fractures and any other fractures. If there is no interaction with other fractures, the position of the new extremity is simply the previous one incremented of a fixed step size in the current direction of propagation or in a direction given by an auxiliary orientation property. In case of interaction between process zones of fractures, the current direction of propagation is deviated. The influence of each fracture contained in the neighborhood is weighted depending on the fractures inter-distance and their relative size. The vector of deviation is computed as the resultant of all vectors, linking the fractures being propagated and those in the neighborhood, with distance-based weights obtained by kriging and size-based weights computed as the ratio of the current fracture size and of those of the neighboring fractures. The sum of the current vector of propagation and of the deviation vector gives the position of the new fracture extremity. This process is repeated until the fractures have been propagated a number of propagation steps defined by the user. Note, however, that fracture arrest my occur before reaching the total number of propagation steps. Fracture termination occurs when a growing fracture intersects another fracture or any other discontinuity such as bedding planes. In case of fracture intersection, the priority is always given to the longest fracture. The smaller fracture is either truncated or its free extremity is projected on the other fracture in the new direction of propagation. The locked extremities are flagged and are no more propagated. The remaining free extremities, however, continue to propagate until intersections occur or the total number of propagation step is reached. In case of intersection with a bedding plane, the fracture may arrest or not depending on a user-given probability. This probability may be also controlled by the rock's mechanical properties. This new method is closer to the mechanics of fractures and admits some variations to integrate all the factors that control the development of fracture networks. Thereby, resulting fracture pattern are geologically more consistent. This opens interesting perspectives for the characterization and modeling of naturally fractured reservoirs, such as, overcoming uncertainties about fracture density, direct integration of observed fracture in the simulation process, better assessment of fracture connectivity.

Download / Links

BibTeX Reference

@INPROCEEDINGS{Henrion2009IPTC,
    author = { Henrion, Vincent and Caumon, Guillaume },
     title = { Stochastic Propagation of Discrete Fracture Networks },
 booktitle = { Proc. International Petroleum Technology Conference },
      year = { 2009 },
       doi = { 10.2523/IPTC-14053-ABSTRACT },
  abstract = { Fractures are ubiquitous structures occuring in a wide variety of rock types and tectonic settings over a broad range of scales. The permeability of these structural heterogeneities may, on average, be a few orders of magnitude higher or lower than those of the surrounding matrix rocks. Consequently, fractures and faults in rocks are essential elements of hydrocarbon migration and entrapment (Aydin, 2000). For that reason, the knowledge and modeling of the geometry and properties of individual fractures and of the fracture pattern are critical to better assess oil recovery and to optimize production of naturally fractured reservoirs. Two approaches can be selected to represent fractures : (1) using an homogenization method of the properties of the intact rocks, plus the individual fractures and its pattern, into a stochastic continuum, (2) or by a discrete fracture network (DFN) representation that models the characteristics of each individual fractures. We focus here on the DFN approach. The technique was created in the early 1980s (e.g. Priest and Hudson 1976, Dershowitz 1984, Chiles 1988) and has been continuously developed afterwards with many applications in civil, environmental and reservoir engineering and other geoscience fields. It emerges as a peerless technique to simulate explicitly 3D fracture patterns. However, the DFN approach based on marked-point process suffers from some restrictions and consequently does not "look" realistic. One reason is that it relies on a simplistic description of the geometry of fractures as a planar surface of predefined shape. Indeed, planar fractures or faults are rather the exception than the rule. Mean fracture orientation can change throughout the reservoir depending on the local driving stresses. Mechanical interactions between adjacent fractures, material properties of the rock and local heterogeneities may also influence the trajectory of fractures and, in some complex manner, it is the interplay of all these factors that is responsible of the resultant fracture pattern. Another limitation is the independence between point positions, between points and marks, and among marks trough the simulation process. Variations of the general Poisson point processes can take into account some spatial correlation among points, for instance clustering effect with parent-daughter model or minimal spacing between points with hard-core process (Stoyan et al., 1995, Chiles, 2004), but none of these can alone reproduce the whole complexity of fracture systems as observed on outcrops. This is a significant problem for dynamic reservoir studies, especially because the connectivity of the DFN realizations is not consistent. Jointly to the development of stochastic fracture models at reservoirs scales, the mechanical processes controlling the initiation, propagation and interactions of fractures have been better and better understood. This leads to the development of process-based algorithms (Renshaw and Pollard, 1994) which provide more realistic synthetic images of fracture patterns. They are, however, prohibitively time consuming for reservoirs scale applications and have been so far limited to simple 2D cases. Hybrid approaches can therefore be very useful. 

The present work motivation is to integrate in a stochastic framework insights from fracture mechanics theory. The proposed approach borrows an original idea of Srivastava et al. (2004) which is to limit the contribution of the marked-point process to the simulation of initial fracture seeds. Each one is then propagated to simulate the growth of fractures. We extend this method to full 3D propagation of fractures. We also change the way to determine the orientation of fractures at each propagation step. The procedure begins with the seeding of points, each corresponding to the center of a fracture. As in the DFN approach the point's positions follows a Poisson distribution law that expresses the probability of a number of points occuring in a volume, given the expected number of fractures per unit volume. In a second step, an initial fracture plane is generated for each central point previously simulated. Their geometrical properties are sampled in the probability density functions of fracture parameters (orientation and sizes) formulated according to available data. Once the distribution of initial fracture has been determined, the geometry of all fractures is simulated by sequential growth. For each propagation step, the fracture are visited by decreasing sizes. Considering fracture mechanics, this is consistent because longer fracture have the greatest impact on the final fracture pattern. By propagating first the longer fractures, we reproduce the effect of differential velocity propagation. Consequently, longer fractures having their process zones larger affect much more the propagation of smaller fractures and not the inverse. During propagation, the positions of new fracture extremities are determined. The first step is to examine if there is mechanical interaction between the growing fractures and any other fractures. If there is no interaction with other fractures, the position of the new extremity is simply the previous one incremented of a fixed step size in the current direction of propagation or in a direction given by an auxiliary orientation property. In case of interaction between process zones of fractures, the current direction of propagation is deviated. The influence of each fracture contained in the neighborhood is weighted depending on the fractures inter-distance and their relative size. The vector of deviation is computed as the resultant of all vectors, linking the fractures being propagated and those in the neighborhood, with distance-based weights obtained by kriging and size-based weights computed as the ratio of the current fracture size and of those of the neighboring fractures. The sum of the current vector of propagation and of the deviation vector gives the position of the new fracture extremity. This process is repeated until the fractures have been propagated a number of propagation steps defined by the user. Note, however, that fracture arrest my occur before reaching the total number of propagation steps. Fracture termination occurs when a growing fracture intersects another fracture or any other discontinuity such as bedding planes. In case of fracture intersection, the priority is always given to the longest fracture. The smaller fracture is either truncated or its free extremity is projected on the other fracture in the new direction of propagation. The locked extremities are flagged and are no more propagated. The remaining free extremities, however, continue to propagate until intersections occur or the total number of propagation step is reached. In case of intersection with a bedding plane, the fracture may arrest or not depending on a user-given probability. This probability may be also controlled by the rock's mechanical properties. This new method is closer to the mechanics of fractures and admits some variations to integrate all the factors that control the development of fracture networks. Thereby, resulting fracture pattern are geologically more consistent. This opens interesting perspectives for the characterization and modeling of naturally fractured reservoirs, such as, overcoming uncertainties about fracture density, direct integration of observed fracture in the simulation process, better assessment of fracture connectivity. }
}