## Abstract

The interconnectedness of a subsurface fracture system depends on five fracture and fracture set geometry characteristics: density, orientation, dimensions (height and length), connectivity and aperture. Orientation and density can be determined with some degree of accuracy but the uncertainty ranges associated with the other characteristics will be large. Calibration with dynamic indicators will bring out the hydraulic properties of the network and can help to reduce uncertainty in the geometric characteristics. The fracture network geometry in the inter-well space can be described through geostatistical analysis but more commonly it is derived from seismic data, through edge detection or the analysis and calibration of seismic velocity anisotropy attributes or curvature analysis. Geomechanical analysis is necessary to evaluate the impact on network conductivity of the changing stress state in the reservoir. Using percolation theory, a qualitative assessment can be made of how the conductivity and interconnectedness of the fracture network is developed in terms of network clusters. This integrated analysis of static and dynamic fracture data leads to the formulation of conceptual models and a set of modelling rules. A Digital (Discrete) Fracture Network model, based on these concept(s), rules and inter-well interpolation data, can then be analysed in terms of expected network cluster size, distribution and hydraulics. This provides a control on the expectations embedded in the concepts prior to dynamic modelling. The uncertainties are normally quite large and very seldom will they allow construction of the definitive fracture model. The systematic and integrated analysis of fracture network geometry, constrained with dynamic indicators and subsequent network cluster analysis, are necessary preparatory steps for construction and analysis of dynamic models.

Fractured reservoirs are reservoirs in which production and recovery is influenced to a greater or lesser extent by fractures. They can be subdivided into four different types (cf. Nelson 2001; Allan & Qing Sun 2003) schematically indicated in Figure 1. The variability in fracture network interconnectedness, and in the architecture and properties of the matrix, are the basic reasons that fractured reservoirs show a large variety of behaviours during hydrocarbon production. These large uncertainties make the appraisal, development and management of fractured reservoirs difficult. Failure to assess uncertainties properly leads to missed opportunities and low hydrocarbon recovery.

The special nature of fractured reservoirs lies in the interaction between, the (relatively) high pore volume, low permeability matrix (the storage domain) and the low pore volume, high permeability fracture system (the flow domain; Fig. 2). This interaction is a function of matrix architecture and fracture network geometry, but also of the mechanisms and physical processes that control the transfer of hydrocarbons from the matrix to the fracture network. The initial and developing stress state and the presence or absence of an aquifer also influence performance. For more mature fractured reservoirs the possible impact of water flood or gas injection must be incorporated in the analysis. In Type I, II & III fractured reservoirs (*sensu* Fig. 1) the geometry of the fracture network is often the main control on well location and orientation. In Types II & III, if the network does not interact with matrix in a sufficient manner, the time it takes to produce hydrocarbons may make a field totally uneconomic.

The natural starting point for the analysis of a fractured reservoir will be the collection and processing of fracture data followed by analysis of the geometrical aspects of that data. In tandem with the geometric analysis, dynamic data must be analysed to calibrate the possible ranges in interconnectedness of the network. Additional data (e.g. seismic) must be analysed and again calibrated with static and dynamic data to constrain the inter-well space. This integrated analysis, which leads to the formulation of conceptual models and rules for fracture network model construction, will be the focus of this paper. The discussion will be illustrated by the analysis of image log data. It is stressed here that what is presented does not represent a complete analysis or modelling exercise. The examples are used only to demonstrate the applicability of the methodology. In the context of this paper fractures are defined as ‘approximately planar features along which cohesion has been lost and only small displacement has taken place.’ This includes displacement normal and parallel to fracture walls. Some would call the latter faults, or micro-faults if displacements were small. The inclusion of faults in a fracture network, whatever their displacement, is entirely appropriate where they provide pathways for flow.

## Data sources for fracture network geometry

In most cases, fracture data will be derived from core analysis and image log interpretation. A third potential data source is outcrop analogues, used to fill scale gaps in the subsurface dataset. Seismic data is an indirect source, the interpretation of which will deliver a fault and horizon framework that forms the structural envelope in which the development and distribution of fractures can be analysed. It is also the main source of information on strain related subseismic structures in the inter-well space.

### Core

Cores are often taken in a few wells only and over limited parts of the reservoir section. They provide the only direct observation of fractures in the subsurface, but that is limited to the borehole scale. Information on the following can be obtained:

#### Orientation

Fracture orientation can be measured directly from core if it is oriented (e.g. with a scribe line), or indirectly using bedding plane orientation or palaeomagnetic methods (Davison & Haszeldine 1985). Under-sampling of fractures making a small angle with the borehole is a common limitation (Fig. 3).

#### Connectivity

Fracture systems develop in a variety of styles depending on stress state and the characteristics of the rock sequence. As a result fractures will abut (or terminate) against other fractures, cross each other or simply terminate within rock. This is captured in the number of Y-, X- or I-nodes (Bech *et al.* 2001; Fig. 4a). The second aspect of connectivity is the interaction of fractures with lithostratigraphy (Fig. 4b; Loosveld & Franssen 1992). Stratabound fractures occur where bedding planes, and/or intervening shale or otherwise incompetent layers, inhibit propagation (cf. Anderson 1981; Hanson *et al.* 1982; Odling *et al.* 1999). The vertical persistence (or crossing probability) is defined as the ratio of the number of fractures that cross the boundary to the total number of fractures in the stratum (Bech *et al.* 2001). Together stratabound and non-stratabound fractures may form the total fracture spectrum and it requires some filtering to separate them into different sets. Quantification of both types of connectivity relationships will help to (a) decide on the fracture mechanism(s), and (b) integrate core observations with the regional structural geological history.

#### Mechanical stratigraphy

Combining (litho-) stratigraphy and fracture characteristics, e.g. connectivity relationships and density, can reveal major contrasts between different parts of the reservoir sequences. These contrasts are commonly observed in a vertical direction related to mechanical contrasts between various units (Wennberg *et al.* 2006). Mechanical contrasts can also be observed in a horizontal direction where they may indicate partitioning of deformation related to folding or faulting (Stephenson *et al.* 2007).

#### Relative age

This can potentially be derived from abutting and cross-cutting relationships (Fig. 5). Consistency in such relationships is a reliable criterion (Ortega *et al.* 2006). The geometric relationships give only relative ages but with more detailed study of fracture fills (cf. Laubach *et al.* 2004), ages can be bracketed against the burial history and form an important input into the overall geological history.

#### Aperture

This can be estimated by measuring the width of cemented or partially cemented fractures in core. The cementation history may show crack–seal episodes of progressive opening and aperture estimates must take that history into account (Laubach *et al.* 2004). The aperture estimate may or may not be correct. Aperture of open fractures will be affected by the *in situ* stress state, which may be different from the stress state that controlled the formation of currently mineralized fractures. Estimates of open fracture aperture from a core are always very uncertain because the width of open fractures at the surface may not be representative of the width at depth. Unloading effects and associated stress release, as well as damage to the fracture during coring and subsequent core handling, will affect the measurements.

#### Dimensions

The dimensions of fractures are captured by length, height and aspect ratio (= length/height). Vertical core may allow heights to be sampled, whereas core from horizontal wells allow for the sampling of length. Both require the observation of fracture terminations and those are usually few and far between.

#### Induced fractures

The process of coring and subsequent handling of cores leads to stress changes, which can cause the formation of coring-induced fractures. They may use and enhance suitably oriented natural fractures when they propagate along them. There are few definitive criteria for distinguishing coring-induced fractures, but some fracture orientations and geometries are typical and when present should alert the observer (Kulander *et al.* 1990).

#### Fracture type

To distinguish between tensile and shear fractures is really only possible in core (Kulander *et al.* 1990; Loosveld & Franssen 1992). Fracture types and the spatial relationship between fractures, especially that between en-echelon fractures, allow an assessment to be made of the remote stress state during fracture formation (Olsen & Pollard 1989). The distinction between fracture types is important because their spatial distribution and dimension characteristics may be significantly different (Odling *et al.* 1999; Cacas *et al.* 2001). The information will help to subdivide fractures in specific sets, assist in determining age relationships, and in general, supply important pointers to the structural history (Hancock 1985).

#### Morphology

Characterizing the fractures in terms of morphology and content is important in the analysis of fractured reservoirs. There is no definitive scheme that gives all the definitions and qualifiers that should be used. It is a fit for purpose argument that determines the choice for a descriptive scheme. Some important characteristics and contrasts (Kulander *et al.* 1990) are captured schematically in Figure 5.

### Image log

Borehole image logs (BHI) are usually collected over a greater length interval than core, and thus provide a more continuous image of the subsurface. There are two main types of BHI tools (Chueng 1999). Micro-resistivity tools create electrical images using pad-mounted resistivity button arrays. Ultrasonic tools create acoustic images using a rotating transducer. In addition, lower resolution BHI tools, based on density-log technology, can be run on drill pipe to provide near real-time information on large fractures and faults (Prensky 1999). BHIs can be used to gather fracture data in the same categories as core data. In general, the same restrictions and possibilities apply but there are differences. There is the obvious difference in resolution and quality between a BHI and a core. In general BHIs will under-sample the fracture system and the resolution restrictions will mean that certain details cannot be observed (Rawnsley *et al.* 2004).

#### Resistive v. conductive

Image log features with a resistive response are generally interpreted as being cemented or closed. Features that are filled with a fluid (formation water or drilling mud) that is more conductive than the host rock, are labelled as conductive and interpreted as fractures that are open and thus potential conduits for fluid flow (Lofts & Bourke 1999). This need not always be true. Certain types of cement or fracture fill could produce a resistivity contrast with the matrix equivalent to that between water and matrix, and open fractures drilled with oil based mud can be more resistive than a high water saturation reservoir.

#### Aperture

For micro-resistivity image logs, the measured signal can be considered to be a function of fracture aperture and the conductivity of the fluid in the fracture (Luthi & Souhaité 1990). Resolution depends very much on the resistivity contrast between rock and fluids. Similarly, in ultrasonic tools the amplitude response can be related to aperture but is sensitive to the roughness at the intersection with the borehole wall (Chueng 1999). Verga *et al.* (2002), using apertures of cemented fractures in thin sections, evaluated the electric fracture aperture using neural network techniques. This process yielded different distributions for core and image log data. In principle, aperture determination from image logs is a possibility, but resolution issues will limit accuracy (Luthi & Souhaité 1990). The obvious advantage is that aperture derived from image log data is observed under reservoir conditions.

#### Drilling induced fractures and borehole breakouts

These are the result of redistribution of the stresses around the borehole while drilling, which may cause the rock to fail (Lofts & Bourke 1999). Depending on borehole orientation, stress orientation and magnitude and mud conditions in the borehole, the rock may potentially fail in tension (tensile fractures in the direction of maximum horizontal stress) or in compression (borehole breakouts in the direction of the minimum horizontal stress). Drilling induced fractures will show up on image logs as open fractures oriented parallel to the borehole axis. Borehole breakouts are usually indicated as broad bands separated by 180°, again parallel to the borehole axis (Lofts & Bourke 1999). The importance of induced fractures and borehole breakouts lies in the fact that they can provide information about the present-day stress field (Trice 1999; Zoback *et al.* 2003).

#### Core v. image log calibration

Qualitative and quantitative calibration of fracture characteristics between BHI and core is essential to reduce uncertainty in concepts and models. The calibration helps to estimate the under-sampling of BHI relative to core and to filter out induced fractures. The data from both cores and image logs is presented in itemized lists which contain information, logged against depth, on fracture type, azimuth and dip, and as much as can be determined on the other characteristics discussed in this section (Fig. 5). Important characteristics are: mechanical fracture type; whether fractures are open, closed or partially closed and show wall to wall contacts or cement; solution effects; type and character of fracture fill and relative ages of fill; aperture and an estimate of the topography of the fracture walls, usually called roughness and sometimes expressed in the Joint Roughness Coefficient (cf. Kulander *et al.* 1990).

### Outcrop analogues

Outcrops provide an opportunity to sample fracture characteristics, which, due to the limited sampling by cores and image logs, are difficult to capture in the subsurface. This applies above all to fracture dimensions and connectivity relationships. It can be used to guide the assessment of specific network aspects in a general or qualitative sense (Mercadier & Mäkel 1991; cf. Van Dijk 1998). Even if all the aspects of the geological history are comparable, the fact that the analogue sequences are found at the surface will introduce a difference. The change in stress state related to the rise of the sequences from the reservoir level (assuming that the outcrop sequence came from that level) to the surface can lead to a considerable change in the fracture network characteristics, such as density. Outcrop data is a useful source of data to construct a mechanical stratigraphy (Trice 1999; Wennberg *et al.* 2006; Stephenson *et al.* 2007), even if only in relative terms. In using outcrop analogue data, a comparison of all aspects of the geological history must be made and a compensation for the potential impact of the difference in depth and history applied before the data is used to fill the subsurface data gaps (cf. Cacas *et al.* 2001). Given that the problems are considerable an element of doubt about the appropriateness of the analogue will remain.

## Analysis of fracture network geometry

The interconnectedness of a fracture network is determined by five geometrical fracture characteristics. To build a Discrete Fracture Network, knowledge of the distribution of these characteristics (geographically, vertically and in value ranges) is a prerequisite. The geometrical fracture characteristics are:

orientation of specific fractures or fracture sets;

density expressed as spacing between fractures in general or in specific sets and clustering describing specifics of spacing distributions;

connectivity, which describes the interaction and overprinting relations between fracture sets and reservoir sequences;

aperture or width and specific properties, such as roughness and cement fill, which impact on hydraulic properties;

length and height, which captures the horizontal and vertical extent of fracture and fractures sets.

The core and/or image log analysis and interpretation will result in an overview of specific fracture sets and relationships between them, and possibly a first pass mechanical stratigraphy scheme. Grouping in fracture sets is often based on orientation (both dip and/or strike direction) and fracture type, but any meaningful attribute or attribute combinations can be used. It must be consistent with the geological history of the reservoir and over- and under-burden (Cacas *et al.* 2001). The fault and horizon interpretation and the geological history provide the overall framework in which the development and distribution of fractures and fracture sets, and their relationships in time and space are analysed. The geological history should comprise an assessment of the various structural phases and stress states the reservoirs experienced over geological time and provide information on the kind of fracturing regimes (regional, fold or fault related) that were involved in fracture network formation. It also needs to detail the kinematics of the fault–horizon framework (cf. Gauthier *et al.* 2000) and should include the effects that diagenesis and hydrocarbon charge had on matrix and fractures.

The analysis of data has to take into account its source in terms of dimensions of the object and the derived parameter, and the scale at which data is sampled (cf. Van Dijk 1998). Data derived from boreholes and scan lines is essentially one-dimensional, whereas data derived from maps and outcrop surfaces is two-dimensional. Derived parameters from borehole data can have no dimensions (number of fractures), be one-dimensional (e.g. number per unit length) or three-dimensional (e.g. orientation of a fracture plane). The inherent resolution of the data source means that mixing data sets from various sources, at various scales, is a hazardous affair if the impact of the resolution and scale differences has not been analysed.

### Orientation

Plotting fracture data in stereographic projections and rose diagrams gives a powerful visual overview of the data (Fig. 6). Rose diagrams or circular histograms show the contrasts in fracture strike orientation. Poles plotted on stereonets present a view of the fracture distribution in space and of angular differences. They can be constructed as equal angle nets, which honour angular relationships, or equal area nets in which spatial distribution is better represented (Phillips 1971). Contour diagrams provide focus on the centres of gravity but lose the contribution of individual measurements. The distribution of poles and/or strike orientations is used to make a first differentiation in fracture sets for specific reservoir or mechanical units. If the description of the fracture network is intended to be at a high level only, it may be sufficient to read off characteristic directions for datasets straight from pole plots or contour diagrams.

Various distributions can be used to characterize a fracture dataset statistically. Commonly used are parameters based on the Fisher distribution, the equivalent for three-dimensional data to the normal distribution (Davis 2002). The data set is characterized by azimuth and dip of a mean vector and a concentration parameter *K*. There are tests to determine if a data set is Fisher distributed, but as a simple rule of thumb a *K* value of > 5 is often good enough to use the Fisher parameters. If that is not the case, another way of characterizing the data is through eigenvector analysis. Eigenvectors and their eigenvalues are calculated from a 3×3 matrix of the sums of squares and products of the direction cosines of a directional data set. They represent the greatest, intermediate and least moment of inertia of a data set, and can be used to locate the principal direction and describe the dispersion in the data set (Mardia 1972; Davis 2002). It makes little sense to use these methods on a multimodal set with clearly separated clusters. There are methods to tackle even that but splitting the set may be just as simple and effective. Splitting the data may also create problems. Defining orientation groups may work well for some fields but be difficult in others when the orientation of otherwise related sets changes in a geographical sense or with position in the mechanical stratigraphy.

### Density

#### Spacing

Density analysis involves the quantification of the spacing between fractures or of its reciprocal, frequency, i.e. the number of fractures per unit length (Ortega *et al.* 2006). Fracture spacing is influenced by lithology, porosity and grain size, bed thickness and stress history (Narr & Suppe 1991). Apart from the probable link to bed thickness, there are no hard and fast rules to determine spacing other than by measuring it in core or image log. Strictly speaking spacing is the orthogonal distance between two fracture planes. The standard measure of spacing in the subsurface is the downhole spacing i.e. the distance between two fractures measured along a borehole (Fig. 7). That is seldom the orthogonal distance since the value depends on the relative orientation of fractures and borehole. Spacing in a fracture set with little dispersion in orientation, and a reasonable angle with the scan line, can be corrected by using the normal downhole spacing or applying the Terzaghi correction (Fig. 7). The latter can be applied in various ways depending on orientation differences between fracture set and scan line (cf. Davy *et al.* 2006). In its original form the correction consists of multiplying downhole (or scan line) spacing by the cosine of the angle between the mean pole of the fractures and the scan line (Terzaghi 1965; Priest 1993). It is not an exact solution to the problem, but the errors introduced by using a scan line of relatively uniform orientation can be reduced to a tolerable level (cf. Terzaghi 1965). Whatever correction is applied, spurious results will impact the summary statistics, especially when data sets are small. To avoid applying very high correction factors to fractures subparallel to the borehole, a truncation limit is sometimes applied or values or value ranges are weighted. The impact of sampling blind spots (Fig. 3) could be compensated for by comparison with data from wells with different orientations. But regional differences may mean that the absence of fractures in the blind spot is genuine. If compensation is applied it will increase uncertainty in the spacing calculations.

Commonly used summary spacing statistics are arithmetic average, standard deviation and coefficient of variation (Cv). The latter is defined as standard deviation divided by the arithmetic mean and is a measure that is used to analyse clustering (Gillespie *et al.* 1993). In cumulative frequency plots the data sets can also be characterized by curve fitting techniques. The average (corrected) spacing is subjective when data is subdivided in various sets and it can only be relative when derived from image logs and can therefore be best used qualitatively rather than quantitatively. Possibilities are simple position in range schemes or breaking the range of values, for related sets, into quartiles.

Calculating fracture intensity measures that are related to the dimension of the measuring region can reduce the problems with spacing parameters. P_{32} is a very useful measure for the three-dimensional intensity of fractures and is an important parameter for fracture network modelling. It is defined as the area of fractures in a volume (m^{2}/m^{3}). Assuming the wellbore to be the volume, then the fracture area for any fracture is determined by the intersection plane with the wellbore (Fig. 8). This carries the implicit assumption that it will cross the whole wellbore. P_{32} is scale independent and not influenced by orientation effects and fracture size (Dershowitz & Herda 1992; cf. van Dijk 1998).

#### Clustering

A first pass assessment of the clustering in a fracture data set is made through the coefficient of variation Cv. Where standard deviation is large with respect to mean, large spacing values are mixed with small ones. The Cv will be >1 and the system is clustered. When Cv is *c*. 1 the large spacing values are balanced by the low values and the system is randomly spaced. Sets with Cv values <1 are uniformly or regularly spaced (cf. Gillespie *et al.* 1993). The scaling characteristics of different parts of the fracture network (geographically as well as fracture-type related) guide interpretation and extrapolation of the data beyond the borehole. In the 1980s and 1990s many publications appeared using fractal analytical methods to address these scaling issues (e.g. Gillespie *et al.* 1993; Hewett 1994). At the simplest level fractal geometry can be defined as being scale independent (or self-similar) at all scales between defined limits. In other words any portion of the system is a scaled version of the whole system (cf. van Dijk 1998).

The standard way of analysing fracture data for fractal geometry is by plotting the cumulative frequency of a fault or fracture attribute in a log–log plot (Walsh *et al.* 1994). Commonly used attributes are spacing, throw, height and length. If the data adheres to a power-law distribution, i.e. a straight line on a log–log plot, the data can be said to have fractal geometry. The distribution is described by the fitted curve and the function exponent (*D*) is taken as the fractal dimension (Gillespie *et al.* 1993). Data in such plots generally show the effect of truncation and censoring (Priest 1993; Bech *et al.* 2001). The former is indicative of a loss of data in the small value realm. For image log data, that could be a consequence of resolution (Ortega *et al.* 2006). Censoring at the other end of the curve occurs when the values of the plotted measure approach the limit of the area of observation. It is important to recognize and analyse these effects because they can lead to incorrect rejection of fractal geometry when that is masked by data problems (Gillespie *et al.* 1993; Ortega *et al.* 2006). A fundamental requirement for the analysis of populations in a log–log frequency plot is that the value range in the data set is large enough. Walsh *et al.* (1994) state that a range of at least one order of magnitude must be sampled.

Another method to analyse the geometry of fracture patterns is the box-counting technique (Walsh & Watterson 1993; Hewett 1994). The method consists of superimposing square grids of given side length *d* on a fracture pattern and recording the number of boxes containing fractures *n*. After repeating with boxes of different sizes the curve log *d* over log *n* is plotted. For a fractal pattern the central segment of the curve is a straight line of which the slope *D* is by analogy called the fractal dimension. Walsh & Watterson (1993) tested the methodology and formulated guidelines for its use with outcrop data. The use of this method in wells, for what is effectively one-dimensional fracture data from the subsurface, has been challenged (Harris *et al.* 1991), because it is sensitive to data resolution and the shape of the area analysed (Cowie *et al.* 1996).

If fractal geometry can be proven for a data set, it means that fracture properties can be extrapolated beyond the domain of the wellbore. Similarly, if it can be demonstrated at the seismic scale, it can be extrapolated to the smaller scale and be used as a proxy for fracture characteristics (Van Dijk 1998). But relating millimetre-sized discontinuities, measured in core or image log, with kilometre-sized structures is not without problems (cf. Gillespie *et al.* 1993). The correspondence between, for example, small and large scale faults showing similar geometries over many orders of magnitude of scale may have been demonstrated in some cases (Tchalenko 1970; Hewett 1994; Pickering *et al.* 1997), but to assume that fault and fracture data sets have, by definition, fractal geometry for the various attributes is far too simplistic (cf. Harris *et al.* 1991; Loosveld & Franssen 1992; Walsh & Watterson 1993; Cowie *et al.* 1996; Nicol *et al.* 1996). Strictly speaking it must be demonstrated that fractal geometry exists at both scales and that they fit the same distribution. Loosveld & Franssen (1992), considering networks of stratabound extensional features in a particular bed, concluded that these are characterized by regular spacing, i.e. are not scale invariant. They postulated that fracture networks consisting mainly of such stratabound systems would not be fractal (cf. Odling *et al.* 1999). Walsh *et al.* (1994) summarized the problem when density of subsurface data is based on total observations of fractures with and without some offset. A non-differentiated subsurface data set will overestimate density when compared with a population of seismic faults. It implies that fault data can only be used to estimate small-scale fault density and that density of fractures without shear displacement cannot be derived in this way.

### Connectivity

The relative number of I-nodes to X- and Y-nodes, as well as the ratio of number of X-nodes to number of Y-nodes, can be analysed by plotting the nodal proportions in a ternary connectivity diagram (Fig. 4). These ratios are considered as first pass indicators of the interconnectedness of a fracture network. Systems dominated by I-nodes will show very low interconnectedness but with a larger proportion of X- and/or Y-nodes that will improve (Fig. 4). The interconnectedness will increase more with proportionally more X-nodes because there is the chance that they will cross yet another fracture (Manzocchi 2002). Systems that are developed in orientation clusters (i.e. relatively isolated, independent clusters of uniformly oriented fractures; Fig. 9) will show a low number of X- and/or Y-nodes but will have high Cv values. Density-clustered systems, which combine fractures of various orientations in localized clusters (Fig. 9), will show large number of X- and/or Y-nodes and also have Cv values that will indicate clustering (Manzocchi 2002). As depicted, these systems are end-members and natural systems can be expected to show more heterogeneity in cluster development. Note also that when systems are not similar, the nature of clustering may differ with scale.

Fracture systems that are characterized by stratabound fractures, i.e. have a low vertical persistence, are usually regularly spaced, e.g. have Cv<1 (Loosveld & Franssen 1992; Manzocchi 2002), with spacing related to thickness of the strata (Narr & Suppe 1991). Non-stratabound systems are typically clustered and will have high vertical persistence values (Bech *et al.* 2001; Fig. 4). The consequence of abutting against bedding is that fracture aspect ratios, especially in thin-bedded sequences, could be very high for stratabound systems. For non-stratabound systems the aspect ratio tends to be much lower (Loosveld & Franssen 1992; Odling *et al.* 1999). Connectivity relationships are difficult to establish from core and virtually impossible from image log (Trice 1999). From core it may be possible to collect a small data set for both connectivity aspects, but the uncertainty that must be attached to the connectivity of the network will be large. The connectivity picture is more often than not completely derived from outcrop analogue data.

### Aperture

If aperture data is available from core or image log, the first-pass evaluation is straightforward. Basic statistics can be used to characterize the distribution of aperture for various sets. There are other methods to at least estimate what the fracture aperture could be; drilling mud-loss analysis (Lietard *et al.* 1996; Sanfillippo *et al.* 1997), combinations of various wireline logs (Trice 1999) or combining porosity data with fracture density (Bona *et al.* 2003) are some of them. In general, any method gives a rough estimate of aperture only. It is difficult to differentiate between matrix and fracture porosity from cores and logs combined, and density data, especially from image logs is not an absolute number. Measured (=mechanical) aperture would be equivalent to the hydraulic aperture if the fracture walls formed perfect parallel plates. This is virtually never the case and deviation from that ideal situation causes a reduction in flow (Witherspoon *et al.* 1980). The more pronounced the roughness of the fracture walls, the more flow is reduced. In engineering geology a number of measures are used to describe roughness (Priest 1993). commonly used in geological problems is the joint roughness coefficient (JRC) empirically defined to estimate the rugosity of fracture asperities (Barton *et al.* 1985; Bakhtar *et al.* 1985). The JRC may be subjectively determined by comparison with standard roughness profiles (Priest 1993).

### Height and length

Cores and image logs do not sample enough of the subsurface, relative to fracture dimensions, to give a data set large enough to determine height and length distributions with any degree of accuracy. Additional data may come from two methods to estimate fracture dimensions. The first is data from outcrop analogues and the second is subsurface faults. The problems with the use of outcrop data have been outlined above. Even if the conditions for their use are met, a considerable uncertainty must be attached to outcrop data that fills the gaps in subsurface data. The implied assumption when seismic fault data is used is that the same power-law describes the length and/or height distributions of faults and fractures (Loosveld & Franssen 1992; Pickering *et al.* 1997). Applied to the problem of length and height, plotting the fault data as derived from seismic and fitting a power-law through the data would then give an indication of the overall length or height distribution of the fault–fracture continuum. One problem is that the total fault set is the result of the total geological history. It may be related to the total fracture set but it is questionable whether it always forms a continuum with the subset of open fractures, and it should not be assumed that fault or fracture data, always adhere to power-law relationships (Nicol *et al.* 1996).

### Network geometry analysis example

A dataset from three wells is used to illustrate some of the points on fracture network analysis made above. These wells are widely spaced (Fig. 10) so the problem addressed is akin to one with a few appraisal wells where fracture data is available, but could equally well represent a situation where image logs are scarce. The wells penetrated three stratigraphic units, Units 1, 2 & 3, which show relatively small differences in matrix permeability development. Fractures are interpreted from micro-resistivity image logs. Only fractures that were interpreted as conductive are considered here. The fracture orientations in Unit 2 are detailed in Figure 11. Well A, drilled in a 350° orientation, shows some effect of undersampling but well B, orientated at 290°, is much less affected. The part of well D that was logged, changes orientation from 070° to 020° and may show some undersampling in Unit 3. Fractures are in general developed in consistent directional sets over the units. A broadly NE–SW direction and a north–south direction are the main orientations. There seems to be variation in direction of fracture sets in a regional sense. It will therefore not be simple to define uniform ranges to capture the main orientations in the fracture network. The four ranges indicated capture some of it but in several instances (e.g. Well D, Unit 2) the range boundaries will split what are directionally consistent trends. Variability can also be seen in dip, from near vertical to as low as 30°. Overall dispersion seems to be the rule because only occasionally are poles really clustered.

The average P_{32} values for each unit shows that well *D* is the most intensely fractured (Table 1). Wells A and B have lower but more-or-less similar values. The normal downhole spacing average for well A correlates with the P_{32} but less so for well B, showing that, in this case, it is not an optimal estimator of fracture density. The calculated fractal dimensions are low, in seven cases lower than 1.0, which indicates a relatively low number of small spacing values. In most cases the cumulative frequency distribution shows a better logarithmic then a power-law fit (Fig. 12; Table 1), but that could be deceptive because a logarithmic fit is less sensitive to the effects of truncation and censoring. The Cv values are indicative of clustered development only in a few cases, for sets that have low *D* values. Taken together this casts doubt on the fractal nature of the total set of fractures per unit. The spatial distribution and spacing data can be interpreted to indicate superposition of several orientation sets. The latter point is reinforced by density measures calculated for data sorted into the four orientation ranges (Table 2). These statistics must be used with care because application of uniform orientation ranges splits what seem to be consistent orientation sets in some units. The results show a large scatter in the spacing averages that are influenced by spurious results in small data sets. In general a considerable increase in Cv values is shown. This points to a considerable degree of clustering and potentially fractal geometry, which would be concealed if multiple sets were inappropriately grouped together. The uniformity in spacing averages across units in a specific well shows that there is an apparent absence of mechanical stratigraphy.

### Summary

In their simplest form, the end products of the geometric analyses are presented as high–medium–low ranges for the various parameters or fixed orientations for specific orientation sets. With adequate data coverage and quality, rules can be more sophisticated, e.g. a distribution function based on curve fitting of available data. Obviously there are differences in quality and size of the fracture data set from one fractured reservoir to another but some problems will be universal. This means that we can generalize the relative value derived from the analysis of the geometrical aspects for the subsequent modelling of the fracture network. The rules determined for fracture orientation and density from core data can be considered to be hard. Density derived from image logs indicates relative differences only, due to under-sampling. Connectivity and aperture are not adequately captured due to sampling limitations. Core and image logs do not provide adequate information on fracture dimensions. That gap is often filled by extrapolation from subsurface fault data or with data from an outcrop analogue. Both methods are not without problems and fracture dimensions come associated with a considerable uncertainty. The consequence of the above is that some of the parameters, notably aperture, are often used as empirical matching parameters in dynamic modelling (Gilman 2003).

## Dynamic calibration

Before a fracture network model is constructed, the fracture network geometry and geometric modelling rules must be calibrated against data that indicates flow into the borehole. Some of these so-called dynamic indicators can pinpoint flow at a specific point or short interval along the wellbore. Examples are drilling mud losses, production logging tools (PLT), temperature logs and mini drill stem tests (mini-DST). Well tests and production data give information on flow for a larger borehole section and a larger reservoir volume. With reservoir depletion significant changes in dynamic behaviour may be introduced. Such changes could be related to changes in effective stress with depletion. The effectiveness of a fracture network can be described in terms of percolation theory (Manzocchi 2002). The use of this theory, linking geometry to flow, helps to formulate constraints for those network characteristics that are difficult to capture by geometry alone. In general, dynamic indicators are used to narrow the possible ranges rather than to generate new ideas about the fracture network.

### Fracture networks and percolation theory

The basic premise of percolation theory is that flow is largely controlled by the continuity of the permeability contrasts. In a network with low fracture density, fractures will be isolated and few, relatively small clusters of intersecting fractures will be present. With increasing density the number and size of network clusters will grow until a point is reached where one, the so-called spanning or percolation cluster, grows rapidly at the expense of others. When the density of fractures is expressed as the probability of a fracture occurring in a grid cell, a critical probability determines when the percolation cluster is formed. This is called the percolation threshold (Stauffer & Aharony 1994; King *et al.* 2001). Analysis of fracture network models shows that with increasing interconnectedness, fracture permeability (*K*_{F}) increases slowly until, on the scale of the analysed block, a level of fracture saturation is reached where *K*_{F} increases rapidly (Manzocchi 2002). After that point, with increasing interconnectedness, *K*_{F} will still increase but far less dramatically (Fig. 13).

Bour & Davy (1997, 1998) analysed both two- and three-dimensional fault network models. Variation of density and length distribution with power-law relationships allowed definition of what they labelled the topology of connectivity. The system connectivity is defined for specific ranges of fault length fractal dimension *D*. However, in a fracture network the percolation threshold depends on the geometry of that network (Odling *et al.* 1999; Fig. 13). Orientation, length, height and aperture are individual fracture characteristics and represent the first level of characterization of network geometry (Manzocchi 2002). Connectivity, density and clustering are determined by the specifics of fracture sets and define the second level of characterization (Bech *et al.* 2001). Connectivity and clustering are the primary influences on the flow potential of a fracture network (Odling *et al.* 1999; Manzocchi 2002). Percolation threshold distributions have been determined, as a function of these two parameters, by Manzocchi (2002) for simple, two-dimensional systems with relative uniform length distributions. The nodal proportions in a ternary connectivity diagram provide a first-pass indicator of the interconnectedness of a fracture network (Fig. 4). Where data plots close to the *X*–*Y* axis, high interconnectedness is possible but its translation to hydraulic performance is dependent on specific aspects of clustering (i.e. density clustering v. orientation clustering) and the dimensions of fractures.

Simplistically, determining the point of rapid change, i.e. the percolation threshold, against the geometric parameters on the scale of investigation, would be a big step towards optimal fracture network modelling. But without a solid mathematical formulation of the relationship between the components (cf. Manzocchi 2002) and given the uncertainties that are attached to especially the second level characterization, this can only be used at best for a qualitative assessment of the flow in the fracture system when that is calibrated with dynamic indicators.

### Mud losses and wireline methods

Drilling mud losses, PLT inflow zones and temperature logs can potentially link fluid inflow to relative short borehole intervals. They can be used to investigate which part of the fracture network (orientation set, mechanical stratigraphy) is driving production performance. Losses in general give only a qualitative assessment of fracture porosity. Only with careful recording and analysis can a more quantitative assessment of effective aperture of the fracture system be made (Lietard *et al.* 1996; Sanfillippo *et al.* 1997). Monitoring of inflow and outflow mud rates with high-resolution flow meters permits the detection of conductive fractures and potentially an evaluation of the hydraulic fracture width (Verga *et al.* 2000). Loss volumes can provide an indication of the extent of a fracture network beyond the wellbore. Linking losses to fractures still requires integration with other data, since they could also be linked to local matrix properties, e.g. high porosity (cf. Beda & Carugo 2001).

PLTs can be used to analyse the proportional contribution to the total flow of specific wellbore intervals and when matrix conductivity is well constrained, to calculate an estimate of apparent fracture permeability or conductivity (Sullivan *et al.* 2006; Barr *et al.* 2007). The result is generally an interpretation of the flow characteristics of a composite fracture zone, rather than of individual fractures. With thin, high permeability beds present, the distinction between fractures or matrix may be complicated. Modelling the PLT response in terms of matrix–fracture interaction for a variety of fracture scales and fracture and matrix properties (Barr *et al.* 2007) will improve the interpretation. Correlation with mud loss data may improve the interpretation of otherwise ambiguous results (Clifford *et al.* 2005). Temperature logs may perform a similar function to PLTs, where points of different temperature in the drilling fluid along the wellbore could represent open and flowing fractures. Typically, inflowing oil is warmer than the drilling mud (cf. Serra 1986). However, gas expands when it enters the borehole and causes a temperature drop.

Wireline formation testing tools can be used for mini-DSTs. Using the tool straddled by two isolating packers, the formation can be tested in various configurations for a variety of objectives (Whittle *et al.* 2003; Coelho *et al.* 2005). With good isolation, detailed information can be obtained on the flow characteristics of a small interval of the formation (Verga *et al.* 2002). Interpretation methodologies are those used for well tests. The formation pressure data acquired through formation testing tools in the standard way may also provide important information in fractured reservoirs. Observed depletion or increased pressures may indicate the network interconnectedness during production or injection.

### Well tests and production data

Well tests, interference and tracer tests and production data, characterize the flow from a considerable sector of the borehole which can be turned into information on the extent of the fracture network beyond the borehole.

Reservoir performance can be evaluated through pressure transient tests by measuring flow rates and pressures under changing flowing conditions. The interpretation and modelling of well test data in a fractured reservoir is a complex subject in itself (Aguilera 1995) and outside the scope of this paper. The objective of well test analysis is to ascertain how reservoir properties change with distance from the tested well (Doe 1991). The shape of the curve of pressure derivative data plotted against time allows an assessment of that spatial development (Wei 2000). The data is analysed for reservoir permeability thickness (*Kh*), the flow capacity or productivity index (PI) and drilling/completion formation damage. The interpretation of pressure transient tests is non-unique and requires a good knowledge of the underlying geological context such as matrix quality and distribution, the hydrocarbon specifics and some conceptual picture of the fracture network (Wei *et al.* 1998). A simple starting point is the ratio of well test permeability (*Kh*) over matrix permeability (Reiss 1980), termed the fracture productivity Index (FPI). This gives a first indication of the extent of the fracture network beyond the wellbore. Increasingly, the tendency is to construct dedicated dynamic models that incorporate detailed and realistic fracture–matrix geology to model and analyse pressure transient test data (Rawnsley & Wei 2001; Leckenby *et al.* 2007). The main aim is then to establish the key geometric elements and configuration and/or physical processes.

An interference test is a multi-well test where pressure changes v. time are measured in one or more observation wells around a producer or injector. These measurements can be analysed for permeability and storage capacity of the fracture system connecting the wells (Araujo *et al.* 1998; Baker *et al.* 2000). Pulse testing is an extension of interference testing where short production or injection periods alternate with shut-in periods. A tracer test has a similar objective using the travel time of an injected tracer. The timing of arrival of pulses or a tracer in observation wells is indicative of the permeability between wells. The changes in production rates of oil, gas and water can provide information on the extent of the fracture network in the vicinity of a well or in between wells. The fracture network extent around a well can be analysed with techniques based on decline curve analysis and associated derivatives (Agarwal *et al.* 1998; Li & Horne 2005). The extent between wells can be analysed with simple rank correlation techniques but more elaborate techniques using principal component analysis are also applied (cf. Main *et al.* 2007).

### Stress and reservoir performance

The conductivity of fracture systems in the subsurface will be sensitive to the magnitude and orientation of the *in-situ* stress state and possible changes thereof. Jolly *et al.* (2000) refer to experimental work suggesting that fractures normal to the minimum principal stress are most conductive. Alternatively there is also evidence that conductivity might be related to the amount of shear stress acting on the fracture plane (Heffer *et al.* 1999). This is formulated in the concept of critical stress, which postulates that during slip, porosity and permeability are increased in response to dilation. The latter is the consequence of the roughness of the fracture surfaces (Barton *et al.* 1995). Both mechanisms are to some extent interdependent. More pronounced roughness might translate into a higher resistance to shear along the plane as a result of interlocking asperities thereby preventing or reducing fracture closure even if normal stress increases. Jolly *et al.* (2000) studied both mechanisms in generic fracture network models and concluded that in both cases their models were sensitive to stress orientation and changes in magnitude and that those changes would be apparent in pressure transient tests.

In most fractured reservoirs, the (local) *in-situ* stress state will change with depletion and with water injection (Jolly *et al.* 2000). In essence the pressure drawdown leads to an increase in effective stress. This in turn may cause a reduction in effective permeability especially if it is primarily related to the fracture network (Buchsteiner *et al.* 1993). When anisotropy in stress state and in fracture network orientation combine in certain ways it may result in anisotropic changes in the permeability field (Heffer *et al.* 1999). With water injection, reservoir pressures will increase locally and that may lead to the formation or reactivation of fractures and the creation of new flow paths. Stress-related change in dynamic behaviour may be indicated by changes in well test results i.e. a reduction in *Kh* and PI over time and by production rate changes. These effects listed above are not exclusively related to stress changes but may be linked to flow in the reservoir such that water or gas block the flow of oil through relative permeabililty effects, or to mechanical problems in a well.

With depletion, the pressure drop may increase the normal stress acting on a fracture plane and cause elastic closure of the fracture. The effects of a considerable pressure drop may be irreversibly plastic (McDermott & Kolditz 2004). Santarelli *et al.* (1998) presented data to show that the stress path, i.e. the change of effective stress magnitude and of horizontal to vertical stress ratios, is an important factor in explaining reservoir performance. If plastic deformation is a consequence of this stress path, it can have important consequences during re-pressurization since it is essentially irreversible. Various empirical models have been developed to assess the effect of depletion on fracture permeability (Walsh 1981; Barton *et al.* 1985; Buchsteiner *et al.* 1993). Such models have to be idealized. They require quantification of parameters (e.g. fracture width, height, asperity characteristics and rock mechanical parameters) that in most cases are not known in detail. The potential importance of stress changes reinforces the need to analyse the initial and changing stress state in a reservoir.

### Dynamic calibration example

Co-visualization diagrams (Fig. 14; cf. Rawnsley *et al.* 2004) of the three wells discussed earlier show the diversity of the fracture network and its flow response. In well A, mud losses start at 3760 m and keep occurring at regular intervals along the well bore. Loss rates start at 2–4 m^{3}/h and remain relatively stable. Near the top of Unit 3 losses and loss rate increase at a point that coincides with a spike on the temperature log and a peak in the NE fracture intensity log. However, earlier temperature log peaks are not associated with losses. In well B losses start at around 4370 m (3 m^{3}/h) and show a similar pattern as in well A. Interpretation of PLT data shows only one inflow interval just before the point where losses started. The fact that this interval coincides with the only major interval of NE fracture density seems significant. The NE and the NS fractures are developed in well defined separate intervals and do not seem to be part of the same fracture spectrum. In well D losses are similar to wells A & B. Losses start near top Unit 2 (3.5 m^{3}/h) but thereafter rates diminish gradually even while new loss points are recorded. The PLT is interpreted to show inflow in four intervals, all contributing roughly the same proportion to total flow. The first interval is not associated with losses but coincides with the most intense NE-fractures interval. The point where losses started coincides with a break in drilling and a small increase in mud weight when drilling was resumed. The bottom PLT interval, contributing slightly more than the others, is also related to a NE-fracture intensity peak.

Well D has a higher degree of fracture intensity as shown for instance by P_{32} data (Table 1). The inflow is developed over considerable parts of the wellbore. This may be related to the matrix permeability in well D combined with relatively high fracture density. In the other two wells the flow is linked to specific parts of the fracture network, i.e. fracture corridors occurring where intensity peaks. The mud loss record supports such an interpretation. Loss rates of a few m^{3}/h and the frequent occurrence of dynamic losses could be indicative of a fracture network where occasional densely fractured zones, i.e. fracture corridors, link up the wellbore to the larger scale fracture network. The fracture intensity curves in the co-visualization diagrams are always very spiky showing the clustered nature of the fracture sets. Major spikes on the different curves do not as a rule coincide but are located at different points along the wellbore. A system that is density clustered (Fig. 9) should show a high degree of coincidence of the fracture clusters in the orientation ranges. The fact that this is not shown here leads to the conclusion that the system is more orientation clustered.

### Summary

Calibration of fracture network geometry with dynamic indicators aims to define conductive properties and calibrate the spatial development of fracture interconnectedness. The latter can only be qualitatively assessed beyond the wellbore, both vertically and horizontally. The calibration should also be used to investigate if aperture assessments, e.g. those derived from cores, image logs and outcrop, are realistic (Gilman 2003). This should take into account the differences between mechanical and hydraulic aperture, and subsurface stress and stress release when core is taken to the surface. Methodologies to evaluate dynamic indicators are numerous and data acquisition conditions vary. During drilling the conditions (e.g. pressure) in the borehole can be significantly different from those when a PLT is taken under flowing conditions after well clean up. The realization that conditions play a role is important for proper interpretation since different indicators need not show a comparable response over the same wellbore interval. Dynamic calibration is made easy by the integration of information and interpretation in co-visualization displays. Where only orientation and density are known in detail, the analysis helps to formulate constraints for dimensions, aperture and connectivity for the network as a whole as well as specific fracture sets.

## Inter-well interpolation

The geometric analysis and dynamic calibration of a fracture network is a prerequisite for the construction of a model of the fractured reservoir on any scale. In addition, some idea must be formed about fracture network characteristics in the inter-well space. Geostatistics are by now a standard approach for conventional reservoirs (Davis 2002), but it does not seem to be in common use for fracture networks. Curvature of geological surfaces derived from seismic data has been used to assess the degree of fracturing. Seismic anisotropy is used to evaluate spatial attributes of fractures and fracture systems. Geomechanical methods can be employed to derive fracture distributions from existing fault geometries and stress states.

### Geostatistics

Interpolation into the inter-well space can be done through the analysis of fracture network characteristics with geostatistical methods, e.g. variograms and subsequent kriging of the data (Olarewaju *et al.* 1997) or neural network techniques (Ouenes & Hartley 2000). Gauthier *et al.* (2002) use standard statistical methods to derive fracture frequency from well data. Several reservoir characteristics, which potentially can be correlated with fracture frequency, are then analysed statistically to determine a possible field-wide fracture frequency distribution. This approach seems to be the general one. Geological drivers (e.g. structure, lithology, porosity) are linked through standard statistics and/or geostatistics to fracture density derived from well data (Bourbiaux *et al.* 2002), whereas dynamic indicators are used as additional constraints (cf. Ouenes & Hartley 2000). For conventional reservoirs the geostatistical analysis is based on well data combined with a formulation of reservoir architecture based on facies analysis (Davis 2002). Often the reservoir architecture is dominated by one direction and shows relatively long correlation distances. A fracture network is normally defined by multiple directions, certainly in density but also in orientation, and will have much shorter correlation distances. In conventional reservoirs, the horizontal dimensions dominate the problem, whereas for fracture networks it is a three-dimensional problem. To validate the geostatistical approach will therefore be more difficult in the fracture case.

### Curvature

In it simplest form curvature of geological surfaces is defined as the reciprocal of the radius of curvature. The computation involves the calculation of the second derivatives of the surface and the specification of various other geometry and scale-related measures (Roberts 1998; Bergbauer & Pollard 2003). Various calculated measures are used to quantify the degree of deformation (or strain) and subsequently to predict fracture orientation and density. Curvature has been used to assess of the degree of fracturing, e.g. through comparison with well fracture data (Lisle 1994) or production data (Belfield 2000). However, the numerical techniques are not universally accepted and without filtering for different scales of surface irregularities, questionable results may be produced (Bergbauer *et al.* 2003). A significant problem with the technique is that total curvature is related to the total structural history whereas open fractures may only be related to part of that history (Heffer *et al.* 1999; Belfield 2000).

### Seismic anisotropy

In principle, seismic velocity anisotropy can be used to determine spatial attributes of fractures smaller than the seismic measurement resolution (Hudson 1981; Crampin 1981; Armstrong *et al.* 1994; Lynn 2004). In rock with vertically aligned fractures, elastic properties vary in the direction crossing the fracture but not along the fracture plane. The effect is that waves travelling along the fracture travel faster than waves crossing the fracture and waves thus affected arrive at the recording point at different times, giving rise to azimuthal velocity anisotropy. Models built with synthetic, controlled data sets show the potential value of the methodology (Will *et al.* 2005; Freudenreich *et al.* 2006). When real data is used the results are more variable. Gray & Head (2000) found good correspondence between Amplitude Versus Offset (AVO) and Amplitude Versus Azimuth (AVAZ) analysis results and well fracture orientation data, but predicted fracture density was more difficult to link to well performance. They explain this by postulating that in their example, well performance is governed by fracture orientation rather than density. The velocity anisotropy approach has been applied to, among others, the Clair field on the UKCS (Smith & McGarrity 2001; Barr *et al.* 2007).

Fracture density and/or orientation estimates derived from seismic anisotropy must be considered to be relative and qualitative until calibrated properly. For instance the basic response of a dense system of short fractures may not differ much from that of a less dense system of longer fractures. The anisotropy response of a multidirectional fracture system, where several directions and associated apertures are combined, will be very difficult to analyse and verify (Maultzsch *et al.* 2003; Liu *et al.* 2003). The spatial averaging implicit in the seismic method is insensitive to factors, such as connectivity, that are critical to flow. In the use of seismic fracture predictors results must be generated with the key geometric aspects of the fracture network in mind. Since the generation of the anisotropy attribute is creating another model, it must be calibrated against the data, both static and dynamic, that was used to describe the fracture network geometry (Boerner *et al.* 2003; Will *et al.* 2005).

### Geomechanical methods

Geomechanical (e.g. finite element and finite difference) techniques use the estimated current or inferred historic stress state to evaluate the orientation and density of fractures with respect to a fault framework (Heffer *et al.* 1999; Bourne *et al.* 2000). With the stress state or displacement boundary conditions, a strain tensor is generated that serves as an indicator of fracture orientation and density. Multiple deformation events can be modelled sequentially by applying variable boundary conditions, e.g. displacements changing in orientation and magnitude with time. The results depend very much on the accuracy of the input data. Both stress state and rock strength parameters are generally poorly constrained. The resolution of the fault framework is also an important factor that determines output quality. Geomechanical methods can be used to couple stress state changes due to depletion with fluid flow (Heffer *et al.* 1999; Bagheri & Settari 2005). This allows an assessment of conductivity changes in the whole fracture network or a specific fracture orientation.

### Summary

The results obtained using all these techniques to evaluate the inter-well space are variable. Special care should be taken when they are based, implicitly or explicitly, on the use of scale invariant or fractal mathematical relationships. Continuity between small and large scale parameters such as fracture and fault dimensions and density, is possible but should not be assumed. Most methodologies present an indicator of fracture density and possibly orientation, often with limited vertical resolution. Whatever technique is used, results must be calibrated against representative well data. Since, as was discussed earlier, density and orientation are not necessarily the most influential parameters on fracture network conductivity, the inter-well interpolation results should come with a considerable error margin.

## Conceptual models and modelling rules

Given the complexity of fractured reservoirs, it is not realistic to assume that addressing the problem in the sequential order of geometric analysis – dynamic calibration – inter-well interpolation will constitute a fully integrated, optimal analysis. In practice, an iterative process is preferred in which specific aspects are studied in conjunction and used to reinforce sub-processes. The process should focus initially on evaluating those geometric elements (orientation and density) that can deliver hard rules. When combined with dynamic indicators, it is possible to constrain the more elusive rules for fracture dimensions, connectivity and aperture. The process involves the formulation of conceptual models and modelling rules. A conceptual model describes the various components of the fracture network in terms of geometry and dynamic behaviour (Rogers & Al Ansari 2004). The starting point for conceptualizing could be the widely used subdivision of fractured reservoirs (Nelson 2001; cf. Allan & Qing Sun 2003; Fig. 1). However, it is not so simple because most fractured reservoirs can be classified vertically and horizontally in two or even more classes. The standard scheme differentiates primarily on the basis of the matrix characteristics and much less on fracture network characteristics (Fig. 1). The discussion given above on percolation theory shows the importance of the variability in that network.

Classification and conceptual model must take both matrix architecture and fracture network interconnectedness into account. The first is a function of the distributions of net reservoir rock, porosity and permeability. The simple scheme shown in Figure 15 cannot be considered to present a rigorous classification but it allows an evaluation of the fractured reservoir in terms of type and in terms of well performance. Such schemes can be used, after fracture network geometry and dynamic indicators have been integrated, to plot wells relative to each other. Note that it does not address the complexity of the transfer process between matrix and fractures. Nevertheless plots like this show what is observed and hence can be used to make a first assessment of the effectiveness of that transfer process. Modelling rules follow from the analysis of the geometrical aspects of the fracture network. Such rules may be uniform over the whole field but in general the set of rules will show differences laterally and vertically. Together with inter-well interpolation results, the conceptual model and the formulated rules form the basis for subsequent modelling of the fractured reservoir.

### Conceptual model and modelling rules example

The fracture data (P_{32}) shows that well D is the most intensely fractured whereas wells A and B are less fractured (Table 1). The data is interpreted as a superposition of several orientation sets with a considerable degree of orientation clustering (Fig.14). In wells A and B the flow is linked to the relatively high density peaks, i.e. fracture corridors, but in well D the network possibly also contributes as it is linked directly to better quality matrix. The analysis of the wells shows that interconnectedness of the fracture network is variable but it seems unlikely that an effective spanning cluster exists at the scale of the three wells. Rather network clusters will be developed at a scale not exceeding that of the inter-well spacing. Fracture orientation is dominated by NE-oriented sets. Purely from a numerical point of view, and indicated by the dynamic calibration, this orientation seems to be the one primarily associated with flow. Conceptually (Fig. 16), the reservoir can be represented by a relatively dense system of dominant NE-oriented fracture clusters with additional clusters in north–south, east–west and NW orientations. The latter are probably similar in density and dimensions. The length of the NE-oriented clusters and the abutting relationships between NE and the other orientations are uncertainties that must be addressed in fracture network modelling.

Potentially well D is the most promising well in terms of performance, based on fracture density and matrix development. Assuming that fracture corridors control the performance, it can be expected that well D will be connected to a larger portion of the total network compared to the other wells. Well B would be the lowest in that respect. The relative positions of the wells on the typing chart of Figure 15, assume a percolation threshold at the well scale and a matrix that will only be capable of low flow rates in the absence of fractures. From the analysis of the example data the following fracture network modelling rules can be specified:

the fracture network consists of network clusters at a maximum size of the well scale;

the network is dominated by NE fractures (53%); fracture clusters in the north–south (20%), east–west (15%) and NW (12%) orientation ranges will be modelled in lesser proportions;

fracture density (P

_{32}) is structurally controlled, i.e. it is characterized by regional differences;connectivity relationships between fractures have not been established and these are used as a modelling uncertainty parameter. Considering that the dynamic indicator calibration of the three wells favoured the NE fracture orientation for flow, the option that the NE orientation represents through going fractures and other orientations abut against it in variable proportions will be considered;

mechanical stratigraphy is not indicated by the fracture density, which seems to be uniform over the units and in the modelling exercise unit boundaries will not inhibit fracture propagation. Note that this is not necessarily correct since similar densities can also occur when boundaries do inhibit propagation;

as a consequence of rule (e) the vertical dimension of fractures is controlled by the thickness of the modelled units; the horizontal dimension will be controlled by the inferred network cluster size, i.e. the maximum length will not exceed well length and will probably be less; analysis of network clusters in constructed models will have to be used to establish suitable ranges.

No inter-well interpolation results are available other than the large-scale permeability in Unit 2 (Fig. 10). This will be used as a proxy for density distributions for the orientation ranges in the model. This map was derived from reservoir architecture modelling work not discussed further here.

## Fracture network modelling

The objective of fracture network modelling is to study the interconnectedness and hydraulic potential of the network and, when it meets certain criteria, pass the essential elements on to a reservoir simulation model. The latter can be a coarse full field model that is used to manage field development and evaluate reserves. Alternatively, it can be a detailed sector model to analyse specific issues concerning one or a few wells. The choice for simulation methodology ranges from relatively simple matrix models with a permeability multiplier to cater for the fracture network, via dual medium (dual porosity or dual permeability) models to explicitly modelled fractures. What determines the choice must be the capacity of the chosen method to represent the complexity of the real fracture to matrix connection while preserving the essence of the real reservoir (Gilman 2003) as has been captured in the conceptual models and modelling rules. Irrespective of method, size and grid dimensions, a good model is one that captures the important flow mechanisms and geometrical aspects and gives reasonable estimates of past and future performance of the area modelled. Given that objectives may be different and that conceptual models, modelling rules and inter-well constraints carry large uncertainty ranges, there are no simple schemes to guide the modelling of fracture networks.

### Discrete Fracture Network modelling

Discrete Fracture Network (DFN) models are constructed by generating, through a combination of deterministic and stochastic methods, fracture sets with variable dimensions, orientations, density and interaction relationships (Rawnsley *et al.* 2004). The modelling of a fracture network usually takes place in the same geocellular framework as has been used for the matrix (Cacas *et al.* 2001; Bourbiaux *et al.* 2002). Considering that lithology contrasts commonly play a role in development of a natural fracture network, this is a logical choice to begin with.

The first phase of constructing a DFN model deals with translating geometry rules and inter-well interpolation results into a fracture network, which is calibrated against the data derived from wells (Bech *et al.* 2001). The second phase involves analysis of the interconnectedness in the fracture network. The numerical proportions of network clusters (i.e. groups of fractures connected with each other) and number of fractures in a network cluster are indicators of the level of interconnectedness at the scale of the model. With that data, an assessment can be made of how the model relates to the percolation threshold as that was estimated in conceptual terms. The last phase deals with assigning hydraulic properties to the fracture network. The most direct way is to assign an (hydraulic) aperture value to the modelled fractures. In simple models that could be a fixed value for all fractures but can be made more complex by using probability density functions that are sampled differently for different fracture sets. An alternative is to work with fracture porosity maps that are interpolated per grid cell, with aperture assigned proportionally to fractures in the cell (cf. Dershowitz *et al.* 2000). A deterministic approach in which conductive features are identified and assigned values or value ranges (cf. Barr *et al.* 2007) is at the other end of the scale.

Fracture porosity, which is a function of fracture spacing and aperture, is difficult to determine and modelling it properly is a challenge. When a fractured reservoir has reasonable matrix permeability and the matrix provides the storage, the value of the fracture porosity is relatively unimportant in terms of volume. The fracture volumes become important in economic terms when matrix storage is relative low. However, fracture porosity is always important as it determines fracture permeability (Witherspoon *et al.* 1980).

### DFN modelling example

The network modelling method used for the example presented here is based on fracture growth on various scales and orientation, using controlled relationships with unit boundaries and between fracture sets. The constraints are combined in a series of probability maps or grids that contain a numerical definition of specific parameters. The method is fully explained by Swaby & Rawnsley (1997) and Rawnsley *et al.* (2004). For the example, fracture sets have been constructed that can be combined in various ways to form a model. Four sets of NE clusters, with maximum lengths 2000 m, 1500 m, 1000 m and 500 m respectively, are constructed using the permeability distribution map (Fig. 10) to control density. High permeability is assumed to equate to high density. Fracture sets for the north–south, NW and east–west orientations were constructed with length distributions constrained by an impedance factor attached to the NE-2000 set. The impedance value indicates the fraction of intersecting fractures for which growth is stopped against the impeding set or unit boundary. For the model these constraints were: (1) fully constrained with impedance 1.0 and (2), (3) and (4) with impedance values of respectively 0.75, 0.5 & 0.25. Density constraints are chosen to ensure that, rather than a full fracture network, all constructed sets model fracture clusters. In this way sixteen fracture sets were generated: four NE sets and twelve (4×3) orientation fracture sets. For the latter the same permeability map was used as a constraint but to account for the observed lower densities, the probability of growing a fracture was half that of the NE sets.

Models created by combining specific fracture sets were investigated for network cluster development. In Table 3 the results of such an analysis for the combinations of NE sets with the other orientation sets are shown for specific impedance values. The number of network clusters (i.e. 25 or more fractures connected) and the percentage of the total number of fractures in the largest cluster, characterize a model. For low NE impedance (0.25), one single cluster makes up almost the total network in all models. Such models will be well above a percolation threshold on the scale of the model and are not considered to reflect the well results. In contrast, for the total impedance case (1.0), the models, with the exception of that with the NE-2000 fracture set, show none or only one cluster that contains more than 25 fractures. These models are well below the percolation threshold on the scale of the model and are also considered invalid. For two of the models the network cluster data is shown in Figure 17. Model #01 (NE-2000 combined with impedance 1.0 orientation sets) shows clusters associated with the well where the NE fractures dominate totally. Proportionally the well A cluster is much larger than that linked to well D. In model #09 (NE-1000 combined with impedance 0.5 orientation sets) the NE fractures are still dominant but the length distribution contrast with the other clusters is much smaller. Proportionally, network clusters linked to wells A and D are balanced in this model and this reflects better what was expected in conceptual terms.

The size of the well A cluster in both models is a consequence of the relatively high matrix permeability in the vicinity of well A (Fig. 10). Density in the generated fracture sets is linked to the size of this area. On the basis of dynamic well analysis a larger network cluster was expected for well D. This could mean that the modelled permeability distribution is a poor proxy for fracture density. On the other hand, differences in regional overprinting relationships, between one or more of the orientation sets with respect to the NE set, may also be a potential factor. Such regional differences could be related to the orientation-clustered nature of the network. It would require more complex models with local impedance differences to investigate this further. The two models are also compared in terms of fracture porosity, which was calculated here in a 50 m×50 m grid model (Fig. 18). For the same aperture assigned to both NE fracture sets and the orientation-clustered fracture sets, the differences between the two models are considerable. Model #01 again shows the dominance of the NE direction, but in model #09 the other orientations are also important. The effects on flow of the various configurations can be analysed dynamically with these fracture network models, e.g. using well test results as matching criteria (cf. Rawnsley & Wei 2001).

The fracture sets generated are very simple, i.e. they have a single orientation and only contain vertical fracture clusters. Given that the models are non-unique, the choice of simplicity was deliberate in order to generate representative sets of models that capture the identified uncertainties. Here, where the focus is on the length of the dominant fracture orientation and its interaction with the other orientations that can be adequately represented with unidirectional sets of vertical fractures. Compared with multidirectional, variable dip sets, these simple sets will elucidate greater contrasts between the different models. Once a preferred model has been chosen, multidirectional, variable dip sets can be generated to achieve greater refinement (cf. Leckenby *et al.* 2007).

### The step to simulation modelling

The basic equation that describes flow in fractures is the cubic law which implies that fracture conductivity for an ideal open fracture, i.e. for a fracture between parallel plates that are not in contact with each other, and assuming laminar flow, is uniquely defined by aperture (Witherspoon *et al.* 1980). The natural deviation from a situation of parallel plates causes a reduction in flow, which for example, can be accounted for by the Joint Roughness Coefficient, when the effective fracture permeability is calculated (Bakhtar *et al.* 1985). Provided that a reasonable estimate of JRC can be made, the impact of the adjusted aperture on fracture permeability *K*_{F} can be evaluated. Using the cubic law without aperture correction implies that the calculated fracture conductivity will represent an absolute upper limit.

The effective permeability of a fracture network depends on the geometry of that network and the conductivity of individual fractures. Using the cubic law the (intrinsic) permeability of a single fracture is related to the square of the aperture (*L*^{2}); fracture conductivity is equal to the permeability multiplied by aperture (*L*^{3}). Schemes to compute the effective fracture permeability range from relative simple vectorial decomposition schemes of the fracture permeability in *X*, *Y* and *Z* directions of the chosen grid, to full-scale simulation methods. Oda's method (Dershowitz *et al.* 2000) determines an effective permeability tensor by integration of fracture conductivity of fracture sets of specific orientation or orientation ranges. It is a quick method but it does not incorporate the effects of fracture size and connection to other fractures. In principle, therefore, it can only be used for well-connected fracture networks (Dershowitz *et al.* 2000; Will *et al.* 2005). Fracture permeability simulation methods derive a directional permeability by applying flow and no-flow boundary conditions to the grid and solving for specific orientations (Dershowitz *et al.* 2000). Bourbiaux *et al.* (1999) describe methods based on the discrete modelling of fractures and matrix in what is called the joint element method. The simulation methods incorporate the connectivity aspects that are neglected in the grid cell solution schemes, but are sensitive to the assumptions made for boundary conditions. The simpler schemes are often chosen because explicit simulation methods are computationally intensive.

The above methodologies involve upscaling from the geological scale, i.e. the geometric and hydraulic characteristics of the fracture network, to a representative formulation of fracture permeability and the fracture network geometry on the reservoir simulation grid scale (Bourbiaux *et al.* 1999). The characteristics of a fracture network are potentially measured at a metre scale or lower whereas the grid dimensions are tens of metres or higher. Transferring the DFN to a grid expresses the problem in terms of averages of the small-scale quantities (Gilman 2003) and effectively the reservoir model becomes a simplified representation of the actual fractured reservoir (Gorell & Bassett 2001). To represent the detailed geology adequately requires that the consequences of upscaling are properly thought through. The dimensions and the interconnectedness of the fracture network and the flow interaction between the fracture network and the matrix must be captured at an appropriate scale, e.g. the representative elementary volume (REV) or the minimum volume for which a physical property can be measured and still be applied as an average value over a larger area (Trice 1999; Dershowitz *et al.* 2000; Bourbiaux *et al.* 2002). The REV approach assumes that a representative elemental volume exists at a scale of the chosen grid cell size. In a fracture network where the scale of interconnectedness varies widely that may not be achieved easily (Dershowitz *et al.* 2000).

The considerations that determine the choice of grid dimensions are related to those that determine the choice of simulation model type and will be discussed below. The objective could be to model the full field, represented by relatively large-scale features such as faults and fracture corridors, or at the other end of the scale to analyse local well test data (Rawnsley & Wei 2001; Gorell & Bassett 2001; Bourbiaux *et al.* 2002). This is not an either–or decision. Small-scale sector models may be indispensable in reducing the uncertainty in concepts and models. Detailed analysis can bring out the key components of the system that must be represented in the larger scale reservoir simulation model (Rawnsley & Wei 2001; Araujo *et al.* 2004). The capacity of computers to handle models of a certain size will be one limitation. If, because of resource constraints, grid dimensions become larger than fracture dimensions, a different modelling strategy or sector modelling to derive ‘pseudo-properties’, which compensate for deficiencies introduced during upscaling, may be required (Gilman 2003).

In fractured reservoirs, the fluid transfer from matrix to fracture commonly occurs at rates that are considerably lower than the flow within the fractures (Kazemi *et al.* 1992; Daly & Müller 2004). When that is the case the permeability of the fracture network and of the matrix must be represented separately and additional terms to govern the flow from matrix to fracture must be included. Bourbiaux *et al.* (2002) summarized the requirements for various simulation methodologies. A single medium (or single porosity) model can be used for dense, well-connected fracture networks where the exchange between the small matrix blocks and the dense fracture system is virtually instantaneous. At the scale of the grid model, the problem has become a homogeneous one (cf. Waldren & Corrigan 1985). A dual porosity–single permeability model is required where the kinematics of matrix to fracture transfer are such that the matrix is only the source of fluids. If matrix to matrix flow has to be modelled, dual porosity–dual permeability representation is required. In the dual porosity–single permeability model flow takes place in the fractures and between matrix and fractures. Matrix blocks are assumed to be completely surrounded by fractures. In dual porosity–dual permeability models, matrix to matrix and fracture to matrix flow is also possible. This representation is more geologically realistic (Gilman 2003) but computationally intensive. In the subsurface, matrix blocks will not be floating but must support stresses through matrix block contacts (Baker & Kuppe 2000). Finally, when the dimensions of the fracture network make them major conduits for fluid flow at a scale well beyond that of the matrix blocks, explicit modelling, where fractures are represented directly by specific grid cells, may be required. Computational requirements mean that the simple representations are often preferred.

The separate representation of fracture flow means that the model must account for the different contribution of matrix and fractures and their interaction. The basic idealized formulation to describe this has been proposed by Warren & Root (1983). The matrix to fracture flow is described by a transfer function in which the shape factor accounts for the shape and dimensions of matrix blocks. This formulation has been extended in numerous ways to account for multi-phase flow and for different types of recovery mechanisms (cf. Waldren & Corrigan 1985; Kazemi *et al.* 1992; Bourbiaux *et al.* 1999 & 2002; If & Frykman 2005). The fluid exchange between matrix and fractures depends on the pressure differences between matrix and fractures and viscous, capillary and gravity effects. Which one of these plays the greater role depends on various factors (Waldren & Corrigan 1985) and the formulation of the shape factor is not without problems. For the various flow mechanisms, different values of the shape factor are required (Bourbiaux *et al.* 1999). Given that there is an obvious sensitivity in the shape factor computation with respect to the grid dimension, it is considered to be a history matching parameter (Gilman 2003). In simple terms, dynamic modelling requires that the permeability distribution of the fracture network and of the matrix will be combined in a meaningful way. The rules and concepts can be used to construct a digital fracture network model, which in turn can be used to evaluate and calculate the effective fracture permeability field. The results of the integrated analysis of static and dynamic data captured in rules and conceptual model are the basis for a considered choice with respect to grid dimensions, simulation model type and transfer function, to ascertain that the constructed model is appropriate for the intended analysis.

## Concluding remarks

Proper modelling of fractured reservoirs to analyse and predict performance is essential to implement the right development options thereby assuring optimal recovery. Available data from cores and image logs, even when available in abundance, cannot provide a description of the geometry of a fracture network that is enough to construct the definitive fracture model. Cores and image logs simply cannot provide sufficient information on fracture dimensions, aperture and connectivity. This means that considerable uncertainty has to attach to these parameters when the data is used for fracture network modelling. The construction of a fracture network model is therefore not a simple matter of combining concepts and rules but requires analysis in its own right. Even with a set of relatively restrictive rules, as used in the example described here, several equi-probable DFN models can be generated. Dynamic modelling can be carried out by capturing the total range of possibilities through scenario modelling (i.e. using the same model with changes in parameter values), or through the building of several models addressing different static and/or dynamic concepts. Analysis of the fracture network models against the background of concepts and rules, and revisiting those when appropriate, is essential to constrain parameter ranges and scenarios.

Complex digital models of fractured reservoirs are computationally time consuming and thus expensive to run and difficult to interpret and adjust. There is thus an obvious incentive to construct simple models that are based on the dominant processes and geometries that control performance. There are no hard rules to govern the choice of modelling strategy and methodology, other than a requirement that it reflects the contrast between fractures and matrix as well as serving the reservoir simulation objectives. In that choice fit-for-purposes-simplicity is generally better then complex modelling for the sake of complexity. The fact that we can build such models, having powerful computers and versatile modelling software, does not mean that we must.

Fractured reservoirs are difficult to analyse, model and manage and there are no universally applicable ‘rules of thumb’. Methodologies that work for one fractured reservoir may fail totally in others. But by using workflows in which the analysis of data is integrated correctly, the understanding of a fractured reservoir can be greatly improved and useful conceptual, static and dynamic digital models can be constructed. The analysis of the fracture network in terms of geometry and dynamics and its translation into concepts and modelling rules, which are then tried and tested in DFN mode, is a key element in this workflow. New methodologies to improve data acquisition and analysis, some of which are discussed in this volume, will continue to improve this workflow.

## Acknowledgments

This paper is published with the permission of Shell International Exploration & Production. The critical comments of D. Barr, P. Frykman and R. Jolly are gratefully acknowledged. The views of numerous Shell colleagues I have worked with over the years influenced my own views on fractured reservoirs and are reflected in this paper. The opinions expressed, however, are the sole responsibility of the author.

- © The Geological Society of London 2007