## ABSTRACT

To optimise recovery in naturally fractured reservoirs, the field-scale distribution of fracture properties must be understood and quantified. We present a method to systematically predict the spatial distribution of natural fractures related to faulting and their effect on flow simulations. This approach yields field-scale models for the geometry and permeability of connected fracture networks. These are calibrated by geological, well test and field production data to constrain the distributions of fractures within the inter-well space. First, we calculate the stress distribution at the time of fracturing using the present-day structural reservoir geometry. This calculation is based on a geomechanical model of rock deformation that represents faults as frictionless surfaces within an isotropic homogeneous linear elastic medium. Second, the calculated stress field is used to govern the simulated growth of fracture networks. Finally, the fractures are upscaled dynamically by simulating flow through the discrete fracture network per grid block, enabling field-scale multi-phase reservoir simulation. Uncertainties associated with these predictions are considerably reduced as the model is constrained and validated by seismic, borehole, well test and production data. This approach is able to predict physically and geologically realistic fracture networks. Its successful application to outcrops and reservoirs demonstrates that there is a high degree of predictability in the properties of natural fracture networks. In cases of limited data, field-wide heterogeneity in fracture permeability can be modelled without the need for field-wide well coverage.

## INTRODUCTION

Natural fracture systems can have a dramatic impact on reservoir performance—they may act as highly permeable flow conduits or as baffles and seals. The complexity of a fracture network typically leads to an extremely heterogeneous and anisotropic permeability distribution within the reservoir. Successful management of these reservoirs is impossible without substantial knowledge of the natural tensile and shear fracture systems. For instance, it is essential to know their spatial distribution and hydraulic properties on an inter-well scale to properly simulate the field-wide recovery processes.

This paper presents a new method for predicting natural fracture distributions related to faulting and their effect on reservoir simulations (Figure 1). We use geomechanical models of rock deformation to compute the field-scale distribution of stress responsible for fracturing based on the observed structural geometry of the field. In so doing, we neglect any fractures formed prior to the evolution of field-scale structure, which is reasonable only if these fractures were subsequently cemented or mechanically closed. Fracture network geometries are obtained by simulating the initiation, growth, and termination of individual fractures within the calculated stress field. The resulting network geometries are partially constrained and validated by the available data. Predicted fracture orientations are compared statistically to those observed from core, borehole image and outcrop data. Predicted fracture connectivity is checked against mud loss and production logging tool data. Thereafter, multi-phase, well-scale or field-scale flow simulations of the fracture model are compared to well test and production data. Fault geometries, rock strength, and the remote stress load are all revised until no significant discrepancies between the observations and the model remain.

In the next section, we explain in more detail the geomechanical method employed to predict stress and fractures, the simplifications made, and the assumptions required. We will focus on hydraulically conductive fractures formed in relation to faulting. Faulting is typically the primary mechanism for strain accommodation in brittle rocks and therefore a simple mechanical model of faulting is expected to be sufficient for fracture prediction in many structural settings.

Hydraulically conductive fractures tend to affect fluid flow severely when they connect to form percolating networks. Once conductive fractures form a single flow path through the matrix, they change the flow pattern significantly—generally more so than sealing fractures, which must first close all holes in the barrier formed by their network, before they have a similar impact on permeability. As such, we follow Thompson’s (2000) concept of effective fracture systems.

The following section describes our method for upscaling and simulating flow through fractured matrix using the predicted fracture networks. In this section, we also discuss how flow simulation and fracture prediction are combined to constrain each other. We finish by presenting one example of a recent application of this method to a real carbonate reservoir where the dominant mode of fracturing is hydraulically conductive. This example aims to demonstrate the value of close integration between geomechanical fracture prediction and flow simulation. We will attempt to show how the use of both static and flow data to constrain a single field-scale model of fracture geometry and permeability, significantly reduces the uncertainty in both in the near- and inter-well regions of the reservoir.

## FRACTURE PREDICTION

It is neither possible nor generally necessary to accurately predict individual fractures within a reservoir. Rather, we restrict our attention to predicting the properties of just those natural fracture networks that are hydraulically conductive. This is achieved by describing the continuum stress field responsible for reservoir fracturing using geomechanics.

### Continuum Mechanical Models

We consider the deformation to be distributed smoothly on length-scales of interest. In this way, we can combine knowledge of the regional stress history and the mean rock mechanical properties to compute a smooth field of stress. Large discontinuities in the deformation field, such as seismically visible faults must be represented in the model explicitly rather than as part of the continuum.

The simplest approach to determine the stress field within a faulted reservoir is to assume that the rocks behave as a homogeneous, isotropic, and linear-elastic material, and the faults as surfaces free of shear stress. As a consequence, layers of different mechanical properties, and mechanical anisotropy due to pre-existing weaknesses from earlier deformations are not represented. In this model, large-scale mechanical heterogeneity is represented only in the form of three-dimensional fault geometries.

The distribution of elastic stress related to faulting is governed by the distribution of slip over the fault network. We choose not to rely on fault throws observed from seismic data since they represent only one component of the fault displacement vector. Instead, we calculate the distribution of slip vectors by loading the network of frictionless faults according to the remote stress that caused the faults to slip. The orientation of this remote stress is estimated from the regional geological history, and the magnitudes according to the mean rock strength prior to faulting. A unique pattern of fault slip exists which releases all shear stress resolved onto the faults by the remote stress and other nearby faults. Such a slip distribution represents the static equilibrium of faults that possess zero shear strength. This is a reasonable approximation given that the residual shear strength of faults is typically at least ten times less than the strength of rock between the faults.

Numerical solutions for the three-dimensional elastic stress tensor field are obtained using a boundary element analysis called Poly3D (Thomas, 1993). Poly3D solves the equations of linear elasticity by representing fault surfaces with a series of triangular elements, each of constant slip (e.g. Comninou and Dunders, 1975; Jeyakumaran et al., 1992). This approximation permits solutions for complicated three-dimensional fault geometries and slip distributions by using a large number of triangles to represent the often irregular fault geometries found in seismic data. This capability is obtained at the expense of simplifying the mechanical properties of the intervening rock, which must be modelled as isotropic and homogeneous. Other numerical methods that use a discretisation of both faults and rocks, (e.g. finite element) can avoid this simplification, but are at present unable to readily represent the required degree of complexity in fault geometry. We proceed under the hypothesis that faults, and not the surrounding rocks, represent the dominant source of mechanical heterogeneity that governs field-scale fracture distribution.

Previous applications of elastic faulting models to geological problems include: fault interaction (e.g. Bürgmann et al., 1994; Willemse et al., 1997), fault growth (e.g. Segal and Pollard, 1980; Cowie and Scholtz, 1992; Scholtz et al., 1993; Cowie and Shipton, 1998), fault linkage (Aydin and Scholtz, 1990; Crider and Pollard, 1998), and fault-related fracturing (e.g. Pollard and Segal, 1987; Martel and Boger, 1998; S.J. Bourne and E.J.M. Willemse, 2000, written communication).

Figure 2 illustrates the distribution of elastic stress calculated for a simple hypothetical network of right-lateral strike-slip faults. Notice the characteristic mean stress reduction in areas just to one side of the fault tips, for example, at (–400, 1,400), where rocks have been stretched by moving along the fault plane away from the fault tip. This causes tensile fractures that propagate at high angles away from the fault in the direction of greatest compression (compare with Figure 3). Elsewhere the stress field is more complex and defies simple intuition. This is a consequence of mechanical interaction between faults so closely spaced that their stress fields overlap, for example, at (300, 1,100).

### Brittle-elastic Failure Models

The fault-related stress fields computed by the boundary element method are wholly elastic, and therefore no allowance is made for stress relaxation that follows secondary fracture formation around faults. This approximation remains reasonable if the secondary fractures are small in size compared to the primary faults. In this case, individual fractures will induce only small perturbations on the elastic stress field, which can be neglected. No attempt is made to predict large-scale fault growth since this would significantly effect the computed field-scale stress field and present a significantly non-linear problem. The elasticity assumption ultimately breaks down once significant portions of the reservoir have undergone bulk brittle failure (e.g. > ~50%) since this introduces variations in effective elastic stiffness over length scales comparable with the faults.

We determine brittle failure states within the reservoir according to the Griffith stress criterion for tensile failure (Griffith, 1921), and a Coulomb stress criterion for shear failure (Coulomb, 1773). These criteria characterise rock strength in terms of two material properties, cohesion and internal friction, both of which we assume to be constant across the field within the mechanical layer in question. Whichever criterion is met first as the stress field changes, determines the mode of first failure (Figure 4a). Figure 4 is an example of the distribution and type of brittle failure associated with the elastic stress field around the right-lateral strike-slip faults in Figure 2. Shear failure is predicted to localise around the fault tips whereas tensile failure occurs within the dilatational regions alongside faults.

The lateral extent of fractured regions depends on the rock strength relative to the remote driving stress: increased rock strength results in fractured areas shrinking towards the faults, and vice versa. We include fracture formation mechanisms such as fluid pressure increase, diagenesis, cooling, and erosion of the overburden as body forces, which produce an isotropic reduction in compressive stress within the calculations. This has the effect of shifting all the stress states plotted in Figure 4(a) towards the left and consequently increases the lateral extent of both tensile and shear failure.

The lateral extent of fractured areas is governed by: (1) rock strength relative to the remote stress, and (2) the magnitude of body forces, both of which are uncertain. Nonetheless, the overall shapes of fracture clusters are far less uncertain as they correspond to contours of the elastic stress field that depend only on fault geometry and the orientation of the principal remote stress, both of which are known. Consequently, a small number of wells distributed across the field can provide sufficient data to validate the elastic stress field and constrain the lateral extent of fracturing.

### Discrete Fracture Networks

Once the distribution, mode, and orientation of failure has been determined, the actual geometry of fracture networks responsible for permeability can be simulated by growing fractures with in the three-dimensional, brittle-elastic stress field. During growth, a forbidden zone around each fracture controls fracture spacing and interaction—this represents an overall reduction in local stress due to the presence of the growing fracture. As a result, other fractures are less likely to nucleate within these zones, which limits the ultimate spacing of fractures. Furthermore, growth rules are included to simulate fracture hooking that arises when the propagating tips of two fractures approach each other. Their mechanical interaction leads to the fractures turning towards each other and connecting (Cruikshank et al.,1991; Olson and Pollard, 1990; Thomas and Pollard, 1993). This is an important mechanism affecting fluid flow as it produces connectivity between essentially parallel fractures.

Figure 5 shows an example of a discrete tensile fracture network grown within a single mechanical layer according to the stress field around the faults. Tensile fractures initiate in areas of tensile failure (Figure 4) and propagate in directions perpendicular to the orientation of least compressive stress in Figure 2(b). Away from faults, fractures propagate in the direction of regional stress, but close to the faults fractures grow along curved trajectories to become either parallel or perpendicular to the faults (compare with Figure 3). This is a consequence of shear stress release on the faults and gives rise to variation in the direction of maximum fracture permeability close to faults.

The overall method described above for calculating elastic stress fields to predict the type and spatial distribution of brittle failure has been validated against an outcrop at Nash Point, UK (S.J. Bourne and E.J.M. Willemse, 2000, written communication) where both faults and tensile fractures have been mapped in detail over an 80m by 100m ar ea (Petit et al., in press).

## FLOW SIMULATION

Static fracture models can be calibrated and constrained locally at the well bore, with core data, image logs and mud losses. Calibration on a sector scale or field scale must be obtained through dynamic data, such as well tests and history matching, as described below.

### Well-test Simulation

Flow simulation of parts of the network without upscaling allows early selection of fracture models based on a comparison with well-test data. Comparison of simulated well tests with actual results reduces uncertainty in fracture permeability and their lateral extent in the vicinity of the well.

To do this, we cut out regions of the fracture network model and simulate these without upscaling using MaficOil (1998). MaficOil is a finite-element simulator for single-phase flow through a three-dimensional fracture network. The matrix is discretised by three-dimensional tetrahedral elements and a parallel plate model at the faces of these elements represents discrete fractures. The fracture aperture and compressibility are inferred from geological data and then constrained by flow data. Previous applications of finite-element flow simulation to infer natural fracture properties from well tests include Wei et. al. (1998) and Wei (2000).

Figure 6 shows the pressure response over time of a fracture network to a producing well located in the middle of the fractures. The response of the well can be compared to actual data within standard well-test evaluation programs, as in Figure 6(d). The simulated response shows a pronounced reaction to all features of the fracture network; such as, the anisotropy of individual fractures, the transition from fracture network to unfractured rock, and the diffusion of the pressure front into nearby fracture clusters. In real well tests, well-bore storage may obscure the explicit response but, if the character of the simulation and field data is the same, well-test simulations constrain the fracture model effectively on a sector scale.

### Field-scale Simulation

Once the fracture models are constrained locally via explicit, single-phase simulation, they must be upscaled to a dual-porosity simulator to allow field-scale, multi-phase flow simulation. Dual-porosity (or dual-permeability) simulators represent reservoir properties on much coarser continuum grids, typically having lateral grid block dimensions of the order of 100 m. This scale allows them to describe much larger sections of the field than discrete fracture-network simulators. In addition, they model multi-phase flow processes such as water imbibition or gas-oil gravity drainage, using transfer functions to couple the matrix and the fracture grids.

To upscale the hydraulic properties of a fracture network to a dual-porosity grid, we use MaficOil, the program also used to simulate well tests. It computes the equivalent permeability of a fracture network within each grid block by flow simulation through the network inside the block. The boundary conditions it applies depend on the capabilities of the simulator receiving the resulting fracture grid. If this simulator only uses the components parallel to the co-ordinate axes of its finite difference grid, the faces of the grid block whose normal vector is perpendicular to the flow direction are no-flow boundaries. For tensor simulators, periodic boundary conditions are applied (Durlofsky, 1991), whereas uniform boundary conditions (Long et al., 1982) serve as an approximate solution for both simulator types. The effects and consistency of these boundary conditions are discussed by, for example, Renard and de Marsily (1997) and Jackson et al. (2000).

In some single-phase reservoirs it is possible to treat the combined system of fractures and matrix as one effective medium, if the fracture spacing is small enough. In that case, fractures enhance matrix permeability locally, but full dual-permeability effects are absent. In general, the permeability of such a combined system is not simply the sum of the effective fracture and matrix permeability (Kamath et al. 1998). Consequently, the dynamic upscaling route accounts for this effective medium by upscaling the combined system of matrix and fractures.

In addition to the effective permeability, several other static parameters are computed. These are the number of matrix blocks in each grid block, the imbibition length, the angle of the permeability eigenvectors with the co-ordinate axes, and the like. These parameters are used to check alignment of the simulator grid and fracture system, and to calculate multi-phase transfer functions between matrix and fractures.

The upscaled fracture grid is then used to obtain a field-scale history match. By simulating the historical production of one fluid from the field, and monitoring how well the simulated associated production data fit historical data, the quality of the fracture model is assessed. If a model based on one particular set of rock properties and remote stress state does not match historical data, even if other reservoir properties are varied within their uncertainty range, then this model is discarded as invalid. The requirement that field-wide models match production data from all wells, makes this route a powerful filter.

In our simple multi-phase example with one well (Figure 7), dual-permeability effects influence the simulated interaction between water and oil. The fracture system connects a flank aquifer with a producer in the middle of the model, but the advancing water front is slowed by water imbibition into the matrix. When the water cone in the fractures finally breaks through to the perforation, oil production must be cut back.

## APPLICATION TO A PRODUCING CARBONATE RESERVOIR

As an example, results are shown for a tight, fractured carbonate oil reservoir in the Shell Group that has an average matrix porosity less than 3 percent and an average matrix permeability less than 0.01 mD. The drive mechanisms for the undersaturated oil in this over-pressured reservoir are compaction and natural depletion. An economically successful well must target a laterally extensive network of connected fractures to increase both the initial production rates and the ultimate recovery (UR). Figure 8 shows a perspective view of the field. To help optimise well planning and reservoir development, the field-scale lateral distribution and hydraulic properties of the fracture systems were predicted using the method previously described (see Figure 1).

### The Stress and Fracture Model

Fracture observations from core and outcrop analogues suggest a fault-related fracturing mechanism for the reservoir. Hence, we interpreted fault geometries from 3-D seismic data and computed elastic stress fields using the boundary element technique described previously. The orientation of remote palaeo-stress responsible for fault slip was determined by reference to regional geological data. Brittle failure analysis of this stress field yields a lateral probability distribution for tensile fracturing as shown in Figure 9. These probabilities are derived under the assumption that rock cohesion is randomly variable within a grid block by a factor of 10 percent. The resulting fracture model is consistent with both fracture orientations and relative intensities as measured from image logs for 14 wells.

As in the example shown in Figures 2 and 4, this fracture map is significantly different to the common view that fault-related fractures are distributed as a simple function of distance to the nearest fault. Instead, the predicted fracture distribution shows large areas adjacent to faults that are essentially not fractured. Areas of significant fracturing are bound by faults but they tend to extend between rather than along them.

### The History Match

The fracture models that were grown according to the tensile fracture probability map, were upscaled dynamically and used in a field-wide flow simulation of the production history of the field. The dynamic validation of fracture models against production data from 10 wells proved to be a strong validation of the approach.

In the simulation, wells were constrained by historical flowing bottom-hole pressure data, whereas stimulations and mechanical integrity problems were modelled as a skin change. Figure 10 shows the simulated oil production rate versus cumulative oil produced, and the corresponding historical data. All plots have the same vertical scale but the horizontal axis is different. The plots show that there are two classes of wells in the field: (i) wells that drain fracture clusters (top row) with high rates and an ultimate recovery that is effectively determined by fracture cluster size; and (ii) wells located in unfractured regions (bottom row) with low rates, whose UR is determined by matrix properties alone. The middle row shows wells in fracture clusters that were stimulated, or that experienced integrity problems.

Our geomechanical fracture model, based on the fault interpretation, matches the (decline) rates in 9 out of the 10 producing wells, as well as pressure readings in observation wells, without any local parameter adjustment. The tenth well can be explained by refining the grid for both stress and flow calculations around this well. The extra resolution appears to be necessary to properly represent small-scale yet large stress perturbations associated with the junction of several nearby faults.

History matching a high proportion of wells greatly reduces the uncertainty in the fracture model away from the currently producing area. Indeed, several re-activations of old, shut-in wells since this work have confirmed the current model. Hence, future well locations can be selected with much more confidence.

### Sensitivity Analysis

The impact of uncertainty in the orientation of remote stress, fault geometry and rock strength were all assessed using several different scenarios for the predicted fracture network geometry. The model described above is the only one, out of the many considered, that fits all static data and gives a dynamic rate and pressure match. Changing either the rock strength or the body force magnitude, which affect the mode and lateral extent of fracturing, or the remote stress, which affects the shape of fractured areas, by more than 10 percent results in no more than 4 out of 10 wells showing a history match.

The results of these sensitivity studies demonstrate that both the location and lateral extent of the predicted fractured regions are required to explain both the static and dynamic data. A random distribution of fractured regions, or one simply related to fault proximity, does not suffice.

### The Value

Being able to predict fracture distribution and hydraulic properties on a field-scale has a clear business impact, as it significantly increases the probability of success (POS) of a new well. Wells located completely at random would have a probability of encountering a connected fracture system of just 0.3; wells at random near faults would have a POS of less than 0.5. However, field-wide fracture prediction increases the probability of locating economic fracture networks to between 0.8 and 0.9.

Well re-completions have confirmed the predictive capabilities of this fracture model. Production rates from all wells are, to date, consistent with the flow forecast from the fracture model.

## DISCUSSION AND CONCLUSIONS

In short, rocks obey physics, and physics can be used to predict the fractures that affect flow on the scale of naturally fractured reservoirs.

For the most part, heterogeneous fault and fracture models have been obtained as discrete fracture networks based on stochastic distributions, such as power-laws (Lapointe, 1988; Gauthier and Lake, 1993). These rely on stochastic realisations of the large numbers of fracture networks calibrated to borehole, and inflow data. This approach is ultimately limited to near-well scales as it lacks information on the variation of fracture statistics away from wells. A field-scale stochastic fracture model would require a large number of evenly distributed wells to allow simple interpolation of fracture statistics between the wells.

Other methods, based on a curvature analysis of seismically observed mechanical layers do offer fracture prediction within the inter-well space (e.g., Harris et al., 1960; Ramsey, 1967; Lisle, 1994; Stewart and Podolski, 1998). Several studies report a relationship between curvature and fracturing (e.g., Gorham et al., 1979; Woodward 1984; Gerla, 1987), and some observed that the best areas of productivity within a hydrocarbon reservoir correspond to areas of maximum curvature of the reservoir strata (e.g., Murray, 1968; Belfield, 2000). Nonetheless, application of curvature methods to faulted reservoirs can be misleading since faulting readily induces large horizontal strains in the surrounding rock that have little or no effect on surface contours yet cause fracturing. Heffer et al., (1999) describe a new geostatistical method that combines strain information derived from surface contour and borehole fracture data to predict fractures. This offers the advantage of being tightly constrained by the known data but, across the inter-well space, predictions will be primarily reliant on the shape of surface contours much like the explicit curvature methods.

The fracture model presented here uses a geomechanical model of faulting to predict the field-scale distribution of fault-related fractures that enhance flow within reservoir simulations. The geometry of seismically observed faults are used to compute the full-stress tensor in three-dimensions responsible for brittle failure within the reservoir using a boundary-element numerical method. Fault network geometry can be represented with arbitrary complexity, but rock mechanical properties must be represented as a homogeneous, isotropic and linear elastic medium or half-space. From the computed field of stress we determine the likely geometry and distribution of fracture networks that can be determined from models of fracture initiation, growth, and interaction. The predicted lateral extent and permeability of fracture clusters are sensitive to uncertainties in both rock strength and hydraulic fracture aperture, but using well-test data and simulation of a multi-phase production history significantly reduces these.

By modelling the physical processes responsible for fractures and flow, more meaningful and realistic fracture systems can be predicted. Moreover, uncertainty can be minimised by integrating all the available static and dynamic data. As the model parameters of rock strength, remote stress and the like are field-scale, information from each well constrains the whole fracture model and not just the areas close to wells. This makes the model suitable for fracture prediction and flow forecasting in all parts of the reservoir and not just those parts around existing wells.

Predictive field-scale fracture models of subsurface natural fracture networks enables improved field development through: (1) better assessment of the recovery mechanism; (2) more reliable production forecasts; (3) well placement for optimal drainage; (4) minimal water-cut; and (5) recognition of drilling hazards associated with fractures.

## ACKNOWLEDGMENTS

This work benefited greatly through the involvement and ideas of Jean Borgomano, Tony Cortis, Arnout-Jan Everts, Karen Foster, Joel Ita, Bettina Kampman, Robin Leinster, Steve Livera, Thomas Mauduit, Dick Nieuwland, Steve Oates, Pascal Richard and Roeland Roeterdink. We thank the Stanford University Rock Fracture Project for access to Poly3D. The paper benefited from thorough reviews by three anonymous referees. None of the above necessarily share the opinions expressed in this paper. The final graphic design was by Gulf PetroLink staff.

### ABOUT THE AUTHORS

**Stephen J. Bourne** received a BA in Natural Sciences from the University of Cambridge, and then completed a PhD and PostDoc. at the University of Oxford on the mechanics of distributed continental deformation. Stephen joined Shell in 1997 and is a geophysicist in the Carbonate Development Team at Shell International Exploration and Production.

**Lex Rijkels** is a reservoir engineer in the Carbonate Development Team at Shell International Exploration and Production. Lex has worked on fractured and carbonate reservoirs since 1996, following a two-year assignment in Shell’s downstream research centre, SRTCA. He received a MSc in theoretical physics from Delft University of Technology.

**Ben J. Stephenson** received a BA in Earth Sciences from Cambridge University and carried out his PhD on Himalayan tectonics and metamorphism at Oxford University. He worked for the Cambridge Arctic Shelf Programme (CASP) research group for two years, specialising in the regional geology of the Arctic and Chinese basins. Ben joined Shell in 1998 and is a structural geologist in the Carbonate Development Team at Shell International Exploration and Production.

**Emanuel J.M. Willemse** received a MSc in Geology from Utrecht University, and a PhD in Applied Earth Sciences from Stanford University, working on rock mechanical and reservoir engineering aspects of faulted and fractured rocks. In his 11 years with Shell, he has worked in various positions, including secondment to Maersk Olie og Gas AS in Denmark. He is currently senior structural geologist in the Carbonate Development Team in Shell International Exploration and Production.