Analytical solution of coupled stress-flow-transport processes in a single rock fracture

Computers & Geosciences 37 (2011) 1437–1449

Contents lists available at ScienceDirect

Computers & Geosciences journal homepage: www.elsevier.com/locate/cageo

Analytical solution of coupled stress-flow-transport processes in a single rock fracture Zhihong Zhao a,n, Lanru Jing a, Ivars Neretnieks b, Luis Moreno b a b

Department of Land and Water Resources Engineering, Royal Institute of Technology, Sweden Department of Chemical Engineering and Technology, Royal Institute of Technology, Sweden

a r t i c l e i n f o

abstract

Article history: Received 9 July 2010 Received in revised form 20 December 2010 Accepted 22 February 2011 Available online 11 March 2011

A closed-form solution is presented for modeling the coupled stress-flow-transport processes along a single fracture embedded in a porous rock matrix. Necessary assumptions were made to simplify the subject into a two-dimensional (2D) problem, considering the changes of fracture aperture and matrix porosity under various stress conditions. The cubic law was assumed to be valid for the fluid flow in the fracture, with an impermeable rock matrix. For transport mechanisms, advective transport along the fracture, longitudinal hydrodynamic dispersion in the flow direction, and the matrix diffusion were considered in three different transport models under constant concentration or constant flux (Danckwerts’) inlet boundary conditions. This analytical solution can be used as a constitutive model, or as an example for validation of similar constitutive models, for modeling the coupled hydro-mechanicalchemical (HMC) processes in fracture networks of crystalline rocks. The influences of stress/deformation processes on different transport mechanisms in a single fracture under different inlet boundary conditions were studied for the first time. The results show that changes of fracture, as controlled by a combination of normal closure and shear dilatancy, have a significant influence on the solute concentration distribution both along the fracture and in the rock matrix, as well as on the solute residence/breakthrough time, especially when shear-induced dilatancy occurs. Under compressions, the decreasing matrix porosity slightly increases the solute concentration along the fracture and in the rock matrix. & 2011 Elsevier Ltd. All rights reserved.

Keywords: Single rock fractures Analytical solution Stress Fluid flow Transport

1. Introduction The fracture networks in crystalline rocks, for example, granites, provide the major pathways for groundwater flow and contaminant migration. The fractures’ dominating influence on flow and transport processes is an important issue in many scientific or technical fields, such as geological radioactive waste repositories, naturally fractured petroleum, and geothermal reservoirs in crystalline rocks (Jing, 2003; Jing and Stephansson, 2007). Due to the in situ stresses, the pattern and flow rate of groundwater in crystalline rocks are influenced by fracture deformation, which, in turn, have an impact on mass transport processes. In the past, the mechanical, hydraulic, and mass transport processes in rock fractures were most often studied separately in different branches of geosciences, and the combined effects of stress/deformation on fluid flow and mass transport have not been investigated using analytical solutions, despite the fact that the individual processes of stress, flow, or transport have been extensively investigated by using analytical or numerical modeling. In this paper, we present a

n

Corresponding author. E-mail address: [email protected] (Z. Zhao).

0098-3004/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.cageo.2011.02.015

closed-form solution for the coupled stress-flow-transport processes in a single rock fracture, for the first time, to the authors’ knowledge. Mass transport through single rock fractures is a complex natural phenomenon including different mechanisms, such as advection and hydrodynamic dispersion, matrix diffusion, and sorption reactions (Bodin et al., 2003a, 2003b). Since the 1970s, a large number of analytical, semianalytical, or numerical models have been proposed to simulate contaminant transport processes in single rock fractures and fracture networks, using random walk particle tracking techniques (Neretnieks, 1980; Tang et al., 1981; Moreno and Rasuson, 1986; Moreno et al. 1988). However, the effects of the coupled stress-flow processes on these complicated transport mechanisms have not been studied, compared with its potential importance in practice for subsurface engineering and environmental protections. Mathematical models of coupled stress/deformation, groundwater flow, and solute transport processes in fractured rocks could be developed by using either equivalent porous continuum or discrete element method (DEM) models. For the equivalent continuum models, the governing equations, derived based on the principles of continuum mechanics, are solved numerically to describe the macroscopic behavior of the fractured rocks, while overlooking the details of the fracture network geometry. It is the

1438

Z. Zhao et al. / Computers & Geosciences 37 (2011) 1437–1449

most commonly applied modeling approach used for large-scale problems of either sparsely or heavily fractured rocks. The DEM models, on the other hand, explicitly represent the geometry and contact mechanisms of blocks (rock matrix) and fracture systems, based on the solutions of the equations of motion of the blocks, and the constitutive models of rock matrix (block materials) and fractures, respectively. A closed-form solution for the coupled stress/deformation, fluid flow, and mass transport processes in a single fracture, therefore, can serve as an effective tool contributing to establish the constitutive models of rock fractures for DEM modeling approaches. The main objective of this paper is to develop an analytical solution for the coupled stress-flow-transport problems in a single fracture bounded by porous rock matrix, considering the change of fracture aperture and rock matrix porosity resulting from stress/deformation, steady fluid flow, and different dispersive and diffusive processes of solute transport. The model represented by this solution can be especially useful for safety/ performance assessments of underground radioactive waste repositories in fractured crystalline rocks, since the simulation of radioactive nuclide transport process depends largely on the validity and reliability of the models of the stress-flow-transport in individual rock fractures.

2. Conceptual model The simplest mathematical model of a single rock fracture is the smooth parallel plate model, which assumes that the two walls of a fracture can be represented by two nominally smooth and parallel surfaces (Snow, 1965), separated by an initial aperture of 2b [L] (Fig. 1). Considering a single fracture embedded in a saturated porous rock matrix, with a contaminant source located at the origin of the fracture (x¼ 0), the solute particles are then transported by the groundwater flow along the fracture (the advection), and may also diffuse in and out of the porous matrix. Although the parallel plate assumption differs from the natural fractures with rough surfaces, it remains, however, the basis of many complicated modeling approaches for flow and transport processes at the local scale (Bodin et al., 2003b). It can provide the basic insights for fluid flow and solute transport through an individual rock fracture under different stress conditions. More importantly, the parallel plate model is the necessary condition for deriving a closed-form solution for such a complex problem of coupled stress-flow-transport processes in a single fracture, because such an analytical solution does not exist for natural fractures of rough surfaces. This, however, does not prevent considerations of some mechanical effects of roughness on the mechanical component of the closed-form solution indirectly, as described in Section 3.

Besides the parallel plate model as noted above, in order to derive an analytical solution, the following basic assumptions were also adopted: (1) The model was defined in a 2D space, with the fracture aperture much smaller than its length. (2) The hydraulic aperture of the fracture was assumed to be equal to its mechanical aperture, and a residual aperture always existed even under high compression stress. This could be regarded as an indirect consideration of the effects of rock fractures surface roughness in reality. (3) The effect of fluid pressure and its variation on fracture deformation in the flow direction was assumed to be negligible, which indicates that the fracture surfaces always remain parallel during the deformation/displacement process. This assumption is reasonable because usually groundwater flows under small hydraulic pressure gradients through a fracture in subsurface rocks in practice. (4) The fracture aperture changes due to normal closure or shear dilatancy were described by commonly adopted mechanical constitutive models of rock fractures in the field of rock mechanics, in order to simplify the derivation of the final analytical solution. This assumption is reasonable, since the scope of the study is the effect of the stress on the flow and transport processes, not the complicated stress/deformation processes in details. (5) Transverse dispersion within the fracture was neglected, compared to longitudinal dispersion. This means that the solute concentration was assumed to be equal at the fixed fracture section. (6) The permeability of the porous matrix was sufficiently low to hinder water flow, but contaminant particles could enter the rock matrix by molecular diffusion. With the above conceptualization and assumptions, we can explicitly describe the individual mechanical, hydraulic, and transport processes in single fractures, separately, and then combine them together to generate a closed-form solution for coupled stress-flow-transport processes in a single rock fracture.

3. Coupled stress-flow-transport processes 3.1. Mechanical process When stress is applied to rocks, both the fractures and the rock matrix will move and/or deform, which results in changes in fracture aperture and matrix porosity. These changes influence groundwater flow and contaminant transport in fractures, as well as the diffusion of solutes into and out of the rock matrices.

z T: tensile strength kn: nonlinear normal stiffness

2b

0

Ss S

ks

T D

kn

D: dilation

x S: slider Ss: shear strength ks: nonlinear shear

Fig. 1. Schematic illustration of parallel plate model for fractures (after Bodin et al., 2003b).

Z. Zhao et al. / Computers & Geosciences 37 (2011) 1437–1449

(1985) suggested that in situ fractures probably behave in a manner similar to the third or fourth loading cycles (Fig. 2a).

3.1.1. Normal closure behavior It has been demonstrated that the normal closure of a rough rock fracture is typically nonlinear under compressive normal stress; i.e., the normal stiffness of a fracture increases with increasing normal stress (Goodman, 1976; Swan, 1983; Cook, 1992; Rutqvist and Stephansson, 2003), and finally levels off to some asymptotic value at high values of the compressive normal stress (Jaeger et al., 2007). The most commonly applied fracture normal closure model is Bandis’ hyperbolic function (Fig. 2b; Bandis et al. (1983)):

sn ¼ kn0

 un , 1un =unmax

3.1.2. Shear behavior Fracture displacement during shear is inherently a coupled process, in which both normal and shear displacements occur (Jaeger et al., 2007). Fig. 3a shows the typical, although much idealized, shear stress–shear displacement behavior of a fresh rock fracture under constant normal stress conditions, during direct shear tests under laboratory test conditions. During this process, the most important mechanism is the so-called shear dilatancy, i.e., the increase of fracture aperture during shear displacement due to roughness of the fracture surfaces. Various empirical models have been suggested to simulate the dilatancy phenomenon in the past, with varying degrees of success. The simplest approach is to use a dilation angle relating the normal and tangential displacements of a fracture:



ð1Þ

where sn [M/LT2] is the normal stress (compression positive), kn0 [M/L2T2] is the initial normal stiffness at a prescribed stress state, un [L] is the normal displacement of the fracture, and unmax [L] is the maximum normal displacement, approached asymptotically with increasing normal stress. The basic parameters kn0 and unmax are determined by experiments (Barton et al., 1985). The fracture’s normal closure can be expressed by rearranging Eq. (1) as un ¼

sn kn0 þ sn =unmax

:

1439

ud ¼ us tan fd ,

ð3Þ

where ud [L] is the shear dilatancy, us [L] is the shear displacement and fd is the dilation angle. Note that we did not consider a dilation angle that varies with normal stress and shear displacement in this study, but assumed that the dilation angle was a constant for deriving a closed-form solution. A simplified model, as shown in Fig. 3b, was adopted to approximate the behavior of rock fractures during shear (Fig. 3a), so that the final analytical model can capture the main

ð2Þ

It is noted here that those parameters in the above equation depend on the stress paths of normal compressive loading– unloading cycles (Jing and Stephansson, 2007). Barton et al.

σn

σn

un

un

Fig. 2. Fracture closure behavior under normal compression (after Bandis et al., 1983). (a) Compressive loading-unloading cycles and (b) normal stress — normal displacement curve.





p

p

r

r

us

us (I) (II) ud

ud udmax

(III)

udmax

φd us

us up

ucs

Fig. 3. Mechanical behavior of fractures during direct shear testing under constant normal stress (after Jing and Stephansson, 2007). (a) Conceptual model and (b) simplified model.

1440

Z. Zhao et al. / Computers & Geosciences 37 (2011) 1437–1449

mechanical features without losing necessary generality. The peak shear stress, tp [M/LT2], is determined by the Mohr– Coulomb criterion, characterized by a cohesion (C) and a frictional angle (f): 9t9r C þ sn tan f ¼ tp :

ð4Þ

Before reaching the peak shear stress, the shear stiffness, ks, is assumed to be constant. In the case of jtj Z tp ,

ð5Þ

Therefore, the stress-dependent porosity of rock matrix can be written as a3 sn

y ¼ 1ð1yi Þea1 ðsn sni Þa2 =a3 ðe

ea3 sni Þ

,

ð12Þ

where yi is the initial matrix porosity at zero stress. This relation links stress with rock porosity that, in turn, is related to matrix diffusion process. 3.2. Fluid flow process

one has

t ¼ signðus Þtp ,

ð6Þ

where sign( ) is the sign function, which extracts the sign of the shear displacement, i.e., the shear direction. Dilatancy starts to occur at the onset of slip of the fracture, governed by the dilation angle. The accumulated dilation is limited by a critical value of shear displacement, ucs [L], corresponding to the observation that crushing of asperities at large shearing would eventually prevent continued increase of dilatancy. In this way, three functions are employed to describe the dilatancy during three different stages (I, II, and III in Fig. 3b) of the shear displacement: 8 0 us r up with up ¼ tp =ks ðIÞ > < ðIIÞ : ud ¼ us tan fd up rus r ucs ð7Þ > : u tan f u Z u ðIIIÞ cs s cs d

3.1.3. Hydraulic aperture Experimental results typically show that an apparent residual aperture, bres [L], exists in rock fractures even under high compressive stresses, when fracture appears to be nominally closed by compression (Witherspoon et al., 1979; Raven and Gale, 1985; Pyrack-Nolte et al., 1987; Cook, 1992; Renshaw, 1995). Additionally, the residual aperture may vary with sizes of the tested fracture specimen, which was not considered in this paper. Therefore, the resultant hydraulic aperture can be expressed as 2b ¼ bres þ ðunmax un Þ þ ud ,

ð8Þ

where un and ud are determined by Eqs. (2) and (7), respectively. 3.1.4. Rock matrix porosity Rock matrix porosity, y, changes under the in situ stresses and can be defined as a function of volume change of micropores in the rock matrix (Han and Dusseault, 2003). In this study, Eq. (9) relates the increment in matrix porosity to the increment in compressive stress and was adopted because of the explicit inclusion of stress (Zimmerman, 1991): dy ¼ ½Cbc ð1yÞCm dsn , 2 2

ð9Þ 2 2

where Cbc [L T /M] and Cm [L T /M] are effective bulk compressibility and rock matrix compressibility, respectively. Note that sn is the compressive stress acting on the rock matrix. In this study we assumed that it was equal to the stress applied on the fracture surfaces. Because Cm is usually small enough to be negligible, Eq. (9) can be simplified into dy ¼ Cbc ð1yÞdsn :

ð10Þ

In reality, Cbc is stress dependent and can be measured by laboratory tests. Based on experimental data, Zimmerman (1991) derived an empirical relation for stress-dependent Cbc as Cbc ¼ a1 þa2 ea3 sn ,

ð11Þ

where a1, a2 and a3 are constants determined from curve fitting of the test data; the unit of sn is MPa here.

The conceptual model shown in Fig. 1 a fracture model for fluid flow in rock fractures with a possible analytical solution (Zimmerman and Bodvarsson, 1996), represented by the wellknown cubic law. With a horizontal fracture of constant crosssection in the x direction, one can argue that the fluid velocity is negligible in the y and z directions. With these simplifications and the ‘‘no-slip’’ boundary condition, the relation between the flow rate and the hydraulic pressure gradient can be represented by the cubic law (Bear, 1972): Qx ¼ w

2b3 dp , 3m dx

ð13Þ

where Qx [L3/T] is the volumetric flow rate, w [L] is the fracture width, m (Pa s) is the dynamic viscosity, x [L] is the coordinate along the fracture axis, and dp/dx [MT2/L2] is the hydraulic pressure gradient. The average flow velocity of a steady flow state flow can be obtained by dividing the flux, Qx, by the cross-sectional area, 2wb [L2]: v¼

Qx b2 dp : ¼ 2wb 3m dx

ð14Þ

This equation links the change of flow velocity in fracture with the stress-dependent aperture. Furthermore, it transfers the effects of stress to the advective transport in fracture, as well as other transport mechanisms. 3.3. Transport process Generally speaking, contaminant transport and spreading in the fractured porous media are influenced by many causes, such as advective transport along the connected fractures, longitudinal hydrodynamic dispersion within the fractures in the flow direction, matrix diffusion (solute particles going into and out of the micropores in rock matrices), adsorption onto the fracture walls, adsorption within the matrix, or decay for radioactive nuclides. It is expected that some of them are affected significantly by the in situ stresses. Many more details on stress effects on particle transport in single fracture are presented in Section 4. In the past three decades, the analytical solutions for representing solute transport in single fractures for different initial and boundary conditions have been extensively studied and developed. Among them, Tang et al. (1981) derived a general analytical solution for contaminant transport along a discrete fracture seated in a porous rock matrix, considering all the transport mechanisms noted above. Additionally, various solutions in semianalytical forms for different boundary conditions at the inlet were also proposed (Neretnieks, 1980; Rasmuson and Neretnieks, 1981; Moreno and Rasuson, 1986; Berkowitz and Zhou, 1996; Sun and Buschec, 2003; Shih, 2007). Although these models considered many processes, mechanisms, and factors, they did not consider the effect of the stress/displacement process of the fractures and the effects of matrix pore volume changes due to matrix deformation. It is also difficult to identify or explore the contributions of different individual transport mechanisms in a fundamental study. Therefore, in the following section, three transport models and two inlet boundary conditions

Z. Zhao et al. / Computers & Geosciences 37 (2011) 1437–1449

1441

are considered in this research, in order to examine the stress effects on different transport mechanisms. To represent the contaminant transport processes in the conceptual fracture model as shown in Fig. 1, two coupled, onedimensional (1D) governing equations are usually needed: one for the fracture and the other for the porous matrix. The coupling conditions are the continuity of fluxes and concentrations along the fracture surface. In this study, we first assume that the continuous source concentration is constant c0 at the inlet, the same as defined by Tang et al. (1981). The boundary conditions for the concentration along the fracture are

the constrictivity and tortuosity are difficult to measure and to be separated experimentally, an empirical relation of Archie’s law (Archie, 1942) is often used to relate the effective diffusion coefficient to the porosity:

cf ð0,tÞ ¼ c0 ,

ð15aÞ

Da ¼

cf ð1,tÞ ¼ 0,

ð15bÞ

cf ðx,0Þ ¼ 0:

ð15cÞ

where rm [M/L3] is the bulk density of the matrix, and Km is the distribution coefficient that was defined by Freeze and Cherry (1979) as the mass of solute adsorbed per unit volume of matrix divided by the concentration of solute in solution. Neretnieks (1980) clarified the distinction between effective and apparent diffusivities. Da is usually used for describing instationary diffusion (Eq. (18)), whereas De is used in the literature for flux calculation (Eq. (17)). The term (y þ rmKm) is called rock capacity, a. The volume sorption coefficient rmKm is equal to zero if no sorption or retardation occurs. The value of rock capacity a is then equal to that of porosity y; i.e., Da ¼Dm for the conservative contaminant transport. Without detailed derivation, we present only the final solutions of Eqs. (17) and (18) below, under the boundary conditions (Eqs. (15) and (16)) (Carslaw and Jaeger, 1959; Neretnieks, 1980):

The boundary conditions for the concentration of matrix are cm ðx,b,tÞ ¼ cf ðx,tÞ,

ð16aÞ

cm ð1,z,tÞ ¼ 0,

ð16bÞ

cm ðx,z,0Þ ¼ 0,

ð16cÞ

where cf [M/L3] is the volumetric concentration of solute in the fracture, written as a function of location and time cf(x,t); cm [M/L3] is the volume concentration of solute in the porous matrix, written as cm(x,z,t); N represents the infinite. The boundary condition of constant inlet concentration is the most used one in the modeling, which may corresponds to a contaminant which is dissolved in the fluid to a concentration determined by its solubility in practice (Moreno and Rasuson, 1986). However, the contaminant is transported by advection and hydrodynamic dispersion in the fractures, and the concentration at the locations close to the inlet can be significantly influenced by the boundary conditions when the value of the Pe´clet number is small. In other words, much more contaminant mass may be introduced into the fractures due to the large dispersive flux at the inlet boundary when dispersion is strong. In order to meet the mass conservation condition, the constant flux (Danckwerts’) inlet boundary conditions were proposed (Danckwerts, 1953), and the corresponding solutions were derived by Rasmuson (1986) and Moreno and Rasuson (1986). They are presented as transport model III in later section. 3.3.1. Transport model I The simplest particle transport model in a single fracture is the one only considering the advective transport in the flow direction and matrix diffusion, without dispersion or adsorption. Considering the mass balance of solute in the fracture and in the pores of the matrix, respectively, one can obtain the following differential equations of solute transport (Neretnieks, 1982):  @cf @cf De @cm  þv  ¼ 0 ðfor fractureÞ, ð17Þ @t @x b @z z ¼ b @cm @2 cm ¼ Da @t @z2

ðfor rock matrixÞ,

ð18Þ

where v [L/T] is the mean fluid velocity in the fracture (Eq. (14)); De [L2/T] is the effective diffusion coefficient, which is related to the pore diffusion coefficient Dm [L2/T] and the rock matrix porosity in the following form: De ¼ yDm :

ð19Þ

Neretnieks (1980) defined the pore diffusion coefficient as Dm ¼ D ðd=o2 Þ, where D* [L2/T] is the molecular diffusion coefficient in the fluid, d is the pore constrictivity ( o1), and o is the tortuosity, which is a function of rock matrix porosity. Because

m

De ¼ D y :

ð20Þ

The exponent m can take the values between 1.3 and 4 for different rocks but is often taken to be 1.6 for granites (Ruffet 1:6 0:6 et al., 1995). Therefore, one has De ¼ y D and Dm ¼ y D . 2 Da [L /T] is the apparent diffusion coefficient, which is defined as De , ðy þ rm Km Þ

tw De

cf b Da0:5 ¼ erfc c0 2ðttw Þ0:5

ð21Þ

! ðfor fractureÞ,

0:5 cm ðtw =bÞðDe =D0:5 a Þ þ zb=Da ¼ erfc c0 2ðttw Þ0:5

ð22Þ

! ðfor matrixÞ,

ð23Þ

where tw is water residence time, and is equal to tw ¼ x=v.

3.3.2. Transport model II Based on the transport model I presented above, we added the hydrodynamic dispersion term in the fracture. The governing equation for the solute transport process of a fracture then becomes  @cf @cf @2 cf De @cm  þv Df  ¼ 0, @t @x b @z z ¼ b @x2

ð24Þ

where Df [L2/T] is the hydrodynamic dispersion coefficient, which describes the spreading of a tracer pulse by local variations in velocity and molecular diffusion. It can be expressed as Df ¼ aL vþ D , where aL [L] is the dispersivity in the direction of the fracture axis (Bear, 1972). Detwiler et al. (2000) proposed a complex theoretical model of Df dependent on the Pe´clet number (Pe), but we use the first one for simplicity in this study. Together with 2 Eq. (18), one obtains the following solutions for t 4 x2 =4x DL : !  Z 1 cf 2 v x2 v2 2 ¼ pffiffiffiffi exp x exp x  2 2Df c0 p l 16x Df 2 ! x2 De 1 pffiffiffiffiffiffi erfc dx ðfor fractureÞ, 2D 2 0:5 2 D 4bx f a 2ðtðx =4x Df ÞÞ  Z cm 2 v ¼ pffiffiffiffi exp x 2Df c0 p l

1

2

exp x 

x2 v2 2

16x Df 2

!

ð25Þ

1442

Z. Zhao et al. / Computers & Geosciences 37 (2011) 1437–1449

0

2 D B4bxx2 D pe ffiffiffiffi Da B f

erfcB  @ 2 t

1 þ Dzb 0:5 C a C dx 0:5 C A

ðfor matrixÞ:

ð26Þ

x2 2 4x D f

3.3.3. Transport model III As noted before, constant inlet concentration boundary condition is an approximation only, if the hydrodynamic dispersion is present (Rasmuson, 1986). The correct one ensuring the mass conservation is the constant flux (Danckwerts’) boundary at inlet, which is expressed as (Danckwerts, 1953) vc0 ¼ vcf ð0,tÞDf

@ c ð0,tÞ: @x f

ð27Þ

It is corresponding to a constant flux vc0 at x ¼0 and there is no hydrodynamic dispersion for xo0 (Moreno and Rasuson, 1986). The solutions for the governing equations of Eqs. (24) and (18) under the constant inlet flux boundary condition (Eq. (27)) were derived by Rasmuson (1986) and Moreno and Rasuson (1986):  Z 1   cf v v v ¼ exp x exp  x D Df Df c0 x " f !  Z 1 2 v x2 v2 2 pffiffiffiffi exp x exp x  2 2Df p l 16x Df 2 0 1 3 x2 De 1 pffiffiffiffiffiffi B 7 0:5 C  erfc@4bx2 Df Da  Adx5dy ðfor fractureÞ, 2 2 t x2 4x Df

ð28Þ  Z 1   cm v v v ¼ exp x exp  x Df Df Df c0 x

"

2  pffiffiffiffi exp

p

0



v x 2Df

Z

1 l

1

2

x2 v2

2

exp x 

!

2

16x Df 2

3

De ffiffiffiffi x p þ zb 0:5 2 B4bx Df Da Da C 7  0:5 Cdx7dy erfcB @ A 5 2 2 t x2

ðfor matrixÞ:

ð29Þ

4x Df

3.4. Solutions for coupled stress-flow-transport process Combining Eqs. (8) and (14), we obtain the average velocity of groundwater for a steady-state flow in a single fracture as 2 dp 1  : v¼ bres þðunmax un Þ þ ud 12m dx

determine the actual integration range. The equations can also be solved directly by Mathmatica or other commercial computation software.

4. Model behavior verification Because there are no available experimental results published in the existing literatures and no laboratory experiments considering coupled stress-flow and transport processes with matrix diffusion and dispersion were ever performed, we apply the above solutions to an example to demonstrate the effect of stress on the contaminant transport of a rock fracture in a porous rock matrix. We test the sensitivity of the different transport mechanisms, such as advection, matrix diffusion, and longitudinal hydrodynamic dispersion, to the normal closure, shear dilatancy, and stress-dependent matrix porosity, respectively. The material properties are given in Table 1. The three transport models with different longitudinal dispersivities under the initial condition of zero stress applied after 10 years are compared in Fig. 4. Without normal stress applied, the initial fracture aperture was 30 mm, and the groundwater velocity in the fracture was 7.5  10  7 m/s for a hydraulic pressure gradient of 10 Pa/m (or water head gradient of 0.001 m/m). If only considering the pure advection in the fracture, the contaminant would be expected to travel far enough after 10 years, to make contaminant concentration equal to c0 in the fracture. However, the relative concentration in the fracture decreased to a very small value after x¼8 m, if the matrix diffusion is considered, indicating that diffusion of the contaminant (for example, radioactive nuclides) into the rock matrix probably is one of the more significant transport mechanisms, especially for final safety assessment of a potential repository of radioactive waste (Neretnieks, 1980). For larger longitudinal dispersivity of aL ¼1 m, the longitudinal hydrodynamic dispersion increases the contaminant concentration in the fracture or in the rock matrix at the initial state under the constant concentration inlet boundary (Fig. 4a and c), which is in line with the results shown by Tang et al. (1981). However, under constant flux boundary conditions, the relative solute concentration estimated by transport model III is smaller than that from transport model I in the locations closer to the inlet (x o3 m), due to the influence of longitudinal hydrodynamic dispersion. After that part, longitudinal hydrodynamic dispersion

Table 1 Material properties of groundwater, rock matrix, and fracture.

ð30Þ

Then substituting the updated half aperture b of the fracture (Eq. (8)), the mean fluid velocity v obtained from Eq. (30), and the effective and apparent diffusion coefficients (Eqs. (20) and (21)) due to the changed matrix porosity represented by Eq. (12) into the three transport models presented above (Eqs. (22), (23), (25), (26), (28), and (29)), one can have the solutions for the coupled stress-flow-transport process in the fracture–matrix system as shown in Fig. 1, and examine the effect of stress on solute transport in the single fracture. We did not present the complete expression here to save space. To evaluate the integrals in the transport models II and III (Eqs. (25), (26), (28), and (29)), we used the Gauss–Legendre quadrature technique with 100 Gauss points. It was found that the numerically significant portion generally extends over a smaller range, even though the integrand theoretically extends from l to infinity (Tang et al., 1981). In this way, before the integration, a prior scanning procedure was carried out to

Material parameters

Value

Groundwater Hydraulic pressure gradient, dp/dx (Pa/m) Dynamic viscosity, m (Pa s)

10 1.0  10  3

Intact rock Initial matrix porosity, y Empirical constants, a1, a2, a3 (  10  4,  10  4, 1) Cementation factor, n

0.01 0.82, 5.35, 0.12 1.6

Fractures Normal stiffness at zero stress state, kn0 (GPa/m) Shear stiffness, ks (GPa/m) Dispersivity along fracture, aL (m) Diffusivity in water, Dn (m2/s) Frictional angle, | (deg) Cohesion, c (MPa) Critical shear displacement for dilation, Ucs (mm) Shear dilation angle, fd (deg) Residual aperture, bres (mm) Maximum normal displacement, Dunmax (mm)

434 434 0.5 1.6  10  9 24.9 0 3 5.0 5 25

Z. Zhao et al. / Computers & Geosciences 37 (2011) 1437–1449

1443

αL=1m

αL=0.1m 1.0

Model I Model II Model III

0.8

R e la tive c o n c e n tra tio n c /c0

For fracture

R e la tive c o n c e n tra tio n c /c0

1.0

0.6

0.4

0.2

0.0

Model I Model II Model III

0.8

0.6

0.4

0.2

0.0

0

2

4

6

8

0

2

Distance along fracture (m)

6

8

1.0

Model I Model II Model III

0.8

R e la tive c o n c e n tra tio n c /c0

For matrix

1.0

R e la tive c o n c e n tra tio n c /c0

4

Distance along fracture (m)

0.6

0.4

0.2

0.0

Model I Model II Model III

0.8

0.6

0.4

0.2

0.0

0.0

0.1

0.2

0.3

0.4

Distance into matrix (m)

0.5

0.0

0.1

0.2

0.3

0.4

0.5

Distance into matrix (m)

Fig. 4. Comparison among the different transport models under the initial conditions without stress applied. Two longitudinal dispersion coefficients (aL ¼1 and 0.1 m) were assumed to study the sensitivity of the transport model to the inlet boundary conditions.

increases the contaminant concentration in the fracture (Fig. 4a). The same phenomenon exists in the matrix for the constant flux inlet boundary also (Fig. 4c). Three transport models give similar concentration profiles for smaller longitudinal dispersivity of aL ¼0.1 m (Fig. 4b and d), which indicates that the longitudinal dispersion plays a negligible role in this case. The concentration discontinuity at the inlet for transport model III is also illustrated in Fig. 4a. The effects of constant concentration or flux inlet boundary conditions are demonstrated here, and the same trend can be found in the following results with stress applied. (1) Normal closure effect by pure normal compression: In this case, we mainly study the effect of the fracture aperture closure under an increasing normal stress on the mass transport along the fracture and into the rock matrix, without considering the shear dilatancy and matrix porosity changes. The issues of the shear dilatancy and matrix porosity change are treated in the later sections. Fig. 5 shows the concentration profiles along the fracture and in the rock matrix (x¼2 m) at a short time scale (t¼10 years), under two normal stresses of 5 and 10 MPa, respectively. The results illustrate that the contaminant concentrations both along the fracture and in the rock matrix decrease with the increasing normal stress value, regardless of what transport models are used. For transport model III, the concentration at the inlet decreases with increasing stress too. The fracture’s normal closure reduces the transport activity in the fractured rocks. However, the effect of hydrodynamic dispersion becomes more significant during the

process of increasing normal stress, and it increases the contaminant concentration along the fracture and especially in the rock matrix. In Fig. 5c and d, the relative concentration from the transport model III is higher than that from transport model I, indicating indirectly an important role of dispersion in the fractures. Without considering the hydrodynamic dispersion, it is almost impossible for particles to diffuse into rock matrix under compressive stress of 10 MPa at this short time scale (t¼10 years). Fig. 6 presents the same information as in Fig. 5, except that it is at a long time scale (t ¼50 years). With the longer time, the contaminant concentrations both along the fracture and in the rock matrix further increase (compared with the results at the time scale of 10 years), and the effects of stress become more obvious. However, the general trend remains the same. Fig. 6c shows that the relative concentration estimated by transport model I without dispersion becomes larger than that of transport model III after a longer time, but the relative concentration of transport model I is still smaller than that of transport model III in Fig. 6d. Considering the relative concentration in the fracture at x¼2 m, transport model I gives a larger value in Fig. 6a than that of transport model III, and inverse in Fig. 6b. It is concluded that larger concentrations in the fracture induce larger concentrations in the rock matrix, and after a long enough time the concentration of transport model I will be larger than that of transport model III, both in fracture and in matrix. Figs. 7 and 8 show the influence of the normal stress on breakthrough curves for the fracture (x¼2 m) and the matrix

1444

Z. Zhao et al. / Computers & Geosciences 37 (2011) 1437–1449

5 MPa

10MPa 1.0

Model I Model II Model III

0.8

R e la tive c o n c e n tra tio n c /c 0

For fracture

R e la tive c o n c e n tra tio n c /c 0

1.0

0.6

0.4

0.2

0.0

Model I Model II Model III

0.8

0.6

0.4

0.2

0.0

0

2

4

6

8

0

2

Distance along fracture (m)

Model I Model II Model III

0.8

R e la tive c o n c e n tra tio n c /c 0

For matrix

6

8

1.0

1.0

R e la tive c o n c e n tra tio n c /c 0

4

Distance along fracture (m)

0.6

0.4

0.2

Model I Model II Model III

0.8

0.6

0.4

0.2

0.0

0.0 0.0

0.1

0.2

0.3

0.4

Distance into matrix (m)

0.5

0.0

0.1

0.2

0.3

0.4

0.5

Distance into matrix (m)

Fig. 5. Concentration profiles along the fracture and into the matrix (x ¼2 m) at the short time scale (t¼ 10 years), under normal stress of 5 and 10 MPa.

(x¼2 m, z¼0.1 m), respectively. Under the initial condition of no stress, the breakthrough curves for models I and II overlap very much after longer times, indicating that the influence of hydrodynamic dispersion is negligible then. However, at the beginning period of transport, relative concentration increases more quickly with longitudinal dispersion under constant concentration inlet conditions. Due to less contaminant mass introduced under constant flux inlet conditions, the relative concentrations of transport model III are lower than the others after a long time. The normal stress decreases the fracture aperture and flow rate, and in turn the rate of contaminant concentration in the fracture and the rock matrix. When a normal stress of 5 MPa is applied, the hydrodynamic dispersion in transport model II increases the rate of contaminant concentration in the fracture and the rock matrix compared with transport model I. The same phenomena happen for transport model III after the transport start immediately, but its relative concentration is smaller after a long time. From the above results, it can be seen that the fracture aperture closure induced by normal stress plays an important role in the flow and transport in fractured rocks. It is expected that the extent of fracture closure would become smaller with further increase of the normal stress, because the normal stiffness increases with the increasing normal stress. In this way, the impact of normal stress would decrease under high normal stresses, especially when fracture aperture is close to its residual aperture value. The influence of hydrodynamic dispersion

becomes more significant, when the fluid velocity (or fracture aperture) becomes smaller. (2) Shear dilatancy effect: To examine the effect of shear displacement of the fracture, we fix the normal stress at 5 MPa, and increase the shear displacement gradually to produce shear dilatancy. From Eqs. (4) and (7), the fracture aperture would start to increase by shear dilatancy, when shear displacement us becomes larger than 5.3 mm. Fig. 9 shows the concentration profiles along the fracture and in the rock matrix (x ¼2 m) at the short time scale (t ¼10 years), with the shear displacement of 0.1 and 0.2 mm, respectively. The influence of increasing shear displacement is just the inverse of that of increasing normal stress. The shear dilatancy increases fracture aperture, then flow rate, and finally the rate of contaminant concentration in the fracture and the rock matrix. The influence of inlet boundary conditions becomes less significant when large shear dilatancy happens, because the advection becomes more significantly induced from increasing fluid velocity. (3) Effect of matrix porosity change: Under normal stress conditions, the influences of matrix porosity change on transport mechanisms are investigated by comparison between the results with and without considering matrix porosity changes, under the compressive normal stress of 5 and 10 MPa, respectively (Fig. 10). For the concentration profiles along the fracture, if we consider the matrix porosity change under normal stresses, the concentration along the fracture increases significantly, compared with the

Z. Zhao et al. / Computers & Geosciences 37 (2011) 1437–1449

1445

5 MPa

10 MPa 1.0

Model I Model II Model III

0.8

R e la tive c o n c e n tra tio n c /c 0

For fracture

R e la tive c o n c e n tra tio n c /c 0

1.0

0.6

0.4

0.2

Model I Model II Model III

0.8

0.6

0.4

0.2

0.0

0.0 0

2

4

6

0

8

2

6

8

1.0

Model I Model II Model III

0.8

R e la tive c o n c e n tra tio n c /c 0

For matrix

1.0

R e la tive c o n c e n tra tio n c /c 0

4

Distance along fracture (m)

Distance along fracture (m)

0.6

0.4

0.2

Model I Model II Model III

0.8

0.6

0.4

0.2

0.0

0.0 0.0

0.1

0.2

0.3

0.4

0.0

0.5

0.1

0.2

0.3

0.4

0.5

Distance into matrix (m)

Distance into matrix (m)

Fig. 6. Concentration profiles along the fracture and into the matrix (x¼ 2 m) at the long time scale (t ¼50 years), under normal stress of 5 and 10 MPa.

R e la tive c o n c e n tra tio n c /c0

1.0

0.8

0.6

0.4

Model I Model II Model III

0.2

0.0 0

20

40

60

80

100

120

140

160

Time (years) Fig. 7. Comparison of breakthrough curves for the fracture at x ¼ 2 m, between the initial case with zero normal stress (solid lines) and the case under the normal stress of 5 MPa (dash lines).

1446

Z. Zhao et al. / Computers & Geosciences 37 (2011) 1437–1449

1.0

R e la tive c o n c e n tra tio n c /c0

0.8

0.6

0.4

Model I

0.2

Model II Model III

0.0 0

20

40

60

80

100

120

140

160

Time (years) Fig. 8. Comparison of breakthrough curves for the matrix (x¼ 2 m, z¼ 0.1 m), between the initial case with zero normal stress (solid lines) and the case under the normal stress of 5 MPa (dash lines).

0.1mm

0.2mm 1.0

Model I Model II Model III

0.8

Model I Model II Model III

0.8

R e la tive c o n c e n tra tio n c /c0

For fracture

R e la tive c o n c e n tra tio n c /c0

1.0

0.6

0.4

0.2

0.0

0.6

0.4

0.2

0.0

0

2

4

6

8

0

2

Distance along fracture (m)

6

8

1.0

Model I Model II Model III

0.8

R e la tive c o n c e n tra tio n c /c0

For matrix

1.0

R e la tive c o n c e n tra tio n c /c0

4

Distance along fracture (m)

0.6

0.4

0.2

0.0

Model I Model II Model III

0.8

0.6

0.4

0.2

0.0

0.0

0.1

0.2

0.3

Distance into matrix (m)

0.4

0.5

0.0

0.1

0.2

0.3

0.4

0.5

Distance into matrix (m)

Fig. 9. Concentration profiles along the fracture and into the matrix (x ¼2 m) at the short time scale (t¼ 10 years), under shear displacements of 0.1 and 0.2 mm.

Z. Zhao et al. / Computers & Geosciences 37 (2011) 1437–1449

1447

5 MPa

10 MPa 1.0

Model I Model II Model III

0.8

R e la tive c o n c e n tra tio n c /c 0

For fracture

R e la tive c o n c e n tra tio n c /c 0

1.0

0.6

0.4

0.2

Model I Model II Model III

0.8

0.6

0.4

0.2

0.0

0.0 0

2

4

6

0

8

2

6

8

1.0

Model I Model II Model III

0.8

R e la tive c o n c e n tra tio n c /c 0

For matrix

1.0

R e la tive c o n c e n tra tio n c /c 0

4

Distance along fracture (m)

Distance along fracture (m)

0.6

0.4

0.2

0.0

Model I Model II Model III

0.8

0.6

0.4

0.2

0.0

0.0

0.1

0.2

0.3

0.4

Distance into matrix (m)

0.5

0.0

0.1

0.2

0.3

0.4

0.5

Distance into matrix (m)

Fig. 10. Comparison of concentration profiles for the fracture and matrix (x¼ 2 m) between the cases with (dash lines) and without (solid lines) considering the matrix porosity change under normal stress conditions, at the short time scale (t ¼10 years).

case of constant matrix porosity. Under higher normal stress, this difference becomes more obvious (Fig. 10a and b). The decreasing matrix porosity caused by increasing the normal stress, inducing decreasing diffusion coefficients, makes the particles more difficult to diffuse into the rock matrix. The concentration in the rock matrix also increases with decreasing matrix porosity, due to the increasing normal stress, compared with the case with constant matrix porosity (Fig. 10c and d). The possible reason is that larger contaminant concentration gradient at the fracture wall tends to push more particles into the matrix pores.

5. Discussion and conclusions An analytical solution for coupled stress-flow-transport processes in a single rock fracture was developed to consider the influence of stress/displacement on the contaminant transport through individual rock fractures. This model explicitly integrates the mechanical, flow, and transport processes in a single rock fracture into a compact mathematical form, by linking the stress and transport processes through changes of aperture and matrix porosity, based on the cubic law and Archie’s law, respectively. Therefore, it can be directly used as a simple constitutive model for modeling the coupled hydromechanical-transport processes in fracture networks in crystalline rocks. However, it is noted here that the validity of this model is

based on the assumptions listed in Section 2. Some important outstanding issues are addressed below. The 2D geometrical model was used mainly for its simplicity, with homogeneous material properties for the rock matrix, in order that a closed-form solution is possible. The stress/deformation relation of the matrix, the groundwater flow velocity, and contaminant transport along fracture or into rock matrix are assumed to be 1D. Their validities depend on the configuration of the particular rocks and fractures. For example, if the rock matrix is heterogeneous with large porosity, the contaminant distribution in the rock matrix may become a 2D problem, and the groundwater flow through porous matrix should not be neglected. Therefore, it is emphasized here that the present model is only appropriate for rocks with small matrix porosity. If there are significant groundwater exchanges between the micropores in matrix and fractures, the model presented in this paper should not be directly used. Across the fracture section, the groundwater velocity and contaminant concentration are assumed to be uniform, following the parallel plate model of a constant aperture. For fractures with variable apertures, the above assumption is not valid. However, for fractures with aperture of micrometer scale, under relatively lower hydraulic pressure gradient and at the time scale of performance/safety assessments for underground radioactive waste repositories, the above assumption may still be applicable.

1448

Z. Zhao et al. / Computers & Geosciences 37 (2011) 1437–1449

The model in this paper was derived based on some empirical mechanical models, so its validity largely depends on the range of validity of those empirical models and associated parameters. Actually, the main contribution of this paper is to highlight the importance of the coupled stress-flow-transport processes in fractured crystalline rocks. Any other mechanical models can replace those models used here, as long as they can consider the changes of fracture aperture and matrix porosity under stress, explicitly in compact mathematical function forms, to meet the requirements of special rocks or sites, and the resultant equations can be integrated analytically. There are a large number of material parameters and empirical constants in the present model. How to obtain the reasonable values for them from laboratory or in situ experiments remains to be a difficult problem and no attempt has been made to test coupled stress-flow-transport processes in a single rock fractures with quantitative representations of transport mechanisms, variables, and parameters, so far, due to experimental difficulty. The uncertainties in measuring the material properties have an important influence on the final assessments of radioactive waste repositories. Moreover, these properties and parameters always vary among different sites. Constant concentration inlet boundary results in large dispersive flux at the inlet boundary for large hydrodynamic dispersion coefficient at the beginning of transport, and more solute mass is injected in the fracture and rock matrix after a long time. Constant flux boundary conditions, however, can eliminate the extra dispersive flux at the inlet, which is more proper physically, but induces an artificial discontinuity at the inlet boundary. If the longitudinal hydrodynamic dispersion is neglected, the constant flux boundary condition becomes the same as the constant concentration boundary condition. The different results between these two boundary conditions indicate that it is a crucial factor for modeling. The concentrations noted in this paper are resident concentrations, and more information on the transformation and to distinguish between flux concentrations and resident concentrations can be found in Kreft and Zuber, (1978), Van Genuchten and Parker (1984), Rasmuson (1986), and Moreno and Rasuson (1986). In the transport models presented in this paper, adsorption on the walls of fracture and within the matrix pores, and the radioactive decay were not included, but there are many available transport models considering these issues in the literature (Tang et al., 1981; Sudicky and Frind, 1984; Cormenzana, 2000; Sun and Buschec, 2003). The basic mechanical and flow models presented here can be directly coupled with other transport models considering adsorption and decay without major difficulties. The present model can also be further extended to a system with multiple parallel fractures or other initial boundary conditions, with the known analytical transport solutions (Sudicky and Frind, 1982; West et al., 2004). From the present study, stresses do have a significant influence on the groundwater flow and contaminant transport in rock fractures, through changing the fracture aperture and matrix porosity. The fracture closure under normal stress and the dilatancy during shear are the two major mechanisms controlling the fluid flow and transport in the fractures, and consequently influence the groundwater velocity in fractures, contaminant concentration distribution, and breakthrough time. The results show that matrix porosity changes corresponding to actual stress conditions should not be neglected in practice, even though we assumed that the porosity was really small in this study and only matrix diffusion was considered. Some basic concluding remarks are presented below. (1) Fracture’s normal closure decreases the contaminant concentrations both along the fracture and into the rock matrix,

but the dilatancy increases the corresponding contaminant concentrations. (2) The importance of the hydrodynamic dispersion depends on the fluid velocity, inlet boundary conditions, and time scales. (3) With decreasing matrix porosity, the concentrations both along the fracture and in the rock matrix increase. (4) The transport models should be chosen according to practical conditions, after the nonnegligible transport mechanisms are identified first. The above outstanding issues and limitations are identified for future research. Among them, the validity range of the parallel plate model of the fractures, the effect of fracture roughness, and the influence of water pressure in fracture and pore pressure in micropores of matrix will be considered in future studies.

Acknowledgment We acknowledge the financial support from the Swedish Nuclear Fuel and Waste Management Co. (SKB) through the international DECOVALEX-2011 project.

References Archie, G.E., 1942. The electrical resistivity log as an aid in determining some reservoir characteristics. Petroleum Technology 146, 54–62. Bandis, S., Lumsden, A.C., Barton, N.R., 1983. Fundamentals of rock joint information. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts 20, 249–268. Barton, N.R., Bandis, S., Bakhtar, K., 1985. Strength, deformation and conductivity coupling of rock joints. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts 22, 121–140. Bear, J., 1972. Dynamics of Fluids in Porous Media. Elsevier, New York, 764 pp. Berkowitz, B., Zhou, Z., 1996. Reactive solute transport in a single fracture. Water Resource Research 32 (4), 901–913. Bodin, J., Delay, F., de Marsily, G., 2003a. Solute transport in a single fracture with negligible matrix permeability: 1. Fundamental mechanisms. Hydrogeology Journal 11, 418–433. Bodin, J., Delay, F., de Marsily, G., 2003b. Solute transport in a single fracture with negligible matrix permeability: 2. Mathematical formalism. Hydrogeology Journal 11, 434–454. Carslaw, H.S., Jaeger, J.C., 1959. Conduction of Heat in Solids, New York, 2nd ed. Oxford University Press, 510 pp. Cook, N.G.W., 1992. Natural joints in rock: mechanical, hydraulic and seismic behaivor and properties under normal stress. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts 29 (3), 198–223. Cormenzana, I., 2000. Transport of a two-member decay chain in a single fracture: simplified analytical solution for two radionuclides with the same transport properties. Water Resource Research 36 (5), 1339–1346. Danckwerts, P.V., 1953. Continuous flow systems: distribution of residence times. Chemical Engineering Science 2, 1–13. Detwiler, R.L., Rajaram, H., Glass, R.J., 2000. Solute transport in variable-aperture fractures: an investigation of the relative importance of Taylor dispersion and macrodispersion. Water Resource Research 36 (7), 1611–1625. Freeze, R.A., Cherry, J.A., 1979. Groundwater. Prentice-Hall, Englewood Cliffs, NJ 604 pp. Goodman, R.E., 1976. Methods of Geological Engineering in Discontinuous Rocks. West Publishing Company, San Francisco, 472 pp. Han, G., Dusseault, B.M., 2003. Description of fluid flow around a wellbore with stress-dependent porosity and permeability. Journal of Petroleum Science and Engineering 40, 1–16. Jaeger, J.C., Cook, N.G.W., Zimmerman, R.W., 2007. Fundamentals of Rock Mechanics, 4th ed. Blackwell Publishing, 488 pp. Jing, L., 2003. A review of techniques, advances and outstanding issues in numerical modelling for rock mechanics and rock engineering. International Journal of Rock Mechanics and Mining Sciences 40, 283–353. Jing, L., Stephansson, O., 2007. Fundamentals of Discrete Element Methods for Rock Engineering—Theory and Application. Elsevier Science Publishers, Rotterdam 562 pp. Kreft, A., Zuber, A., 1978. On the physical meaning of the dispersion equation and its solutions for different initial and boundary conditions. Chemical Engineering Science 33, 1471–1480. Moreno, L., Rasuson, A., 1986. Contaminant transport through a fissured porous rock: impact of the inlet boundary condition on the concentration profile in the frock matrix. Water Resource Research 22 (12), 1728–1730.

Z. Zhao et al. / Computers & Geosciences 37 (2011) 1437–1449

Moreno, L., Tsang, Y.W., Tsang, C.F., Hale, F.V., Neretnieks, I., 1988. Flow and tracer transport in a single fracture: a stochastic model and its relation to some field observations. Water Resource Research 24 (12), 2033–2048. Neretnieks, I., 1980. Diffusion in the rock matrix: an important factor in radionuclide retardation? Journal of Geophysical Research 85 (B8), 4379–4397. Neretnieks, I., 1982. Tracer movement in a single fissure in granitic rock: some experimental results and their interpretation. Water Resource Research 18 (4), 849–858. Pyrack-Nolte, L.J., Myer, L.R., Cook, N.G.W., Witherspoon, P.A., 1987. Hydraulic and mechanical properties of natural fractures in low-permeability rock. In: Herget, G., Vongpaisal, S. (Eds.), Proceedings 6th International Congress on Rock Mechanics. Montreal, Canada, pp. 225–232. Rasmuson, A., 1986. Exact solution of some models for the dynamics of fixed beds using Danckwerts’ inlet condition. Chemical Engineering Science 41 (3), 599–600. Rasmuson, A., Neretnieks, I., 1981. Migration of radionuclides in fissured rock: the influence of micropore diffusion and longitudinal dispersion. Journal of Geophysical Research 86 (B5), 3749–3758. Raven, K.G., Gale, J.E., 1985. Water flow in a natural rock fracture as a function of stress and sample. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts 22 (4), 251–261. Renshaw, C.E., 1995. On the relationship between mechanical and hydraulic apertures in rough-walled fractures. Journal of Geophysical Research 100 (B12), 24629–24636. Ruffet, C., Dardot, M., Gueguen, Y., 1995. Surface conductivity in rocks: a review. Surveys in Geophysics 16 (1), 83–105. Rutqvist, J., Stephansson, O., 2003. The role of hydromechanical coupling in fractured rock engineering. Hydrogeology Journal 11, 7–40. Shih, D.C.F., 2007. Contaminant transport in one-dimensional single fractured media: semi-analytical solution for three-member decay chain with pulse and Heaviside input sources. Hydrological Processes 21, 2135–2143.

1449

Snow, D.T., 1965. A parallel model of fractured permeable media. Ph.D. dissertation. University of California, 331 pp. Sudicky, E.Q., Frind, E.O., 1982. Contaminant transport in fractured porous media: analytical solutions for a system of parallel fractures. Water Resource Research 18 (6), 1634–1642. Sudicky, E.A., Frind, E.O., 1984. Contaminant transport in fractured porous media: analytical solutions for a two-member decay chain in a single fracture. Water Resource Research 20 (7), 1021–1029. Sun, Y., Buschec, T.A., 2003. Analytical solutions for reactive transport of N-member radionuclide chains in a single fracture. Journal of Contaminant Hydrology 62–63, 695–712. Swan, G., 1983. Determination of stiffness and other joint properties from roughness measurements. Rock Mechanics and Rock Engineering 16, 19–38. Tang, D.H., Frind, E.O., Sudicky, E.A., 1981. Contaminant transport in fractured porous media: analytical solution for a single fracture. Water Resource Research 17 (3), 555–564. Van Genuchten, M.Th., Parker, J.C., 1984. Boundary conditions for displacement experiments through short laboratory soil columns. Soil Science Society of America Journal 48, 703–708. West, M.R., Kueper, B.H., Novakowski, K.S., 2004. Semi-analytical solutions for solute transport in fractured porous media using a strip source of finite width. Advances in Water Resources 27, 1045–1059. Witherspoon, P.A., Amick, C.H., Gale, J.E., Iwai, K., 1979. Observations of a potential size effect in experimental determination of the hydraulic properties of fractures. Water Resource Research 15, 1142–1146. Zimmerman, R.W., 1991. Compressibility of Sandstones, Developments in Petroleum Science, vol. 29. Elsevier, New York, 181 pp. Zimmerman, R.W., Bodvarsson, G.S., 1996. Hydraulic conductivity of rock fractures. Transport in Porous Media 23, 1–30.