Fluid flow through a two-dimensional fracture network has been simulated using a discrete fracture model. The computed field-scale permeabilities were then compared to those obtained using an equivalent continuum approach in which the permeability of each grid block is first obtained by performing fine-scale simulations of flow through the fracture network within that region. In the equivalent continuum simulations, different grid-sizes were used, corresponding to N by N grids with N = 10, 40, 100 and 400. The field-scale permeabilities found from the equivalent continuum simulations were generally within 10% of the values found from the discrete fracture simulations. The discrepancies between the two approaches seemed to be randomly related to the grid size, as no convergence was observed as N increased. An interesting finding was that the equivalent continuum approach gave accurate results in cases where the grid block size was clearly smaller than the ‘representative elementary volume’.


Fractured rock hydrology and nuclear waste disposal

Understanding the permeability of a fractured rock system is of interest in many different fields of earth science and engineering. One such field is the assessment of the long-term safety of a geological repository for radioactive wastes. Many countries, for example the UK, are considering using deep geological disposal for the isolation of radioactive waste (Nuclear Decommissioning Authority, 2010). In one illustrative disposal concept, based on the Swedish KBS–3V concept (SKB, 2009), vitrified high level radioactive wastes or spent nuclear fuels will be placed in appropriate containers that will be placed in excavated tunnels that are typically between 250 and 1000 m below the surface. The tunnels will then be backfilled with clay-like materials. Generally, the canisters, backfill material and the surrounding rock mass (together usually referred to as a ‘multi-barrier concept’) are intended to provide levels of isolation to prevent the radioactive waste from entering the human environment.

A number of countries are considering developing underground geological disposal facilities in crystalline rock. Consequently, the associated national research programs have been tasked to investigate the transport properties of crystalline rock systems (Neretnieks, 1993). The disposal facilities will be designed to function for an extended period of time, but some of the radionuclides in high-level radioactive wastes have very long half-lives, and therefore an understanding must be established as to what would happen to the radionuclides if the containers were breached. Although crystalline rocks generally have low intrinsic matrix permeabilities, on a macroscopic scale they often contain networks of interconnected fractures that may provide hydraulic pathways for the migration of radionuclides. Consequently, the macroscopic scale transport properties of fracture networks can be of crucial importance.


In 1992 the Swedish Nuclear Power Inspectorate recognized the need to enhance the theoretical background and to develop models capable of simulating coupled thermo-hydro-mechanical-chemical (THMC) processes in fractured rock, and initiated the DECOVALEX project (development of coupled models and their validation by experiments). The DECOVALEX project is an international cooperative project, with research teams funded by nuclear waste organizations in various countries (Tsang, 2009). DECOVALEX 2011 is the fifth phase of the project, and the present work was produced as a part of this phase of the project. DECOVALEX 2011 was divided into three tasks. Each task focussed on models and experiments involving different geologies and coupling behaviours. Task C focussed on hydro-mechanical coupling in fractured crystalline rocks. There were two parts to this task: (1) the benchmark test, where the computational models of each of the teams were compared as a form of verification of the methods (i.e. a demonstration that the models were solving the mathematical equations correctly); and (2) the test case, where the computational models were compared against experiments for validation. The results presented here are from the benchmark test of Task C.

Simulation approaches for fractured rock masses

Zhang and Sanderson (2002) divided mathematical models of fluid flow into two broad groups: discrete fracture networks (DFN) and equivalent continuum (EC) models. The DFN models are, in principle, the most rigorous. However, application of this method currently is limited, partly because of its computational intensity and partly because of the difficulty in acquiring detailed knowledge of the fractures at a given site. The EC models are both less demanding computationally and conceptually simpler than DFN models. The purpose of this paper is to compare the results obtained from simulations of groundwater flow when different approaches are used to model the fractures. In a discrete fracture model, each fracture is modelled explicitly, commonly as either a planar feature (in 3D models) or a linear feature (in 2D models), where it is assumed to have a specified hydraulic transmissivity. Fractures need not be planar or linear, though assuming them to be so is a routinely used simplification. Although semi-analytical methods can be used to estimate the flow through such a network (Zimmerman and Bodvarsson, 1996a), in most cases the flow through the fracture network is simulated by discretizing each fracture segment, and using a finite element or finite difference formulation to solve the equations of fluid flow. In essence, each fracture segment is represented as an element through which the volumetric flow rate is proportional to the pressure drop along the length of that segment. Discrete fracture network modelling is widely used for both site-specific simulations in engineering applications, as well as a tool to evaluate conceptual models. In the present work, the discrete fracture modelling has been conducted using NAPSAC, a simulation code that is a part of the CONNECTFLOW package that was developed by Serco (Serco, 2011a).

The equivalent continuum model approach assumes that continuum-based numerical methods can be used to model inherently discrete and discontinuous systems. In such models, individual fractures are not explicitly represented. Rather, their effect is averaged out over the continuous region, which is then treated as a porous medium. Fluid flow through the porous medium is assumed to be governed by Darcy's law (Jaeger et al., 2007):  

Where q (m s−1) is the volumetric flow rate per unit area, k (m2) is the permeability tensor, μ (Pa s) is the fluid viscosity, and P (Pa) is the fluid pressure. Use of Darcy's law is usually assumed to be valid when the characteristic volume considered in the problem (such as the size of an element in the finite element mesh) is greater than the so-called representative element volume (REV) of the medium (Bear, 1972). The REV is supposed to be a critical volume such that constitutive properties can be assumed to be constant if defined over regions larger than the REV. Long et al. (1982) concluded that a fractured system would behave more like a continuum if it has a high fracture density and a narrow aperture distribution.

For modelling underground water movement, the hydraulic properties for the ‘equivalent continuum’ can either be obtained from hydrogeological material properties (from interpretations of measurements from either the field or the laboratory), or estimated using some other simulation techniques, such as from discrete fracture network models on a smaller scale. In the present work, the equivalent continuum modelling has been conducted using NAMMU (Serco, 2011b), a finite element code for simulating fluid flow in porous media.

The foregoing discussion applies to rock masses in which the matrix rock has very low permeability, in which case the fluid flow takes place almost exclusively through the fracture network. If the matrix permeability cannot be neglected, the problem can again be addressed either by representing each fracture explicitly, and modelling fluid transfer between each individual fracture and the matrix rock, or by smoothing out the hydraulic properties and using a dual-porosity continuum approach (Zimmerman et al., 1996). Porous media in which flow within the matrix rock cannot be neglected are not addressed in the present work.

Fluid flow simulations

Benchmark test fracture network

A 2D fracture network is generated stochastically. Fractures are modelled as straight lines, and the aperture of each fracture does not vary along the fracture length. Although fractures in the field are unlikely to be straight lines with non-varying apertures, this is a common assumption to make for simplifying the modelling process. The transmissivity T (m4) of each fracture is assumed to follow the cubic law: T = wh3/12, where h (m) is the aperture and w (m) is the width within the fracture plane normal to the direction of flow (Zimmerman and Bodvarsson, 1996b)1. The volumetric flow Q (m3 s−1) is then given by Q = −(T/μ)gradP. As mentioned above, fluid flow is assumed to take place only through the fractures, and the rock matrix is assumed to be impermeable.

The fracture network occupies a macroscopic region of 20 × 20 m, and contains 7797 fractures, and approximately 75,000 fracture segments, defined as a portion of a fracture lying between two adjacent intersections with other fractures. Fracture lengths vary between 0.25 and 30 m, and are distributed according to a truncated power law distribution. The fracture length and orientation distributions are based on data from the Sellafield area, west Cumbria, England, where formations in the Borrowdale Volcanic Group, a thick sequence of Ordovician volcaniclastic rocks, have been characterized; the statistical parameters can be found in Min et al. (2004). The fracture apertures vary between 0.5 and 200 μm, distributed according to a truncated lognormal distribution, with a lognormal standard deviation of 1. The fracture lengths are correlated to the apertures, using the method described by Baghbanan and Jing (2007). It must be noted that although the Sellafield fracture data are used to create a fracture network having some geological and hydrological realism, the simulation results generated from this set of parameters do not necessarily represent the general hydraulic behaviour of the Sellafield area.

A visual representation of the fracture network is shown in Fig. 1. With 7797 fractures in the space of 400 m2, the fracture density of this fracture network is 19.5 m−2. The ‘horizontal’ direction aligns with the x axis, and the ‘vertical’ direction aligns with the z axis. These are nominal directions, for the purpose of simplifying the discussion of the results; gravity plays no role in the flow processes under consideration.

Fluid flow through the network is induced by a constant pressure gradient of ΔP = 200 kPa imposed across the fracture network, in the vertical direction. The nominal pressure gradient is therefore 10 kPa m−1. The system is assumed to be of unit depth in the third direction.

Discrete fracture network simulations

The flows through the network are first evaluated using NAPSAC, a finite element code that was designed for simulating steady-state constant-density groundwater flow through a fracture network. The basic algorithm is very simple. Assuming the fracture transmissivity is related to the fracture aperture by the ‘cubic law’ (Zimmerman and Bodvarsson, 1996b), the flow in each fracture is determined numerically using a Galerkin finite element approach, using the constraints that water pressure is continuous between any two intersecting fractures, and water mass is conserved at each intersection. The flow through the overall network is obtained by summing up the flows in the different fractures. Full details of this algorithm are given in the NAPSAC Technical Summary (Serco, 2011a). Discrete fracture modelling can be readily applied to fluid flow through the benchmark test network, as the locations, lengths, orientations and transmissivities of each of the fractures are explicitly known.

With a linear pressure gradient of 10 kPa m−1 imposed from the top to the bottom of the network, the flow rate through each fracture segment is calculated by NAPSAC. The outflow at the bottom of the 20 × 20 m region is plotted in Fig. 2. The total volumetric flow rate through this boundary is 1.36 × 10−4 m3 s−1. From Darcy's law, with μ = 0.001 Pa s, and A = 20 m2, the computed kzz component of the macroscopic permeability tensor is found to be to 6.8 × 10−13 m2. The detailed outflow profile and the total flow through the lower outflow boundary are used as the main basis for comparison between the discrete fracture and equivalent continuum models.

The full macroscopic permeability tensor k (m2) can also be calculated using NAPSAC. This is done by solving the same problem as described above, for a range of orientations of the imposed pressure gradient, chosen to give a uniform coverage of directions. The effective permeability tensor is determined from a ‘best fit’ between the fluxes and head gradients. Again, a more detailed explanation of the algorithm is given in the NAPSAC Technical Summary (Serco, 2011a). The permeability tensor for the benchmark test network was computed by NAPSAC to be:  
where the x axis is oriented horizontally in Fig. 1. Note that the value of kzz in the above permeability matrix does not precisely agree with the ‘exact’ value of 6.8 × 10−13 m2, because the permeability matrix is found by fitting fifteen directional permeabilities over the entire range of 360°.

The small, but non-zero, off-diagonal terms indicate that the principal directions of the permeability tensor are not exactly aligned with the horizontal and vertical axes of the test network. The principal permeabilities, which are the eigenvalues of k, are 7.43 × 10−13 m2 and 4.89 × 10−13 m2; the direction of the larger principal permeability is rotated by 7.4° from the z axis. A convenient scalar parameter to use to quantify the permeability tensor is its determinant, det k = keff, which is also equal to the geometric mean of the two principal permeabilities, and is sometimes referred to as the effectivepermeability. This parameter rigorously appears, for example, in the equations for well testing for a vertical well located in an anisotropic formation (de Marsily, 1986).

Equivalent continuum simulations

Fluid flow through the benchmark test network can also be simulated using the finite element continuum code NAMMU (Serco, 2011b). To do this, the network is first divided into an N × N grid of equal-area squares, where N is taken to be 10, 40, 100 and 400 (Fig. 3). The full permeability tensor of each of these small grid blocks is calculated by performing a discrete fracture network simulation using NAPSAC. Each grid block is now treated as a porous continuum having the permeability tensor calculated by NAPSAC, and fluid flow through this heterogeneous, locally anisotropic continuum is simulated using NAMMU.

The network was modelled using a grid with 10 × 10 quadratic finite elements (100 elements in total), each with dimensions 2 × 2 m. The permeability of each of the grid blocks was calculated using NAPSAC, and their effective permeabilities are plotted in Fig. 4. The effective permeabilities are defined as the geometric mean of the two principal permeabilities; they are not used in any flow calculations, but are convenient for illustrating the heterogeneity of the network. Despite the heterogeneity of the fracture network, the permeabilities at this scale vary only within about one order of magnitude. At this scale, the permeability distribution appears random, and there is little obvious macroscopic structure.

The flow results are shown in Fig. 5. The flows through the outlet boundary of the effective continuum model are compared with the flows from the discrete fracture network model. The results from the discrete fracture model are plotted as a histogram with a bin size corresponding to the element size of the continuum model to allow comparison. Note that the flow profile of the NAMMU model does not consist of a series of straight lines: a quadratic profile is fitted across each finite element. The flow profile and the total flow (1.50 × 10−4 m3 s−1, compared to NAPSAC's value of 1.36 × 10−4 m3 s−1) show that the results between the continuum and discrete models are reasonably compatible. The permeability in the vertical direction calculated from the flow is 7.50 × 10−13 m2, which is 10% greater than the value calculated from the discrete fracture model.

Figure 6 shows the distribution of the effective permeability for the 40 × 40 grid, with 1600 grid blocks in total, each with dimensions 0.5 × 0.5 m. As the grid blocks over which the permeability is defined become smaller, the distribution of the permeability becomes wider; approximately 40% of the grid blocks have permeabilities less than 1 × 10−13 m2, and about 3% of the grid blocks have permeabilities greater than 1 × 10−12 m2. At this scale, some structure seems to be appearing in the geometric distribution of permeabilities.

The flow results for the 40 × 40 grid are shown in Fig. 7. The detailed flow profiles follow reasonably closely the trend from the discrete fracture simulations. The total outflow is 1.60 × 10−4 m3 s−1, which corresponds to a permeability of 8.00 × 10−13 m2, only 18% greater than the values given by the discrete fracture simulations.

The network was also modelled using a grid with 100 × 100 finite elements (1 × 104 elements in total), each with dimensions 0.2 × 0.2 m. This gives an average of 7.5 fracture segments per element. The distribution of effective permeabilities is shown in Fig. 8. The majority of the grid blocks have very low permeabilities. With such a small number of fracture segments in each grid block, the permeability tensors are likely to be very anisotropic. The permeability distribution is again plotted in Fig. 9, but now in terms of the maximum permeability. Unlike the permeability maps plotted previously, the high permeability regions in Fig. 9 arranged themselves into lines, which essentially follow the paths of the larger (hence higher permeability) fractures seen in Fig. 1.

The flow results for the 100 × 100 grid are shown in Fig. 10. The detailed flow profile agrees qualitatively with the profile computed by the discrete fracture network, and the total flow of 1.46 × 10−4 m3 s−1 compares well to NAPSAC's value of 1.36 × 10−4 m3 s−1. The permeability in the vertical direction calculated from the continuum model is 7.30 × 10−13 m2, which is 7% higher than the value obtained from the discrete fracture model.

The flow results for the N = 400 model (in which case each element has dimensions 0.05 × 0.05 m, and there are only approximately 0.5 fracture segments per element) are plotted in Fig. 11. At this level of discretization, each grid block would typically have between 0–4 fracture segments. An assembly of such grid blocks into a finite element model can be compared to constructing a discrete fracture model. The flow profile and the total flow (1.25 × 10−4 m3 s−1, compared to NAPSAC's value of 1.36 × 10−4 m3 s−1) show that the results between continuum and discrete models are reasonably compatible. The permeability in the vertical direction calculated from the 400 × 400 continuum model flow is 6.25 × 10−13 m2, which is only 8% lower than the value from the discrete fracture model.

Discussion and conclusions

A series of effective continuum models of the benchmark test fracture network were created, at different levels of discretization, and their flow behaviours were compared against results obtained from discrete fracture network modelling. Although at the various levels of discretization the permeability displayed very different distributions, the overall flow through the region seemed to be robust over a large range of N. Note that the discrete fracture network simulations can be considered to represent the case of N = 1. The value N = 400 represents a sensible maximum value, as at this scale (corresponding to grid blocks of 0.05 m size) most grid blocks contains no more than four fracture segments. This study shows that as long as the full permeability tensors for each of the grid blocks are considered, using continuum models should give reasonable answers, regardless of the level of discretization. This study also suggests the possibility of constructing discrete fracture network-like models using effective continuum type simulation packages. Another interesting conclusion is that the continuum models gave accurate results, even in cases where most of the individual grid blocks were below what would normally be considered to be the REV, as they contained a very small number of fracture segments.

It has generally been held, as far back as the seminal paper by Long et al. (1982) that for length scales below the REV it is not possible to replace a discrete fracture network with an equivalent porous continuum. Specifically, Long et al. (1982) stated that: “the following criteria must be met in order to replace a heterogeneous system of given dimensions with an equivalent homogeneous system for the purposes of analysis: (1) There is an insignificant change in the value of the equivalent permeability with a small addition or subtraction to test volume. (2) An equivalent symmetric permeability tensor exists which predicts the correct flux when the direction of gradient in a REV is changed.”

However, the current results seem to indicate that, even if the chosen computational cells are clearly smaller than the REV, it is still sensible to represent each of these small regions by a porous continuum. When viewed at a larger scale, this equivalent porous medium will appear as a heterogeneous medium with different permeability tensors in each cell. Nevertheless, the macroscopic flow rates calculated through this porous medium will closely match the flow rates calculated through the discrete fracture network, regardless of whether or not an REV scale even exists.

Of course it should be noted that all of the results and conclusions of this work pertain to two-dimensional fracture networks; extension of these ideas to three dimensions will be the subject of future work.


This work was funded by the UK Engineering and Physical Sciences Research Council (EPSRC) and the Radioactive Waste Management Directorate (RWMD) of the UK Nuclear Decommissioning Authority. The authors thank Cherry Tweed and Simon Norris of RWMD for supporting this work.

Other definitions of the transmissivity are also in wide use, such as T = ρgh3/12μ, where ρ (kg m−3) is the fluid density, g (m2 s−1) is the magnitude of the acceleration due to gravity, and μ (Pa s) is the fluid viscosity. The latter definition is based on using hydraulic head rather than pressure, P (Pa), as the driving force in Darcy's law. Both definitions are essentially equivalent.
P1Freely available online through the publisher-supported open access option.