Transp Porous Med (2012) 93:777–797 DOI 10.1007/s11242-012-9982-3

Estimating the Hydraulic Conductivity of Two-Dimensional Fracture Networks Using Network Geometric Properties Colin T. O. Leung · Robert W. Zimmerman

Received: 9 June 2011 / Accepted: 6 March 2012 / Published online: 21 March 2012 © Springer Science+Business Media B.V. 2012

Abstract Fluid flow through random two-dimensional fracture networks is investigated, with the aim of establishing a methodology for estimating the macroscopic effective hydraulic conductivity based on the parameters of the fracture network. A wide range of isotropic networks is examined: the lengths are either uniform, or follow a power law or lognormal distribution; the apertures are either uniform, or proportional to the fracture lengths. A methodology is developed that utilises the fracture density and the aperture distribution, but does not require explicit solution of the flow equations. This method provides an accurate estimate of the macroscopic hydraulic conductivity, for all cases considered, spanning ten orders of magnitude. Keywords Permeability · Hydraulic conductivity · Effective medium theory · Fracture network · Percolation

1 Introduction Most oil and gas reservoirs are naturally fractured, and the fractures have a large influence on hydrocarbon production (van Golf-Racht 1982). In addition, the potential consequences of any natural fracturing must be considered when determining the suitability of a rock mass to serve as the host rock for nuclear waste disposal (Wu et al. 1999). In fracture-bearing geologies, the network of fractures will provide the main path for fluid to flow through the rock mass. In many cases, the fracture density is so high as to make it impractical to model it with a discrete fracture network (DFN) approach. For such rock masses, it would be useful to have recourse to analytical, or semi-analytical, methods to estimate the macroscopic hydraulic conductivity.

C. T. O. Leung · R. W. Zimmerman (B) Department of Earth Science and Engineering, Imperial College London, London SW7 2AZ, UK e-mail: [email protected] C. T. O. Leung e-mail: [email protected]

123

778

C. T. O. Leung, R. W. Zimmerman

The exact geometry of the fractures in fractured rock masses is sometimes unavailable, and modellers would have to generate stochastic approximations of such networks to apply their numerical methods. In this study, the aim is to develop a method to relate a fracture network’s hydraulic conductivity to some simple statistical parameters of the network, such as, for example, the mean and standard deviation of the lengths and apertures. The simple case of fracture networks with identical apertures for each fracture is first studied, and two approaches are presented to estimate the hydraulic conductivities. The first approach uses the connectivity, the average fracture lengths, and the number of fractures, to estimate the hydraulic conductivity. The second approach uses the nodal density of the network to describe the random fracture network as a regular checkerboard-like lattice of conducting segments of equal conductance. The hydraulic conductivity of fracture networks with variable apertures can be estimated using the checkerboard method. In contrast to the uniform aperture case, the conductance of the segments are no longer equal, and some effective medium technique must be applied to the segments before calculating the macroscopic permeability. Throughout this investigation, DFN simulations are conducted to provide the “exact” results against which the models are tested.

2 Review of Previous Work An analytical method was developed by Snow (1969) for calculating the (possibly anisotropic) permeability tensor of an idealised jointed rock. Each joint or fracture was assumed to have smooth parallel plane walls of infinite extent, and an arbitrary but uniform aperture. Snow showed that the permeability tensor of this rock mass, ki j , is given by ki j =

1  b3 (δi j − n i n j ), |n i Di | 12

(2.1)

where δi j is the Kronecker delta, n i is the direction cosine of the normal of the fracture, and b is the aperture of the fracture. The summation is taken over all members of a joint set intersected by a sampling line Di , since the flow through each fracture is assumed to be additive. The location of any specific member of a joint set is not explicitly accounted for. In reality, fractures in rock are not infinitely long. Much work has been done to extend Snow’s study to finite-length fracture networks. Long and Witherspoon (1985) studied twodimensional networks of fractures in which all the fractures have uniform aperture and length, examining the relationship of the degree of interconnection to the macroscopic hydraulic conductivity. Their study was based on numerical simulations of fluid flow in fracture networks. They showed that knowledge of the fracture frequency (the number of fractures in a given set per unit length of sample line) alone is not sufficient to determine the hydraulic conductivity without measurements of lengths and fracture densities. Hestir and Long (1990) used simple two-dimensional Poisson network models with uniform-aperture fractures to study the macroscopic conductivity of partially connected networks as a function of the degree of interconnection and scale of measurement. They proposed a relationship between a random fracture network and an analogous regular lattice network having a fraction p of its bonds present. For a given fracture density, the approximation proposed by Snow, in which all the fractures are assumed to have infinite length, will give an upper bound on the permeability. For a given fracture network, Hestir and Long use this upper bound to be the permeability of the analogous regular lattice network with all bonds present ( p = 1). They then proposed several methods to find a value of p that is representative of the

123

Estimating the Hydraulic Conductivity of 2D Fracture

779

random fracture network, based on its connectivity, which is defined as the average number of intersections per fracture. Once a representative value of p is estimated, the permeability was calculated either by ‘percolation theory’ or ‘effective medium theory’, as these methods are typically used on regular lattices. Kirkpatrick (1973) developed an analytical solution for the effective conductance of a resistor network in which the conductances of each of the resistors are random and uncorrelated. If each of the resistors in the network, with conductance Ci , is replaced with the effective ˆ the overall conductivity of the network should not change. Kirkpatrick found conductance C, ˆ that C can be estimated by solving the following equation: N  i=1

Cˆ − Ci = 0, [(z/2) − 1]Cˆ + Ci

(2.2)

where the co-ordination number z is defined to be the average number of bonds that meet at each node. Effective medium theory in the form used by Kirkpatrick was compared with numerical simulations for randomly generated 2D lattices, representing porous media, in a study by Adler and Berkowitz (2000). The conductivities were distributed by a lognormal distribution, with the lognormal standard deviation taking on values between 0 and 8. They showed that effective medium theory gave good agreement with numerical simulations for the 2D case, if the lognormal standard deviation is less than about 6. Desbarats (1992) proposed that the effective conductance of a porous medium could be ˆ is estimated using a generalised power law average. The power average conductivity, C, given by  1/ω N  1 Cˆ = Ciω for ω  = 0, N i=1   N  1 Cˆ = exp ln Ci for ω = 0, (2.3) N i=1

where ω is the power mean exponent. Desbarats’ expression for Cˆ reduces to the arithmetic, geometric, and harmonic mean of the conductances when ω = 1, 0, or −1, respectively. Desbarats showed that for a porous medium network in which the conductances are distributed lognormally, a value of ω of about 1/3 gave good agreement with numerical results, when the lognormal variance is less than 2. de Dreuzy et al. (2001a,b) carried out studies on random fracture networks, first for cases with uniform fracture apertures, and then for cases with lognormal aperture distribution. The lengths of the fractures for both studies followed a power law distribution, n(l) ∼ l −a , where a network with a power law exponent a < 2 is one that is dominated by long fractures, and a network with mostly small fractures corresponds to a > 3. For the uniform fracture aperture case, they presented, supported up by numerical simulations, several different expressions for obtaining the permeability from a, p (a unique connectivity parameter), L (the system size), and lmin (the minimum fracture length), where the expressions are chosen according to the values of a and p. For the second part of their studies, de Dreuzy et al. (2001b) studied systems where the apertures of the fractures vary according to a lognormal distribution. For the case where the fracture apertures are not correlated to fracture length, they showed that the expression for the hydraulic conductivity of a fracture network can be broken down into two terms: the first

123

780

C. T. O. Leung, R. W. Zimmerman

term corresponding to the conductivity of the network with fractures of identical apertures, and the second term to account for the variation of the conductivity, somewhat similar to an effective medium calculation. They concluded that the expressions used to describe these systems are complicated, and depend on the parameters of the fracture network being modelled, as in uniform aperture cases. Similar conclusions can be draw for the case where length and aperture are perfectly correlated, but the overall conductivity in the latter case tends to be much higher. Zimmerman and Bodvarsson (1996) divided a fracture network into a network of nodes and segments, using the cubic law to calculate the conductance of each segment, and used the effective medium approximation of Kirkpatrick to obtain the effective conductance, assuming that the coordination number for an irregular network is the average coordination number of all the nodes in the network. The segment density of the network was approximated by finding the average number of intersections that the fractures have with each of the boundaries of the network. The overall conductivity of the network was approximated using the segment density and the effective conductance of the segments. Using and extending some of these ideas from previous work, the aim for the present study is to show that it is possible to estimate the hydraulic conductivity of a fracture network from the fracture generating parameters, without running full numerical simulations.

3 Defining Fracture Networks The work presented below proceeds as follows. First, a series of fracture networks is generated, and the conductivities of these networks are calculated numerically. Then, using the parameters of the fracture networks and the numerically calculated conductivities, some empirical or semi-analytical relationships are developed between these parameters. Fracture networks are generated stochastically. Each fracture is modelled as a pair of straight parallel plates. (This simplification is not crucial to the investigation. It merely allows the conductance distribution to be discussed in terms of the aperture distribution.) The orientation of each fracture is uniformly randomly distributed. Hence, only nominally isotropic systems are considered in this study. The main variables for generating fracture networks are the fracture densities and the distribution of lengths. The fracture density, ρ, is the number of fracture centres per unit area. In this study, ρ typically lies between 10–30 m−2 . For each fracture network, fracture centres are seeded over an extended domain that is twice the length of the network boundary. For example, if the desired simulation domain is 20 m × 20 m, the fracture centres are seeded over an extended domain of 40 m × 40m. Fractures that lie wholly outside of the simulation domain are removed, and ones that are partially within it are truncated so that they are wholly within the boundary. The purpose of this truncation method is to eliminate edge effects (Long and Witherspoon 1985). Within the simulation domain, only fractures that belong to a cluster that spans from one boundary to the opposite can contribute to the flow. Such a cluster may not always exist. The relationship between the cluster sizes and geometric parameters of the network can potentially be obtained by using percolation theory. Tóth and Vass (2011) studied the relationship between the size of the percolation cluster with parameters that describes the number, the length distribution, and the dip of the fractures using numerical modelling. Zhang and Sanderson (1998) demonstrated the use of percolation theory for studying the flow behaviour of a fracture network having uniform apertures.

123

Estimating the Hydraulic Conductivity of 2D Fracture

781

For this study, percolation is not considered. Only fractures that contribute hydraulically are included; fractures that are hydraulically isolated are systematically removed from the network. Fractures are said to be isolated when they are not connected to any other fracture in the network, or when they are connected to only one other fracture. The connectivity, ζ , of a fracture is the number of times that it intersects another fracture, and so isolated fractures are those fractures with ζ < 2 (assuming that any isolated clusters have been removed). The removal of hydraulically inactive fractures is performed in three stages. First, fractures that are not connected to any other fractures are removed. Then, since only fracture clusters that touch the boundary contribute to the flow, clusters that cannot be traced back to the boundaries are removed. Finally, fractures with ζ = 1 do not contribute hydraulically, and are therefore removed. The process of removing them is an iterative process: for each fracture with ζ = 1 that is removed, the value of ζ of the adjoining fracture is reduced by 1, and potentially to less than 2, rendering it isolated. Networks that do not contain a spanning cluster after the removal of isolated fractures are not considered in this study. For a network that is below the percolation limit most, if not all, of the fractures can be removed. Fracture lengths in a fracture system are typically distributed according to some form of skewed distribution, for example, a lognormal distribution or power law distribution (Long and Witherspoon 1985). Both of these distributions are considered in this article. A lognormal distribution is a probability distribution of a random variable whose logarithm is normally distributed. For the networks where the fracture lengths follow a lognormal distribution, mean fracture lengths, lmean , are all set to 1 m, and the log-standard deviation, σln = 0, 0.5, 1, 3, 4 and 5. When σln = 0, all the fractures have the same length, and when σln = 5, the fracture networks are dominated by a few long fractures, with many small disconnected micro-fractures. Figure 1 shows the length distributions of the fracture networks, and the effect that removing hydraulically inactive fractures will have on the overall distribution. Figure 2 gives a visualisation of two typical realisations of the fracture network, with σln = 1 or 4, in both cases with ρ = 20 m−2 and lmean = 1 m. The process of removing fractures alters the properties of the fracture networks. For networks with a high spread of fracture lengths, the overall density ρ will decrease and lmean will increase, because of the higher proportion of disconnected shorter fractures. Since the processes of truncating the fractures and the removal of isolated fractures have the effect of changing the fracture density. The notation ρi is introduced to indicate the initial fracture density before any manipulation of the fracture network, while ρ is the actual fracture density of the network being investigated. Power-law length distributions are also included in this study, in order to develop a method for estimating permeability that is not restricted to a specific distribution function. For the power law distribution, the fracture network parameters of Min et al. (2004) are used. The fracture lengths are given by a power law: N = 4l −D ,

(3.1)

where N is the number of fractures with length longer than l per unit area. From Min et al., the exponent D is fitted to be 2.2, giving lmean ∼ 1. A truncated distribution is used, where the minimum and maximum cut-offs (lmin , lmax ) are 0.5 and 250 m, respectively. The fracture lengths are generated by the following equation: −D −D −D −1/D l = [lmin − F(lmin − lmax )] ,

(3.2)

123

782

C. T. O. Leung, R. W. Zimmerman

Fig. 1 Typical distribution of fracture lengths l for the lognormal distribution, before (top) and after (bottom) removal of hydraulically inactive fractures, for σln = 0, 1, 3, 4, 5

where F is a random number that is uniformly distributed between 0 and 1. Values of D = 1.2, and D = 3.2 are also used to investigate the effect of this parameter. A complete list of parameters used is listed in Table 1.

123

Estimating the Hydraulic Conductivity of 2D Fracture

783

Fig. 2 Examples of fracture networks with fracture length distributed according to a lognormal distribution, with σln of 1 and 4 (top and bottom), with ρi of 10 m−2 and 20 m−2 , mean fracture length of 1 m, before and after removing hydraulically inactive fractures (left and right). Fracture lengths are coloured according to their lengths: blue (0–0.5 m), cyan (0.5–1 m), green (1–5 m), and red (more than 5 m)

It was noted earlier that the linear size of the extended domain of fracture centre seeding is twice the simulation domain, so that the contributions of longer fractures that can potentially be generated outside the simulation domain but are partially within it. With the parameters used here for fracture generation, the fracture lengths can be up to 250 m for the power law distribution and limitless for the log normal distribution, such that an extended domain size of twice the simulation domain may not be sufficient to capture all the relevant fractures. When the extended domain is 40 m × 40 m, the minimum length of a fracture that is generated outside the extended domain, but can potentially influence the simulation domain, is 20 m. Figure 1 shows that for fractures with lengths distributed according to a lognormal distribution, the case where σln = 3 gives the highest frequency of fractures with length over 20 m, at 0.64 %. It is reasonable to assume that only a small number of these “long” fractures will have the required orientation and location to reach the simula-

123

784

C. T. O. Leung, R. W. Zimmerman

Table 1 List of scenarios studied Case number

Scenario name

Parameters

Realisations

σln

ρi (m−2 )

Log-normal distributiona 1

M1SD0R10

0

10

10

2

M1SD0R12

0

12

10

3

M1SD0R14

0

14

10

4

M1SD0R16

0

16

10

5

M1SD0R18

0

18

10

6

M1SD0R20

0

20

10

7

M1SD0R25

0

25

10

8

M1SD05R10

0.5

10

10

9

M1SD05R12

0.5

12

10

10

M1SD05R14

0.5

14

10

11

M1SD05R16

0.5

16

10

12

M1SD05R18

0.5

18

10

13

M1SD05R20

0.5

20

10

14

M1SD1R10

1

10

10

15

M1SD1R12

1

12

10

16

M1SD1R14

1

14

10

17

M1SD1R16

1

16

10

18

M1SD1R18

1

18

10

19

M1SD1R20

1

20

10

20

M1SD3R20

3

20

10

21

M1SD4R20

4

20

10

22

M1SD5R20

5

20

10

Power Law Distribution

D

ρi (m−2 )

23

D12R10b

1.2

10

6

24

D12R12b

1.2

12

6

25

D12R14b

1.2

14

6

26

D22R10

2.2

10

6

27

D22R15

2.2

15

6

28

D22R20

2.2

20

6

29

D22R25

2.2

25

6

30

D32R20

3.2

20

6

31

D32R25

3.2

25

6

32

D32R30

3.2

30

6

a The average lengths for all the networks with log-normal fracture distribution is set to 1 m b For these scenarios, the fracture network is 10m × 10m, instead of 20m × 20 m as in the rest of the scenarios

tion domain. The same argument can be applied to the power law distribution networks, but now with reference to Fig. 3. Since the focus of this article is not on the domain size, it seems reasonable to keep the extended domain to be twice the length of the simulation domain.

123

Estimating the Hydraulic Conductivity of 2D Fracture

785

Fig. 3 Typical distribution of fracture lengths l for power law distribution, before hydraulically inactive fractures are removed, with power law exponent d = 1.2, 2.2, 3.2

4 Networks of Uniform Aperture Fractures 4.1 Numerical Simulations Numerical simulations have been carried out using NAPSAC (version 9.8.5), a DFN code developed by Serco Technical Consulting Services, Harwell, UK. The code has been extensively verified and used in several fracture network studies; for example, Jackson et al. (2000). Fractures are divided into segments between intersections, and the flowrate through each segment is modelled by a constitutive relation. The cubic law was used for this study, although modifications of the cubic law can also be implemented. At fracture intersections, water pressure is continuous, and water mass is conserved. Fluxes and pressure drops are calculated using the finite element method. For this model, the matrix is impervious, and water flows only through the fractures. NAPSAC calculates the flux through a given fracture system using a number of differently oriented head gradients, chosen to give a uniform coverage of directions, and the conductivity tensor is then fitted to those fluxes. A more detailed explanation of the algorithm is documented in the NAPSAC Technical Summary (Serco 2008). The networks investigated in this study are nominally isotropic, although the calculated permeability tensors usually exhibit very mild anisotropy. The nominal (isotropic) permeability, K NAPSAC , is then defined to be the geometric mean of the two principal macroscopic permeabilities.

123

786

C. T. O. Leung, R. W. Zimmerman

4.2 Relating K to various Parameters of the Network This study hypothesises that K can somehow be correlated with parameters of the fracture network, for example ρ or lmean . The aim of this part of the study is to identify the key parameters. The apertures in this part of the study, h, are set uniformly to 65μm. (Since the flow through each fracture is modelled by the cubic law, the hydraulic conductivity of a network with the same geometry but different aperture can obtained by simple scaling.) A parameter K 0 is defined against which to normalise the conductivity: K0 =

h3 , 12L

(4.1)

where L is the length of the macroscopic region. K 0 represents the conductivity of a fracture network consisting of an orthogonal pair of fractures, with uniform aperture h, each passing straight through the system. For this study, the macroscopic shape of the fracture network is always an L × L square. As a starting point, ρ is plotted against the conductivity obtained from numerical simulations, K NAPSAC (Fig. 4). It is clear from this graph that K cannot be predicted from ρ alone,

Fig. 4 Normalised conductivity as a function of ρ for uniform aperture networks. For a given distribution of fracture lengths, the conductivity increases linearly with ρ. Cases 1–19 correspond to “narrow” lognormal distributions (0 < logσd < 1), cases 20–22 correspond to “broad” lognormal distributions (3 < log σd < 5), and cases 23–32 correspond to power law distributions

123

Estimating the Hydraulic Conductivity of 2D Fracture

787

but as long as the other parameters of the fracture network are kept constant, K increases linearly with ρ. It is reasonable to believe that the length, l, of the fractures will have an effect on the overall conductivity. Sayers and Kachanov (1990) introduced a dimensionless crack density parameter, ε, to study the elastic properties of cracked solids. It is defined by   1  li 2 ε= , (4.2) A 2 i

where A is the area of the network, and li is the length of the ith fracture. This parameter contains information on both ρ and l. Figure 5 compares ε with K /K 0 ; it is seen that the data points show much less scatter than when K is plotted against ρ alone. As with Fig. 4, data points from fracture networks with similar parameters (σln or D) group themselves into straight lines, but the clear improvement compared to using ρ alone is that all the straight lines now pass roughly through the same point. An alternative parameter similar to ε is now proposed. The crack mean density, εmean , can be defined as  n 2  2 n 1  li n l¯ εmean = = , (4.3) A n 2 A 2 i=1

Fig. 5 Normalised conductivity as a function of crack density ε, which is defined in Eq. 4.2, for uniform aperture networks. Note that the conductivities form a linear relationship with ε for a given fracture length distribution, and that the lines for different distributions intercept at the same point

123

788

C. T. O. Leung, R. W. Zimmerman

where n is the total number of fractures, and l¯ is the arithmetic mean of all the lengths. Whereas ε is “the sum of the squares of the fracture lengths”, the parameter εmean can be interpreted as “the square of the mean of fracture lengths”. Figure 6 shows εmean plotted against K . Figure 6 displays similar trends as Fig. 5, in that data points from fracture networks with similar parameters are grouped into straight lines, and the lines seemed to cross at the same point. However, the groupings of data points are tighter when εmean is used instead of ε, with cases 1–19 now collapsed into a single group. It is therefore concluded that εmean is preferred over ε for relating the effects of l and ρ have on K . It is clear, however, that more information is required to establish the full relationship. The segment density, ρseg , is the number of segments, n seg , in the fracture network within a unit area. A segment is the length of fracture between two nodes (fracture intersections). Here n seg is calculated for each fracture network from the following relationship: n seg = n + 2n node ,

(4.4)

where n node is the number of nodes in each fracture network. Equation 4.4 is written from the observation that two new segments are created from each intersection, and the minimum number of segments is the number of fractures when there are no fracture intersections.

Fig. 6 Normalised conductivity as a function of the parameter εmean , which is defined in Eq. 4.3, for uniform aperture networks. Note that when plotted against εmean , the conductivity of fracture networks with the same mean fracture lengths collapse into a single linear relationship. For example, cases 1–19, as well as cases 26–29 (the left most group of cases 23–32) have mean fracture lengths of approximately 1 m. Cases 1–29 are generated using a distribution that specifies mean fracture lengths of 1 m, but due to high spread of fracture lengths in cases 20–22, the actual mean fracture lengths are much higher. In this study, very short fractures tend to be excluded because they are often hydraulically inactive

123

Estimating the Hydraulic Conductivity of 2D Fracture

789

Fig. 7 Normalised conductivity as a function of segment density, ρseg , for uniform aperture networks. ρseg is used as a parameter as an alternative to εmean , because it contains information on connectivity instead of fracture lengths

Figure 7 shows the relationship between ρseg and K /K 0 . As was seen in Figs. 5 and 6, data points for the conductivity for fracture networks generated using similar parameters group into straight lines. Both εmean and ρseg showed promise in describing K , but neither of them tells the full story. Since each contains information not included in the other, perhaps combining the two parameters would give the best description. Hence, the parameter η is defined as √ η = L εmean ρseg . (4.5) A graph of η against K /K 0 (Fig. 8) shows a straight line through the origin, implying direct proportionality between the parameters, i.e., K l¯ √ = η = B nn seg , K0 2L

(4.6)

where B is a dimensionless constant. The overall connectivity of a fracture network, ζ , can therefore be estimated as follows: n node ζ = , n nl¯ . (4.7) η = 1 + 2ζ 2L Substituting this into Eq. 4.6 gives

123

790

C. T. O. Leung, R. W. Zimmerman

Fig. 8 Normalised conductivity as a function of the parameter η, which is defined in Eq. 4.7, for uniform aperture networks. Note the linear relationship, with a slope of 0.185

K nl¯ . = B 1 + 2ζ K0 2L

(4.8)

Equation 4.8 describes fully the relationship for the wide range of networks investigated in this study. From this study B, the dimensionless constant of proportionality, is found to be 9.26 × 10−2 . 4.3 Semi-analytical Method In this section a different method of approximating K is developed, which is used to compare against the numerical values (K NAPSAC ), and to aid the next step of the study, which is to investigate networks with a distribution of fracture apertures. Zimmerman and Bodvarsson (1996) used effective medium theory to develop a model for the effective hydraulic conductivity of two-dimensional fracture networks. Their model can serve as the basis of a semi-analytical approach for estimating the conductivity of a fracture network. In their approach, a two-dimensional fracture network is modelled as a network of conductive segments, where a “segment” is defined as the region of a fracture between two consecutive intersections of that fracture with another fracture. In network terminology, these intersections serve as nodes in the fracture network, and the segments serve as bonds. The coordination number for each node is assumed to be 4, since a coordination number of three corresponds to a fracture that terminates on another fracture, which rarely occurs in the networks studied in this work. With this restriction, the number of segments in a network, n seg ,

123

Estimating the Hydraulic Conductivity of 2D Fracture

791

can be calculated from the number of fractures, n, and the number of nodes/intersections, n node , by Eq. 4.4. The hydraulic conductance C (units: m3 ) of a fracture segment is defined such that the flowrate through the segment is given by Q=

1 C P, μ

(4.9)

where μ is the fluid viscosity, and P is the pressure drop along the segment. The conductance is calculated using the cubic law (Zimmerman and Bodvarsson 1996): C=

h3w , 12lseg

(4.10)

where h is the aperture of the fracture, w is the width of the fracture in the third dimension (set here to unity), and lseg is the characteristic length of the segments (defined below). When all the fractures have the same aperture, the conductances of all the segments are equal. The segments in these fracture networks can be arranged into an orthogonal “checkerboard”-like grid, with the number of segments along the vertical direction equal to the number in the horizontal. The unit cell of such a grid is a block with dimensions lseg × lseg × w. For a square grid, the macroscopic conductivity is simply the conductivity of each unit cell (Fig. 9b). The conductivity of the unit cells, K CB , is calculated using the cubic law, Eq. 4.10, with w set equal to 1: K CB =

h3 . 12lseg

(4.11)

This value can be called the “checkerboard” conductivity. The characteristic segment length, lseg , is a key parameter for calculating this conductivity. Three ways were investigated for calculating lseg : by taking the mean segment length, by balancing the segment density, ρseg , or by using the nodal density, ρnode . The mean segment length is obtained by dividing the total length of the fractures by the total number of segments. The parameter lseg can be defined from ρseg and ρnode using the following argument: in a perfect square lattice network, a unit cell contains two segments and one node, with edges of length lseg , and these parameters can be related by



2 1 = . (4.12) lseg = ρseg ρnode For a square lattice, lseg calculated using either ρseg or ρnode should give the same value, but for a random fracture network this is not the case. Of these three methods for calculating lseg , the value calculated by matching ρnode seems to give the most accurate results. The rearranged fracture network contains both conducting and non-conducting segments: for each fracture there are two end segments that do not conduct (‘dead end segments’). The fraction of the segments that do not conduct, f , is therefore given by f =

2n 2 = . n seg 1 + 2ζ

(4.13)

Application of Kirkpatrick’s effective medium approximation, Eq. 2.2, to a square network containing a fraction f of non-conducting bonds, and a fraction (1 − f ) of bonds having

123

792

C. T. O. Leung, R. W. Zimmerman

Fig. 9 Abstraction and modelling of a random fracture network: a Truncating fractures outside the boundary and removing isolated fractures; b Arranging the fracture network into a regular square lattice to obtain the ‘checkerboard’ conductivity; c Obtaining an effective conductance from the lattice network that was abstracted from the random fracture network

unit conductance, shows that the effect of the non-conducting bonds can be accounted for by multiplying K CB by the factor (1 − 2 f ), i.e., K CB =

h3 (1 − 2 f ). 12lseg

(4.14)

Figure 10a shows a graph of K CB /K 0 plotted against η. The graph shows direct proportionality between the two parameters, with a dimensionless slope of 2.14 × 10−1 . K CB can be used as an upper bound for K for fracture networks with uniform aperture fractures, since it is reasonable to assume that for a fracture network with a given value of ρnode , the highest value of K that the network can have is when all the nodes are arranged

123

Estimating the Hydraulic Conductivity of 2D Fracture

793

(a)

η

(b)

η Fig. 10 a Normalised checkerboard conductivity as a function of η, which is defined in Eq. 4.3. The checkerboard conductivity, K CB , an analytical expression for conductivity proposed in this study, correlates strongly and linearly with η, with slope 2.14 × 10−1 . b The ratio between the numerically calculated and the checkerboard conductivity as a function of η. The ratio converges to α, where α = 0.864 (dotted line)

123

794

C. T. O. Leung, R. W. Zimmerman

into a regular square lattice. Using Figs. 8 and 10a, both of which can be assumed to show a straight line passing through the origin, the relationship between K CB and K is given by K = α K CB ,

(4.15)

where α = 0.864. Figure 10b shows a graph of K /K CB against η, highlighting the value of α. Though there is scatter when η is low (i.e., close to the percolation threshold), α converges to a steady value at higher η. The parameter α represents the error caused by approximating a geometrically “random” fracture network with a checkerboard-like lattice. Fortunately, it seems to take on a value that is nearly universal.

5 Networks of Variable Aperture Fractures Attention is now turned to cases where the apertures are not identical. In this part of the study, the cases in Table 1 will again be used, but instead of identical fracture apertures, the apertures are now taken to be directly proportional to the fracture lengths. Specifically, a 1 m long fracture is taken to have an aperture of 65 μm. As before, the hydraulic conductivities of these networks are calculated numerically using NAPSAC, and then are compared to the conductivities from the approximations. If the apertures of the fractures in the network are not identical, fracture segments will have different conductances. One of the reasons to develop the “checkerboard grid” in the previous section was to establish a base to apply the “effective conductance” approach. The “effective conductance” of a network is defined as the value Cˆ such that, if the conductance of each segment is replaced by Cˆ the overall conductivity of the network does not change (Fig. 9c). The effective conductance is typically estimated from the geometric mean of the conductance distribution, but can also be calculated by methods such as the effective medium theory of Kirkpatrick (1973), i.e., Eq. 2.2, taking z = 4 for a checkerboard geometry: n ˆ  C − Ci = 0. Cˆ + Ci

(5.1)

i=1

Equation 4.14 was first adapted to variable aperture cases, initially using the geometric mean and Kirkpatrick’s approximation to find the effective conductance Cˆ to use in Eq. 4.10 in place of C. The conductivities calculated using these values of K CB were compared with K NAPSAC (not shown in figures). It was found that using the geometric mean or Kirkpatrick’s approximation gives almost identical values of K , but that both of them under-predict the actual conductivity by one or two orders of magnitude. This observation agrees to some extent with the conclusions of David et al. (1990) on applying the EMT of Kirkpatrick to porous medium networks modelled as regular lattices, with a broad distribution of conductances. They found that in such cases, this method tends to underestimate the effective conductivity by a factor of two. Kirkpatrick’s theory, as well as most effective medium approximations, assumes that there is no correlation between the conductances of the segments entering a given node. This is, however, not true for a fracture network, since at each node there are always two pairs of segments having the same conductance, because they are part of the same fracture. Hence, it is not surprising that the geometric mean does not provide a good estimate of the effective conductivity for fracture networks, whereas it works quite well for porous media such as

123

Estimating the Hydraulic Conductivity of 2D Fracture

795

Fig. 11 Comparison of the conductivities predicted using the proposed semi-analytical method, K a , with results from numerical simulations, K NAPSAC . In general, there is good agreement over networks whose conductivities span ten orders of magnitude

consolidated sandstones (Lock et al. 2002), where the assumption that the conductances of the pores that enter a given node are uncorrelated is more accurate. The geometric mean is a special case of power law averaging, a method proposed by Desbarats (1992) to estimate the effective conductivity of a porous medium. Effective conductivity by power law averaging is given by  n 1/ω  1 ω Cˆ = Ci , (5.2) n i=1

where ω = 0 corresponds to the geometric mean, ω = 1 gives the arithmetic mean, and ω = −1 gives the harmonic mean. The initial intention was to match the conductivity as given by Eq. 4.15 with K NAPSAC , recalling the relationship between C and K CB noted in Eq. 4.10 and Eq. 4.11, by varying ω. However, during preliminary studies it was found that using the arithmetic mean to calculate Cˆ gives a very good match to the numerical results, when Cˆ is used in Eq. 4.10 in place of C (Fig. 11). The fact that the effective conductivity of the fracture segments can be approximated by the arithmetic mean can be explained as follows. If the fracture aperture (conductivity) is

123

796

C. T. O. Leung, R. W. Zimmerman

positively correlated with length, then the network will be dominated by a relatively small number of long fractures having high conductivities. These long, highly conductive fractures will act more as if they are in “parallel” with each other than in “series”, in the sense that each of these fractures can transmit fluid across large distances, without the need for fluid to pass through the less conductive fractures of the network. The effective conductance of a set of conductors arranged in parallel is governed by the arithmetic mean rather than the geometric or harmonic mean. As shown in Fig. 11, the conductivity predictions are accurate for networks having conductivities that span many orders of magnitude. At first, it might be thought that the range exhibited in Fig. 11, from 10−14 to 10−3 m2 , is higher than might be expected to be encountered in the field. However, this is merely an artefact of setting the aperture of a one-meter-long fracture to 65 μm. If the aperture of an l = 1 m fracture were taken to be 6.5μm, for example, this would only have the effect of shifting both axes in Fig. 11 downwards by a factor of 10−3 , so that the range would cover 10−17 to 10−6 m2 . The main point is that the prediction method seems to account for vast differences in fracture densities, and for wide distributions of apertures.

6 Conclusions In this study, the effective conductivity of a random fracture network is estimated by representing it as a network of conducting segments between nodes. The method proposed is essentially composed of two parts. First, a relationship is obtained between a random fracture network and one in which all the segments are arranged into a checkerboard grid, when all the segments have the same aperture and hence same conductance, i.e., Eq. 4.14. It is found that the conductivity ratio between a random fracture network and its equivalent checkerboard grid is equal to some constant α, where α = 0.864, which represents the geometry differences between the random and lattice-like networks. The conductivity of a fracture network in which the apertures are correlated with the length can be estimated using the following procedure. First, the fracture intersections (nodes) are identified, and this information is used to remove isolated fractures from the network. Then, use the nodal density of the remaining fractures to create the equivalent checkerboard network to find lseg , Eq. 4.11. Divide each fracture in the actual fracture network into segments of length lseg , and calculate the conductance of each of those segments using the local cubic law, Eq. 4.10. The arithmetic mean of these individual conductances is then used to calculate K CB , from Eq. 4.14. Finally, Eq. 4.15 is used to convert K CB into the effective conductivity. This approach has worked well for all cases considered, regardless of whether the underlying distribution was lognormal or power law. Finally, it should be noted that the methods developed in this article have been applied only to two-dimensional networks formed by the intersection of linear features. Three-dimensional networks formed by the intersection of planar fractures will have a fundamentally different topological character, and it is by not yet clear to what extent, if any, the present methods can be applied to those networks. Acknowledgments This work was funded by the UK Engineering and Physical Sciences Research Council (EPSRC), the UK Nuclear Decommissioning Authority (NDA), and Serco Technical Consulting Services, through CASE studentship award 08002638. The authors thank Cherry Tweed and Simon Norris of the NDA for supporting this work, and thank Andrew Hoch of Serco for providing the NAPSAC code, and guidance in its use. The authors also thank the two anonymous reviewers for their constructive comments.

123

Estimating the Hydraulic Conductivity of 2D Fracture

797

References Adler, P.M., Berkowitz, B.: Effective medium analysis of random lattices. Transp. Porous Media 40, 145–151 (2000) David, C., Gueguen, Y., Pampoukis, G.: Effective medium theory and network theory applied to the transport properties of rock. J. Geophys. Res. 95, 6993–7005 (1990) de Dreuzy, J.-R., Davy, P., Bour, O.: Hydraulic properties of two-dimensional random fracture networks following a power law length distribution 1. Effective connectivity. Water Resour. Res. 37, 2065–2078 (2001a) de Dreuzy, J.-R., Davy, P., Bour, O.: Hydraulic properties of two-dimensional random fracture networks following a power law length distribution 2. Permeability of networks based on lognormal distribution of apertures. Water Resour. Res. 37, 2079–2095 (2001b) Desbarats, A.J.: Spatial averaging of hydraulic conductivity in three-dimensional heterogeneous porous media. Math. Geol. 24, 249–267 (1992) Hestir, K., Long, J.C.S.: Analytical expressions for the permeability of random two-dimensional Poisson fracture networks based on regular lattice percolation and equivalent media theories. J. Geophys. Res. 95(B13), 21565–21581 (1990) Jackson, C.P., Hoch, A.R., Todman, S.: Self-consistency of a heterogeneous continuum porous medium representation of a fractured medium. Water Resour. Res. 36(1), 189–202 (2000) Kirkpatrick, S.: Percolation and conduction. Rev. Mod. Phys. 45(4), 574–588 (1973) Lock, P.A., Jing, X.D., Zimmerman, R.W., Schlueter, E.M.: Predicting the permeability of sandstone from image analysis of pore structure. J. Appl. Phys. 92, 6311–6319 (2002) Long, J.C.S., Witherspoon, P.A.: The relationship of the degree of interconnection to permeability in fracture networks. J. Geophys. Res. 90, 3087–3098 (1985) Min, K.-B., Jing, L., Stephansson, O.: Determining the equivalent permeability tensor for fractured rock masses suing a stochastic REV approach: method and application to the field data from Sellafield, UK. Hydrogeol. J. 12, 497–510 (2004) Sayers, C.M., Kachanov, M.: A simple technique for finding effective elastic constants of cracked solids for arbitrary crack orientation statistics. Int. J. Solids Struct. 27(6), 671–680 (1990) Serco TAS: NAPSAC Technical Summary, Release 9.6. Serco, Harwell, Didcot (2008) Snow, D.T.: Anisotropic permeability of fractured media. Water Resour. Res. 5, 1273–1289 (1969) Tóth, T.M., Vass, I.: Relationship between the geometric parameters of rock fractures, the size of percolation clusters and REV. Math. Geosci. 43, 75–97 (2011) van Golf-Racht, T.D.: Fundamementals of Fractured Reservoir Engineering. Elsevier, Amsterdam (1982) Wu, Y.S., Haukwa, C., Bodvarsson, G.S.: A site-scale model for fluid and heat flow in the unsaturated zone of Yucca Mountain, Nevada. J. Contam. Hydrol. 38, 185–215 (1999) Zhang, X., Sanderson, D.J.: Numerical study of critical behaviour of deformation and permeability of fractured rock mass. Marine Petrol. Geol. 15, 535–548 (1998) Zimmerman, R.W., Bodvarsson, G.S.: Effective transmissivity of two-dimensional fracture networks. Int. J. Rock Mech. Min. Sci. 33(4), 433–488 (1996)

123