Current simulation technology for naturally fractured reservoirs is based on either the continuum or the discrete-fracture approach. The more commonly used continuum model can simulate complex recovery mechanisms. However, it uses a very simplified representation of the fracture system for calculating effective fracture permeability. The discrete-fracture flow method can handle complex fracture geometry. However, its use has been typically limited to basic flow calculations through a connected fracture system embedded in zero-matrix-permeability rock.We have developed a new technique for estimating the effective permeability of gridblocks used in conventional simulators. The idea behind this technique is to integrate the realism of fracture systems, as captured by discrete-fracture models, with the complexity of the flow calculations offered by continuum models. The end product of developing this technique is an efficient numerical code based on the boundary-element method. This code permits the fracture system to be complex and poorly connected, and it also includes the contribution from flow through the matrix rock. For fluid flow in the matrix rock, the fractures are treated as planar-source distributions. Periodic boundary conditions, for the flow properties, are used for the calculation of the effective permeability of individual gridblocks.We first use a simple fracture system to demonstrate the validity of our method and to evaluate the sensitivity of the results to matrix and fracture properties. We then use fracture statistics data from the Mesaverde sandstone, effective permeability values from our code, and a continuum simulator to calculate tracer-flow patterns for a more realistic system.

Summary In this paper we describe a hierarchical approach for modeling fluid flow in a naturally fractured reservoir with multiple length-scale fractures. Based on fracture length relative to the finite-difference grid size, fractures are classified as short disconnected fractures, medium-length fractures, and long fractures. Effective grid-block permeabilities, associated with the short- and medium-length fractures, are used as input to a finite-difference simulator. We also present a separate transport equation for flow through long fractures to capture effects of large-scale high permeability pathways. Simulation results from the new hierarchical approach are compared with those from the conventional dual porosity/permeability model. Introduction The permeability of a fracture is usually much larger than that of the rock matrix. Consequently, the fluid will flow mostly through the fracture network, if the fractures are well connected. This implies that the fracture connectivities and their distribution will determine fluid transport in a naturally fractured reservoir.1 With this fact in mind, some reservoir simulators have been constructed that only model the flow through fracture systems. However, such simulators cannot be used to compute effective conductivities for fractured gridblocks because of their failure to include the flow through the permeable rock matrix. The loss of fracture connectivity occurs often and has been observed, for example, by Laubach2 and by Lorenz and Hill.3 In the model we propose, flow through both the fractures and the matrix is considered, so that losses in fracture connectivity do not automatically cause the calculated effective permeability to drop to zero. Most reservoir simulators employ dual continuum formulations for naturally fractured reservoirs, in which matrix blocks are divided by very regular fracture patterns.4,5 Part of the primary input into these simulation models is the permeability of the fracture system assigned to the individual gridblocks. Currently, this value can only be reasonably calculated if the fracture systems are regular and well connected. However, field characterization studies have shown that fracture systems are very irregular, often disconnected, and occur in swarms.2,3 Incorporating the statistical description of fractures into the characterization of simulation grids will significantly improve the modeling technology for fractured reservoirs. In our previous papers6,7 we applied periodic boundary conditions in computing the effective conductivities of a gridblock containing fractures. Fluid flow in the fractures was treated as two-dimensional (2D) potential flow governed by Darcy's law. The permeability is proportional to h2/12 where h is the fracture aperture.8 The flow in the rock matrix is a three-dimensional (3D) potential flow which also obeys Darcy's law. We developed an efficient method to solve these coupled potential problems using the boundary element method.6,7 Kamath et al.9 showed that disconnected fractures cannot be ignored because they can significantly contribute to the overall flow through the rock matrix. They also showed that the effective permeability is much larger than a simple sum of fracture and matrix permeabilities, as is commonly assumed. Jensen et al.10 applied this new method to two field cases and obtained some successful simulation results. Most naturally fractured reservoirs include fractures of multiple length scale. The effective gridblock permeability calculated by the boundary element method becomes expensive as the number of fractures increases. The calculated effective properties of gridblocks also underestimates the properties of long fractures whose length scale is much larger than the gridblock size. In this paper, we propose an efficient hierarchical modeling approach that overcomes these two unsatisfactory conditions based on the fracture length relative to the finite difference grid size. We assume that short fractures are randomly distributed and contribute to increasing the effective matrix permeability. An asymptotic solution representing the permeability contribution from short fractures is derived. With the short fracture contribution to permeability, the effective matrix permeability can be expressed in general tensor form. Thus, we also develop a boundary element method for Darcy's equation with tensor permeability. For medium-length fractures in a gridblock, a coupled system of Poisson equations with tensor permeability is solved numerically using a boundary element method. The gridblock effective permeabilities are used with a finite-difference simulator to compute flow through the fracture system. The simulator is enhanced to use a control-volume finite-difference formulation11–14 for general tensor permeability input (i.e., a 9-pt stencil for 2D and a 27-pt stencil for 3D). For simplicity of derivation, both the medium-length and long fractures considered in this paper are perpendicular to the bedding surfaces. Because the effective medium-scale calculation underestimates the permeability contribution from long fractures, we explicitly model them as fluid conduits. Thus, we derive a well-like equation for a finite reservoir simulator to model flow through these discrete long fractures. In the next section, we describe the governing equations, and then we present our new hierarchical approach to modeling flow in porous media with multiple length-scale fractures. To illuminate this new approach, we present two numerical examples: a two-dimensional model and a three-dimensional three-layer model. Using the three-layer model, we also compare our single continuum results with those from a conventional dual porosity/permeability model. Finally we conclude with a summary of our main results. Governing Equations The main premise behind our model is that the fracture aperture is small enough to allow the fractures to be treated as planar high-flow regions embedded in the homogeneous porous matrix. The flow in the matrix is governed by Darcy's law, v m = − 1.127 × k m μ ∇ ϕ m = − k ^ m ∇ ϕ m , ( 1 ) where vm (B/D/ft2) is the velocity, km (darcy) is the matrix permeability, ? is the fluid viscosity (cp), and ?m (psi) is the matrix potential or pressure. Clearly k^m=1.127×km/μ and is more convenient to work with than km.

TX 75083-3836, U.S.A. fax 01-214-952-9435. AbstractWe have developed a throat-bulb network model to simulate multiphase flow experiments. The model has been efficiently implemented on massively parallel computers (MasPar and eM5). This model has allowed us to study the dynamics of displacement processes. Using grids of size 64 x 64 x 128 (half a million cells), we are able to incorporate random variations in the throat and pore size distributions and observe flow patterns such as viscous fingering, stable front displacement, and invasion percolation. Furthermore, we are able to discern certain scaling rules in terms of macroscale variables.Displacement processes in a network model with random distribution of throat and pore size can be characterized by three parameters: M (viscosity ratio), Ca (capillary number), and Cg (the ratio of capillary number and Bond number). Macroscale variables, f{ rand Sor, are strongly dependent on flow patterns, whereas the dynamic capillary pressure, Pc, is independent of flow patterns.Analysis oflaboratory measurements of these macroscopic variables is generally based on steady state approximations and moreover, it offers no clue as to how one might extrapolate the results to the reservoir scale. Using the pore network simulator, we can model dynamic flow behavior and understand the scaling rules for experimental measurements. The simulator will eventually become a practical tool to investigate multi-phase flow in various flow conditions. Furthermore, a better understanding of experimental measurements will help us construct effective grid-block properties for reservoir simulation.

In finite difference reservoir simulation, a well is generally treated as a point source or sink. As a result, a well productivity index must be specified to relate the difference between wellblock pressure and wellbore pressure to the production rate. However, analytical formulas for productivity indices are known only for a vertical or horizontal well in a homogeneous formation. In this paper, the method of reflections and slender body theory are shown to provide a means to obtain accurate measures of the productivity index in complex reservoir formations. Using slender body theory, we have derived an asymptotic solution of single-phase flow around an inclined well that is applied here to the determination of the productivity index. The method of reflections can account for the presence of layer boundaries through the application of image solutions. In combination, these techniques permit the calculation of the productivity index in anisotropic formations with layer boundaries. A comparison of the results with numerical solutions verifies their accuracy. The previous formulas sire shown to be recovered as a special case. Several examples are presented. The effect of layer boundaries on the productivity index is shown in both two and three dimensions. We consider both permeable and impermeable layers. We also investigate the productivity index of an inclined well in a layered anisotropic formation.