NERC Micro-to-Macro Thematic
Grant Number: GST/02/2506
Principal Investigators: R. J. Knipe & D. B. Ingham
Scaling of Fluid Behaviour Associated with Flow Through
Complex Geological Structures
1111122S. D. Harris, R. Pecher, N. E. Odling, R. J. Knipe, J. A. Ellis, L. Elliott and D. B. Ingham 12Rock Deformation Research, School of Earth Sciences & Centre for Computational Fluid Dynamics
University of Leeds, Leeds LS2 9JT, UK
This NERC μ2M grant provides funding for two post-doctoral researchers to develop methods for integrating
the geological structural characterisation of fault zones with the numerical modelling of fluid flow behaviour at different scales. The project is a collaboration between a structural geology research group (the Rock Deformation Research Group), an applied mathematics group (the Centre for Computational Fluid Dynamics) and a series of industrial partners.
2 Key Achievements
; 2D, finite-difference flow model developed to handle flow through rocks with faults as partial
; Application of the flow model to sections through the stochastic model of fault damage zone geometry; ; Comparisons of flow properties predicted by flow model and minimum fault rock thickness pathways
through stochastic model;
; Establishing scaling laws for bulk permeability for the fault damage zone model for a range of fracture
properties (length distribution and spatial characteristics);
; Development of accurate 3D models of flow through complex fault zones using control volume finite-
element and Green element approaches.
3 Project Aims
We focus on creating a detailed understanding of the scaling behaviour of fluid flow through volumes containing complex geological structures, i.e. fault and fracture arrays. The correct scaling in such situations is critical to the understanding of fluid flow in hydrocarbon reservoirs, in groundwater systems, waste disposal and in the development of mineralised ore bodies. The results enable the identification, quantification and evaluation of key variables needed for modelling fluid flow in fault zones.
The objectives of this work are to: (i) develop a 3D flow simulator which models the impact of large numbers of structural heterogeneities (densely clustered arrays of faults / fractures acting as either flow retarders or enhancers) under different stress conditions on fluid flow; (ii) provide a means of capturing the micro- and meso-scale flow through complex fault zones, for input, as effective permeability tensors, into conventional simulators modelling macro-scale flow behaviour; (iii) assess the level of detailed geological characterisation and modelling needed before successful upscaling can be achieved and to evaluate the use of statistical parameters to monitor uncertainties in the upscaling. We are focusing on modelling the impact on
the flow of established and defined critical structural parameters (e.g. the size, spacing and permeabilities of fault / fracture arrays).
4 Modelling Development
Significant progress has been achieved over the past year towards realising the aims and objectives of this project. The development of modelling tools for geological domains has followed several approaches, dealing with different levels of structural complexity and linked such that the largest scale model receives „upscaled‟ input from the smaller scale models of fault damage zones. Some additional work is currently being undertaken to compare modelling results with those obtained using the Eclipse reservoir modelling package. We now briefly discuss the development of each of these tools.
Fig. 1. (a) A three-dimensional stochastically generated fault zone. (b) Contours of fault numbers in slices through a fault zone model showing clustering. (c) A three-dimensional array of faults clustered around a single main fault. (d) Slices through (c).
4.1 Stochastic Modelling
The three-dimensional fault zone flow model developed by the project allows the analysis of tortuosity and connectivity through arrays of faults using percolation techniques and pathway searching algorithms. The spatial clustering of minor faults about major faults is a critical element of the model described here (see Fig. 1), with the principal aim being to account for the „sub-seismic‟ scale damage zones which exist about
seismically mappable faults. Techniques for modelling spatial clustering, as found from core, outcrop and seismic scale data sets, have now been developed which lead to „sub-clusters‟ of faults within the damage
To provide an approach which bridges the gap between the essential details of a complex fault network and geological simplifications needed in a reservoir simulator, we have introduced a new methodology which both derives the minimum value of the fault rock thickness along flow paths traversing the fault zone and predicts areas of reduced fault zone connectivity for matrix host rock (K) and fault rock (K) of varying hf
permeabilities. This is an advance on other techniques used to model the impact of faults on cross fault flow since it maintains the true complexity and true spatial distribution of the fault damage zone, rather than treating the fault as a single fault plane or homogenising the damage zone properties. A critical threshold value of this permeability ratio is observed to exist at which the flow characteristics transfer from longer, tortuous pathways (high K/K) to short, direct pathways (low K/K) which encounter an increased thickness h f h f
of fault rock. Fig. 2 illustrates the observed fault rock thickness variation when pathways are constrained to cross the main fault at that location. If the flow pathways have been restricted to a two-dimensional layer or direct one-dimensional routes, then, when expressed as a proportion of the minimum fault rock thickness of the major fault within the domain, the additional fault rock thickness contributed by the damage zone faults was found to range up to 250% for three-dimensional pathways and 375% for direct (one-dimensional) pathways.
Fig. 2. The variation in the pathway lengths along the most efficient routes across the domain and contoured over the major fault plane for (a) three-dimensional, (b) two-dimensional (vertically exaggerated) and (c) one-dimensional pathways.
4.2 Stochastic Model Input to Discrete Fracture Flow Model
A 2D, finite difference, discrete fracture model (Odling and Webman 1991) for flow in fractured rocks has been modified to simulate flow in porous rocks with fractures as partial to complete flow barriers. In this flow model, both the faults and the rock matrix are discretized on a regular square grid and grid elements assigned permeability values representative of rock matrix and fault rock. Procedures are implemented to ensure that the connectivity of the fracture network is preserved in the process of resolving the fracture pattern onto the model grid and to account for the orientation and thickness of individual fractures.
The model has been applied to 2D slices through the stochastic models of fault damage zones described above. 4Fault thickness to length ratios of 1:10 and a permeability contrast of fault rock to host rock of 4 orders of
magnitude, typical of deformation bands and minor faults in permeable sandstones, is used. Fig.3 (left) shows an example of the simulated flow field demonstrating the preferential paths created by the fault network. Sample sizes of 5, 10, 20 and 50m across, containing between 10 and 2000 fractures, were used as input to the flow model to determine a total of 260 estimates of bulk permeability perpendicular and parallel to the main fault for each sample size. Examples using fault size distribution exponents of 1.8 and 2.2, with clustered and random spatial distributions have been run. An example of the flow modelling results in Fig.3 (right) shows the „efficiency‟ of the fracture network regions (50 by 50 m) by comparing the bulk permeability obtained with that of the same region of host rock containing a single continuous fault, oriented perpendicular to flow, which contains the same proportion of fault rock area as the fault damage zone simulation. The figure shows that 50 by 50 m regions of the fault damage zone are around 50% efficient perpendicular to the fault and between 1% and 10% efficient parallel to the fault. The model results perpendicular to the fault for (1/K-1) show a power law scaling effect with an exponent (slope of the graph) of around 1.2. This indicates that as the proportion of fault rock area increases, the efficiency of the network decreases in a power law fashion.
The flow fields resulting from the model were also used in a contaminant transport code (Odling and Roden 1997) to investigate the geometry of flow paths and the total thickness of fault rock crossed. These results are compared to computations of pathways calculated from purely geometrical considerations (section 3.1, Harris et al. 1999), see Fig.4. Preliminary results indicate that the geometrical method gives path lengths very similar to those determined by the flow model but predicts slightly lower values for the fault rock thickness crossed. These results suggest that, once calibrated, the geometrical method could provide a fast and reliable tool for estimating bulk rock permeability. These relationships will be further investigated for different types of system.
The 2D discrete fracture flow model is also being further developed to 3D. The 3D model will be computationally much more intensive which places restrictions on the number of fractures that can be modelled. However, the 3D model will provide an opportunity to test the validity, and calibrate the results, of the 2D modelling.
Fig.3 Left – flow rates in a fractured region (white- low and black – high flow). Right – model results for bulk
permeability (K) versus proportion of fault rock for 50 by 50 m regions, red circles – perpendicular, and blue
triangles – parallel to the major fault. Permeability (K) is normalized with respect to matrix permeability.
Fig.4 Flow pathways predicted by: left – the geometrical model, and right – the flow model.
4.3 Control Volume Finite-Element Model
Two alternative approaches to modelling fluid flow through complex geological structures are currently being investigated. The first study describes a control volume based finite-element scheme which can be used with unstructured grids in two- and three-dimensional domains. The use of flexible grids enables the geological features in a reservoir to be represented accurately since they can be aligned with reservoir layer boundaries, faults and horizontal/inclined wells. This method is currently being developed to incorporate clustered damage zones, utilising input from the stochastic fault zone and discrete fracture flow models.
Although methods such as Cartesian local grid refinement and related hybrid schemes provide an accurate numerical scheme and allow the incorporation of full and symmetric permeability tensors, they are difficult to align to the grid when varying facies, faults, etc. are provided. Control volume finite-difference methods utilise Voronoi or perpendicular bisector (PEBI) grids which are locally orthogonal with block boundaries normal to lines joining nodes on either side of the boundary and they enable the accurate computation of inter-block transmissibility for heterogeneous, isotropic permeability. In this study we are developing models based upon the control volume finite-element (CVFE) method, where the permeability tensor can be full, asymmetric and anisotropic and can vary from cell to cell.
In the CVFE method, the potential is treated in a finite-element manner and the majority of the domain is composed of triangles (2D) or tetrahedra (3D). The control volume approach allows local mass conservation and upstream weighting of phase mobilities over flexible grids based upon a triangular mesh of connections. In 2D, control volumes are formed around grid nodes by joining the mid-points of triangular edges to the barycentre of the triangle, with a natural extension to 3D. The assumption that the permeability tensor inside each triangular element is uniform may lead to significant numerical errors when the permeability contrasts between the triangles are large. Therefore, we assume uniform porosity and permeability inside each control volume. Since the permeability varies inside each triangle, the potential gradient inside each triangle cannot be estimated by a simple linear combination of potentials at vertices. The CVFE method maintains flux continuity across the grid block interfaces of each triangle (2D) or tetrahedron (3D) (by assuming a linear variation of the potential within each region) and potential continuity at specific points on the interfaces. This leads to a generalised harmonic average transmissibility (instead of permeability) factor.
The accurate representation of physical boundaries due to facies and faults / fractures requires that this construction be also implemented for quadrilaterals in 2D and hexahedra in 3D. In terms of facies, complex layer boundaries can then be described by property discontinuities. However, the incorporation of faults is not so simple and we have been developing methods to implement this within the CVFE approach. Completely sealing faults can simply be treated using no-flow (normal) boundary conditions for each phase. Faults which only allow cross-fault flow (no along-fault flow) can be treated using Darcy's law across the thickness of the fault to maintain a pressure gradient between the nodes, which can be considered to be coincident in the large-scale model. For open fractures, pressure continuity is assumed between nodes either side of the discontinuity. In 3D (2D), the difference in the fluxes at the nodes acts as a source/sink for 2D (1D) flow over the fracture
surface (fracture thickness varying with position), with reference to a coordinate system describing the „surface‟ of the fracture, in a variable permeability field. Faults, in general, will be neither pure flow conduits or flow retarders and the flow will be a combination of cross- and along-fault flow. The pressure gradient will establish a cross-fault flow component and the remaining fluid flux provides a source/sink for the flow over the surface of the fault.
4.4 A Highly-Accurate General Numerical Model
The two key objectives of this approach are:
1. to describe the transport processes in heterogeneous porous media at the macro-scale as generally as
2. to solve the governing equations to the highest degree of accuracy achievable by numerical techniques.
The combination of these two crucial factors has been missing in the available solutions so far. Often the complexity of the descriptive equations prevents the high accuracy of their solution and vice versa. The ambitious objective of our approach relies on several powerful tools, from which the most significant ones are the Green element method (GEM) and the object-oriented numerics (OON).
The transient fluid flow in porous media was described by means of a general partial differential equation of a diffusion type. In this form, the equation covers practically any “advanced” phenomena, such as multiphase
fluid flow, non-Darcy effects, anisotropic heterogeneous media, various internal sources, etc. This equation and the associated initial and boundary conditions represent the basic mathematical model, which may be later extended by considering other physical phenomena occurring simultaneously with the fluid flow. The coupling of the basic model with processes such as heat transfer, rock deformation, or chemical and phase changes appears realizable even within the current project due to the successful design of the computer model (see below).
The numerical technique adopted to solve the differential equation or system of equations is the Green element method. Derived from the boundary integral equation theory over a domain mesh, the GEM combines the high accuracy of boundary element techniques with the versatility of finite element techniques. Two improvements were proposed in order to keep the accuracy of the method high while increasing its efficiency for problems in two- and three-dimensional domains: (i) multi-element flux approximation  and (ii) flux-vector based integral equation . The latter is a conceptually new GEM formulation and exhibits several additional benefits in handling fluid velocity fields.
The actual treatment of the structural heterogeneities in a reservoir domain to be modelled will closely depend on the performance of the GEM-based numerical simulator and the computational resources available. Although it would be desirable to treat each geological fault in a deterministic sense, the large number of such faults in a typical fault zone prevents us to do so. Therefore, only the major faults are suitable to be handled as true discontinuities in the otherwise continuous properties of the surrounding rock, and the rest of the faults must be handled stochastically. To allow more major faults to be modelled deterministically, their governing equations were adapted to the dominant dimensionality of these geological bodies while producing more complex boundary constraints in the overall system of equations. That, in turn, has put more demands on the computer implementation of the numerical model.
The complexity of the final computer model is a result of several extreme preferences made: the singular-integral theory of the GEM, the general form of the governing equation, the complicated constraints at the fault boundaries, and the high demands on accuracy without compromising efficiency and stability. To cope with such a difficult task, the entire model was partitioned into a set of “polymorphic” models according to the principles of the object-oriented numerics. By introducing such an abstraction into the solution design, the mathematical complexity of the system of coupled equations is shifted towards the algorithmic complexity and spread among the self-contained abstract data types. Each polymorphic model is then responsible for registering its unknowns and producing the corresponding coefficients for the global matrix in a process that may be called “run-time coupling”. Although the immediate coding difficulties are no less than before, on the contrary, the ultimate benefit is: once the software is properly built, the main complexity is removed from every subsequent mathematical model. The ongoing process of the software development is proceeding well and it is expected to culminate in the second quarter of year 2001.
Harris S.D., McAllister E., & Knipe R.J. 1999 Scaling of fluid flow associated with flow through fault thdamage zones and networks. Proc. 5 Annual Conf. IAMG‟99, pub. Tapir, Trondheim,711-716.
Odling N.E. and Webman 1991 A conductance mesh approach to the permeability of natural and simulated fracture pattern. Water Resour. Research, 27, 2633-2643.
Odling N.E. and Roden J. 1997 Contaminant transport in fractured rocks with significant matrix permeability, using natural fracture geometries. J. Contam. Hydrology, 27, 263-283.
The following publications are associated with this grant:
 Harris, S.D., Pecher, R., Knipe, R.J., McAllister, E., Elliott, L., and Ingham, D.B., Scaling of fluid flow
associated with flow through complex geological structures, submitted to Modelling Permeable Rocks:
Advances in the Quantification of Geology for Petroleum and Groundwater Applications, Cambridge,
UK (March 2001).
 Harris, S.D., McAllister, E. and Knipe, R.J., 2000, Capturing the complexity of clustered fault damage
zones, submitted to Modelling Permeable Rocks: Advances in the Quantification of Geology for
Petroleum and Groundwater Applications, Cambridge, UK (March 2001).
 Odling, N.E., Harris, S.D., McAllister, E. and Knipe, R.J., 2000, Geometrical and permeability
properties of fault damage zones, submitted to Modelling Permeable Rocks: Advances in the
Quantification of Geology for Petroleum and Groundwater Applications, Cambridge, UK (March 2001).  Pecher, R., Harris, S.D., Knipe, R.J., Elliott, L., and Ingham, D.B., 2000. A more accurate Green
element method in two and three spatial dimensions. 22nd International Conference on the Boundary
Element Method, Cambridge, UK (Sept. 4-6).
 Pecher, R., Harris, S.D., Knipe, R.J., Elliott, L., and Ingham, D.B., 2000. New formulation of the Green
element method to maintain its second-order accuracy in 2D/3D. submitted to the International Journal
of Engineering Analysis with Boundary Elements.
 Harris, S.D., Pecher, R., Knipe, R.J., McAllister, E., Elliott, L., and Ingham, D.B., Scaling of fluid flow
associated with flow through complex geological structures, submitted to the AAPG Annual Meeting,
Denver, Colorado, USA (June 2001).
 Odling, N.E., Harris, S.D., McAllister, E. and Knipe, R.J., 2000, Geometrical and permeability
properties of fault damage zones, submitted to the AAPG Annual Meeting, Denver, Colorado, USA
Descriptions of the research activities of the Rock Deformation Research Group, including this μ2M grant
and other current projects, can be viewed at the following website: