Normal-stress dependence of fracture hydraulic properties including two-phase flow properties Hui-Hai Liu & Ming-Yao Wei & Jonny Rutqvist Abstract A systematic approach has been developed for determining relationships between normal stress and fracture hydraulic properties, including two-phase flow properties. The development of a relationship between stress and fracture permeability (or fracture aperture and fracture closure) is based on a two-part Hooke’s model (TPHM) that captures heterogeneous elastic-deformation processes at a macroscopic scale by conceptualizing the rock mass (or a fracture) into two parts with different mechanical properties. The developed relationship was verified using a number of datasets in the literature for fracture closure versus stress, and satisfactory agreements were obtained. TPHM was previously shown to be able to accurately represent testing data for porous media as well. Based on the consideration that fracture–aperture distributions under different normal stresses can be represented by truncated-Gaussian distributions, closed-form constitutive relationships were developed between capillary pressure, relative permeability and saturation, for deformable horizontal fractures. The usefulness of these relationships was demonstrated by their consistency with a laboratory dataset. Keywords Fractured rocks . Hydraulic properties . Multiphase flow . Hydro-mechanical processes . Fracture apertures

Introduction Mechanical deformation processes in fractured rock and their coupling to hydrological processes are important for many applications (e.g., van Golf-Racht 1982; Bandis et al. 1983; Rutqvist and Stephansson 2003; Tsang et al. 2012), including geothermal energy development, oil and Received: 30 May 2012 / Accepted: 1 October 2012 Published online: 10 November 2012 * Springer-Verlag Berlin Heidelberg 2012 H.-H. Liu ()) : M.-Y. Wei : J. Rutqvist Earth Sciences Division, Lawrence Berkeley National Laboratory, Mail Stop 84-171, Berkeley, CA 94720, USA e-mail: [email protected] Fax: +510-486-5686 Hydrogeology Journal (2013) 21: 371–382

gas extraction, nuclear waste disposal, geological sequestration of carbon dioxide, and deep-well injection of liquid and solid wastes. With the significant advancement of computer technology in recent decades, numerical models have been increasingly employed for evaluating coupled hydro-mechanical processes associated with these applications (e.g., Jing and Hudson 2002). Constitutive relationships are fundamental for modeling coupled hydro-mechanical processes, because they determine how coupled processes are actually coupled in a numerical simulation. An important constitutive relationship for fractured rock is the stress dependence of fracture hydraulic properties. In this study, the focus was on the relationship between normal stress and fracture hydraulic properties; the latter is also a strong function of shearing (e.g., Bandis et al. 1983). Note that fracture permeability and fracture apertures (or closures) are closely related through cubic law and its variations (e.g., Zimmermann and Bodvarsson 1996). The stress dependence of fracture hydraulic properties (e.g., permeability, aperture or closure) has been investigated by a number of researchers. Goodman (1974) proposed an empirical hyperbolic relationship between normal stress and fracture closure that was later modified by Bandis et al. (1983) to better fit observations. Other models, based on Hertzian theory, have also been developed to describe the non-linear stress-deformation behavior (McDermott and Kolditz 2006). These models suggested that the observed non-linear behavior could be attributed to the increasing contact areas as the normal stress increases. For example, Brown and Scholz (1985; 1986) used Hertzian theory to study fracture closure as a function of normal stress. Based on the Hertzian theory of deformation of spheres, Gangi (1975) developed a relationship for stress-dependent fracture permeability. Most recently, Liu et al. (2009) proposed a two-part Hooke’s model (TPHM) for describing the relationship between stress and elastic strain. The main idea of the TPHM is to capture heterogeneous deformation processes at a macroscopic scale, resulting from the existence of heterogeneity of rock mass, by conceptualizing the rock mass (or a fracture) into two parts with different mechanical properties. While the TPHM has been validated for different rock types, evaluation of its validity for fractures is limited (Liu et al. 2009; Zhao and Liu 2012; Liu et al. 2011). The first objective of this study was to DOI 10.1007/s10040-012-0915-6

372

perform a more comprehensive evaluation of the TPHM using data on fracture deformation gathered from the literature. Multiphase flow occurs in a number of practical applications associated with fractured rock such as geological sequestration of CO2 and nuclear waste disposal. However, most of the relevant studies available in the literature have focused on constitutive relationships for single-phase flow properties such as permeability (and fracture aperture or closure). A systematic study of the relationships between multiphase flow properties for fractures (such as relative permeabilities for different phases) and normal stress is still lacking in the literature. Even without considering stress-dependence, development of constitutive relationships for multiphase flow is a challenging task for fractures. There are currently two kinds of approaches for developing the related constitutive relationships. One is a porous-medium approach, in which constitutive relationships (between capillary pressure, relative permeability and saturation) developed for porous media are simply borrowed or modified for either single fractures or fracture networks (e.g., Therrien and Sudicky 1996). However, the physical meaning of relevant parameters in these relationships is not always clear for fractures, because the geometry (and connectivity) of fracture apertures is essentially two dimensional on the fracture plane, whereas that of pores in porous media is three-dimensional. The second approach, called the pipe-flow approach, has been developed based on similarities between the observed multiphase-flow behavior within fractures and pipes, including the fact that both wetting and non-wetting phases are discontinuous and yet mobile (Wong et al. 2008; Weerakone et al. 2011), or stratified owing to density differences within horizontally inclined fractures (Fourar and Lenormand 1998; Indraratna et al. 2002). This approach may be valid when fluid-phase distributions are not strongly dependent on fracture apertures (or capillary force). From laboratory test results for nitrogen/water flow within single horizontal fractures, Diomampo (2001), however, did not observe pipe-flow behavior. It is obvious that different regimes may exist for multiphase flow within fractures. The focus here is on two-phase flow within a horizontal fracture when flow paths are mainly determined by capillarity (or the aperture distribution). The second objective of this study was to develop new and closed-form relationships between two-phase flow properties and normal stress for horizontal fractures. This paper is organized as follows: in section Fracture aperture or closure under different normal stresses, a comprehensive evaluation of the validity of the TPHM for describing the stress-dependence of fracture aperture or closure is performed, using datasets from public literature; in section Two-phase flow properties, closed-form relationships of two-phase flow properties are proposed for deformable horizontal fractures. Their usefulness is partially demonstrated by comparing calculation results and test data. Hydrogeology Journal (2013) 21: 371–382

Fracture aperture or closure under different normal stresses Fracture aperture is an important parameter for both mechanical and hydraulic processes within a fractured rock; for example, flow processes primarily occur in fractures. In this section, a relationship is developed for the dependence of fracture aperture (or closure) on the normal stress, based on the TPHM. Also note that the term stress, as used in this paper, refers to effective stress, to take account of the effect of pore-liquid pressure.

A brief description of the TPHM For the sake of completeness, the TPHM of Liu et al. (2009) is briefly discussed here. Liu et al. (2009) argued that the true strain, rather than the engineering strain, should be used in Hooke’s law for accurately modeling elastic deformation, unless the two strains are essentially identical (as they might be for small mechanical deformations). In terms of volumetric strain, the true strain refers to volume change divided by rock volume at the current stress state, and engineering strain refers to volume change divided by the unstressed rock volume. In the literature of rock mechanics and other related scientific areas, however, engineering strain is now used exclusively. Liu et al. (2009) further argued that natural rocks are inherently heterogeneous and thus their mechanical deformation is better described using a heterogeneous system. To deal with this issue, they conceptualize the heterogeneous rock as having two parts, and hypothesize that one part (a portion of pore volume or fracture apertures) obeys true-strain-based Hooke’s law, and the other part approximately follows the engineering-strainbased Hooke’s law, because its deformation is small (as previously indicated, true strain is practically identical to the engineering strain when the deformation is small). For simplicity, the first part is called the “soft” part and the other called the “hard” part. This conceptualization can be represented by the hypothesized composite spring system shown in Fig. 1. These two springs are subject to the same stress, but follow different variations on Hooke’s law. Mavko and Jizba (1991) also considered rock porosity to consist of a soft part and a stiff part, when studying grainscale fluid effects on velocity dispersion in rocks. The “two-part” conceptualization presented here is generally consistent with this previous study. Based on these two considerations, Liu et al. (2009) and Zhao and Liu (2012) developed stress-strain σn

σn

Hard Spring

Soft Spring

Fig. 1 Schematic diagram of a composite spring system consisting of two springs. The hard and soft springs follow engineering-strainbased and natural-strain-based Hooke’s law, respectively (Liu et al., 2009) DOI 10.1007/s10040-012-0915-6

373

relationships and a number of related constitutive relationships for porous media. Remarkable consistency between theoretical and observed results has been obtained for a variety of rock samples. In this study, the TPHM was applied to fracture deformation. For a given fracture, the “soft” part corresponds to a range of apertures whose relative changes with stress are significant for low normal stress (or small contact areas), and the “hard” part corresponds to the range of apertures associated with stabilized contacting areas as the stress changes. The validity of this approximation is evaluated in this paper.

Relationship between normal stress and fracture aperture or closure A derivation of the relationship between normal stress and fracture aperture, based on the TPHM, was given in Liu et al. (2009). The derivation procedure will now be reviewed and extended to include fracture closure. Herein, following Liu et al. (2009), the subscripts 0, e, and t are used to denote the unstressed state, the hard part (where the engineering-strain-based Hooke’s law applies) and the soft part (where the true-strain-based Hooke’s law applies), respectively. Consider a fracture to be subject to a normal stress σn and divide fracture space into hard and soft parts along the direction normal to the fracture plane. Then, the volumetrically averaged fracture aperture (b) is given by: b0 0b0;e þ b0;t

ð1Þ

under unstressed conditions, and b0be þ bt

ð2Þ

under stressed conditions. Hooke’s law for the two parts can be expressed by ds n 0  KF;e

ds n 0  KF;t

dbe b0;e

ð3Þ

dbt bt

ð4Þ

where K is the fracture modulus and subscript F refers to fracture. Note that the stress in the two equations above refers to far-field normal stress, rather than local stress. Combining Eqs. (1)–(4) gives db0dbe þ dbt 0  b0;e

ds n ds n  bt KF;e KF;t

ð5Þ

Integrating the above equation and using Eq. (1) and the following relationship obtained from Eq. (4):
 sn bt 0b0;t exp  ð6Þ KF;t Hydrogeology Journal (2013) 21: 371–382

the following equation can be obtained:
 sn sn b0b0;e ð1  Þ þ b0;t exp  KF;e KF;t

ð7Þ

Note that in the derivation of Eq. (7), the condition used is that for σn 00, be 0b0,e and bt 0b0,t. In Eq. (7), the stress-dependent behavior of fracture aperture is controlled by the second term at low stress and the first term at high stress. However, in many laboratory tests, fracture closure, rather than aperture, is measured. The fracture closure (χ) as a function of normal stress can be derived from Eq. (7) as follows:
 σn σn c0b0  b0b0;e ð Þ þ b0;t ½1  exp   ð8Þ KF;e KF;t Equations (7) and (8) give relationships between fracture aperture, closure and normal stress. In the aforementioned equations, fracture aperture b is a mechanical aperture, rather than the hydraulic fracture aperture that determines fracture permeability through cubic law (Witherspoon et al. 1980; Olsson and Barton, 2001). These two apertures are generally different, except for the case of smooth fractures. Barton et al. (1985) and Zimmerman and Bodvarsson (1996) proposed mathematical expressions to relate these two apertures that are functions of fracture roughness and contact area. When information on fracture roughness and contact area is available or can be estimated, the stress-dependence of fracture permeability can be estimated through the cubic law and one of these relations between mechanical and hydraulic apertures. A simpler approach is to use the ratio of mechanical apertures at different stresses to approximate the corresponding ratio of hydraulic apertures for a given fracture. In this case, the cubic law leads to:
2 k b 0 k0 b0

ð9Þ

where k is the fracture permeability corresponding to fracture aperture b, and again subscript 0 refers to the zero-stress condition. Note that cubic law means that the flow rate of water through a fracture is proportional to the cube of the hydraulic aperture (that implies Eq. 9 for fracture permeability). This treatment seems to be supported by several studies. For example, Elliott and Brown (1988) experimentally show that flow rate through a fracture is approximately proportional to the cube of mechanical aperture for a range of aperture values (15–21 μm). It is obvious from the cubic law that their results essentially indicate the proportionality of mechanical aperture to hydraulic aperture under the corresponding test conditions. Liu et al. (2009, 2011) successfully match two datasets of fracture permeability as a function of stress using Eq. (9). Therefore, the preceding equation can be used as a first-order approximation for practical applications. DOI 10.1007/s10040-012-0915-6

374

Finally, it is useful to further emphasize that the TPHM has two key elements. Firstly, it hypothesizes that in Hooke’s law, true strain, rather than engineering strain, should be used. However, because the hard part is subject to a small deformation, engineering and true strains are practically identical for the hard part. In this case, the engineering strain, for the mathematical simplicity, is used to describe the mechanical deformation of the hard part. Secondly, the TPHM conceptualizes the rock mass (or a fracture) to consist of two parts, hard and soft. This is largely motivated by the well-known heterogeneity of subsurface material. The existence of the soft part in porous media, in terms of micro-cracks and pore space near the grain contacts, has been discussed in previous investigations (e.g., Mavko and Jizba 1991). As previously discussed, the soft and hard parts for fracture apertures are those corresponding to relatively small and large contact areas within fractures, respectively. Nevertheless, the TPHM, as a macroscopic model that is not derived from micro-mechanics, should be ultimately evaluated by how consistent it is with a large range of observations.

Comparisons with experimental observations and discussion To evaluate the usefulness of Eqs. (7) and (8), it was examined whether these equations can satisfactorily match experimental observations from different sources. Since most laboratory data are for fracture closures, the evaluation focused on Eq. (8). A relatively comprehensive literature survey of the closure-stress data was performed. The selection of the datasets that were to be used for the comparisons was based on several considerations. Firstly, the dataset and its measurement procedure must be well documented, such that there is a good understanding of the dataset and the possible errors involved in the measurements. Secondly, as suggested by Barton et al. (1985), the closure-stress measurements may be representative of in situ fracture behavior from the third or fourth loading cycles. Therefore, the analysis was focused on the data obtained from these cycles, or data corresponding to the largest cycle numbers when measurements from the third or fourth loading cycles are not available. Table 1 presents information on the datasets used in this study. Readers are referred to the original data sources listed in Table 1 for the details of data-collection procedures. The literature survey was by no means exhaustive, but Table 1 should include typical datasets for closurestress measurements available in the public literature. To give some examples, Fig. 2 shows comparisons between theoretical results from Eq. (8) and measured data from two experiments. The comparisons are typical for all the other experiments. To avoid nonuniqueness in the parameter estimation, values for parameters in Eq. (8), for a given dataset, were Hydrogeology Journal (2013) 21: 371–382

determined in the following manner: in a closure-stress plot (Fig. 2), data at high stress can be approximately represented by a linear relationship; based on Eq. (8), the slope of the linear relationship is b0,e /KF,e and its intersection with the closure axis is b0,t; then KF,t can be estimated from a data point at a low stress. In order to estimate b0,e, additional information on either KF,e or total aperture b0 under unstressed conditions are needed. In some cases, the information is not available (Table 1). This, however, does not impact on the evaluation of the TPHM, because Eq. (8) only requires the ratio b0,e/KF,e , rather than b0,e, to calculate a fracture closure. Figure 2 and Table 1 indicate that the TPHM adequately represents experimental observations. As expected, significant portions of the fracture aperture (or closure) are soft and characterized by high nonlinearity, as shown in Fig. 2, and the modulus for the soft part is significantly lower than that of the hard part. Note that for most datasets collected for different rock types and by different researchers (Table 1), the values for the correlation coefficient of curve fitting (R2) are above 0.95. It is important to mention that the TPHM, as previously discussed, is a macroscopic model that deals with mechanisms of micro-mechanics in a phenomenological manner. Its validity is demonstrated by its consistence with a number of datasets (Table 1). The fitted values for KF,t are relatively small, because they characterize fracture deformation at low normal stress corresponding to small contact areas within fractures. The values for KF,e are much larger, because KF,e characterizes fracture deformation corresponding to relatively stabilized contact areas as stress changes. The upper limit of KF,e should be of the same order of magnitude as that for the rock matrix. Note that being a relatively simple model, the TPHM has two unique aspects in describing rock deformations. Firstly, it was developed by extending the stress-strain relationship described by Hooke’s law (the TPHM is reduced to the conventional form of Hooke’s law when the soft-part portion is zero). Thus, it can be used to derive relevant constitutive relationships for hydraulic and mechanical properties using a consistent set of parameters with clear physical meanings (Liu et al. 2009; Zhao and Liu 2012), because the corresponding stress-strain relationship is the foundation for deriving other constitutive relationships. Secondly, comparisons between results calculated from the TPHM and test data show that the TPHM is quite general and can be applied to both fractures and porous media (Liu et al. 2009; Zhao and Liu 2012).

Two-phase flow properties As previously indicated, while multiphase flow occurs in a number of practical applications associated with DOI 10.1007/s10040-012-0915-6

375 Table 1 Values for fitting parameters from experimental data Rock type

Initial aperture

Sample No. (loading cycle)

Fitting parameters (as defined in the text) b0,e (mm)

KF,e (MPa)

b0,e/KF,e (mm/MPa)

b0,t (mm)

KF,t (MPa)

R2

Granite

Slate

0.016 mm 0.006 mm 0.0045 mm 0.38 mm 0.505 mm 0.481 mm 0.683 mm < = 0.1 mm

0.00748 0.00315 0.00033 0.128 0.046 0.042 0.064 0.0750

4,179.3 15,181 16,160 39.635 36.508 39.628 25.106 2,309.814

0.000120 3.29359E-05 3.09406E-05 0.003229 0.001260 0.001060 0.002549 3.24749E-05

0.00852 0.00285 0.00417 0.25200 0.45913 0.43949 0.61854 0.02499

3.81679 3.81679 1.79856 2.31922 2.07031 2.38408 2.00781 7.35294

0.985 0.985 0.947 0.985 0.975 0.998 0.990 0.991

Dolerite

0.15 mm

0.1179

4,678.191

2.52003E-05

0.03211

5.84795

0.987

Limestone

0.2 mm

0.1875

8,386.160

2.23624E-05

0.01247

3.96825

0.999

Slate

0.5 mm

0.4303

1,506.220

0.000286

0.06965

5.52486

0.996

Limestone

0.5 mm

0.4336

2,004.390

0.000216

0.06644

9.70874

0.977

Siltstone

0.6 mm

0.5286

2,374.240

0.000223

0.07137

6.62252

0.996

Fused-silica glass Fused-silica glass Fused-silica glass Cheshire quartzite Cheshire quartzite Marble Carnmenellis granite –

Max=113 um

Sample 1 Sample 2 Sample 3 Dr-1 Grd-1 Grd-2 Grd-3 Fresh-slate cleavage (3) Fresh-dolerite joint (3) Fresh-limestone bedding (3) Weathered-slate cleavage (3) Weathered-limestone joint (3) Weathered-siltstone bedding (3) CGL965 (cycle 4,5)

0.1024

657.89

0.000156

0.01060

0.69686

0.928

Max=37.8 um

CGL0394

0.0304

377.39

8.05533E-05

0.00735

0.87184

0.969

Max=28 um

CGL0234

0.0204

90.054

0.000227

0.00757

0.72727

0.988

Max=35.9 um

CCQ0073

0.028

557.34

5.02386E-05

0.00793

6.36943

0.960

Max=74 um

CCQ0102 (cycle 3)

0.0657

2,204.2

2.98067E-05

0.00828

6.32911

0.927

Max=37.9 um –

CSM0042 (cycle 3) Sample 1

0.0299 –

156.04 –

0.000192 0.001453

0.00798 0.06019

0.94518 2.18341

0.951 0.955

–

–

–

– 0.12

– 148.16

0.000591 0.001979 0.000569 0.000810

0.10327 0.34471 0.10939 0.14202

1.5674 3.63636 0.54975 0.90416

0.997 0.992 0.988 0.903

Diorite Granodiorite rock

Ref.

Indraratna et al. (1999)

Malama and Kulatilake (2003)

Bandis et al. (1983)

Brown and Scholz (1985)

Brown and Scholz (1986)

Elliott and Brown (1988)

Granite Kikuma granodiorite Kikuma granodiorite Kikuma granodiorite Inada granite Chichibu schist Kimachi sandstone Carnmenellis granite

– 262 um

Mated Unmated – NKGD (cycle 2)

627 um

TKGD (cycle 2)

0.4511

272.84

0.001653

0.17588

3.11526

0.959

201 um

SKGD (cycle 2)

0.0842

230.94

0.000365

0.11678

1.35685

0.924

650 um 283 um

TIGN (cycle 2) TCSH (Cycle 2)

0.4464 0.1626

743.5 437.63

0.000600 0.000372

0.20361 0.11945

2.43309 0.96993

0.905 0.984

396 um

TKSS (cycle 2)

0.0837

36.668

0.002283

0.31226

0.80645

0.882

–

–

–

–

0.002053 0.003940 0.003762 0.007062 0.001617 0.011696 0.009166

0.04239 0.04854 0.08679 0.09304 0.0536 0.01032 0.0644

0.38183 0.87336 0.48309 0.77942 1.14943 1.01937 0.8244

0.952 0.967 0.943 0.984 0.993 0.981 0.997

Zhao and Brown (1992)

Artificial tension fractures

Charcoal black granite

–

–

–

0.526 mm 0.589 mm 0.39 mm 0.697 mm – –

0.438 0.5128 0.3386 0.61 – –

219 410.2 176.07 586.54 – –

0.000246 0.000623 0.001226 0.002049 0.002 0.001250 0.001923 0.001040 0.000599 4.76281E-05 4.70057E-05

0.02448 0.06751 0.09352 0.07512 0.088 0.07625 0.05141 0.087 0.085 0.00493 0.00151

4.29185 3.10559 1.14286 1.51976 0.83542 1.80505 1.83486 1.443 1.28866 12.1951 3.77358

0.992 0.923 0.948 0.990 0.982 0.965 0.998 0.949 0.994 0.994 0.993

Raven and Gale (1985)

Inada granite

NJ1 NJ2 (180 °C) NJ2 (200 °C) NJ6 EF1 EF3 EF7 (2nd loading120 °C) Sample 1 (cycle Sample 2 (cycle Sample 3 (cycle Sample 5 (cycle IG10 IG15 IG20 IG25 Sample 1 E30 (cycle 2) E32

Granite block Quartz monzonite (Stripa granite)

Hydrogeology Journal (2013) 21: 371–382

3) 3) 3) 3)

Goodman (1976) Iwai (1976) Iwano (1995)

Matsuki et al. (2001)

Sharifzadeh et al. (2008) Pyrak-Nolte et al. (1987)

DOI 10.1007/s10040-012-0915-6

376 Table 1 (continued) Rock type

Initial aperture

Granite samples Granodiorite

Westerley (Rhode Island) granite Austin chalk.

Sample No. (loading cycle)

b0,e (mm)

KF,e (MPa)

b0,e/KF,e (mm/MPa)

b0,t (mm)

KF,t (MPa)

R2

E35 – GO3 1,001 um Load cycle 2 600 um Uniform flow field 760 um Radial flow field 0.11–0.21 mm Sample 1, mated 0.59 mm Sample 1, offset

– 0.8711 0.4968 0.5961 –

– 403.07 170 168.55 –

0.000287 0.000799 0.002161 0.002922 0.003537 0.000125 0.0005

0.00341 0.0157 0.12993 0.10319 0.16386 0.1 0.27

2.89017 10.6383 2.88184 0.5848 1.54799 17.5439 16.9492

0.994 0.995 0.991 0.974 0.941 0.964 0.968

–

–

–

0.006351 0.006774 0.009991

0.04623 0.12593 0.11868

1.41243 1.54321 1.29534

0.996 0.998 0.996

2nd cycle Offset 1 mm Offset 3 mm

fractured rock, studies on the stress-dependence of multiphase flow properties (such as relative permeabilities for different phases) are still lacking in the

(a) 30

Test data (IG10) Calculated curve (IG10) Test data (IG15) Calculated curve (IG15) Test data (IG20) Calculated curve (IG20) Test data (IG25) Calculated curve (IG25)

Normal stress

(MPa)

25 20

15 10 5 0 0

0.02

0.04

0.06

0.08

0.1

0.12

Closure (mm)

(b)

Normal stress

(MPa)

60

Test data (Grd-1) Calculated curve (Grd-1) Test data (Grd-2) Calculated curve (Grd-2) Test data (Grd-3) Calculated curve (Grd-3) Test data (Dr-1) Calculated curve (Dr-1)

50 40 30

20 10 0 0

Ref.

Fitting parameters (as defined in the text)

0.2

0.4

0.6

0.8

Closure (mm)

Fig. 2 Comparisons between measured data (from a Matsuki et al. (2001) and b Malama and Kulatilake 2003) and results calculated from Eq. (8). Items in brackets refer to sample numbers from Table 1 Hydrogeology Journal (2013) 21: 371–382

Chen et al. (2000) Schrauf and Evans (1986) Durham and Bonner (1994) Olsson and Brown (1993)

literature. In this section, relationships are developed between normal stress and multiphase flow properties. Two-phase flow properties include capillary pressure and relative permeability as functions of liquid saturations. The focus here is on two-phase flow within a horizontal fracture when flow paths are mainly determined by capillarity or fracture-aperture distribution. The standard approach used in porous media is based on the local equilibrium assumption (LEA) that capillary pressure is uniformly distributed within the pore space corresponding to representative elementary volume (REV; e.g., Mualem 1976). The LEA essentially implies that liquid distribution for a given phase is completely controlled by capillary force and is independent of liquid flux. That is why capillary pressure and relative permeability can be expressed as functions of saturation only. However, in many cases, the LEA does not hold, especially when instability (or fingering) occurs. For example, for the upward flow of CO2 in a saline aquifer at a geological CO2 sequestration site, significant CO2 fingering will occur because CO2 is lighter than ambient brine. When a numerical grid block contains a number of CO2 fingers (which is generally the case, because the grid-block size cannot be too small for practical applications), LEA will be violated at the grid-block scale. Liu (2011) recently proposed a new theory for multiphase flow to deal with this challenging issue, based on the optimality principle that unstable liquid-flow patterns are self-organized in such a way that total flow resistance is minimal. In this case, the relative permeability is a function not only of saturation, but also liquid flux. Nevertheless, it is appropriate to apply the LEA to individual horizontal fractures that are not subject to gravitational instability for flow process.

Capillary pressure The capillary pressure Pc, defined as the difference between non-wetting-phase and wetting-phase pressures, DOI 10.1007/s10040-012-0915-6

377

is given by the Young-Laplace equation for a fracture element (with aperture B) at the fluid interface: Pc 0

Ts cos a 2B

ð10Þ

integration to be one. To meet this requirement, the righthand side of Eq. (11) can be multiplied by a normalized factor. However, as demonstrated by derivation procedures to be discussed later, the direct use of Eq. (11) is valid for this study, because normalization processes are always involved when this equation is used. Based on the LEA and Eq. (11), the non-wetting phase saturation (Snw) for a capillary pressure given by Eq. (10) is obtained as:

where Ts is surface tension and α is contact angle. Based on the LEA and above equation, fracture apertures larger than Tscosα/(2Pc), for a given capillary pressure Pc, are generally occupied by non-wetting phase and the rest by wetting phase. This reasoning is the foundation for deriving the capillary pressure-saturation relationship (Mualem 1976). The relationship between capillary pressure and fluid saturation is determined by the probability density function (pdf) of fracture aperture and its dependence on stress. Approximately Gaussian or truncated-Gaussian aperture distributions are commonly observed in both natural and man-made fractures (Walsh et al. 2008). As shown in Fig. 3, with increasing stress, small apertures are reduced to zero and the pdf is close to a truncatedGaussian distribution. A number of researchers reported that the pdf under normal load might be described by lognormal distributions to characterize its skewed behavior (e.g., Pruess and Tsang 1990). Recently, Sharifzadeh et al. (2008) developed a new method for measuring fractureaperture distributions and reported measurement results for several artificial fractures. They found that the observed pdfs are similar to Poisson or log-normal distributions. Thus, different kinds of distribution may be used for describing the pdf for fractures. In this study, for simplicity, it was assumed that the fracture-aperture pdf is described by a truncated-Gaussian distribution:

Bu u tb 0 pffiffiffi ; t00 0 pffiffiffi 2d 2d

ðBuÞ2 1 f ðBÞ0 pffiffiffiffiffi e d2 for B 0 2p d

Pc B 0 P B

ð11Þ

Probability Density

where u and δ are mean and standard deviation, respectively, for the corresponding non-truncated Gaussian distributions. Note that the integration of f(B) within the valid range of B≥0 is less than one, as a result of truncation. A mathematically rigorous pdf requires the

u

u0 Aperture b

Fig. 3 Gaussian (dashed curve) and truncated-Gaussian (solid curve) aperture distributions. The solid curve is subject to a larger normal stress. The parameters u0 and u are means for the Gaussian distribution for the dashed curve and the corresponding Gaussian distribution for the solid curve, respectively. They correspond to the largest probability values Hydrogeology Journal (2013) 21: 371–382

R1 Snw 0

B R1

xf ðxÞdx xf ðxÞdx

0

2 pdffiffiffiffi etb 2p 2 pdffiffiffiffi et00 2p

þ u2 erfcðt b Þ þ u2 erfcðt 00 Þ

ð12aÞ

0

ð12bÞ

where f(x) is the pdf defined by Eq. (11). The corresponding wetting-phase saturation (Sw), by definition for a two-phase flow system, is: Sw 01  Snw

ð13Þ

In Eq. (12), fracture aperture B can be replaced by capillary pressure Pc using Eq. (10). If a capillary pressure Pc 0P* is known for a fracture aperture B0B* that corresponds to fracture elements at the interface between wetting and non-wetting phases, then Eq. (10) can be rewritten as: ð14Þ

Equations (10) (or 14), (12) and (13) together comprise the relationship between Pc, Snw, and Sw. For a given capillary pressure, Eq. (10) or Eq. (14) is used to estimate fracture aperture B associated with the interface between wetting and non-wetting phases. Then, saturations are calculated from Eqs. (12) and (13). To determine the stress dependence of capillary pressure–saturation relationships, it must be known how the fracture-aperture pdf changes with normal stress. The latter corresponds to complex mechanical deformations because of the variability in fracture asperity and contact areas. Recently, Walsh et al. (2008) conducted a detailed numerical simulation of mechanical deformation within a fracture under different normal stresses. They found that the aperture pdf becomes more skewed with increasing normal stress, but is still closer to a truncated-Gaussian distribution than a log-normal distribution. They also noticed that the pdf evolution can be approximately described by a so-called penetration model, in which fracture aperture changes uniformly while aperture values are limited to being zero or positive. The same model was used by other researchers in determining stress-dependent DOI 10.1007/s10040-012-0915-6

378

fracture properties (e.g., Brown 1987; Oron and Berkowitz 1998). In this study, the penetration model is adopted. In other words, parameter δ is treated as being independent of stress, but u will change with stress, as illustrated in Fig. 3, where the solid curve represents a pdf corresponding to a larger stress. The parameter u can be related to average aperture b(Eq. 7) by: R1

Bf ðBÞdB

b0 0R1

f ðBÞdB

0u þ

qffiffiffi 2 2 t00 p de

ð15Þ

erfcðt00 Þ

0

where t00 is defined in Eq. (12b). Note that in Eq. (15), Bf(B)dB represents the volume occupied by fracture segments with apertures from B to B+dB. The stress dependence of b is given in Eq. (7). For a given value of b, u can be derived from Eq. (15), but an iterative numerical procedure is needed. Then the capillary pressure-saturation relationship can be calculated through the dependence of u on stress. Figure 4 shows curves of capillary pressure (as a function of wetting-phase saturation) for two different stresses (or u values). As expected, the capillary pressure curve strongly depends on the normal stress.

Relative permeability The relationship between fracture permeability and stress is given in Eqs. (9) and (7). Modeling of multiphase-flow processes also requires the relative permeability for different phases. As previously indicated, studies on fracture relative permeability are very limited and no systematic study has been found in the literature that investigates stress-dependent fracture relative permeability. 10

2

Pc / P *

101

In this section, relationships are developed between fracture relative permeability, saturation, and stress for a two-phase flow system. The general approach for developing closed-form relativity relationships for porous media involves two key aspects. Firstly, based on the pdf of pore-size distribution, a conceptual model of multiphase flow through the medium is developed. There are currently two commonly used conceptual models. One was proposed by Mualem (1976), in which the spatial distribution of soil pores is completely random. The other one is the model developed by Burdine (1953), in which a perfect spatial correlation of soil pores is assumed, such that flow paths can be conceptualized as a group of parallel capillary tubes. In both models, pore size is determined from capillary pressure–saturation relationship using the Young-Laplace equation that suggests that pore-size is inversely proportional to the capillary pressure. Note that pore size for porous media corresponds to fracture aperture herein. The second aspect of this general approach is to make further corrections to the results, determined from these conceptual models, by multiplying them with a tortuosity factor. The rationale behind this is that in those conceptual models, flow paths are assumed to be straight, while in reality, they are tortuous as a result of the geometric complexity of the pore spaces. The tortuosity factor is generally assumed to be a power function of saturation for a given phase and the exponent value needs to be empirically determined from measurements. In the current study, this general approach is followed to develop relative permeability for a horizontal fracture. Note that the corresponding relationships cannot be directly borrowed from porous media, simply because geometry (and connectivity) of fracture apertures is essentially two-dimensional and is different from that of pores in porous media, which is three-dimensional. Burdine’s (1953) model was employed herein, based on the consideration that flow paths in fractures are close to parallel capillary tubes that have been implied by commonly reported long-range correlations (fractal behavior) of aperture distribution and observed channelized multiphase flow patterns (e.g., Walsh et al. 2008; Brown 1987; Power and Tullis 1992; Chen and Horne 2006). The wetting phase relative permeability, kwr, is expressed as (Burdine, 1953): RB

100

kwr 0

Bmin R1 Bmin

10

-1

0

0.25

0.5

0.75

1

Sw Fig. 4 Normalized capillary pressure (Pc/P*) as a function of wetting-phase saturation (Sw) for two stresses, with the solid curve corresponding to a smaller normal stress. The parameter values are: u=0.145 mm (for the solid curve) and 0.045 mm (for the dashed curve), and δ=0.03 mm (parameters as defined in the text) Hydrogeology Journal (2013) 21: 371–382

1 Pc2

dV

1 Pc2

dV

Swm

ð16Þ

where the power function term for Sw is the tortuosity factor with m as an empirical parameter, and dV is the volume element corresponding to fracture aperture B associated with capillary pressure Pc: dV / Bf ðBÞdB

ð17Þ DOI 10.1007/s10040-012-0915-6

379

In Eq. (16), Bmin refers to a critical fracture aperture. When the wetting phase is confined to fracture apertures smaller than Bmin, the wetting phase becomes discontinuous and thus immobile. Inserting Eqs. (10) and (17) into (16) yields: RB

Zy

u3 t3 f ðtÞdt0 pffiffiffi I1 ðx; yÞ p x rffiffiffiffi 2 6d 2 u 2 I2 ðx; yÞ þ pffiffiffi I3 ðx; yÞ þ3du p p rffiffiffiffi 2 I4 ðx; yÞ þ 2d 3 p

It ðx; yÞ0

B3 f ðBÞdB

B

min kwr 0Swm R1

B3 f ðBÞdB

ð18Þ

Bmin

Similarly, the relative permeability for the non-wetting phase (knwr) can be written as: BRmax m B knwr 0Snw BRmax

where u and δ are parameters associated with the fractureaperture pdf defined in Eq. (11). Inserting Eq. (11) into Eqs. (18) and (19) and using Eqs. (20a)–(20e), the relative-permeability expressions are obtained for wetting and non-wetting phases as:

kwr 0Swm

B3 f ðBÞdB ð19Þ

It ðtmin ; tb Þ It ðtmin ; 1Þ

where Bmax is the other critical fracture aperture. When the non-wetting phase is confined to fracture apertures larger than Bmax, the non-wetting phase becomes immobile as a result of blockage by the wetting phase. Note that Bmin and Bmax correspond to residual saturations for wetting and non-wetting phases, respectively. As previously discussed, the fracture-aperture pdf can be presented by a truncated-Gaussian distribution given in Eq. (11). Before presenting the relative permeability relationships subject to the pdf, it is useful to define several derived integration functions: I1 ðx; yÞ0

e

t 2

pffiffiffi p dt0 ½erfcðxÞ  erfcðyÞ 2

ð20aÞ

x

Zy I2 ðx; yÞ0

tet dt0 2

i 1 h x2 2 e  ey 2

ð20bÞ

x

Zy I 3 ðx; yÞ0

ð21Þ

B3 f ðBÞdB

0

Zy

ð20eÞ

 1 1  x2 2 2 xe  yey t 2 et dt0 I 1 ðx; yÞ þ 2 2

m knwr 0Snw

tmin 0

It ðtb ; tmax Þ It ð0; tmax Þ

ð22Þ

Bmin  u Bmax  u pffiffiffi ; tmax 0 pffiffiffi 2d 2d

ð23Þ

where tb is defined in Eq. (12a). Relative permeability (Eqs. 21 and 22) is a function of the parameter u that is related to stress through fracture aperture (Eq. 15). The relationship between average fracture aperture and stress is described by the TPHM (Eq. 7). Relative permeability is also a function of Bmin or Bmax. For the non-wetting phase, the residual saturation can be determined using the following percolation procedure: consider a horizontal fracture to be initially filled with a wetting phase; then the largest apertures occupied by the wetting phase are replaced by a nonwetting phase step by step until a continuous non-wetting flow path is formed. The corresponding non-wetting phase saturation is the residual saturation. Since fracture aperture changes uniformly (Fig. 3), regions within a fracture that are occupied by a residual non-wetting phase remain the same under different stress conditions. In other words, Bmax follows:

ð20cÞ

x

Bmax  Bmax;ref 0u  uref

Zy

where the subscript ref refers to a reference condition under which related parameter values can be estimated from testing data. Applying a similar consideration to the wetting phase, that smallest apertures initially occupied by the non-wetting phase are replaced by a wetting phase step

I 4 ðx; yÞ0

t 3 et dt0 2

 o 1 n x2  2 2 e x þ 1  ey y2 þ 1 2

x

Hydrogeology Journal (2013) 21: 371–382

ð20dÞ

ð24Þ

DOI 10.1007/s10040-012-0915-6

380

Bmin  Bmin;ref 0u  uref

ð25Þ

where the smallest value for Bmin is limited to being zero. Equation (25) is based on an argument that fracture regions containing a residual wetting phase (including fracture space that is closed as result of increasing stress) remain unchanged with changing stress. Unlike the nonwetting phase, the situation for the wetting phase is more complex than the non-wetting phase, because the residual regions involve opening and closing fracture apertures with changing stress, which can considerably alter wetting-phase connectivity in some cases. Thus, Eq. (25) is considered to be a first-order approximation at this point. It is also useful to note that Bmin 00 and Bmax 0∞ correspond to zero residual saturations for wetting and non-wetting phases, respectively. As an example, Fig. 5 shows relative permeability curves for two stresses (or different u values), indicating considerable sensitivity to the stress change.

even exists, in the literature. This subsection therefore evaluates these relationships for fractures that are not subject to mechanical deformation. Two aspects in this evaluation are of particular interest: (1) whether the relationships can satisfactorily represent experimental observations; and (2) whether a tortuosity factor with the same m value (Eqs. 21 and 22) can be applied to different phases for different tests. The universal m value has been used in porous media (Brooks and Corey 1964). The dataset of Chen and Horne (2006) were used for evaluating the theoretical results, because this dataset, including information on the fracture–aperture pdf

(a)

1

0.9 0.8 0.7

kwr ,knwr

by step until a continuous wetting flow path is formed, gives:

kwr

knwr

0.6 0.5 0.4

Evaluation of the new relationships In the preceding two subsections, closed-form relationships between capillary pressure, saturation, and stress were developed for two-phase flow in a horizontal fracture. The dependence of these relationships on stress is realized though parameter u which changes with stress. Unfortunately, testing data on the stress dependence of multiphase flow properties for fractures is very rare, if it 1

0.3 0.2 0.1 0

(b)

0

0.25

0.5 Sw

0.75

1

1

0.9 0.9

0.8 0.8

kwr

0.6

0.7

knwr kwr ,knwr

kwr , knwr

0.7

0.5 0.4

kwr

knwr

0.6 0.5 0.4

0.3 0.3

0.2 0.2

0.1 0.1

0

0

0.25

0.5

0.75

1

Sw Fig. 5 Relationships between relative permeability and saturation for two different normal stresses with the dashed curve corresponding to a larger stress. The parameter values are: u= 0.145 mm (for the solid curve) and 0.045 mm (for the dashed curve), Bmin =0.13 mm (for the solid curve) and 0.03 mm (for the dashed curve), Bmax =0.165 mm (for the solid curve) and 0.065 mm (for the dashed curve), and δ=0.03 mm (parameters as defined in the text) Hydrogeology Journal (2013) 21: 371–382

0

0

0.25

0.5 Sw

0.75

1

Fig. 6 Comparisons between the relative permeability data of Chen and Horne (2006) (black squares and blue circles) and calculated results (solid curves) from Eqs. (21) and (22) for two fractures: a homogeneously rough; and b randomly rough. Parameter values are: a u=0.145 mm, Bmin =0.13 mm, Bmax = 0.165 mm, and δ=0.03 mm; b u=0.24 mm, Bmin =0.22 mm, Bmax = 0.28 mm, and δ=0.05 mm (parameters as defined in the text) DOI 10.1007/s10040-012-0915-6

381

(Gaussian distribution), and the corresponding experimental procedures, are very well documented. They conducted concurrent, steady state, air–water-flow experiments in horizontal fractures at room temperature. The fractures for the experiments, with dimensions of around 35 cm × 15 cm, were created by mounting a rough glass plate on top of a smooth aluminum plate, confined by a metal frame bolted to the bottom plate. All measurements were digitized automatically with a high-speed data acquisition system. Two fractures with distinct surface roughness, homogeneously rough (HR) and randomly rough (RR), were used in the tests; the former does not have correlation in spatial-aperture distribution and the latter does. Their test results show strong channelized flow behavior and they also reported relative permeability data as a function of water saturation for both the wetting phase (water) and the non-wetting phase (air). Details of the dataset and its collection procedure can be found in Chen and Horne (2006). Data for a smooth fracture were also reported in Chen and Horne (2006), but were not used in this study because all natural fractures are rough. Figure 6 shows comparisons between the relative permeability data of Chen and Horne (2006) and results calculated from Eqs. (21) and (22) for two (RR and HR) fractures. Parameter values for u and δ were taken directly from Chen and Horne (2006). Values for parameters Bmin, Bmax, and m in Eqs. (21–23) were adjusted to achieve curve fitting. In general, matches between calculated results and data are reasonable, suggesting that the theoretical relationships can satisfactorily represent the test data. It is especially of interest to note that a single m value of 1.0 was used for both wetting and non-wetting phases for the two fractures. Note that m02 in the widely used Brooks and Corey (1964) model that was developed for soils from Burdine’s (1953) model, which means that for a given saturation, soils have smaller tortuosity-factor values than fractures. This makes sense physically, because flow paths are more tortuous in three-dimensional pore space for soils than in essentially two-dimensional fracture space. However, it is acknowledged that further evaluation of the relationships with more datasets, when they become available, is necessary to further confirm the result of m01 for fractures.

Summary and conclusions Mechanical deformation processes in fractured rock and their coupling to hydrological processes are important for many applications. Constitutive relationships are fundamental for modeling coupled hydro-mechanical processes, because they determine how coupled processes are actually coupled in a numerical simulation. An important constitutive relationship for fractured rock is the stress-dependence of fracture hydraulic properties. A systematic approach has been presented for determining relationships between normal stress and fracture properties. The relationship between stress and fracture permeability (or related aperture and closure) is based on a Hydrogeology Journal (2013) 21: 371–382

two-part Hooke’s model (TPHM) that captures heterogeneous elastic-deformation processes at a macroscopic scale by conceptualizing the rock mass (or a fracture) as having two parts with different mechanical properties. The developed relationship was verified using a number of datasets obtained from the literature for fracture closure versus stress and satisfactory agreements were obtained. The TPHM was shown to be able to accurately represent testing data for porous media as well. Multiphase flow occurs in a number of practical applications associated with fractured rock. Based on the consideration that fracture-aperture distributions under different normal stresses can be captured by truncatedGaussian distributions, closed-form constitutive relationships were developed between two-phase flow properties (capillary pressure, relative permeability and saturations) for deformable horizontal fractures. The usefulness of these relationships was partially demonstrated by their consistency with a laboratory dataset that was collected without involving mechanical deformation. Acknowledgements We thank Drs. Jim Houseworth, Daisuke Asahina, and Daniel Hawkes at Lawrence Berkeley National Laboratory for reviewing the initial version of the paper. We also appreciate the constructive comments from the associate editor (Dr. Philipp Blum) and two anonymous reviewers. This work was supported by the Assistant Secretary for Fossil Energy, Office of Sequestration, Hydrogen, and Clean Coal Fuels of the US Department of Energy under Contract No. DE-AC02-05CH11231. In particular, we would like to acknowledge In Salah JIP and their partners BP, Statoil, and Sonatrach for providing valuable discussions on the subject.

References Bandis SC, Lumsden AC, Barton NR (1983) Fundamentals of rock joint deformation. Int J Rock Mech Min Sci Geomech Abstr 20 (6):249–268 Barton N, Bandish S, Bakhtar K (1985) Strength, deformation and conductivity coupling of rock joints. Int J Rock Mech Min Sci Geomech Abstr 22(3):121–140 Brooks RH, Corey AT (1964) Hydraulic properties of porous media. Hydrology paper no. 3, Civil Engineering Dept., Colorado State Univ., Fort Collins, CO Brown SR (1987) Fluid flow through rock joints: the effect of surface roughness. J Geophys Res Solid Earth 92: 1337–1347 Brown SR, Scholz CH (1985) Closure of random surfaces in contact. J Geophys Res 90(B7):5531–5545 Brown SR, Scholz CH (1986) Closure of rock joints. J Geophys Res 91(B5):4939–4948 Burdine NT (1953) Relative permeability calculations from poresize distribution data. Petr Trans Am Inst Mining Metall Eng 198:71–77 Chen CY, Horne RN (2006) Two-phase flow in rough-walled fractures: experiments and a flow structure model. Water Resour Res 42:W03430. doi:10.1029/2004WR003837 Chen Z, Narayan SP, Yang Z, Rahman SS (2000) An experimental investigation of hydraulic behaviour of fractures and joints in granitic rock. Int J Rock Mech Min Sci 37:1061–1071 Diomampo GP (2001) Relative permeability through fractures. Rep. SGP-TR-170. Stanford University, Stanford, CA DOI 10.1007/s10040-012-0915-6

382 Durham WB, Bonner BP (1994) Self-propping and fluid flow in slightly offset joints at high effective pressures. J Geophys Res 99(B5):9391–9399 Elliot GM, Brown ET (1988) Laboratory measurement of the thermo-hydro-mechanical properties of rock. Q J Eng Geol Hydrogeol 21:299–314 Fourar M, Lenormand R (1998) A viscous coupling model for relative permeabilities in fractures. SPE 49006, paper presented at the SPE Annual Technical Conference and Exhibition, New Orleans, LA Gangi AF (1975) Variation of whole and fractured porous rock permeability with confining pressure. Int J Rock Mech Min Sci Geomech Abstr 15:249–257 Goodman RE (1974) The mechanical properties of joints. Proc. 3rd Congr ISRM, vol 1, Denver, CO, September 1974, pp 127–140 Goodman RE (1976) Method of geological engineering in discontinuous rocks. West, New York Indraratna B, Ranjith PG, Gale W (1999) Single phase water flow through rock fractures. Geotech Geol Eng 17:211–240 Indraratna B, Price J, Ranjith P, Gale W (2002) Some aspects of unsaturated flow in jointed rock. Int J Rock Mech Min Sci 39:555–568 Iwai K (1976) Fundamental studies of fluid flow through a single fracture. PhD Thesis, University of California, Berkeley, USA Iwano M (1995) Hydromechanical characteristics of a single rock joint. PhD Thesis, Massachusetts Institute of Technology, USA. Jing L, Hudson JA (2002) Numerical methods in rock mechanics. Int J Rock Mech Mining Sci 39:409–427. doi:10.1016/ j.ijrmms.2008.04.005 Liu HH (2011) A conductivity relationship for steady-state unsaturated flow processes under optimal flow conditions. Vadose Zone J 10(2):736–740. doi:10.2136/vzj2010.0118 Liu HH, Rutqvist J, Berryman JG (2009) On the relationship between stress and elastic strain for porous and fractured rock. Int J Rock Mech Min Sci 46(2):289–296 Liu HH, Rutqvist J, Birkholzer JT (2011) Constitutive relationships for elastic deformation of clay rock: data analysis. Rock Mech Rock Eng 44(4):463–468. doi:10.1007/s00603-010-0131-4 Malama B, Kulatilake PHSW (2003) Models for normal fracture deformation under compressive loading. Int J Rock Mech Min Sci 40:893–901 Matsuki K, Wang EQ, Sakaguchi K, Okumura K (2001) Timedependent closure of a fracture with rough surfaces under constant normal stress. Int J Rock Mech Min Sci 38: 607–619 Mavko G, Jizba D (1991) Estimating grain-scale fluid effects on velocity dispersion in rocks. Geophysics 56:1940–1949 McDermott CI, Kolditz O (2006) Geomechanical model for fracture deformation under hydraulic, mechanical and thermal loads. Hydrogeol J 14(4):487–498 Mualem Y (1976) A new model of predicting the hydraulic conductivity of unsaturated porous media. Water Resour Res 12:513–522 Olsson R, Barton N (2001) An improved model for hydromechanical coupling during shearing of rock joints. Int J Rock Mech Min Sci 39:555–568

Hydrogeology Journal (2013) 21: 371–382

Olsson WA, Brown SR (1993) Hydromechanical response of a fracture undergoing compression and shear. Int J Rock Mech Min Sci Geomech Abstr 30(7):845–851 Oron AP, Berkowitz B (1998) Flow in rock fractures: the local cubic law assumption reexamined. Water Resour Res 34:2811–2825 Power WL, Tullis TE (1992) Euclidean and fractal models for the description of rock surface roughness. J Geophys Res Solid Earth 96:415–424 Pruess K, Tsang YW (1990) On two-phase relative permeability and capillary pressure of rough-walled rock fractures. Water Resour Res 26(9):1915–1926 Pyrak-Nolte LJ, Myer LR, Cook NGW (1987) Witherspoon PA. Hydraulic and mechanical properties of natural fractures in lowpermeability rock. 6th ISRM Congress, Montreal, QB, September 1991 Raven KG, Gale JE (1985) Water flow in a natural rock fracture as a function of stress and sample size. Int J Rock Mech Min Sci Geomech Abstr 22(4):251–261 Rutqvist J, Stephansson O (2003) The role of hydromechanical coupling in fractured rock engineering. Hydrogeol J 11:7–40 Schrauf TW, Evans DD (1986) Laboratory studies of gas flow through a single natural fracture. Water Resour Res 22(7):1038–1050 Sharifzadeh M, Mitani Y, Esaki T (2008) Rock joint surface measurement and analysis of aperture distribution under different normal and shear loading using GIS. Rock Mech Rock Eng 41(2):299–323 Therrien R, Sudicky EA (1996) Three-dimensional analysis of variably-saturated flow and solute transport in discretelyfractured porous media. J Contam Hydrol 23:1–44 Tsang CF, Barnichon JD, Birkholzer J, Li XL, Liu HH, Sillen X (2012) Coupled thermo-hydro-mechanical processes in the near field of a high-level radioactive waste repository in clay formations. Int J Rock Mech Min Sci 49:31–44. doi:10.1016/ j.ijrmms.2011.09.015 Van Golf-Racht TD (1982) Fundamentals of fractured reservoir engineering. Elsevier, New York Walsh R, McDermott C, Kolditz O (2008) Numerical modeling of stress-permeability coupling in rough fractures. Hydrogeol J 16:613–627 Weerakone WMSB, Wong RCK, Mehrotra (2011) Single-phase (Brine) and two-phase (DNAPL-Brine) flow in induced fractures. Transp Porous Med 89:75–95 Witherspoon PA, Wang JSW, Iwai K, Gale JE (1980) Validity of cubic law for fluid flow in deformable rock fracture. Water Resour Res 16(6):1016–1024 Wong RCK, Pan X, Maini BB (2008) Correlation between pressure gradient and phase saturation for oil-water flow in smooth- and rough-walled parallel plate models. Water Resour Res 44: W02418. doi:10.1029/2007/WR006043 Zhao J, Brown ET (1992) Hydro-thermo-mechanical properties of joints in the Carnmenellis granite. Q J Eng Geol 25:279–290 Zhao Y, Liu HH (2012) An elastic strain-stress relationship for porous rock under anisotropic stress conditions. Rock Mechanics and Rock Engineering. doi:10.1007/s00603-011-0193-y Zimmerman RW, Bodvarsson GS (1996) Hydraulic conductivity of rock fractures. Transp Porous Media 23:1–30

DOI 10.1007/s10040-012-0915-6