Influence of inherent and stress-induced anisotropy of hydraulic conductivity on the
ground water flow around a rock tunnel
Huang-Kuei Chu, Po-Sung Lai, Jia-Jyun Dong
Graduate Institute of Applied Geology, National Central University, Taiwan (R.O.C.)(No. 300, Zhongda Rd., Zhongli District, Taoyuan City 32001, Taiwan (R.O.C.),
Key words; inherent anisotropy, stress-induced anisotropy, tunnel excavating, Oda model, JRC-JCS model

Abstract
The hydraulic conductivity around the disposal tunnel is one of the key parameters for the safety assessment ofradioactive waste disposal. This study aims to explore the inherent anisotropy (orientation of the discontinuities)
and stress induced anisotropy of the hydraulic conductivity around a rock tunnel. JRC-JCS model is used toestimate the aperture of discontinuities under stress. Based on the calculated stress field via Kirsch solution andthe equivalent continuum model, the hydraulic conductivities around a circular tunnel can be calculated. The
groundwater inflow of the tunnel is further evaluated via finite difference method. The result shows that thehydraulic conductivity on the tunnel wall is about 1 ~ 2 orders of magnitude larger than the one away from thetunnel (or the one of rock mass under boundary stress). The major principal hydraulic conductivity on the tunnel
wall can be 6 ~ 9 times larger than the minor principal value. The principle directions of the hydraulicconductivity near the tunnel wall are also significantly deviated from the tangential and radial directions when
the inherent anisotropy is considered. Groundwater flow analysis shows that the total head and the flow velocity
are dominated by the inherent and stress induced anisotropy of hydraulic conductivity. Surprisingly, the inflowof the tunnel is insignificantly influenced by the spatial variation of hydraulic conductivity around the tunnel
wall.

1. Introduction
Radioactive waste is usually deposited in rocktunnels. The potential of nuclear species migrating
through the rock and the impacts on the biosphere areimportant issues for tunnel disposal. It is important toassess the permeability of the rock mass around thetunnel for evaluating the groundwater flow field. Dueto the presents of primary and secondary
discontinuities, rock masses frequently demonstrateinherent anisotropy of hydraulic conductivity. In
addition, the stress also affects the aperture of the
fissure, which in turn affects the hydraulicconductivity of the rock mass. Therefore, thehydraulic conductivity of rock mass could beanisotropic if the stress tensor is not isotropic.
Oda (1986) proposed a quasi-continuum model toevaluate the hydraulic conductivity of rock mass. The
equivalent hydraulic conductivity of rock mass isfunction of fracture density, persistency, aperture,and orientation of the fractures. Since the aperture offractures can be evaluated using the stiffness andstresses acting on the fractures, the influences ofinherent anisotropy and stress induced anisotropy canbe considered (Cheng, 2006).
In this study, we consider the inherent and stressinduced anisotropy of rock mass around the tunnel
wall, where the stress redistribution is significant.
Finite difference method (FDM) is used to analyze
the flow field around the tunnel. To simplify the
discussion, the boundary stress is assumed as
isotropic. That is, the stress induced anisotropy is
only due to stress redistribution caused by tunnel
excavation.
2. Hydraulic conductivity around tunnel
The hydraulic conductivity of rock mass isevaluated via Oda model and JRC-JCS model, which
are briefly introduced below:(1) Quasi-continuum model
Oda (1985) proposed the concept of fracture tensor,
where the fractures of rock mass is assumed as disk
shape and the intact rocks are impermeable. The
fracture tensor can be expressed as follow Eq. (1):

6th East Asia Forum on Radwaste Management Conference
November 27-29, 2017, Osaka, Japan

Pij 

 tm rm4

0 0  r t nˆi nˆ j E (nˆ ) f (r ) g (t )ddrdt
2 3

(1)

where  is the fracture density, r is the fractureradius, t is the fracture aperture, n̂ is the normal
vector of fracture. The E (nˆ ) , f (r ) , g (t ) are thedensity function of normal vector, radius and apertureof fractures. The hydraulic conductivity of thefractured rock mass can be expressed as function of
fracture density:kij 1PKK  ij  Pij 
12(2)

The inherent anisotropy of hydraulic conductivity
is mainly controlled by the density function of normal
vector of fractures, and the density function can beapproximated by fabric tensor (Kanatani, 1984). In
this study, only the second order fabric tensor wasconsidered. The stress induced anisotropy is mainly
dominated by the different aperture of fracture withdifferent orientation, which influenced by the stresses
acting on the fractures.
(2) JRC-JCS model
JRC0 

 0.2 c  0.1
5 
JCS0


(3)

 c for the intact rock uniaxial compresive strength.
If the discontinuity is fresh, then JCS   c . The
maximum closer value of aperture t max comes
empirical formula for the regression of test data
(Bandis et al., 1983) as follow Eq. (4).

 JCS 0 
t max  0.3  0.006 JRC0  2.24
 t0 

0.25(4)Bandis et al. (1983) consider the relationshipbetween normal stress  n and the aperture closure
t , which can be expressed by hyperbolic curve as
follows Eq. (5):
1
n
k0
t 
1
n
1
k 0 t maxt  t0  t

(6)

The fracture aperture could be varied under
different shear displacement. Based on the fracture
shear testing results, Barton (1982) proposed anormalization procedure to evaluate the aperture
variation due to shear stress. The JRCmob is roughness
coefficient of the discontinuous surface under theaction of shear stress, then shear displacement
induced normal displacement un (shear dilatancy)
can expressed by the following Eq. (7) JCS n   u n   tan  JRC mob log   n 
(7)

JRCmob log JCSn  n  is the dilation angle driven by
shear stress can be estimated by the single pressure
intensity JCS and the normal stress acting on the
discontinuous surface.
(3) Mechanical and hydraulic apertures

Bandis (1982) proposed a JRC-JCS model toevaluate the deformability of discontinuities. The
fracture aperture can be obtained if the roughness
coefficient JRC0 and the uniaxial compressive
strength JCS0 (subscript 0 indicates that the test sizeis 10 cm) are available. The initial aperture (mm) can
be calculated using Eq. (3).
t0 

stiffness, p is length of discontinuous), normal
stress  n , roughness coefficient JRC0 and uniaxial
compressive strength JCS0.

(5)

Finally, we can get the aperture (Eq. (6)) cause bynormal stress which is function of initial normal
stiffness k 0 (  h p , h is constant of normal

The aperture t in Eq. (6) is mechanical aperture.
However, the aperture in Eq. (1) should be hydraulicaperture, where two parallel plates as the fracture
walls were assumed. Based on extensively
experimental data, Barton et al. (1985) proposed theempirical relation between the hydraulic aperture e(m) and mechanical aperture E (m) as follows:E e2e  JRC02.5(8)

JRC0 is the roughness coefficient for the test size of
100 mm of the discontinuities. Eq. (8) is only
applicable to the condition of E > e and this empiricalequation is used in this study.
(4) Kirsch’s solution
To investigate the influence of the stressanisotropy due to the redistribution of the stress after
excavation of the tunnel, Kirsch’s solution (Kirsch,1898) is adopted. Two dimensional stresses around atunnel can be calculated using Eqs. (9) - (11).
r h v2 

 (1 h v

 r  2

 (1 h v
2 v
a2
4  a2 3  a4) h
 (1  2  4 )  cos 2 (9)
2
r
2
r
r

 v
a2
3  a4) h (1  4 )  cos 2 (10)2r2r

 (1 

2  a2 3  a4 4 )  sin 2
r2
r

(11)

In the above equations,  r is the radial stress, is the tangential stress,  r is the shear stress,  v
and  h are the boundary vertical stress and6th East Asia Forum on Radwaste Management Conference
November 27-29, 2017, Osaka, Japanhorizonatl stress, a is the tunnel radius and r is
the distance from the tunnel center (Fig. 1).

Fig. 2 Schematic diagram of coordinate system
Fig. 1 Schematic diagram of stress after tunnel
excavation
The parameters of this study are shown in Table 1.
The boundary stress of the tunnel is assumed to be
isotropic (the vertical stress equals to the horizontal
stress). Radius of the tunnel is 5 meters. Vertical andhorizontal effective stresses are 7.5 MPa (about thestress at a depth of 500 m). Two fabric tensors wereassumed: (1) isotropic distribution of normal vector(D11 = D22 = D33 = 0) and (2) transversely isotropic
distribution of normal vector (D11 = D22 = -1, D33 =2).

3. Results I- Hydraulic conductivity around tunnel
Fig. 3 shows the principal stresses and principal
directions (black cross) of hydraulic conductivity
tensor around the tunnel (at distances from the centerof the tunnel equal to 1.0, 1.5 and 2.0 times of radius)
under an isotropic boundary stress of 7.5 MPa.

Table 1 Parameters and conditions used for
calculating the hydraulic conductivity
Stressstate

Vertical stress  v (MPa)

7.5

Horizontal stress  h (MPa)7.5

Length of discontinuous (m)1

JRC of discontinuous (-)

15

JCS of discontinuous (MPa)100

Volume density (strip/m3)

5(a) isotropic (D33 = 0)

Discontinuous surface normal vector distribution
density function E (nˆ )  1  Dij ni n j(a) D11 = D22 = 0,
D33 = 0(b) D11 = D22 = -1,
D33 = 2Fracturefeatures(b) anisotropic (D33 = 2)
Fig. 3 Distribution of hydraulic conductivity of rockaround tunnel
The coordinate system adopting by this study is asfollows: the axis 1 represents the east, the axis 2
represents the north, the axis 3 represents the upperside and the direction of the tunnel axis (Fig. 2). The
tunnel excavation surface of this study will fall on theaxis 2 axis 3 plane. The vector can be represented by
the space cosine coordinate system (Lee and Farmer,
1993) or in the direction cosine.

Fig. 3 (a) shows the hydraulic conductivity tensorof an isotropic normal vector distribution (D33 = 0).
The principal directions are in the radial andtangential directions, which is accordance to theprincipal stress directions. It is indicated that the eventhe boundary stress and fracture distribution are
isotropic, the hydraulic conductivity is anisotropic for

6th East Asia Forum on Radwaste Management Conference
November 27-29, 2017, Osaka, Japan

the stress redistribution due to tunnel excavation. The
maximum principal direction is in the tangential
direction for the fracture parallel to the radial
direction tends to close more for larger normal stressthan that of the fracture perpendicular the radial
direction. The anisotropy reduced when the location
away from the tunnel wall for the stress tensor
gradually became isotropic. From Fig. 3 (b), theprincipal directions are no longer in the radial andtangential direction for the anisotropic distribution of
the fracture orientation.
Fig. 4 shows the hydraulic conductivity anisotropy
ratio (maximum principal value / minimum principal
value). When D33 = 0 (Fig. 4 (a)), the hydraulicconductivity anisotropy ratio equals to 6 times. Thatis the hydraulic conductivity in tangential direction is
6 times larger than that in radial direction. With
increasing distance from the center of the tunnel, thehydraulic conductivity anisotropy reduced rapidly.
When D33 = 2 (Fig. 4 (b)), the anisotropy ratio at theroof of tunnel can be goes up to about 9 (Fig. 4 (a))
while the anisotropic ratio reduced to slightly larger
than 2 times at the side wall.

4. Results II- Flow field around tunnel
The flow field around the tunnel was evaluated.
The model covers an area of 175 m × 175 m. The
tunnel is located in the center of the square and theradius of the tunnel is 5 m. The govern equation ofsteady-state, two-dimensional groundwater flow canbe expressed as:
K xx h
 2 h K yy h
 2 h K xy h K xx 2  K yy 2 
x x
x
y y
y
x y

(11)
 K xyK xy hhh
 K xy
0
xy
y x
yx
2

2

where the K xx , K yy , K xy , K yx are the four
components of hydraulic conductivity tensor. The h
is the total head. The finite difference method (FDM)
is used to calculate the equal potentials, flow velocity,and flow rate flow into the tunnel. Dirichlet boundary
was applied on the outer boundaries of the square. A
zero total head was assumed on the inner boundary of
the tunnel wall.
The grid is square in shape with a size and 0.5 m ×
0.5 m. The stress tensor at each grid was calculatedusing the Kirsch’s solution (Eqs. (9) - (11)) and thehydraulic conductivity can be calculated by Oda’s
model (Eq. (1)).
Fig. 5 is the flow velocity around the tunnel which
the influence of shear dilatancy effect was neglected.
The red line segment represents the region withhigher flow velocity. The flow direction is in theradial rate around the tunnel is very consistent, andthe total flow rate into the tunnel wall is 0.032 m3 / s.

(a) roof of tunnel

Fig. 5 The flow velocity around the tunnel neglecting
the shear dilatancy effect (isotropic (D33 = 0))(b) side wall of tunnel
Fig. 4 The relationship between the anisotropy ratioof hydraulic conductivity near the tunnel wallFig. 6 (a) is the flow velocity around the tunnel
which the shear dilatancy effect was considered. The
red line segment represents the higher flow rate,
which is larger than the average velocity value in this

6th East Asia Forum on Radwaste Management Conference
November 27-29, 2017, Osaka, Japanregion. The total flow rate flowing into the tunnel
wall is 0.057 m3 / s, which is about two times of thatwhen the shear dilatancy was neglected.
Fig. 6 (b) shows the flow velocity around thetunnel for the rock mass with inherent anisotropy.
The fabric tensor D33 = 2. It can be observed that theflow directions are no longer toward to the center of
the tunnel for the presence of inherent anisotropy.
Generally, the flow velocity is larger at the roof thanthat at the side wall for the combine effect of normal
vector distribution and stress redistribution. The total
flow rate flowing into tunnel is 0.054 m3/s, which isalmost identical to the case for D33 = 0.

the tunnel excavation caused by the adjacent tunnel
wall hydraulic conductivity coefficient increased
significantly. This can be seen from the tunnel wallabout 0.5 times the radius of the hydraulic gradientslow down.
Fig. 7 (b) shows the equipotential lines (contour
lines of total head) around the tunnel with anisotropic
normal vector distribution (D33 = 2). It can beobserved that the equipotential lines are elliptical(long axis on horizontal direction). The long axis ispassing through the side wall where the hydraulicconductivity is larger than the one on the roof.

(a) isotropic (D33 = 0)(a) isotropic (D33 = 0)(b) anisotropic (D33 = 2)(b) anisotropic (D33 = 2)
Fig. 7 The equipotential lines around the tunnel

Fig. 6 The flow velocity around the tunnel consider
shear dilatancy effect

5. Conclusion

Fig. 7 (a) shows the equipotential lines (contour
lines of total head) around the tunnel with isotropic
normal vector distribution (D33= 0). The equipotential
line concentrates the boundary to the tunnel, because

Based on the Oda model, the stress-induced
anisotropy of hydraulic conductivity caused by thestress redistribution around the tunnel was evaluated.
Finite difference method was used to assess theground water flow around the tunnel. The flow

6th East Asia Forum on Radwaste Management Conference
November 27-29, 2017, Osaka, Japanvelocity, equipotential lines, and inflow rate were
calculated and the influence of stress-induced andinherent anisotropy of hydraulic conductivity on theflow field was evaluated. Important findings arelisted below:
It is found that the redistribution of tunnel
excavation stress will obviously change the tensor(principal value, main direction and anisotropy ratio)around the wall of the tunnel and affect the
groundwater flow behavior.
1. The hydraulic conductivity around the tunnel
increased 1 to 2 orders of magnitude before andafter the tunnel excavation, where only the
influence of stress redistribution was considered.
The influence zone is about one radius of thetunnel.
2. Neglecting shear dilatancy effect significantly
underestimates the apertures of fractures and thehydraulic conductivity of rock mass around thetunnel.
3. The hydraulic conductivity tensor (principal values,
principal directions and anisotropy ratio) of rockmass around the tunnel is strongly dominated by
the anisotropic stress around the tunnel andanisotropic distribution of normal vectors offractures. If the normal vector distribution isisotropic, the principal directions will be in theradial and tangential directions of the tunnel.
However, the principal directions will deviate fromthe radial and tangential directions if the
distribution of normal vector of fractures isanisotropic. When D33 = 2, the maximum principal
value of hydraulic conductivity is 9 times larger
than the minimum principal value at tunnel roof.
4. The inflow rate is not significantly influence by thestress redistribution due to tunnel excavation andthe normal vector distribution. However, the flow
field around the tunnel is significantly affected by
the anisotropic hydraulic conductivity, which is
influenced by the anisotropic stress field inducedby tunnel excavation and the anisotropic
distribution of the normal vector of fractures. If thenormal vector distribution is isotropic, the flow
vectors are toward to the center of the tunnel andthe flow velocity is symmetrical. For the case of
D33 = 2, the flow vectors are no longer along theradial directions and the flow velocities at tunnelroof are larger than the ones at side wall.

Mining Science & Geomechanics Abstracts, 20 (6),249-268, 1983.
[3] Barton, N., Bandis, S. and Bakhtar, K., Strength,
deformation and conductivity coupling of rockjoints, International Journal of Rock Mechanics
and Mining Science & Geomechanics Abstracts,
22 (3), 121-140, 1985.
[4] Knatani, K., Distribution of directional data andfabric tensors. International Journal of Engineering
Science, 22 (2), 149-164, 1984.
[5] Lee, C. H., Farmer, I., Fluid Flow in
Discontinuous Rocks.. Chapman and Hall London,
169 pp, 1993.
[6] Oda, M., Permeability tensor for discontinuousrock masses, Geotechnique, 35 (4), 483-495, 1985.
[7] Oda, M., An equivalent continuum model forcoupled stress and fluid flow analysis in jointed
rock masses, Water Resources Research, 22 (13),1845-1856, 1986.
[8] Oda, M., Yamabe, T., Ishizuka, Y., Kumasaka,
H., Tada H. and Kimura, K., Elastic stress andstrain in jointed rock masses by means of crack
tensor analysis, Rock Mechanics and RockEngineering, 26, 89-112, 1993.
[9] Po-Sung Lai, Inherent and stress-induced
anisotropy of hydraulic conductivity around a rocktunnel - equivalent continuum (master’s thesis
Chinese version), graduate institute of appliedgeology of National Central University, 2017.
[10] Yun-Chia Cheng, Inherent and stress-dependent
anisotropy of permeability for jointed rock masses(master’s thesis Chinese version), graduate
institute of applied geology of National Central
University, 2006.

Reference[1] Bandis, S., Lumsden, A. C. and Barton, N. R.,
Experimental studies of scale effects on the shear
behavior of rock joints, International Journal ofRock Mechanics and Mining Science &Geomechanics Abstracts, 18 (1), 1-21, 1981.
[2] Bandis, S., Lumsden, A. C. and Barton, N. R.,
Fundamentals of rock joint deformation,
International Journal of Rock Mechanics and6th East Asia Forum on Radwaste Management Conference
November 27-29, 2017, Osaka, Japan