The use of discrete fracture network simulations in the design of horizontal hillslope drainage networks in fractured rock

Engineering Geology 163 (2013) 132–143

Contents lists available at ScienceDirect

Engineering Geology journal homepage: www.elsevier.com/locate/enggeo

The use of discrete fracture network simulations in the design of horizontal hillslope drainage networks in fractured rock Donald M. Reeves a,⁎, Rishi Parashar a, Greg Pohll a, Rosemary Carroll a, Tom Badger b, Kim Willoughby b a b

Division of Hydrologic Sciences, Desert Research Institute, 2215 Raggio Parkway, Reno, NV 89512, United States Washington State Department of Transportation, P.O. Box 47365, Olympia, WA 98504, United States

a r t i c l e

i n f o

Article history: Received 7 September 2012 Received in revised form 17 May 2013 Accepted 26 May 2013 Available online 7 June 2013 Keywords: Hillslope drainage Hillslope stability Rock fractures Drainage networks Fracture networks

a b s t r a c t Characteristics of fracture networks are explored in a discrete fracture network framework to provide design guidelines for horizontal drainage networks in fractured rock. Central to the study is defining how fracture attributes relate to fracture network structure and network-scale fluid flow, and in turn, how flow characteristics of fracture networks influence horizontal drain length and orientation. Multiple realizations of stochastic fracture networks, generated from both synthetic and field-specific data sets, serve as a basis for understanding physical fracture network structure and resultant global flow and for performing intersection analyses of hillslope drains with flowing fractures. Study results indicate that the logarithm of the standard deviation of fracture transmissivity, log(σT), is the single most important attribute for drainage network design, as higher values of log(σT) describe heterogeneous flow patterns where only a small portion of the network conducts a significant quantity of fluid. Thus recommended drain lengths for intersecting significantly conductive fractures increase with increases in log(σT). Fracture trace length, orientation, and density also play a role, albeit secondary to the distribution of transmissivity, in defining drain length as a function of drain orientation relative to the mean fracture set orientation. The spatially discontinuous nature of fracture networks and the wide range in transmissivity values found in natural fracture networks tend to produce high degrees of variability in computed intersection distances between drains and fractures conducting significant quantities of fluid. To account for this variability, a conservative approach is recommended where horizontal drain lengths along a pre-defined orientation are scaled by discrete fracture network computed intersection distances equal to the upper 95th confidence interval. © 2013 Elsevier B.V. All rights reserved.

1. Introduction The presence of water is one of the most critical factors contributing to hillslope instability. The installation of horizontal drains and drain networks is common practice used to decrease the elevation of the water table surface. Lowering of the watertable dries a large portion of the hillslope, which increases the shear strength of the soil and decreases the probability of slope failure. Effectiveness of hillslope drainage is typically described in terms of the increase in the factor of safety, defined as the ratio of shear strength to shear stress, once horizontal drains are installed. The need for design guidelines for horizontal drains used to promote hillslope stability has been noted by several researchers (Choi, 1974; Kenney et al., 1977; Prellwitz, 1978; Nonveiller, 1981; Forrester, 2001). Drainage system design is most developed for irrigated areas (Donnan, 1946; Israelsen, 1950; Maasland, 1956; Kirkham, 1958; Talsma and Haskew, 1959; U.S. Department of the Interior, 1978) with analytic solutions focusing on simplified conditions including idealized slopes, ⁎ Corresponding author. Tel.: +1 775 673 7605; fax: +1 775 673 7363. E-mail address: [email protected] (D.M. Reeves). 0013-7952/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.enggeo.2013.05.013

constant pre-drain groundwater levels and recharge, homogeneous hydraulic conductivity, and regularly-spaced drains that are typically parallel to the slope (Hooghoudt, 1940; Bouwer, 1955; Schmid and Luthin, 1964; Wooding and Chapman, 1966; Childs, 1971; Towner, 1975; U.S. Department of the Interior, 1978; Lesaffre, 1987; Ram and Chauhan, 1987; Fipps and Skaggs, 1989). Many of the limitations of the analytic solutions, such as complex hillslope geometry, transient water levels and recharge, heterogeneous hydraulic conductivity fields and complex drain configurations, can be relaxed through the use of numerical models of hillslope drainage and stability (e.g., Angeli et al., 1998; Cai et al., 1998; Rahardjo et al., 2003; Pathmanathan, 2009). Numerical models are also useful for quantifying the sensitivity of design parameters to overall hillslope drain performance and resultant factor of safety. Previous research on hillslope drainage and stability has primarily focused on soil and unconsolidated sediment where flow occurs through the interconnected pore space of the medium. Many watersheds are underlain by fractured rock that are either directly involved in slope stability or are host aquifers that can adversely influence the stability of overlying soils (Al Homoud and Tal, 1997; Montgomery et al., 1997; Jiao et al., 2005; Furuya et al., 2006; Gerscovich et al., 2006; Ghosh et al., 2009). In these cases, drainage networks installed within the

D.M. Reeves et al. / Engineering Geology 163 (2013) 132–143

overlying soil are likely to be ineffective to sufficiently lower pore pressures in or reduce recharge from the underlying fractured rock mass. Installation of a drainage network within the fractured bedrock must then be considered. Fractured rock presents a very different challenge to the design of hillslope drainage networks. This is because fractured rock typically has little or negligible primary porosity and permeability in the rock matrix itself, and connected networks of discontinuous fractures impart secondary porosity and permeability that govern the flow of groundwater. Unlike porous media where flow occurs at the porescale, flow in fractured rock systems occurs through complex patterns of interconnected, conductive fractures (e.g., Long et al., 1982; Smith and Schwartz, 1984; Renshaw, 1999; de Dreuzy et al., 2001; Berkowitz, 2002; Neuman, 2005; Reeves et al., 2008a,b; Klimczak et al., 2010; Reeves et al., 2010). Hence, the design of horizontal hillslope drainage networks in fractured rock must take into account the characteristics of the fracture networks to maximize the probability for drains to intersect fractures conducting significant amount of flow that, in turn, will sufficiently reduce pore pressures. In this paper, we provide guidelines for the design of horizontal hillslope drainage networks for fractured rock by utilizing a discrete fracture network framework to explore how specific fracture attributes relate to physical network structure, how networks conduct flow given various geometric configurations and degrees of heterogeneity in hydraulic properties, and how network-scale flow characteristics influence horizontal drain orientation and length. Stochastic fracture networks generated from both synthetic and field-specific data sets serve as a basis for understanding physical network structure and global flow, and for performing intersection analyses of hillslope drains with fractures conducting a significant amount of flow. Emphasis is placed on linking the most relevant fracture attributes with horizontal hillslope drainage network design. 2. Properties of fractured rock Bedrock typically has little or no primary porosity and permeability, and networks of fractures serve as primary conduits for fluid flow. These networks are spatially discontinuous and highly irregular in geometry and hydraulic properties. The variability in geometric and hydraulic properties is the result of the complex interplay between current and past stress fields, inhomogeneous rock mechanical properties (e.g., Young's modulus, Poisson's ratio), mechanical fracture interaction and distributions of flaws or weakness in a rock mass. Full characterization of fractured rock masses is not possible since known fracture locations and their attributes consist of an extremely small sample of the overall fracture network, i.e., any fracture characterization effort grossly undersamples a field site due to limited accessibility to the fractures themselves. Fractured rock masses are typically characterized during field campaigns that measure fracture attributes from a number of sources including boreholes (e.g., electrical resistivity, ultrasonic, optical televiewers, etc.), rock outcrops, road cuts, tunnel complexes, seismic images and hydraulic tests. These fracture data can then be used to generate representative, site-specific discrete fracture networks (DFN) through the derivation of probabilistic descriptions of fracture location, orientation, spacing, length, aperture, hydraulic conductivity/transmissivity and values of network density (e.g., Figure 1). Mapped fractures are often included into DFN models as deterministic features, though the addition of stochastic fractures honoring site-specific fracture data is always necessary to maintain proper fracture density and network connectivity. The limited accessibility to the network leaves an incomplete understanding of the patterns of fracturing within a rock mass that can often be improved through visual inspection of representative networks generated according to site-specific statistics. Relationships between fracture statistics, network structure, and network-scale flow are explored in this section. Fracture data are only discussed in the

133

context of common probability distributions used to describe variability of specific fracture attributes found in natural fractured rock systems. Readers are encouraged to refer to Munier (2004) and Reeves et al. (2012) for additional detail on rock fracture statistical analysis. 2.1. Network structure Network structure is defined as the complex geometrical and hydraulic configuration of interconnected fracture segments as influenced by natural variability in fracture orientation, length, transmissivity and density. A defining characteristic of natural fracture networks is trace length which is most often a power-law distributed variable in natural fracture networks: P ðL > lÞ ¼ CL

−a

ð1Þ

where the power law exponent a lies between 1 and 3 (e.g., Davy, 1993; Bour and Davy, 1997; Renshaw, 1999; Bonnet et al., 2001). Note that average fracture length decreases as values of a increase. Networks typically consist of at least two or more fracture sets grouped by orientation, and fracture length for each set may or may not have the same value of a. The distribution of orientation about a mean set orientation is typically modeled using a Fisher distribution (Fisher, 1953): pðθÞ ¼

κ⋅sinθ⋅eκ⋅cosθ eθ −e−θ

ð2Þ

where θ (degrees) is symmetrically distributed ð−π2 ≤θ≤π2Þ according to a constant dispersion parameter, κ. The extent to which individual fractures cluster around the mean orientation is described by κ where higher values of κ describe higher degrees of clustering. The spacing between fractures is typically found to have a negative exponential trend which can be replicated using a Poisson point process (Ross, 1985). Hydraulic properties of fractures, whether aperture, hydraulic conductivity, or transmissivity, generally follow either log-normal trends (Stigsson et al., 2001; Andersson et al., 2002a,b) or power law trends (Gustafson and Fransson, 2005; Reeves et al., 2008a). Here we use a log-normal distribution to assign fracture transmissivity independently to individual fracture segments, T: " # 2 1 −ðlogT−logðμ T ÞÞ   pðT Þ ¼ qffiffiffiffiffiffiffiffiffiffiffiffi exp 2log σ 2T T 2πσ 2T

ð3Þ

where μT is the mean and σT is the standard deviation. Values of log(σT) are commonly around 1 for fractured media (Stigsson et al., 2001; Andersson et al., 2002a,b), which describes variability in transmissivity for individual fractures encompassing 5 to 6 orders of magnitude. Distributions of hydraulic conductivity and transmissivity can also be inferred from distributions of aperture values using the well-known cubic law (Snow, 1965). The generation of fracture networks involves adding individual fractures with random location, orientation, length, and transmissivity into a finite domain until a density criterion is satisfied. This density criterion can vary from one-dimensional to three-dimensional measures, and a two-dimensional fracture density criterion, defined as: ρ2D = A−1 ∑ ni = 1 li, is used in this work to relate the sum of fracture lengths, li, to two-dimensional domain area, A. The ρ2D criterion is equivalent to the P21 density criterion used in some discrete fracture network studies (e.g., Munier, 2004). Synthetic data consisting of two fracture sets with power law distribution of lengths with exponent values in the range 1 ≤ a ≤ 3, moderate fracture density, mean orientations of ± 45° with variability described by a Fisher distribution with κ = 20, and a log-normal transmissivity distribution with log(σT) = 1 are used to generate three different network types (Figure 2). Once a network is generated, the hydraulic backbone is identified by eliminating dead-end

134

D.M. Reeves et al. / Engineering Geology 163 (2013) 132–143

100

100

90

90

80

80

70

70

60

60

50

50

40

40

30

30

20

20

10

10

0

0

10

20

30

40

50

60

70

80

90 100

0 0

10

20

30

40

50

60

70

80

90 100

Fig. 1. Three-dimensional discrete fracture network (DFN) generated according to two fracture sets with significant variability about mean fracture orientations, a power-law length distribution exponent of a = 2, and a relatively sparse density. The two dimensional network at the bottom left is computed by projecting all fractures onto the horizontal slice located in the center of the three-dimensional DFN (i.e., the two dimensional networks are tracemaps of the three dimensional DFN). The two-dimensional network on the bottom right is the result of identifying the hydraulic backbone by eliminating all dead-end fracture segments and non-connected clusters. From Reeves et al. (2012).

segments and isolated clusters. This is accomplished in our model using both geometric and flow techniques (Parashar and Reeves, 2012). The hydraulic backbone represents the interconnected subset of a fracture network that is responsible for conducting all flow and 100

100

a=1

90

transport across a domain. Hence, analysis of backbone characteristics can provide insight into these processes. The generated networks in this section do not explore the full parameter space for fractured media. The wide range of fracture length

a=2

90

100 90

80

80

80

70

70

70

60

60

60

50

50

50

40

40

40

30

30

30

20

20

20

10

10

10

0 0

10 20 30 40 50 60 70 80 90 100

0

0 0

10 20 30 40 50 60 70 80 90 100

a=3

0

10

20

30

40

50

60

70

80

90 100

Fig. 2. Two-dimensional hydraulic backbones for networks with a = 1, a = 2 and a = 3, respectively, from left to right. Backbone densities, defined as the sum of backbone fracture length normalized by domain area, are very similar for direct comparison of network structure as a function of a. From Reeves et al. (2012).

D.M. Reeves et al. / Engineering Geology 163 (2013) 132–143

exponents, however, provide sufficient variability and produce three distinct types of hydraulic backbones. Networks generated with a = 1 produce backbone structures dominated by long fractures, and networks with a = 3 produce backbone structures dominated by short fractures (Figure 2). Backbones with a mixture of short and long fractures are produced for networks generated with a = 2 (Figure 2). Another feature of these networks is that density of the network increases from 1.0 m/m2 to 2.0 m/m2 as the value of a increases from 1 to 3. These fracture network densities are arbitrarily assigned; however, the increase in density with increasing values of a is necessary to maintain a percolating backbone that promotes fluid flow from one side of the domain to the other. For example, the density values assigned to a = 1 and a = 2 (ρ2D = 1.0 m/m2 and ρ2D = 1.5 m/m2, respectively) result in non-percolating networks if used with a = 3 (Table 1). Conversely, networks generated with a = 1 and ρ2D = 2.0 m/m2 (density assigned to a = 3) produce unrealistically dense networks (Table 1). The relationship between network density, power-law exponent for trace length, and backbone percolation will be further discussed shortly. 2.2. Network-scale flow To our knowledge, trends between fracture length, density, and the total amount of fluid flow a network conducts have not been studied. There has been some work, however, that relates these same fracture parameters to contaminant transport behavior (Reeves et al., 2008a,b,c; Zhang et al., 2010). This lack of knowledge on fluid flow can, in part, be explained by the fact that each natural fracture network exhibits unique combinations of fracture attributes, making it difficult to compare total flow through networks with differing distributions of length and transmissivity and values of density. Nonetheless, understanding how geometric and hydraulic properties of fracture sets influence overall network scale flow is helpful in identifying network types that may or may not be suitable for installation of drainage networks. We conducted a numerical investigation to understand how total network-scale flow may be related to backbone structure, i.e., primarily fracture trace length and density. A total of 50 network realizations were generated on a 100 m by 100 m domain and solved for flow for each of the network types shown in Fig. 2 (Table 1). Fracture networks are highly heterogeneous and the number of realizations needed to explore the full range in network-scale flow responses may quite possibly be in the range of tens of thousands; however, we limited the analysis to 50 realizations to provide a general comparison of total flow for different network types. Computation of flow in twodimensional discrete fracture networks involves solving for hydraulic head at all internal nodes (intersection point of two or more fractures) inside the domain according to Darcy's law (Priest, 1993; Klimczak et al., 2010; Parashar and Reeves, 2012). Constant gradient boundary conditions are applied to the networks to induce flow from top to bottom, and fluid flow is solved using iterative techniques for large and sparse linear systems to obtain the steady-state flow field (Parashar and Reeves, 2012). This process is repeated for each statistically equivalent realization of a given network type. Each network type was solved for flow with two different distributions of transmissivity. The first transmissivity distribution assigns all fractures the same value of T = 1.0 × 10−5 m2/s to facilitate

135

direct comparisons between fracture trace length distribution and overall network-scale flow, Q [m3/s]. The second distribution of transmissivity assigns a log-normal distribution with log(σT) equal to 1 to individual fractures, which provides a good representation of variability found in natural fracture networks. Flow through these networks is denoted by Q σ T ¼0 and Q σ T ¼1 respectively in Table 1. Analysis of the DFN data yields several insights. As previously discussed, network percolation requires greater fracture densities as values of a increase. This can be observed in Table 1, where values of ρ2D for a = 1 and a = 2 are 1.0 and 1.5 m/m2, respectively, yet produce nearly the same backbone density of 0.70 m/m2. A value ρ2D = 2.0 m/m2 for a = 3 networks supports this trend, albeit with a slightly higher backbone density of 0.83 m/m2. The explanation for this trend is simple: mean fracture length decreases with increases in fracture length exponent values, requiring more fractures to reach a sufficient degree of network connectivity. The use of a constant T value for all fracture segments along with comparable backbone densities indicates that fracture network flow decreases as values of a increase (Figure 3, log σT = 0 case). Implementation of a log-normal T distribution with log(σT) = 1 increases network scale flow for a = 1 networks by approximately an order of magnitude (Figure 3, σT = 1 case). This dramatic increase in network-scale flow results from domain-spanning fractures in a = 1 networks, which provide continuous pathways for fluid flow from one end of the domain to the other, being assigned T values greater than the mean. Network-scale flow for a = 2 and a = 3 is similar. This reflects the higher densities of these networks where T values lower than the mean serve as “bottlenecks” which have a restrictive influence on total flow. Variability in network-scale flow for the constant T simulations was evaluated according to minimum, maximum and mean flow to study overall network connectivity and backbone structure for individual realizations. Networks with a = 2 showed the greatest variability with up to 90% differences in flow from the mean listed in Table 1. This indicates that these networks exhibit the greatest variability in backbone structure and connectivity, and this is most likely caused by the influence of both short and long fractures in producing complex backbone structures (Figure 2). Networks with a = 3 have up to 45% differences in flow from the mean listed in Table 1. These networks are dominated by short fractures, and differences in flow are attributed to the degree of connectivity between individual realizations. Networks with a = 1 exhibited the lowest degree of variability with only up to 20% differences in flow from the mean listed in Table 1. Networks with a = 1 are dominated by very long, often domain-spanning fractures that lead to similar degrees of connectivity for individual realizations. Variability in network-scale flow for the log-normal T (log σT = 1) shows similar trends where a = 1 networks have the lowest variability (up to a factor of 2 from the mean) and a = 2 and a = 3 networks have higher variability (up to a factor of 3 from the mean). Results from the above analysis assume that fracture length and transmissivity are uncorrelated properties. This is common practice for hydrological investigations (e.g., Smith and Schwartz, 1984; de Dreuzy et al., 2001; Reeves et al., 2008a). There is some evidence from the rock fracture mechanics literature that length and aperture, and hence transmissivity, may be correlated for mechanically isolated joints (Schultz et al., 2008; Klimczak et al., 2010), although the correlation of length and aperture for mechanically interacting fractures

Table 1 Discrete fracture network parameters and results from 50 network realizations. a

ρ2D [m/m2]

Fractures

Backbone ρ2D [m/m2]

Q σ T ¼0 [m2/s]

STDEV Q σ ¼0 T [m2/s]

Q σ T ¼1 [m2/s]

STDEV Q σ ¼1 T [m2/s]

1 2 3

1.0 1.5 2.0

1116 3898 6757

0.72 0.70 0.83

2.81 × 10−3 1.69 × 10−3 1.32 × 10−3

2.60 × 10−4 2.25 × 10−4 1.95 × 10−4

2.39 × 10−2 4.10 × 10−3 1.85 × 10−3

1.71 × 10−2 3.80 × 10−3 2.11 × 10−3

136

D.M. Reeves et al. / Engineering Geology 163 (2013) 132–143

1.00E-01

and it is recommended that specific drainage projects perform these analyses to maximize the probability for success.

1.00E-02

3.1. Differentiating fracture types

Q [m2/s] 1.00E-03

1.00E-04 0.5

1

1.5

2

2.5

3

3.5

a Fig. 3. Network simulations with comparable backbone densities indicated that total fracture network flow decreases with increasing values of a. This trend linearly decreases when all fracture segments have the same values of T (log(σT) = 0). Application of a realistically parameterized log-normal distribution with log(σT) = 1 further increases network-scale flow for all cases, with the largest increase in flow for the a = 1 networks.

(joints and faults) in a network is less clear. The transmissivity of mechanically-interacting fractures within a stress field can be altered through the combined effects of shear dilation and mechanical aperture dilation and closure (e.g., Min et al., 2004; Chen et al., 2007; Zhang et al., 2007; Nemoto et al., 2009). The topic of how stress fields precisely influence the correlation between fracture length and aperture at the network-scale is beyond the scope of this work and is left to future research. 3. Fracture intersection analysis Fracture networks created by tectonic stress, particularly in extensional regimes, consist mostly of fractures that are oriented vertical to sub-vertical (e.g., Zoback, 2010), and drain networks are oriented horizontal to sub-horizontal. The two-dimensional DFN framework introduced in the previous section is used to evaluate the intersection of these mutually orthogonal features. Note that the intersection of a three-dimensional DFN with, in this case, a horizontal plane results in a two-dimensional network (Figure 1). Here we assume that the drainage networks are perfectly horizontal, and hence, the intersection analyses for two-dimensional and three-dimensional networks are equivalent. Intersection of horizontal hillslope drainage networks is determined by projecting domain spanning scanlines across the networks and recording intersections of three fracture types: all fractures, backbone fractures and dominant fractures. These three fracture types will be defined shortly. The number of intersections is then normalized by the scanline length to provide an average distance of intersection r as a function of angle θ for each of the three fracture types. Fracture networks are spatially discontinuous features with non-random orientations, and the scanline orientation θ is rotated from 0 to 360° in 5 degree increments to fully capture the directional nature of intersection distances. This process is repeated for 100 DFN flow realizations and ensemble (the average of all realizations) trends are computed. Distances of scanline intersection with dominant fractures represent horizontal drain lengths necessary for successful drainage of a rock mass. These distances may or may not exhibit directional dependence as illustrated in the following examples. A systematic analysis using both synthetic and site-specific fracture data is conducted to understand the influence of fracture orientation, length, density and transmissivity on ensemble intersection trends. Differences in intersection trends between single fracture network realizations and ensemble-averaged trends are both qualitatively and quantitatively defined. The following sections provide information on links between fracture network properties and hillslope drainage network design. It must be emphasized that all rock masses are unique,

Fracture networks are highly heterogeneous and standard deviations of transmissivity on the order of log(σT) = 1 are typical for natural fractures. This level of heterogeneity in transmissivity corresponds to 5–6 orders of magnitude. The discrete fracture network simulations show that a subset of all the fractures belong to the hydraulic backbone and conduct flow, and an even smaller subset of these hydraulic backbone fractures conduct the majority of flow across a network for log(σT) = 1. This finding is consistent with field observations that indicate only 10% (or less) of the total fracture population contributes significantly to flow (e.g., Long et al., 1991), and necessitates the grouping of fractures into different categories to compute intersection distances according to fracture type. Discrete fracture networks are analyzed for intersection with three fracture types: all fractures, hydraulic backbone fractures, and dominant fractures. “All fractures” refer to all fractures present in a rock mass, whereas “hydraulic backbone” fractures refer only to the interconnected fractures of the hydraulic backbone. These two fracture types – all fractures and backbone fractures – are physical properties of a network assuming that all fractures within the DFN are permeable to flow. “Dominant” fractures, as defined in this subsection, refer to fractures of the hydraulic backbone that conduct flow exceeding a specified value: FD≥

GF NB =2

ð4Þ

where flow through a dominant fracture, FD, must be equal to or exceed the global flow through the corresponding DFN realization, GF, normalized by the number of boundary fractures, NB. Global flow is defined as the total flow that passes through the fracture network for a given domain size. Boundary fractures intersect the DFN model boundary and are responsible for all flow in and out of the model domain. Half of the boundary fractures are used because global DFN flow represents the net flow through a network, and on average approximately half of the boundary fractures conduct flow into the network and approximately half of the boundary fractures conduct flow out from the network. The use of a robust mathematical definition of dominant fractures in Eq. (4) is necessary for comparing networks with different degrees of heterogeneity in transmissivity, defined by the log variance of transmissivity, log(σT) (Figure 4). In the following intersection analysis, “successful” hillslope drainage attributed to lowering of pore pressures depends on the intersection of dominant fractures with horizontal hillslope drains. The criterion of a hillslope drain intersecting a dominant fracture is conservative. Dominant fractures consist of a very small subset of the total fracture population that conducts the majority of flow through the network. Assuming that the conductance of the drain network is sufficiently high, intersection of a dominant fracture with a hillslope drain will, by definition, always dewater the surrounding fracture network and sufficiently reduce pore pressures. However, there may be some cases where a drain network may not intersect a dominant fracture, yet enough conductive fractures of the hydraulic backbone are intersected by the drain network such that pore pressures are sufficiently lowered. These cases would occur infrequently at a fractured rock site and are intentionally excluded from the analysis. 3.2. Influence of fracture length and transmissivity Discrete fracture networks are generated according to synthetic parameter sets to allow for systematic study of specific fracture attributes on intersection distances. Since network connectivity and global flow

D.M. Reeves et al. / Engineering Geology 163 (2013) 132–143

137

Fig. 4. Ranked plots of flow through individual fracture segments of the fracture backbone along with the identification of dominant fractures using Eq. (4) (left). Flow through individual fracture segments of the hydraulic backbone with line thickness proportional to segment flow (right). The plots are made from networks with identical geometric properties (a = 2, θ = ± 25°, κ = 25, ρ2D = 1.5 m/m2) with log(σT) = 0 (top) and log(σT) = 1 (bottom). Note that networks with homogeneous distributions of T (log(σT) = 0) contain many dominant fractures (top), whereas networks with heterogeneous distributions of T (log(σT) = 1) exhibit focused flow through a small subset of fractures of the hydraulic backbone (bottom).

are highly dependent on distributions of fracture length and transmissivity, the influence of these parameters on intersection distances is explored first. Six sets of discrete fracture networks consisting of a = 1, a = 2 and a = 3, with standard deviation of transmissivity described by log(σT) = 0 (homogeneous) and log(σT) = 1 (heterogeneous) are studied. All of these networks consist of two fracture sets with orientations of ±45° and deviations in orientation described by κ = 20. Radial plots of distance r as a function of angle θ describe the distance of intersection for all fractures, backbone fractures, and dominant fractures (Figure 5). The top half of Fig. 5 from left to right describe the influence of these fracture trends with increasing values of a given log(σT) = 0. For the most part, trends of intersection for all fracture types are similar in distance, with the exception that the intersection distance of all fracture and hydraulic backbone groups are very similar for the a = 1 networks. This indicates that even though the a = 1 networks have a larger fracture spacing in these examples, many of the fractures are connected to the hydraulic backbone, because the a = 1 length distribution favors longer, domain-spanning fractures. The effect of deviations in transmissivity about the mean value is readily observed during visual inspection of the top (log(σT) = 0) and bottom (log(σT) = 1) plots. Heterogeneity in transmissivity greatly extends intersection distances for the dominant fractures from approximately 3–4 m to 6–10 m along the Cartesian coordinate axes. This reflects the trends shown in Fig. 5 where increasing heterogeneity effectively concentrates flow through a smaller subset of

fractures, increasing the spacing between dominant fractures. Larger intersection distances with increases in log(σT) indicate that the number of total fractures intersected by a horizontal drain on average needs to effectively double (approximately 3 to 6) to ensure that a dominant fracture is intersected. As a final comment, the square shape of average intersection distances signifies a directional dependence on hillslope drain length. If the radial distance to dominant fractures were equal, then the overall shape of these intersection trends would be described by a circle of fixed radius as shown in the next example. Rather, the square shape indicates that the distance to intersection with all types rangespfrom ffiffiffi a minimum along the Cartesian coordinate axes to a factor of 2 at the diagonals. This pattern reflects the mean fracture set orientations of ±45°, where the minimum distance represents the angles of maximum intersection by the two fracture sets (0°, 90°, 180°, 270°), and the maximum distance represents the diagonals, where intersection of the scanline is reliant on only one fracture set as the other fracture set is parallel to the scanline. 3.3. Influence of fracture set orientation and density Previous intersection analyses have focused on networks containing two fracture sets with mean orientation of ±45°. In this section, networks with a = 2 and ρ2D = 1.5 m/m2 are utilized to study the influence of fracture set orientation and density assigned to individual fracture sets on resultant intersection distances.

138

D.M. Reeves et al. / Engineering Geology 163 (2013) 132–143

a

d

b

c

e

f

Fig. 5. Ensemble-averaged directional distances of scanline intersection with all fractures (red), backbone fractures (green) and dominant fractures (blue) for a = 1 (a,d), a = 2 (b,e) and a = 3 (c,f). Values of r denote radial distance of average fracture intersection distance (meters) to each fracture type according to angle θ. All networks contain two fracture sets with orientations of ±45° with each set having equal density. Plot axes all have the same scale for consistent visual comparison. Note how intersection distances for all fractures and hydraulic backbone fractures are more or less consistent for values of a, and that intersection distances for dominant fractures are highly dependent on values of log(σT) = 0 (top, a–c) and log(σT) = 1 (bottom, d–f).

Intersection trends with fracture set orientation of ±45° are square with the largest distances to fracture intersections coinciding with mean fracture orientation (Figure 5). Changing mean fracture set orientation from ± 45° to ± 90° effectively rotates the square shape to a diamond where the largest intersection distances of fracture intersection coincide with mean fracture orientation (Figure 6a). Changing values of κ used to describe variability in fracture orientation about the mean was observed to only slightly alter overall intersection trends (not shown). The use of a random fracture orientation, however, results in circular intersection trends devoid of

a

b

directional dependence (Figure 6b). In this case, the design of a hillslope drainage network would not need to preferentially scale drain length as a function of θ. Trends of intersection plots, whether square, circular or diamond in shape, all exhibit symmetry with respect to fracture orientation. This is specifically caused by the use of two fracture sets with orthogonal orientations (or random in the case of Figure 6b) and assigning an equal number of fractures to each set. These fracture network parameters were intentionally held constant to investigate the influence of fracture length and transmissivity. However, prior probability, defined as the ratio of fractures applied

c

Fig. 6. Influences of orientation and density on ensemble-averaged intersection trends with all fractures (red), backbone fractures (green) and dominant fractures (blue) for networks with two fracture sets oriented at ± =90°. The previous intersection analyses in Fig. 5 contained networks with two fracture sets oriented at ±45°. For these networks, intersection patterns were square and the largest distances to fracture intersections coincide with mean fracture orientation. The use of two fracture sets oriented at ±90° (left) effectively rotates the square shape to a diamond where the greatest intersection distances coincide with mean fracture orientation. Use of random fracture orientation (center) results in circular intersection plots absent of directional dependence. Unequal weights assigned to each fracture set can greatly skew directional intersection trends (right), where 20% and 80% of fractures are assigned to the set with orientation −45° and 45°, respectively.

D.M. Reeves et al. / Engineering Geology 163 (2013) 132–143

to each fracture set, is commonly unequal for natural fracture networks. The effect of unequal distribution of fractures among each fracture set is evaluated by assigning 20% and 80% of fractures to the −45° and 45° sets, respectively (Figure 6c). This effectively yields a fracture spacing that is 4 times lower along the 45° fracture set than the −45° fracture set. This type of unequal fracture distribution results in strongly skewed and asymmetric plots with the smallest intersection distances along −45° (due to intersecting the 45° set with higher density), and largest intersection distances along 45° (due to intersecting −45° set with a more sparse density). 3.4. Deviations between individual realizations and the ensemble Thus far, the analysis has focused on ensemble-averaged trends from a total of 100 individual fracture network realizations. Spatial configuration and structure of fracture networks are highly uncertain, and this is reflected in the fracture generation methodology itself. To study deviations of individual realizations with the ensemble trends, individual realizations (light gray) are plotted against the ensemble (dark blue) in Fig. 7 for networks with a = 1, θ = ± 45°, κ = 20, and log(σT) = 1. This is the same network shown with ensemble trends in Fig. 5d. Variability of intersection distances to dominant fractures for 100 individual realizations about the ensemble is large, with the majority

a

b

c

d

139

of realizations clustered around the ensemble mean and a small subset of individual realizations showing extreme variability (Figure 7a, b). For example, some scanlines for specific network realizations intersect only 1–2 dominant fractures, which generate intersection distances in excess of 40–50 m. No dominant fractures were intersected along a given angle θ for a very small subset of realizations. This finding reflects that only a very small number of dominant fractures are present in the simulated fracture networks which is consistent with observations from natural fractured rock flow systems. These dominant fractures emerge in the flow simulations through a complex interaction between fracture connectivity, length, and transmissivity. Thus, the observed large or undefined intersection distances to dominant fractures for some of the fracture network realizations reflects what one may expect at a field site. To restrict variability to a quantitative limit, 95%-confidence intervals are computed from the 100 intersection distance trends (Figure 7c,d). The most important confidence interval for the purpose of hillslope drainage network design is the 95th percentile, as this metric defines the distance necessary to intersect at least one dominant fracture for 95% of the simulations. This is the distance that will be used later on to guide directional length of hillslope drains for the site-specific networks. It is apparent from examination that the 95th percentile is not fully formed and its shape is inconsistent with the square shape exhibited by the 5th percentile and ensemble mean (Figure 7c). The

Fig. 7. Variability of intersection distances to dominant fractures for 100 individual realizations about the ensemble-averaged trends for a = 1, κ = 20, θ = ± 45°, and logðσ T Þ ¼ 1. This is the same network simulations shown in Fig. 5d. These plots show a high degree of variability about the mean (a,b). To provide quantitative estimates of uncertainty about the ensemble average, 95% confidence intervals are computed from the individual realizations (bottom). The plots indicate that the 5th percentile and mean trends are well-developed and follow a square trend, yet the 95th percentile is irregular due to incomplete sampling of the fracture network parameter space. The modified 95th percentile shown in (d) was developed by taking the mean of the irregular 95th percentile over all θ, and rescaling the mean trend. This method worked well in preserving the amplitude and frequency of the 95th percentile, as shown by the good match to the theoretical trend for a perfect square (dotted line).

140

D.M. Reeves et al. / Engineering Geology 163 (2013) 132–143

shape of the 95th percentile should be approximately square with a larger length scale than the mean. The 95th percentile trend in Fig. 7c is a result of undersampling from 100 realizations, and quite likely that several thousand realizations would be required to fully form the 95th percentile trend. This high number of realizations is not computationally feasible for this analysis. Instead, we compute a ‘modified’ 95th percentile trend by combining properties of the raw 95th percentile and the ensemble curve. First, the intersection distance of the 95th percentile is averaged over all θ. This average distance is then multiplied by the normalized ensemble trend. This process generates the modified 95th percentile

Table 2 Statistics of fractures sets for Climax granite stock.

Prior probability Mean strike Mean dip Dispersion (κ)

Set 1

Set 2

Set 3

Set 4

Set 5

Set 6

Set 7

0.03 125 19 65

0.13 317 25 37

0.10 360 85 33

0.14 321 83 24

0.13 289 82 23

0.32 48 80 18

0.15 N/A N/A N/A

a a

b

b

c

Fig. 8. Intersection distances for the LCA3 fracture networks for all fracture types (a) along with trends of distance to dominant fractures with ensemble mean and 95% confidence intervals (b,c). Note how the non-orthogonal fracture sets and non-equal number of fractures in each set leads to skewed trends.

c

Fig. 9. Intersection distances for the Climax granite stock fracture networks for all fracture types (a) along with trends of distance to dominant fractures for ensemble mean and 95% confidence intervals (b,c). Note that the presence of multiple fracture sets with several unique orientations creates intersection trends that are nearly circular.

D.M. Reeves et al. / Engineering Geology 163 (2013) 132–143

a

b

c

d

141

Fig. 10. Discrete fracture network realizations and resultant hillslope drainage network for LCA3 (a,b) and Climax (a,b) field sites. Drain length as a function of θ is based on the 95th percentile trend. Note that the drainage network is more symmetric for Climax than for the LCA3, and that hillslope drain lengths are approximately three times larger for the Climax site (~30 m) than for LCA3 site (~10 m).

curve (red) in Fig. 7d that is remarkably similar to the trend of a perfect square (dotted purple). This check validates that the constructed 95th percentile trend correctly reproduces the square shape observed in the 5th percentile and ensemble-mean plots.

intersection distances for all fracture types are nearly circular (Figure 9), indicating that the result of multiple fracture sets at the Climax site is a near uniform distribution of orientation. 4. Summary of guidelines for horizontal drainage network design

3.5. Site specific fracture networks The directional-intersection analysis method is applied to two sets of site-specific fracture data in this subsection. The first field site is the Lower Carbonate Aquifer (LCA3) located on the Nevada National Security Site (NNSS), and the second is the Climax Fractured Granite Stock also located on the NNSS. To briefly summarize, the LCA3 fracture networks are very dense (ρ2D = 6.9 m/m2), with two sets of non-orthogonal fractures with strikes of θ1 = 227° and θ2 = 292° and prior probabilities equal to 0.64 and 0.36, respectively. This translates to networks with two non-orthogonal fracture sets separated by 65° and an unequal distribution of fractures between the sets, where one set contains approximately two-thirds of the fractures. These factors result in skewed intersection distance trends (Figure 8). The Climax stock is an intrusive, densely-fractured (ρ2D = 4.3 m/m2) granite body on the NNSS that was characterized and analyzed for flow and transport properties by Reeves et al. (2010). This rock mass was subjected to multiple fracturing events under various stress-field configurations, which reflects the multiple groups of fracture orientation (Table 2). All fracture sets were assigned fracture lengths according to a power-law distribution with exponent a = 1.6, and transmissivity according to a lognormal distribution with log(σT) = 1. The resultant

The information contained in this paper is intended to assist in the design of hillslope drainage networks in fractured rock. The first step in designing a hillslope drainage network in fractured rock is to collect site-specific fracture data. These field data can then be analyzed to determine distributional properties of orientation, location, length, aperture and hydraulic conductivity/transmissivity and values of fracture density. Generation of discrete fracture networks from this data provides valuable insight into the structure and flow characteristics of the underlying fracture networks. An analysis of the distance of intersection with three defined fracture types systematically identified the influences of fracture length, transmissivity, orientation and density on directional intersection trends. The overall shape of intersection trends are highly sensitive to mean fracture orientation, where minimum and maximum distances to two sets of orthogonal fractures coincide with the directions most perpendicular and parallel to the mean fracture set orientations, respectively. All examples, regardless of fracture network attributes, indicate that many fractures need to be intersected on average before a dominant fracture is crossed. The spatially discontinuous nature of fracture networks and large variability in transmissivity indicates that a high degree of variability

142

D.M. Reeves et al. / Engineering Geology 163 (2013) 132–143

exists between individual realizations and ensemble averages. For example, the sampling scanline for a very small subset of individual fracture network realizations did not intersect a dominant fracture. To account for this variability, a conservative approach where horizontal drain length is scaled to the distance defined by the 95th percentile is recommended. This concept is illustrated by Fig. 10, where the length of individual drains is scaled to the directional dependence of the 95th percentile trend for the LCA3 and Climax granite stock fracture statistics. It is important to note that scaling drain lengths based on the 95% percentile of dominant fracture intersection distance does not guarantee that a dominant fracture will be intersected at all drain installations. We emphasize however that this guideline serves as both a realistic and conservative metric to scale drain lengths. A result of this analysis not yet discussed is the possibility that some sparsely fractured networks may require drain lengths (or perhaps, drain densities) that are not economically feasible for successful dewatering. As an example, average drain length for the Climax site is approximately 30 m. This calls into question the practice of installing horizontal hillslope drains prior to a rigorous analysis of fracture network statistics. It is highly recommended that the intersection analysis contained in this section be conducted prior to the design and installation of drainage networks in fractured media. Acknowledgments DMR, RP, GP, and RC would like to thank Washington State Department of Transportation (WSDOT) for project funding. This manuscript is the result of a collaborative study with WSDOT as well as Department of Transportations from 10 additional states. This work is currently being compiled into a manual that presents guidelines for horizontal drainage networks in both porous and fractured media utilizing principles of hydrology, site characterization, and numerical models of hillslope drainage and stability. The authors would also like to thank the helpful comments from two anonymous reviewers which improved the quality of the manuscript. References Al Homoud, A.S., Tal, A.B., 1997. Engineering geology of the Aqaba Ras El Naqab highway re-alignment, Jordan. Engineering Geology 47 (1–2), 107–148. Andersson, P., Byegärd, J., Dershowitz, B., Doe, T., Hermanson, J., Meier, P., Tullborg, E.L., Winberg, A., 2002a. Final report on the TRUE Block Scale Project: 1. Characterization and model development. Technical Report TR-02-13. Swedish Nuclear Fuel and Waste Management Company (SKB), Stockholm, Sweden. Andersson, P., Byegärd, J., Winberg, A., 2002b. Final report on the TRUE Block Scale Project: 2. Tracer tests in the block scale. Technical Report TR-02-14. Swedish Nuclear Fuel and Waste Management Co. (SKB), Stockholm, Sweden. Angeli, M.G., Buma, J., Gasparetto, P., Pasuto, A., 1998. A combined hillslope hydrology/ stability model for low-gradient clay slopes in the Italian Dolomites. Engineering Geology 49 (1), 1–13. Berkowitz, B., 2002. Characterizing flow and transport in fractured geological media: a review. Advances in Water Resources 25 (8–12), 861–884. Bonnet, E.O., Bour, O., Odling, N., Davy, P., Main, I., Cowie, P., Berkowitz, B., 2001. Scaling of fracture systems in geologic media. Reviews of Geophysics 39 (3), 347–383. Bour, O., Davy, P., 1997. Connectivity of random fault networks following a power law fault length distribution. Water Resources Research 33 (7), 1567–1583. Bouwer, H., 1955. Tile drainage of sloping fields. Agricultural Engineering 36, 400–403. Cai, F., Ugai, K., Wakai, A., Li, Q., 1998. Effects of horizontal drains on slope stability by three-dimensional finite element analysis. Computers and Geotechnics 23 (4), 255–275. Chen, Y., Zhou, C., Sheng, Y., 2007. Formulation of strain-dependent hydraulic conductivity for a fractured rock mass. International Journal of Rock Mechanics and Mining Sciences 44, 981–996. Childs, E.C., 1971. Drainage of groundwater resting on a sloping bed. Water Resources Research 7 (5), 1256–1263. Choi, Y.L., 1974. Design of horizontal drains. Journal of Engineering Society 12, 37–49. Davy, P., 1993. On the frequency-length distribution of the San Andreas fault system. Journal of Geophysical Research 98 (B7), 12,141–12,151. de Dreuzy, J.R., Davy, P., Bour, O., 2001. Hydraulic properties of two-dimensional random fracture networks following a power-law length distribution: 1. Effective connectivity. Water Resources Research 37 (8), 2065–2078. Donnan, W.W., 1946. Model tests of a tile-spacing formula. Proceedings of the Soil Science Society of America 11, 131–136.

Fipps, G., Skaggs, W., 1989. Influence of slope on subsurface drainage of hillsides. Water Resources Research 25 (7), 1717–1726. Fisher, R., 1953. Dispersion on a sphere. Proceedings of the Royal Society of London Series A 217, 295–305. Forrester, K., 2001. Subsurface Drainage for Slope Stabilization. ASCE Press, Reston. Furuya, G., Suemine, A., Sassa, K., Komatsubara, T., Watanabe, N., Marui, H., 2006. Relationship between groundwater flow estimated by soil temperature and slope failures caused by heavy rainfall, Shikoku Island, southwestern Japan. Engineering Geology 85 (3–4), 332–346. Gerscovich, D.M.S., Vargas, E.A., de Campos, T.M.P., 2006. On the evaluation of unsaturated flow in a natural slope in Rio de Janeiro, Brazil. Engineering Geology 88 (1–2), 23–40. Ghosh, A., Sarkar, S., Kanungo, D.P., Jain, S.K., 2009. Slope instability and risk assessment of an unstable slope at Agrakhal, Uttarakhand. Proceedings of the India Geotechnical Conference, Guntur, India, pp. 783–786. Gustafson, G., Fransson, A., 2005. The use of the Pareto distribution for fracture transmissivity assessment. Hydrogeology Journal 14, 15–20. http://dx.doi.org/10.1007/ s10040-005-0440-y. Hooghoudt, S.B., 1940. General consideration of the problem of field drainage by parallel drains, ditches, watercourses, and channels. Publ. No.7 in the Series Contribution to the Knowledge of Some Physical Parameters of the Soil (Titles Translated from Dutch). Bodemkundig Instituut, Groningen, The Netherlands. Israelsen, O.W., 1950. Irrigation Principles and Practices. John Wiley and Sons, New York. Jiao, J.J., Wang, X.-S., Nandy, S., 2005. Confined groundwater zone and slope stability in weathered igneous rocks in Hong Kong. Engineering Geology 80 (1–2), 71–92. Kenney, T.C., Pazin, M., Choi, W.W., 1977. Design of horizontal drains for soil slopes. Journal of Geotechnical Engineering 103 (11), 1311–1323. Kirkham, D., 1958. Seepage of steady rainfall through soil into drains. Transactions of the American Geophysical Union 39, 892–908. Klimczak, C., Schultz, R.A., Parashar, R., Reeves, D.M., 2010. Cubic law with correlated aperture to length and implications for network scale fluid flow. Hydrogeology Journal. http://dx.doi.org/10.1017/s10040-009-0572. Lesaffre, B., 1987. Analytical formulae for traverse drainage of sloping lands with constant rainfall. Irrigation and Drainage Systems 1, 105–121. Long, J.C.S., Remer, J.S., Wilson, C.R., Witherspoon, P.A., 1982. Porous media equivalents for networks of discontinuous fractures. Water Resources Research 18 (3), 645–658. Long, J.C.S., Karasaki, K., Davey, A., Peterson, J., Landsfeld, M., Kemeny, J., Martel, S., 1991. An inverse approach to the construction of fracture hydrology models condition by geophysical data: an example from the validation exercises at the Stripa Mine. International Journal of Rock Mechanics and Mining Science and Geomechanics Abstracts 28 (2–3), 121–142. Maasland, M., 1956. The relationship between permeability and the discharge, depth, and spacing of tile drains. Bulletin No. 1, Groundwater and Drainage Series. Water Conservation and Irrigation Commission, New South Wales, Australia. Min, K.-B., Rutqvist, J., Tsang, C.-F., Jing, L., 2004. Stress-dependent permeability of fractured rock masses: a numerical study. International Journal of Rock Mechanics and Mining Sciences 41, 1191–1210. Montgomery, D.R., Dietrich, W.E., Torres, R., Anderson, S.D., Heffner, J.T., Loague, K., 1997. Hydrologic response of a steep, unchanneled valley to natural and applied rainfall. Water Resources Research 33 (1), 91–109. Munier, R., 2004. Statistical analysis of fracture data adopted for modeling discrete fracture networks — version 2. Rep. R. 04-66. Swedish Nuclear Fuel and Waste Management Company (SKB), Stockholm, Sweden. Nemoto, K., Watanabe, N., Hirano, N., Tsuchiya, N., 2009. Direct measurement of contact area and stress dependence of anisotropic flow through rock fracture with heterogeneous aperture distribution. Earth and Planetary Science Letters 281, 81–87. http://dx.doi.org/10.1016/j.epsl.2009.02.005. Neuman, S.P., 2005. Trends, prospects and challenges in quantifying flow and transport through fractured rocks. Hydrogeology Journal 13, 124–147. http://dx.doi.org/ 10.1007/s10040-004-0397-2. Nonveiller, E., 1981. Efficiency of horizontal drains on slope stability. Proceed. 10th Int. Conf. Soil Mech. Found. Eng. 3, Stockholm, Sweden, pp. 495–500. Parashar, R., Reeves, D.M., 2012. On iterative techniques for solving flow in large twodimensional discrete fracture networks. Journal of Computational and Applied Mathematics 236 (18), 4712–4724. Pathmanathan, M.L., 2009. Numerical Simulation of the Performance of Horizontal Drains for Subsurface Slope Stabilization. M.S. Thesis Washington State University. Prellwitz, R.W., 1978. Analysis of parallel drains for highway cut-slope stabilization. Proceed. 16th Ann. Idaho Eng. Geol. Soils Eng. Symp, Boise, Idaho, pp. 153–180. Priest, S.D., 1993. Discontinuity Analysis of Rock Engineering. Chapman and Hall, London. Rahardjo, H., Hritzuk, K.J., Leong, E.C., Rezaur, R.B., 2003. Effectiveness of horizontal drains for slope stability. Engineering Geology 69, 295–308. Ram, S., Chauhan, H.S., 1987. Drainage of sloping lands with constant replenishment. Journal of Irrigation and Drainage Engineering 113 (2), 213–223. Reeves, D.M., Benson, D.A., Meerschaert, M.M., 2008a. Transport of conservative solutes in simulated fracture networks: 1. Synthetic data generation. Water Resources Research 44, W05401. http://dx.doi.org/10.1029/2007WR006069. Reeves, D.M., Benson, D.A., Meerschaert, M.M., 2008b. Transport of conservative solutes in simulated fracture networks: 2. Ensemble solute transport and the correspondence to operator-stable limit distributions. Water Resources Research 44, W05410. http://dx.doi.org/10.1029/2008WR006858. Reeves, D.M., Benson, D.A., Meerschaert, M.M., 2008c. Influence of fracture statistics on advective transport and implications for geologic repositories. Water Resources Research 44, W08405. http://dx.doi.org/10.1029/2007WR006092.

D.M. Reeves et al. / Engineering Geology 163 (2013) 132–143 Reeves, D.M., Pohlmann, K., Pohll, G., Ye, M., Chapman, J., 2010. Incorporation of conceptual and parametric uncertainty into radionuclide flux estimates from a fractured granite rock mass. Stochastic Environmental Research and Risk Assessment. http:// dx.doi.org/10.1007/s00477-010-0385-0. Reeves, D.M., Parashar, R., Zhang, Y., 2012. Hydrogeologic characterization of fractured rock masses intended for disposal of radioactive waste. In: Rahman, R.O.A. (Ed.), Radioactive Waste. InTech Publishing (ISBN 979-953-51-0551-0). Renshaw, C.E., 1999. Connectivity of joint networks with power law length distributions. Water Resources Research 20 (9), 2661–2670. Ross, S.M., 1985. Introduction to Probability Models, 3rd ed. Academic Press, Orlando, FL. Schmid, P., Luthin, N.J., 1964. The drainage of sloping lands. Journal of Geophysical Research 69 (8), 1525–1529. Schultz, R.A., Soliva, R., Fossen, H., Okubo, C., Reeves, D.M., 2008. Dependence of displacement–length scaling relations for fractures and deformation bands on the volumetric changes across them. Journal of Structural Geology 30, 1405–1411. http://dx.doi.org/10.1016/j.jsg.2008.08.001. Smith, P.A., Schwartz, F.W., 1984. An analysis of the influence of fracture geometry on mass transport in fractured media. Water Resources Research 20 (9), 335–348. Snow, D.T., 1965. A Parallel Plate Model of Fractured Media. Ph.D. Dissertation University of California, Berkeley.

143

Stigsson, M., Outters, N., Hermanson, J., 2001. Äspö Hard Rock Laboratory, Prototype Repository Hydraulic DFN Model no. 2. IPR-01–39. Swedish Nuclear Fuel and Waste Management Company (SKB), Stockholm, Sweden. Talsma, T., Haskew, H.C., 1959. Investigation of water-table response to tile drains in comparison with theory. Journal of Geophysical Research 64 (11), 1933–1944. Towner, G.D., 1975. Drainage of groundwater resting on a sloping bed with uniform rainfall. Water Resources Research 11 (1), 144–147. U.S. Department of the Interior, 1978. Drainage manual. Water Resources Technical Publication (Washington, D.C., 286 pp.). Wooding, R.A., Chapman, T.J., 1966. Groundwater flow over a sloping impermeable layer. 1: application of the Dupuit–Forchheimer assumption. Journal of Geophysical Research 71 (12), 2895–2902. Zhang, J., Standifirid, W.B., Roegiers, J.-C., Zhang, Y., 2007. Stress-dependent fluid flow and permeability in fractured media: from lab experiments to engineering applications. Rock Mechanics and Rock Engineering 40 (1), 32-21. http://dx.doi.org/ 10.1007/s00603-006-0103-x. Zhang, Y., Baeumer, B., Reeves, D.M., 2010. A tempered multiscaling stable model to simulate transport regional-scale fractured media. Geophysical Research Letters 37, L11405. http://dx.doi.org/10.1029/2010GL043609. Zoback, M.D., 2010. Reservoir Geomechanics. Cambridge University Press, Cambridge, New York.