# A Model of Anisotropic Property of Seepage and Stress for Jointed Rock Mass

Authors: Pei-tao Wang; Tian-hong Yang; Tao Xu; Qing-lei Yu; Hong-lei Liu
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

Joints often have important effects on seepage and elastic properties of jointed rock mass and therefore on the rock slope stability. In the present paper, a model for discrete jointed network is established using contact-free measurement technique and geometrical statistic method. A coupled mathematical model for characterizing anisotropic permeability tensor and stress tensor was presented and finally introduced to a finite element model. A case study of roadway stability at the Heishan Metal Mine in Hebei Province, China, was performed to investigate the influence of joints orientation on the anisotropic properties of seepage and elasticity of the surrounding rock mass around roadways in underground mining. In this work, the influence of the principal direction of the mechanical properties of the rock mass on associated stress field, seepage field, and damage zone of the surrounding rock mass was numerically studied. The numerical simulations indicate that flow velocity, water pressure, and stress field are greatly dependent on the principal direction of joint planes. It is found that the principal direction of joints is the most important factor controlling the failure mode of the surrounding rock mass around roadways.


---

## Body

## 1. Introduction

Underground mining has been considered a high-risk activity worldwide. Violent roof failure or rock burst induced by mining has always been a serious threat to the safety and efficiency of mines in China. Accurate and detailed characterization for rock masses can control stable excavation spans, support requirements, cavability, and subsidence characteristics, and thus influence the design of mining layouts and safety of mines. Rock mass is a geologic body composing of the discontinuities which have a critical influence on deformational behavior of blocky rock systems [1]. The mechanical behavior of this material depends principally on the state of intact rock whose mechanical properties could be determined by laboratory tests and existing discontinuities containing bedding planes, faults, joints, and other structural features. The distributions and strength of these discontinuities are both the key influencing factors for characterizing the discontinuous and anisotropic materials. Amadei [2] pointed out the importance of anisotropy of jointed rock mass and discussed the interaction existing between rock anisotropy and rock stress. Because of computational complexity and the difficulty of determining the necessary elastic constants, it is usual for only the simplest form of anisotropy, transverse isotropy, to be used in design analysis [3]. The key work is to study the principal direction of elasticity or permeability and then assume the jointed rock mass as transversely isotropic geomaterial.Extensive efforts have been made to investigate the mechanical response of transversely anisotropic rock material. Zhang and Sanderson [4] used the fractal dimension to describe the connectivity and compactness of fracture network and found that the deformability and the overall permeability of fractured rock masses increase greatly with increasing fracture density. By using the artificial transversely isotropic rock blocks, the mechanical properties with different dip angles were obtained by Tien and Tsao [5]. Brosch et al.  [6] evaluated the fabric-dependent anisotropy of a particular gneiss by studying the strength values and elastic parameters along different directions. Exadaktylos and Kaklis [7] presented the explicit representations of stresses and strains at any point of the anisotropic circular disc compressed diametrically. The stress concentration in uniaxial compression for a plane strain model was investigated, and a formulation of failure criteria for elastic-damage and elastic-plastic anisotropic geomaterials is formed by Exadaktylos [8]. A methodology to determine the equivalent elastic properties of fractured rock masses was established by explicit representations of stochastic fracture systems [9], and the conditions for the application of the equivalent continuum approach for representing mechanical behavior of the fractured rock masses were investigated. Such anisotropy has an effect on the interpretation of stress distributions. Amadei’s results indicate that an anisotropy ratio (damaged elastic modulus to intact elastic modulus) of between 1.14 and 1.33 will have a definite effect on the interpreted in situ state of stress [10]. Hakala et al.  [11] concluded that the anisotropy ratio of Finland rock specimens is about 1.4 and suggested to be taken into account in the interpretation of stress measurement results.Moreover, in rock engineering like mining and tunnel engineering, interactions between in situ stress and seepage pressure of groundwater have an important role. Groundwater under pressure in the joints defining rock blocks reduces the normal effective stress between the rock surfaces and therefore reduces the potential shear resistance which can be mobilized by friction. Since rock behavior may be determined by its geohydrological environment, it may be essential in some cases to maintain close control of groundwater conditions in the mine area. Therefore, accurate description of joints is an important topic for estimating and evaluating the deformability and seepage properties of rock masses.Many researchers have done a lot of studies on the seepage and anisotropic mechanical behaviors of the jointed rock mass. Using homogeneous samples like granite or basalt, Witherspoon et al.  [12] investigated and defined the permeability by fracture aperture in a closed fracture. Based on geometrical statistics, Oda [13] has studied and determined the crack tensor of moderately jointed granite by treating statistically the crack orientation data via a stereographic projection. Using Oda’s method, Sun and Zhao [14] determined the anisotropy in permeability using the fracture orientation and the in situ stress information from the field survey. For fractured rock system, attempts have been made by Jing et al. [15–18] and Bao et al. [19] to investigate the permeability. Using UDEC code, Jing et al. studied the permeability of discrete fracture network, such as the existence of REV in [15], relations between fracture length and aperture in [16] or stress effect on permeability in [17], or solute transport in [18]. Bao et al. [19] have discussed the mesh effect on effective permeability for a fractured system using the upscaled permeability field. It would help the workers to spend less computational effort and memory requirement to investigate effective permeability. The key for determining deformation modulus and hydraulic parameters is to study representative elementary volume (REV) and scale effects of fractured rock mass [9, 20–22]. However, the difficulty of testing jointed rock specimens, at scales sufficient to represent the equivalent continuum, indicated that it is necessary to postulate and verify methods of synthesizing rock mass properties from those of the constituent elements like intact rock and fractures.Although great progress has been made, it is difficult to study the anisotropic deformation of large-scale rock mass, ground water permeability change, and interaction between stress and seepage due to the geological complexity of discontinuities. Qiao et al.  [23] have pointed out that the rock mass which is not highly fractured and has only few sets of joint system usually behaves anisotropically. However, the mechanical properties directionality of highly jointed rock mass is usually ignored. In particular, since the limitation of mesh generation, the numerical method can hardly deal with the highly jointed specimens. The discontinuity like fault could be treated as specific boundary. An effective solution should be found to characterize the influence of joints on the rock mass.In this paper, based on equivalent continuum theory and theoretical analysis, a mathematical model for anisotropic property of seepage and elasticity of jointed rock mass is described. A permeability analysis code is developed to evaluate the anisotropic permeability for DFN model based on VC++ 6.0 in this paper. The DFN model could be analyzed to evaluate the REV size and anisotropic property of permeability, which would provide important evidence for the finite element model. The outline is as follows. First, 3D images involving detailed geometrical properties of rock mass, such as trace lengths, outcrop areas, joint orientations, and joint spacing, were captured using ShapMetriX3D system. According to the statistical parameters for each set of discontinuities, fracture network is generated using Monte Carlo method, and jointed rock samples in different sections could be captured. Next, permeability tensors and elasticity tensors of rock mass of the sample are calculated by discrete medium seepage method and geometrical damage theory, respectively. Finally, using the finite element method (FEM), an anisotropic mechanical model for rock mass is built. An engineering practice of roadway stability at the Heishan Metal Mine, in Hebei Iron and Steel Group Mining Company, is described. Then the stress and seepage fields surrounding the roadway were numerically simulated. On the basis of the modeled results, the influences of joint planes on stress, seepage, and damage zone were analyzed. It is expected throughout this study to gain an insight into the influences of discontinuities on the mechanical behavior of rock mass and offer some scientific evidence for the design of mining layouts or support requirements.


## 2. Generation of a Fracture System Model by ShapeMetriX3D

Traditional methods for rock mass structural parameters in mining engineering applications include scanline surveying [24] and drilling core method [25]. The process generally requires physical contact with rock mass exposure and therefore is hazardous. Additionally, taking manual measurements is time consuming and prone to errors due to sampling difficulties or instrument errors. ShapeMetriX3D is a tool for the geological and geotechnical data collection and assessment for rock masses [26]. The equipment is used for the metric acquisition of rock and terrain surfaces and for the contact-free measurement of geological/geotechnical parameters by metric 3D images. The images are captured using Nikon D80 camera with 22.3 megapixel. As shown in Figure 1, two images are acquired to reproduce the 3D rock mass face. Details information related to the measurements is reported in [27] by Gaich et al. The measurements could be taken at any required number and extent, even in regions that are not accessible. During measuring process, two digital images taken by a calibrated camera serve for a 3D reconstruction of the rock face geometry which is represented on the computer by a photorealistic spatial characterization as shown in Figure 1. From it, measurements are taken by marking visible rock mass features, such as spatial orientations of joint surfaces and traces, as well as areas, lengths, or positions. Finally, the probability statistical models of discontinuities are established. It is generally applied in the typical situations such as long rock faces at small height and rock slopes with complex geometries [28]. Almost any rock face can be reconstructed at its optimum resolution by using this equipment and its matching software. By using this system we can increase working safety, reduce mapping time, and improve data quality.Figure 1
Imaging principle of the 3D surface measurement.In order to accurately represent rock discontinuities, ShapeMetriX3D (3GSM) in this paper is used for the metric acquisition of rock mass exposure and for the contact-free measurement of geological parameters by metric 3D images. Stereoscopic photogrammetry deals with the measurement of three-dimensional information from two images showing the same object or surface but taken from two different angles, just as shown in Figure1.From the determined orientation between the two images and a pair of corresponding image pointsP1(u,v) and P2(u,v), imaging rays (colored in red) are reconstructed whose intersection leads to a 3D surface point P(X,Y,Z). By automatic identification of corresponding points within the image pair, the result of the acquisition is a metric 3D image that covers the geometry of the rock exposure. Once the image of a rock wall is ready, geometric measurements can be taken as shown in Figure 2. There are a total of three groups of discontinuities in this bench face.Figure 2
Geological mapping and geometric measurements of stereoscopic restructuring model.Figure3 is facilitated to show the 3-dimensional distribution of the trace along the section of Profile I-I from the stereoscopic model shown in Figure 2. The height of this rock face is approximately 5.0 m with a distance of 2.0 m perpendicular to the trace. The measured orientation in a hemispherical plot could also be captured as in Figure 4(a). Figure 4(b) shows the results of the distribution of joints.Figure 3
Length and sketch of Profile I-I by stereoscopic model.Stereographic plot and the distribution of total three sets of joints.
(a)
(b)The generated 3D rock face is about 5 m × 5 m with a high resolution enough to distinguish the fractures. According to the survey results, the geologic data like joint density, dip angle, trace length, and spacing can be acquired. Baghbanan and Jing [16] generated the DFN models whose orientations of fracture sets followed the Fisher distribution. In this paper, four different probability statistical models are used to generate the DFN models. The dip angle, dip direction, trace length, and spacing all follow one particular probability statistical model. Table 1 shows the basic information about the fracture system parameters. Type I of the probability statistical model in Table 1 stands for negative exponential distribution, Type II for normal distribution, Type III for logarithmic normal distribution, and Type IV for uniform distribution. Based on the probability models, the particular fracture network could be generated using Monte Carlo method [29, 30].Table 1
Characteristics of 3 sets of discontinuities for the slope exposure.
Set
Density  (m−1)
Fracture characteristics
Dip direction (°)
Dip angle (°)
Trace length (m)
Spacing (m)
Type
Mean value
Standard  deviation
Type
Mean value
Standard  deviation
Type
Mean value
Standard  deviation
Type
Mean value
Standard  deviation
1#
1.2
II
90.73
13.88
III
71.24
15.05
II
0.86
0.23
IV
0.83
0.75
2#
1.3
II
289.96
22.89
IV
59.85
8.4
II
0.89
0.28
IV
0.76
0.73
3#
3.3
III
189.05
14.55
III
55.29
13.35
II
0.85
0.25
IV
0.3
0.31


## 3. Constitutive Relation of Anisotropic Rock Mass

Due to the existence of joints and cracks, the mechanical properties (Young’s modulus, Poisson’s ratio, strength, etc.) of rock masses are generally heterogeneous and anisotropic. Three preconditions should be confirmed in this section: (i) anisotropy of rock mass is mainly caused by IV or V-class structure or rock masses containing a large number of discontinuities; (ii) rock masses according to (i) could be treated as homogeneous and anisotropic elastic material; (iii) seepage tensor and damage tensor of rock mass with multiset of joints could be captured by the scale of representative elementary volume (REV).


### 3.1. Stress Analysis

Elasticity represents the most common constitutive behavior of engineering materials, including many rocks, and it forms a useful basis for the description of more complex behavior. The most general statement of linear elastic constitutive behavior is a generalized form of Hooke’ Law, in which any strain component is a linear function of all the stress components; that is,(1)εij=[S]σij,
where [S] is the flexibility matrix, εij is the strain, and σij is the stress.Many underground excavation design analyses involving openings where the length to cross-section dimension ratio is high are facilitated considerably by the relative simplicity of the excavation geometry. In this section, the roadway is uniform cross section along the length and could be properly analyzed by assuming that the stress distribution is the same in all planes perpendicular to the long axis of the excavation (Figure5). Thus, this problem could be analyzed in terms of plane geometry.Figure 5
A diagrammatic sketch for underground excavation.The state of stress at any point can be defined in terms of the plane components of stress (σ11,σ22, and σ12) and the components (σ33,σ23, and σ31). In this research, the Z direction is assumed to be a principal axis and the antiplane shear stress components would vanish. The plane geometric problem could then be analyzed in terms of the plane components of stress since the σ33 component is frequently neglected. Equation (1), in this case, may be recast in the form
(2)[ε11ε22ε12]=[S][σ]=[1E1-ν312E3-(ν12E1+ν312E3)0-(ν12E1+ν312E3)1E1-ν312E30002(1+v12)E1]×[σ11σ22σ12],
where [S] is the flexibility matrix of the material under plane strain conditions. The inverse matrix [S]-1 (or [E]) could be expressed in the form(3)[σ11σ22σ12]=[E][ε]=[E12(E3-E1v312)Δ-E12(E3v21+E1v312)Δ0-E12(E3v21+E1v312)ΔE12(E1v312-E3)Δ000E12(1+v12)][ε11ε22ε12],where Δ is expressed as follows:
(4)Δ=-E1E3+E1E3v212+E12v312+E12v312+2E12v21v312.The moduliE1 and E3 and Poisson’s ratios v31 and v12 could be provided by uniaxial strength compression or tension in 1 (or 2) and 3 directions.The mechanical tests including laboratory and in situ tests for rock masses in large scale can hardly capture the elastic properties directly. Hoek-Brown criterion [31–33] could calculate the mechanical properties of weak rocks masses by introducing the Geological Strength Index (GSI). Nevertheless, the anisotropic properties cannot be captured using this criterion, and thus the method for analyzing anisotropy of jointed rock mass needs a further study.In this paper, the original joint damage in rock mass is considered as macro damage field. In elastic damage mechanics, the elastic modulus of the jointed material may degrade and the Young’s modulus of the damaged element is defined as follows [34]:
(5)E=E0(1-D),
where D represents the damage variable and E and E0 are the elastic moduli of the damaged and intact rock samples, respectively. In this equation, all parameters are scalar.With the geometric information of the fracture sample, the damage tensor [35] could be defined as
(6)Dij=lV∑k=1Nα(k)(n(k)⊗n(k)),(i,j=1,2,3),
where N is the number of joints, l is the minimum spacing between joints, V is the volume of rock mass, n(k) is the normal vector of the kth joint, and a(k) is the trace length of the kth joint (for 2 dimensions).According to the principal of energy equivalence [36], the flexibility matrix for the jointed rock sample can be obtained as
(7)Sij′=(1-Di)-1Sij(1-Dj)-1,
where Sij is the flexibility matrix for intact rock; Di and Dj are the principal damage values in i and j directions, respectively. For plane strain geometric problem, the constitutive relation, where the coordinates and principal damage have the same direction, could be expressed as follows:(8)[σ11σ22σ12]=[E0(v0-1)(D1-1)22v02+v0-1-E0v0(D1-1)(D2-1)2v02+v0-10-E0v0(D1-1)(D2-1)2v02+v0-1E0(v0-1)(D2-1)22v02+v0-1000E0(D1-1)(D2-1)2(1+v0)][ε11ε22ε12],where E0 is the Young’s modulus of intact rock and v0 is Poisson’s ratio for intact rock.The parameters could be easily captured by laboratory tests. Equation (7) gives the principal damage values D1 and D2. Based on geometrical damage mechanics, all elements in this matrix could be obtained by the method mentioned previously, and the anisotropic constitutive relation of jointed rock sample could finally be confirmed.


### 3.2. Seepage Analysis

The seepage parameters of rock mass are quantized form of permeability and also are the basis to solve seepage field of equivalent continuous medium. Based on the attribute of fractured rock mass and by taking engineering design into account, fractured rock mass is often considered to be anisotropic continuous medium. In the fracture network shown in Figure6, a total number of N cross points or water heads and M line elements are contained. Parameters can be obtained with the model such as water head, the related line elements, equivalent mechanical fissure width, seepage coefficient, and so on. The corresponding coordinates of each point could be acquired. For a fluid flow analysis based on the law of mass conservation, the fluid equations on a certain water head take the form [37]
(9)(∑j=1N′qj)i+Qi=0,(i=1,2,…,N),
where qj is the quantity of flow from line element j to water head i,N′ is the total number of line elements intersect at i, and Qi is the fluid source term. In the joint network, each line element would be assigned a length lj and fissure width bj to investigate the permeability.Figure 6
Fracture network schematic diagram in seepage area.MN and M′N′ are the constant head boundaries with the hydraulic head of hb1 and hb2 respectively; MM′ and NN′ are the impervious boundaries with flow V of zero.According to the hydraulic theory [38], for a single joint seepage, the flow quantity qj of line element j can be expressed as
(10)qj=ρbj312μ·Δhjlj,
where Δhj is the hydraulic gradient, μ is the coefficient of flow viscosity, ρ is the density of water, and bj is the fissure width of joint. On the basis of (9) and (10), the governing equations can be represented as (11) for seepage in fracture network
(11)(∑j=1N′ρbj312μ·Δhjlj)i+Qi=0.Based on the discrete fracture network method [39], an equivalent continuum model for seepage has been established [28, 40, 41]. The hydraulic conductivity can be acquired based on Darcy’s Law on the basis of water quantity in the network (see (12)).Boundary conditions in the model shown in Figure6 are as follows:(a)
MN and M′N′ as the constant head boundary;(b)
MM′ and NN′ as the impervious boundary with flow V of 0.Then the hydraulic conductivity can be defined as(12)K=Δq·M′MΔH·MN,
where Δq is the total quantity of water in the model region (m2·s−1), ΔH is the water pressure difference between inflow and outflow boundaries (m), K is the equivalent hydraulic conductivity coefficient in the MM′ direction (m·s−1), and MN and M′Mare the side lengths of the region (m).Based on Biot equations [42], the steady flow model is given by
(13)Kij∇2p=0,
where Kij is the hydraulic conductivity and p is the hydraulic pressure. For plane problems, the dominating equation of seepage flow is as follows in (14). The direction of joint planes is considered to be the principal direction of hydraulic conductivity
(14)Kij=[K11K12K21K22].


### 3.3. Coupling Mechanism of Seepage and Stress

#### 3.3.1. Seepage Inducted by Stress

The coupling action between seepage and stress makes the failure mechanism of rock complex. The investigations on this problem have pervasive theoretical meaning and practical value. The principal directions associated with the symmetric crack tensor are coaxial with those of the permeability tensor. The first invariant of the crack tensor is proportional to the mean permeability, while the deviatoric part is related to the anisotropic permeability [13]. Generally, the change of stress which is perpendicular to the joints plane is the main factor leading to the increase or decrease of the ground water permeability. In the numerical model, seepage is coupled to stress describing the permeability change induced by the change of the stress field. The coupling function can be described as follows as given by Louis [43]:
(15)Kf=K0e-βσ,
where Kf is the current groundwater hydraulic conductivity, K0 is the initial hydraulic conductivity, σ is the stress perpendicular to the joints plane, and β is the coupling parameter (stress sensitive factor to be measured by experiment) that reflects the influence of stress. The larger β is, the greater the range of stress induced the permeability [44].


#### 3.3.2. Stress Inducted by Seepage

On the basis of generalized Terzaghi’s effective stress principle [45], the stress equilibrium equation could be expressed as follows for the water-bearing jointed specimen:
(16)σij=Eijklεkl-αijPδij,
where σij is the total stress tensor, Eijkl is the elastic tensor of the solid phase, εkl is the strain tensor, αij is a positive constant which is equal to 1 when individual grains are much more incompressible than the grain skeleton, P is the hydraulic pressure, and δij is the Kronecker delta function.


### 3.4. Coordinate Transformation

Generally, planes of joints are inclined at an angle to the major principal stress direction as shown in Figure7. In establishing these equations, the X,Y, and 1, 3 axes are taken to have the same Z(2) axis, and the angel θ is measured from the x to the 1 axis.Figure 7
The relationship between planes of joints and the principal stress direction.Considering the hydraulic pressure needs no coordinate transformation, and (16) will be expressed as
(17)[σ11σ22τ12]=[Tσ]-1[E][Tε][ε11ε22γ12]-αij[PP0]δij,
where [Tσ]-1 is the reverse matrix for stress coordinates transformation and [Tε] is the strain coordinates transformation matrix, and they can be expressed as follows:
(18)[Tσ]-1=[cos2θsin2θ-2cosθsinθsin2θcos2θ2sinθcosθsinθcosθ-sinθcosθcos2θ-sin2θ],[Tε]=[cos2θsin2θcosθsinθsin2θcos2θ-sinθcosθ-2sinθcosθ2sinθcosθcos2θ-sin2θ].Similarly, the coordinate transformation form can be expressed as(19)K11=K11′+K22′2+K11′-K22′2cos2θ,K12=K21=-K11′-K22′2sin2θ,K22=K11′+K22′2-K11′-K22′2cos2θ,
where K11′ is the hydraulic conductivity coefficient along the distribution direction of joint planes, K22′ is perpendicular to the distribution direction, and θ is measured from the coordinate system x to the optimal direction of permeability.


## 3.1. Stress Analysis

Elasticity represents the most common constitutive behavior of engineering materials, including many rocks, and it forms a useful basis for the description of more complex behavior. The most general statement of linear elastic constitutive behavior is a generalized form of Hooke’ Law, in which any strain component is a linear function of all the stress components; that is,(1)εij=[S]σij,
where [S] is the flexibility matrix, εij is the strain, and σij is the stress.Many underground excavation design analyses involving openings where the length to cross-section dimension ratio is high are facilitated considerably by the relative simplicity of the excavation geometry. In this section, the roadway is uniform cross section along the length and could be properly analyzed by assuming that the stress distribution is the same in all planes perpendicular to the long axis of the excavation (Figure5). Thus, this problem could be analyzed in terms of plane geometry.Figure 5
A diagrammatic sketch for underground excavation.The state of stress at any point can be defined in terms of the plane components of stress (σ11,σ22, and σ12) and the components (σ33,σ23, and σ31). In this research, the Z direction is assumed to be a principal axis and the antiplane shear stress components would vanish. The plane geometric problem could then be analyzed in terms of the plane components of stress since the σ33 component is frequently neglected. Equation (1), in this case, may be recast in the form
(2)[ε11ε22ε12]=[S][σ]=[1E1-ν312E3-(ν12E1+ν312E3)0-(ν12E1+ν312E3)1E1-ν312E30002(1+v12)E1]×[σ11σ22σ12],
where [S] is the flexibility matrix of the material under plane strain conditions. The inverse matrix [S]-1 (or [E]) could be expressed in the form(3)[σ11σ22σ12]=[E][ε]=[E12(E3-E1v312)Δ-E12(E3v21+E1v312)Δ0-E12(E3v21+E1v312)ΔE12(E1v312-E3)Δ000E12(1+v12)][ε11ε22ε12],where Δ is expressed as follows:
(4)Δ=-E1E3+E1E3v212+E12v312+E12v312+2E12v21v312.The moduliE1 and E3 and Poisson’s ratios v31 and v12 could be provided by uniaxial strength compression or tension in 1 (or 2) and 3 directions.The mechanical tests including laboratory and in situ tests for rock masses in large scale can hardly capture the elastic properties directly. Hoek-Brown criterion [31–33] could calculate the mechanical properties of weak rocks masses by introducing the Geological Strength Index (GSI). Nevertheless, the anisotropic properties cannot be captured using this criterion, and thus the method for analyzing anisotropy of jointed rock mass needs a further study.In this paper, the original joint damage in rock mass is considered as macro damage field. In elastic damage mechanics, the elastic modulus of the jointed material may degrade and the Young’s modulus of the damaged element is defined as follows [34]:
(5)E=E0(1-D),
where D represents the damage variable and E and E0 are the elastic moduli of the damaged and intact rock samples, respectively. In this equation, all parameters are scalar.With the geometric information of the fracture sample, the damage tensor [35] could be defined as
(6)Dij=lV∑k=1Nα(k)(n(k)⊗n(k)),(i,j=1,2,3),
where N is the number of joints, l is the minimum spacing between joints, V is the volume of rock mass, n(k) is the normal vector of the kth joint, and a(k) is the trace length of the kth joint (for 2 dimensions).According to the principal of energy equivalence [36], the flexibility matrix for the jointed rock sample can be obtained as
(7)Sij′=(1-Di)-1Sij(1-Dj)-1,
where Sij is the flexibility matrix for intact rock; Di and Dj are the principal damage values in i and j directions, respectively. For plane strain geometric problem, the constitutive relation, where the coordinates and principal damage have the same direction, could be expressed as follows:(8)[σ11σ22σ12]=[E0(v0-1)(D1-1)22v02+v0-1-E0v0(D1-1)(D2-1)2v02+v0-10-E0v0(D1-1)(D2-1)2v02+v0-1E0(v0-1)(D2-1)22v02+v0-1000E0(D1-1)(D2-1)2(1+v0)][ε11ε22ε12],where E0 is the Young’s modulus of intact rock and v0 is Poisson’s ratio for intact rock.The parameters could be easily captured by laboratory tests. Equation (7) gives the principal damage values D1 and D2. Based on geometrical damage mechanics, all elements in this matrix could be obtained by the method mentioned previously, and the anisotropic constitutive relation of jointed rock sample could finally be confirmed.


## 3.2. Seepage Analysis

The seepage parameters of rock mass are quantized form of permeability and also are the basis to solve seepage field of equivalent continuous medium. Based on the attribute of fractured rock mass and by taking engineering design into account, fractured rock mass is often considered to be anisotropic continuous medium. In the fracture network shown in Figure6, a total number of N cross points or water heads and M line elements are contained. Parameters can be obtained with the model such as water head, the related line elements, equivalent mechanical fissure width, seepage coefficient, and so on. The corresponding coordinates of each point could be acquired. For a fluid flow analysis based on the law of mass conservation, the fluid equations on a certain water head take the form [37]
(9)(∑j=1N′qj)i+Qi=0,(i=1,2,…,N),
where qj is the quantity of flow from line element j to water head i,N′ is the total number of line elements intersect at i, and Qi is the fluid source term. In the joint network, each line element would be assigned a length lj and fissure width bj to investigate the permeability.Figure 6
Fracture network schematic diagram in seepage area.MN and M′N′ are the constant head boundaries with the hydraulic head of hb1 and hb2 respectively; MM′ and NN′ are the impervious boundaries with flow V of zero.According to the hydraulic theory [38], for a single joint seepage, the flow quantity qj of line element j can be expressed as
(10)qj=ρbj312μ·Δhjlj,
where Δhj is the hydraulic gradient, μ is the coefficient of flow viscosity, ρ is the density of water, and bj is the fissure width of joint. On the basis of (9) and (10), the governing equations can be represented as (11) for seepage in fracture network
(11)(∑j=1N′ρbj312μ·Δhjlj)i+Qi=0.Based on the discrete fracture network method [39], an equivalent continuum model for seepage has been established [28, 40, 41]. The hydraulic conductivity can be acquired based on Darcy’s Law on the basis of water quantity in the network (see (12)).Boundary conditions in the model shown in Figure6 are as follows:(a)
MN and M′N′ as the constant head boundary;(b)
MM′ and NN′ as the impervious boundary with flow V of 0.Then the hydraulic conductivity can be defined as(12)K=Δq·M′MΔH·MN,
where Δq is the total quantity of water in the model region (m2·s−1), ΔH is the water pressure difference between inflow and outflow boundaries (m), K is the equivalent hydraulic conductivity coefficient in the MM′ direction (m·s−1), and MN and M′Mare the side lengths of the region (m).Based on Biot equations [42], the steady flow model is given by
(13)Kij∇2p=0,
where Kij is the hydraulic conductivity and p is the hydraulic pressure. For plane problems, the dominating equation of seepage flow is as follows in (14). The direction of joint planes is considered to be the principal direction of hydraulic conductivity
(14)Kij=[K11K12K21K22].


## 3.3. Coupling Mechanism of Seepage and Stress

### 3.3.1. Seepage Inducted by Stress

The coupling action between seepage and stress makes the failure mechanism of rock complex. The investigations on this problem have pervasive theoretical meaning and practical value. The principal directions associated with the symmetric crack tensor are coaxial with those of the permeability tensor. The first invariant of the crack tensor is proportional to the mean permeability, while the deviatoric part is related to the anisotropic permeability [13]. Generally, the change of stress which is perpendicular to the joints plane is the main factor leading to the increase or decrease of the ground water permeability. In the numerical model, seepage is coupled to stress describing the permeability change induced by the change of the stress field. The coupling function can be described as follows as given by Louis [43]:
(15)Kf=K0e-βσ,
where Kf is the current groundwater hydraulic conductivity, K0 is the initial hydraulic conductivity, σ is the stress perpendicular to the joints plane, and β is the coupling parameter (stress sensitive factor to be measured by experiment) that reflects the influence of stress. The larger β is, the greater the range of stress induced the permeability [44].


### 3.3.2. Stress Inducted by Seepage

On the basis of generalized Terzaghi’s effective stress principle [45], the stress equilibrium equation could be expressed as follows for the water-bearing jointed specimen:
(16)σij=Eijklεkl-αijPδij,
where σij is the total stress tensor, Eijkl is the elastic tensor of the solid phase, εkl is the strain tensor, αij is a positive constant which is equal to 1 when individual grains are much more incompressible than the grain skeleton, P is the hydraulic pressure, and δij is the Kronecker delta function.


## 3.3.1. Seepage Inducted by Stress

The coupling action between seepage and stress makes the failure mechanism of rock complex. The investigations on this problem have pervasive theoretical meaning and practical value. The principal directions associated with the symmetric crack tensor are coaxial with those of the permeability tensor. The first invariant of the crack tensor is proportional to the mean permeability, while the deviatoric part is related to the anisotropic permeability [13]. Generally, the change of stress which is perpendicular to the joints plane is the main factor leading to the increase or decrease of the ground water permeability. In the numerical model, seepage is coupled to stress describing the permeability change induced by the change of the stress field. The coupling function can be described as follows as given by Louis [43]:
(15)Kf=K0e-βσ,
where Kf is the current groundwater hydraulic conductivity, K0 is the initial hydraulic conductivity, σ is the stress perpendicular to the joints plane, and β is the coupling parameter (stress sensitive factor to be measured by experiment) that reflects the influence of stress. The larger β is, the greater the range of stress induced the permeability [44].


## 3.3.2. Stress Inducted by Seepage

On the basis of generalized Terzaghi’s effective stress principle [45], the stress equilibrium equation could be expressed as follows for the water-bearing jointed specimen:
(16)σij=Eijklεkl-αijPδij,
where σij is the total stress tensor, Eijkl is the elastic tensor of the solid phase, εkl is the strain tensor, αij is a positive constant which is equal to 1 when individual grains are much more incompressible than the grain skeleton, P is the hydraulic pressure, and δij is the Kronecker delta function.


## 3.4. Coordinate Transformation

Generally, planes of joints are inclined at an angle to the major principal stress direction as shown in Figure7. In establishing these equations, the X,Y, and 1, 3 axes are taken to have the same Z(2) axis, and the angel θ is measured from the x to the 1 axis.Figure 7
The relationship between planes of joints and the principal stress direction.Considering the hydraulic pressure needs no coordinate transformation, and (16) will be expressed as
(17)[σ11σ22τ12]=[Tσ]-1[E][Tε][ε11ε22γ12]-αij[PP0]δij,
where [Tσ]-1 is the reverse matrix for stress coordinates transformation and [Tε] is the strain coordinates transformation matrix, and they can be expressed as follows:
(18)[Tσ]-1=[cos2θsin2θ-2cosθsinθsin2θcos2θ2sinθcosθsinθcosθ-sinθcosθcos2θ-sin2θ],[Tε]=[cos2θsin2θcosθsinθsin2θcos2θ-sinθcosθ-2sinθcosθ2sinθcosθcos2θ-sin2θ].Similarly, the coordinate transformation form can be expressed as(19)K11=K11′+K22′2+K11′-K22′2cos2θ,K12=K21=-K11′-K22′2sin2θ,K22=K11′+K22′2-K11′-K22′2cos2θ,
where K11′ is the hydraulic conductivity coefficient along the distribution direction of joint planes, K22′ is perpendicular to the distribution direction, and θ is measured from the coordinate system x to the optimal direction of permeability.


## 4. A Case Study

### 4.1. Description of the Area under Study

The model described above is applied to evaluate the seepage field and stress field of jointed rock roadway in the Heishan Metal Mine (Figure8). The mine is located in Chengde city, Hebei province, in northern China. Heishan Metal Mine has transferred from open pit mining to underground mining since 2009. The deformation and stability of the roadways for mining the hanging-wall ore become the key technique issue. The elevation of the crest of the slope is 920 m. The roadway in the research is in 674 m level of the high northern slope. The rock mass is mainly composed of anorthosite and norite in the northern slope area. The anisotropic properties of seepage and stress field will be discussed in detail.Figure 8
The rock mass exposure of northern slope in Heishan open pit mine.


### 4.2. Capture of Joint Network

Depending on the 3D contact-free measuring system, discontinuities on the exposure could be easily captured. The collection process for the geology information of northern slope in Heishan Metal Mine has been discussed previously. A 3D fracture network of rock mass could be obtained using the Monte Carlo method by analyzing the statistical parameters (Figure9(a)). Finally, the fracture network is generated and expanded to the roadway, and the fracture network perpendicular to the center line of roadway, could be easily captured (Figure 9(b)).Joint network based on the statistical data using Monte Carlo method. (a) 3D joint network in rock mass; (b) profile with 10 m edge length perpendicular to the roadway.
(a)
(b)


### 4.3. Permeability Investigation

According to the symmetry of geometry, the method for solving the hydraulic conductivity coefficient in different directions of fracture network is presented to save computational time [41].Rotated every 15° clockwise and exerted certain water pressure on the boundary shown in Figure10, the hydraulic conductivity coefficient of verticaldirection K[90], 75° direction K[75], 60° direction K[60], 45° direction K[45], 30° direction K[30], and 15° direction K[15] could be acquired, respectively. Then the permeability tensors by fracture parameters can be determined. By enlarging the geometric size of the fracture network, the permeability scale effect can be investigated finally.Figure 10
Generated fracture network by Monte Carlo method (side length of internal squares is 1 m, 4 m, 7 m, 8 m, and 9 m, resp.).Based on the algorithm in Section3.2, the size effect of rock mass in different sizes of statistics window is studied. The region of fracture network is enlarged from 1 m to 9 m with different step lengths until the equivalent parameters in different directions achieve a stable value. With a fixed center of this region, the equivalent permeability is calculated rotating every 15° clockwise as shown in Figure 11.Figure 11
Fracture network in the study area with edge length of 7 m.The variation of hydraulic conductivity coefficients with the increase of direction angles and sample size is depicted in Figure12. The main direction angle is approximated to 15°. It can be seen that the anisotropy of seepage property is apparent with the distribution of joints. Seven hydraulic conductivity coefficients of different sizes (1 m, 3 m, 4 m, 5 m, 6 m, 8 m, and 9 m) are chosen to be compared with those of the sample whose size is 7 m. The coefficient deviation between one particular size and 7 m would be found and presented in Figure 13. Values of permeability decrease with the increase of sample size and tend to stabilize when the sample size comes to 7 m.Figure 12
Permeability tensor of the rock mass (unit: m/s).Figure 13
Deviation of the permeability values to that of   7 m sample under different directions of samples in different sizes.The principal values of permeability in maximum and minimum in a rock sample with joint plane angle being 15° are listed in Table2. The principal permeability values decrease from 4.23 (×10−6m/s) to 2.66 (×10−6m/s) as the sample size increases from 1 m to 9 m along the main direction. According to the results discussed previously, the REV is 7 m × 7 m for this fracture sample, and the related hydraulic conductivity is 2.77 (×10−6m/s) for maximum and 0.87 (×10−6m/s) for minimum. The ratio for the maximum to minimum hydraulic conductivity value is 3.18.Table 2
Principal valuesK[15] of fracture hydraulic conductivity at different sizes.
Sample size (m)
1
4
7
8
9
Principal value of  hydraulic conductivity (×10−6 m/s)
Maximum
4.23
3.28
2.77
2.69
2.66
Minimum
0.80
0.97
0.87
0.88
0.85
Ratio
5.29
3.38
3.18
3.06
3.13


### 4.4. Damage Tensor

Similar to the algorithm in Section3.1, the size effect of sample damage in different sizes of statistics window is studied. The region of fracture network is enlarged from 3 m to 12 m with different step lengths until the damage values in different directions stabilize. With a fixed center of this region, the principal damage values are calculated rotating every 15° clockwise. Figure 14 shows the sample damage in different sizes and directions. Similar to the deviation analysis of permeability, the damage values also fluctuatingly tend to stability according to the deviation of the damage tensor in different sizes and directions shown in Figure 15. In this study, the size of representative elementary volume is also 7 m in length, and the initial damage tensor of REV of jointed rock mass could also be obtained. The principal damage values D1 and D2 for fracture sample are 0.17 for minimum and 0.50 for maximum, and the principal direction angle θ of damage tensor is approximately 105° which is perpendicular to the direction of principal permeability. The related rock parameters, including undamaged Young’s modulus E0, Poisson’s ratio v0, and uniaxial tensile strength, et al, are listed in Table 3. It should be noted that the elasticity and strength of rock could be determined in the laboratory tests. Finally, the constitutive relation can be found according to (8). Equations (17) and (19) can be introduced into the FEM simulations.Table 3
Parameters used in the model to validate the model in simulating the anisotropic properties of stress and seepage.
Material parameters
Values
Undamaged Young’s modulusE0 (GPa)
65
Undamaged Poisson’s ratiov0
0.23
Density (kg/m−3)
2700
Principal direction angle of joint planeθ (°)
15
Coupling parameterβ (Pa−1)
0.5Figure 14
Damage tensor (×10−1).Figure 15
Deviation of the damage tensor to 7 m sample under different directions of samples in different sizes.


### 4.5. Geometry and Boundary Conditions

A numerical model of Heishan Metal Mine is established in order to simulate the mechanism of stress and seepage in jointed rock mass, taking into account the anisotropic property, as shown in Figure16. Table 3 shows the basic mechanical properties. Parameters listed in Table 4 are the hydraulic conductivity coefficients Kij, damage tensor Dii and Young’s modulus Eij, and shear modulus Gij of the fracture sample along joint plane direction used in the finite element code. The model contains two roadways with a 3.5 m × 3 m three-centered arch section within a 50 m × 50 m domain. The bottom boundary of the domain is fixed in all directions, and the left and the right boundaries are fixed in the horizontal direction. In this regard, a pressure (σs) of 5.97 MPa is applied on the top boundary of the model to represent the 221 m deep overburden strata. Under the steady-state groundwater flow condition, a hydrostatic pressure, pw, of 2.21 MPa is applied upon top boundary. No-flow conditions are imposed on the three boundaries of the rectangular domain. The initial water pressure on the roadways boundaries is 0 MPa. All the governing equations described previously are implemented into COMSOL Multiphysics, a powerful PDE-based multiphysics modeling environment. The model is assumed to be in a state of plane strain (with no change in elastic strain in the vertical direction) and static mechanical equilibrium.Table 4
The hydraulic conductivity coefficientsKij, damage tensor Dii and Young’s modulus Eij, and shear modulus Gij of the fracture sample along joint plane direction.
Subscriptsij
K
i
j  (×10-6 m/s)
D
i
i
E
i
j (GPa)
G
i
j (GPa)
11
2.77
0.50
18.83
—
22
0.87
0.17
51.91
—
12
—
—
9.34
10.97Plane strain roadway model: (a) boundary conditions; (b) the finite element mesh.
(a)
(b)


### 4.6. Results and Discussion

#### 4.6.1. Stress Distribution

Adverse performance of the rock mass in the postexcavation stress field may be caused by either failure of the anisotropic medium or slip on the weakness planes [3]. The elastic stress distribution around the roadways directly influences the deformation of rock and thus determines the design process.Figure17 shows the contour of first principal stress coupled with the seepage process. The orientation and magnitude of maximum principal stress controlled the distribution of the stress concentration in the heterogeneous media. The simulation result shows that the principal stress concentration zones appear mainly in rock surrounding the roadways. There exist maximum stress concentration areas in the arch foot and floor. Measures should be taken to control the deformation and assure the construction safety.Figure 17
Contours of the first principal stress forθ=15° case.To characterize the response of the stress to the hydraulic mechanics, a comparison of two scenarios is also presented as shown in Figure18. The first principal stress in the stress-seepage coupled model along the horizontal section A-A′ where y=27 is compared with a decoupled models. The result shows that the first principal stress increases when the seepage process is considered. Figure 19 shows a plot of normal stress in joint plane direction and hydraulic conductivity along the horizontal section A-A′ where y=27 as is shown in Figure 16. The coupled and decoupled model could be analyzed using Comsol Multiphysics code. In the coupled case, the governing equations for solid and fluid phase are solved in weakly coupled sense. For the sake of convenient contrast, the normal stress distribution as well as hydraulic conductivity is plotted when no seepage-coupled process is considered. The result shows that the normal stress will be underestimated when no coupling of seepage process is included. Moreover, the compressive stress leads to the decrease of permeability, and the hydraulic conductivity increases when seepage process is included.Figure 18
Contrast of first principal stress between coupled and decoupled model.Figure 19
Observed changes of normal stress in joint plane direction and hydraulic conductivity between coupled and decoupled models.


#### 4.6.2. Seepage Distribution

Flow velocity with the angle of joint plane being 15° is shown in Figure20. All the process is not considered time dependent. Therefore, the velocity change is induced by the permeability variation caused by the compressive stress upon the joint plane. In this case study, water flows along the principal direction of joint plane where the permeability is largest in the model. An asymmetric seepage field is observed, and maximum of flow velocity is distributed on right top of the roadways roof. Moreover, the seepage pressure is also asymmetrically distributed as shown in Figure 21.Figure 20
Darcy’s velocity magnitude (m/s).Figure 21
The seepage pressure distribution (Pa).


#### 4.6.3. Damage Zone

For underground mining, two principal engineering properties of the joint planes should be considered. They are low tensile strength in the direction perpendicular to the joint plane and the relatively low shear strength of the surfaces. The anisotropic strength parameters as tensile and compressive strength parameters for jointed rock mass are discussed by Chen et al. [46], Claesson and Bohloli [47], Nasseri et al. [48], Gonzaga et al. [49], and Cho et al. [50]. As discussed previously, the fluid pressure and velocity are sensitive to the joint plane angles θ. In this section, the Hoffman anisotropic strength criterion is used to assess the damage zone in this numerical model as shown in
(20)σ12XtXc-σ1σ2XtXc+σ22YtYc+Xc-XtXtXcσ1+Yc-YtYtYcσ2+τ122S2=1.The propertiesXt and Yt represent the tensile strength along joint plane direction and perpendicular to joint plane direction, respectively. Xc and Yc represent the compressive strength along joint plane direction and perpendicular to joint plane direction, respectively. S is the shear strength of the material along the joint direction. σ1 represents the normal stress along the principle direction of elasticity and σ2 represents the normal stress perpendicular to the principle direction of elasticity. τ12 represents the shear stress.The shear strength of the joints can be described by the simple Coulomb law(21)S=c+σ2tanϕ,
where c is cohesive strength and ϕ is the effective angle of friction of the joint surfaces.It should be noted that the mechanical parameters be acquired from laboratory or in situ tests. However, it is difficult to directly employ the strength parameters for jointed rock mass due to the inaccessibility of the tests for huge rock mass. Based on the in situ and laboratory tests of Heishan Metal Mine, the tensile, compressive, and shear strength parameters are listed in Table5.Table 5
Mechanical properties of rock mass in the direction of joint plane.
Tensile strength (MPa)
X
t
16.0
Y
t
4.0
Compressive strength (MPa)
X
c
120
Y
c
96
Shear strength (MPa)
c
1.5
ϕ
50The direction of joint plane is 15°, and thus the tensile strength in the direction perpendicular to joint plane is relatively low, compared with that in other directions. Figure22 shows the damage zone in this case study and the failure area mainly distributes in the direction perpendicular to the weakness plane, which gives an illustration of the roadway failure mode in tabular orebodies. Moreover, the failure of covered rock mass and the rock pillar in and between the two roadways does not influence each other in this case study, and thus the choice of roadway’s interval is proper from the aspect of the mechanics analysis.Figure 22
Damage zone of the simulating model (the direction of damage zone is approximately perpendicular to that of the joint planes whose angle is 15°).


### 4.7. Further Discussion

These simulations were performed to develop an understanding of the mechanics of joints and influence on stress and seepage fields and to gauge the ability of the proposed transversely anisotropic model to capture the response of jointed rock mass. For this purpose, a total of six scenarios with the joint plane anglesθ ranging from 0° to 150° with an interval of 30° are simulated in order to examine the effect of joint plane directions; see Figure 7 for the definition of θ. And the anisotropic properties of seepage field and damage zones are examined in these simulations.


#### 4.7.1. Seepage Field

Figure23 presents the fluid pressure distribution and flow field arrows with different joint plane directions. It can be seen that the maximum pressure, which is located on the top boundary, equals the initial fluid pressure (2.21 MPa). The fluid pressure distribution on top of the roadways differs with the increase of directions.Contour of the fluid pressure and flow velocity vectors in different joint plane directions.
(a)
(b)
(c)
(d)
(e)
(f)When the joint plane is parallel or perpendicular to the floors of roadways, the fluid pressure distributions and flow arrows are all axially symmetric, which agrees with expectations. The fluid pressure in the 30° or 60° case is distributed unsymmetrically as shown in Figures23(b) and 23(c).The flow velocity iny direction along the horizontal section A-A′ where y=27 m as shown in Figure 16(a) curves are plotted in Figure 24 for θ=0°, 30°, 60°, and 90°, respectively. In all cases the absolute value of velocity increases approximately exponentially above the roadways. As shown in Figure 24, the flow velocity in scenario where θ=0° is in the lowest level. The reason is that the joint plane in this scenario is horizontally distributed and the compressive stress leads to the decrease of permeability. When angle of joint plane increases to 30°, flow velocities above both roadways increase and Roadway I has a more higher velocity than Roadway II. The scenario where θ=60° has similar performance. However, when the joint plane is vertically distributed, flow velocity above Roadway I is lower than that where θ=60°. The flow velocity in the case where joints are vertically or horizontally distributed is symmetric.Figure 24
The flow velocity iny direction along the horizontal section A-A′ for different numerical models, at θ = 0°, 30°, 60°, and 90°.


#### 4.7.2. Damage Zone

The effect of joint plane angle on the damage zone is illustrated by using the Hoffman anisotropic strength criterion as shown in (23). The damage zones in different angles of joint planes are shown in Figure25 and could be an index to visualize the potential failure mode of the roadway.The damaged zone under different principal elastic directions (red arrows represent the flow velocity vector).
(a)
(b)
(c)
(d)
(e)
(f)It is clear from Figures25(a) to 25(f) that an increase of the joint plane angle has significantly influenced the shape and size of damage zone. When the joints are horizontally distributed (θ=0°), the tensile strength in the direction perpendicular to joint planes is much lower and thus the damage zone mainly concentrates within roof and bottom of roadways. This failure mainly manifests as roof falling and floor heave. Similarly, when the angle increases to 30°, the damage zones mainly concentrate in the rock mass of left roof and right floor within the roadways which is also in the direction perpendicular to joint planes. Compared with the result in Figure 25(b), the direction of damage zone with θ=60° shown in Figure 25(c) rotated significantly and the area also increases. When the joints are vertically distributed (θ=90°), the principal direction of damage zone is nearly horizontal. Lateral rock mass surrounding the roadways stabilizes with the increase of joint plane angle in this study, and the pillar between two roadways is also stable. The scenarios where θ=120° or θ=150° have the similar response to the joint plane direction with that of 60° or 30°. Qualitatively, the simulation results are in good agreement with the results by Jia et al. [51] and Huang et al. [52]. The model in this study could to some extent be effectively used to analyze the anisotropic property for jointed rock mass.


## 4.1. Description of the Area under Study

The model described above is applied to evaluate the seepage field and stress field of jointed rock roadway in the Heishan Metal Mine (Figure8). The mine is located in Chengde city, Hebei province, in northern China. Heishan Metal Mine has transferred from open pit mining to underground mining since 2009. The deformation and stability of the roadways for mining the hanging-wall ore become the key technique issue. The elevation of the crest of the slope is 920 m. The roadway in the research is in 674 m level of the high northern slope. The rock mass is mainly composed of anorthosite and norite in the northern slope area. The anisotropic properties of seepage and stress field will be discussed in detail.Figure 8
The rock mass exposure of northern slope in Heishan open pit mine.


## 4.2. Capture of Joint Network

Depending on the 3D contact-free measuring system, discontinuities on the exposure could be easily captured. The collection process for the geology information of northern slope in Heishan Metal Mine has been discussed previously. A 3D fracture network of rock mass could be obtained using the Monte Carlo method by analyzing the statistical parameters (Figure9(a)). Finally, the fracture network is generated and expanded to the roadway, and the fracture network perpendicular to the center line of roadway, could be easily captured (Figure 9(b)).Joint network based on the statistical data using Monte Carlo method. (a) 3D joint network in rock mass; (b) profile with 10 m edge length perpendicular to the roadway.
(a)
(b)


## 4.3. Permeability Investigation

According to the symmetry of geometry, the method for solving the hydraulic conductivity coefficient in different directions of fracture network is presented to save computational time [41].Rotated every 15° clockwise and exerted certain water pressure on the boundary shown in Figure10, the hydraulic conductivity coefficient of verticaldirection K[90], 75° direction K[75], 60° direction K[60], 45° direction K[45], 30° direction K[30], and 15° direction K[15] could be acquired, respectively. Then the permeability tensors by fracture parameters can be determined. By enlarging the geometric size of the fracture network, the permeability scale effect can be investigated finally.Figure 10
Generated fracture network by Monte Carlo method (side length of internal squares is 1 m, 4 m, 7 m, 8 m, and 9 m, resp.).Based on the algorithm in Section3.2, the size effect of rock mass in different sizes of statistics window is studied. The region of fracture network is enlarged from 1 m to 9 m with different step lengths until the equivalent parameters in different directions achieve a stable value. With a fixed center of this region, the equivalent permeability is calculated rotating every 15° clockwise as shown in Figure 11.Figure 11
Fracture network in the study area with edge length of 7 m.The variation of hydraulic conductivity coefficients with the increase of direction angles and sample size is depicted in Figure12. The main direction angle is approximated to 15°. It can be seen that the anisotropy of seepage property is apparent with the distribution of joints. Seven hydraulic conductivity coefficients of different sizes (1 m, 3 m, 4 m, 5 m, 6 m, 8 m, and 9 m) are chosen to be compared with those of the sample whose size is 7 m. The coefficient deviation between one particular size and 7 m would be found and presented in Figure 13. Values of permeability decrease with the increase of sample size and tend to stabilize when the sample size comes to 7 m.Figure 12
Permeability tensor of the rock mass (unit: m/s).Figure 13
Deviation of the permeability values to that of   7 m sample under different directions of samples in different sizes.The principal values of permeability in maximum and minimum in a rock sample with joint plane angle being 15° are listed in Table2. The principal permeability values decrease from 4.23 (×10−6m/s) to 2.66 (×10−6m/s) as the sample size increases from 1 m to 9 m along the main direction. According to the results discussed previously, the REV is 7 m × 7 m for this fracture sample, and the related hydraulic conductivity is 2.77 (×10−6m/s) for maximum and 0.87 (×10−6m/s) for minimum. The ratio for the maximum to minimum hydraulic conductivity value is 3.18.Table 2
Principal valuesK[15] of fracture hydraulic conductivity at different sizes.
Sample size (m)
1
4
7
8
9
Principal value of  hydraulic conductivity (×10−6 m/s)
Maximum
4.23
3.28
2.77
2.69
2.66
Minimum
0.80
0.97
0.87
0.88
0.85
Ratio
5.29
3.38
3.18
3.06
3.13


## 4.4. Damage Tensor

Similar to the algorithm in Section3.1, the size effect of sample damage in different sizes of statistics window is studied. The region of fracture network is enlarged from 3 m to 12 m with different step lengths until the damage values in different directions stabilize. With a fixed center of this region, the principal damage values are calculated rotating every 15° clockwise. Figure 14 shows the sample damage in different sizes and directions. Similar to the deviation analysis of permeability, the damage values also fluctuatingly tend to stability according to the deviation of the damage tensor in different sizes and directions shown in Figure 15. In this study, the size of representative elementary volume is also 7 m in length, and the initial damage tensor of REV of jointed rock mass could also be obtained. The principal damage values D1 and D2 for fracture sample are 0.17 for minimum and 0.50 for maximum, and the principal direction angle θ of damage tensor is approximately 105° which is perpendicular to the direction of principal permeability. The related rock parameters, including undamaged Young’s modulus E0, Poisson’s ratio v0, and uniaxial tensile strength, et al, are listed in Table 3. It should be noted that the elasticity and strength of rock could be determined in the laboratory tests. Finally, the constitutive relation can be found according to (8). Equations (17) and (19) can be introduced into the FEM simulations.Table 3
Parameters used in the model to validate the model in simulating the anisotropic properties of stress and seepage.
Material parameters
Values
Undamaged Young’s modulusE0 (GPa)
65
Undamaged Poisson’s ratiov0
0.23
Density (kg/m−3)
2700
Principal direction angle of joint planeθ (°)
15
Coupling parameterβ (Pa−1)
0.5Figure 14
Damage tensor (×10−1).Figure 15
Deviation of the damage tensor to 7 m sample under different directions of samples in different sizes.


## 4.5. Geometry and Boundary Conditions

A numerical model of Heishan Metal Mine is established in order to simulate the mechanism of stress and seepage in jointed rock mass, taking into account the anisotropic property, as shown in Figure16. Table 3 shows the basic mechanical properties. Parameters listed in Table 4 are the hydraulic conductivity coefficients Kij, damage tensor Dii and Young’s modulus Eij, and shear modulus Gij of the fracture sample along joint plane direction used in the finite element code. The model contains two roadways with a 3.5 m × 3 m three-centered arch section within a 50 m × 50 m domain. The bottom boundary of the domain is fixed in all directions, and the left and the right boundaries are fixed in the horizontal direction. In this regard, a pressure (σs) of 5.97 MPa is applied on the top boundary of the model to represent the 221 m deep overburden strata. Under the steady-state groundwater flow condition, a hydrostatic pressure, pw, of 2.21 MPa is applied upon top boundary. No-flow conditions are imposed on the three boundaries of the rectangular domain. The initial water pressure on the roadways boundaries is 0 MPa. All the governing equations described previously are implemented into COMSOL Multiphysics, a powerful PDE-based multiphysics modeling environment. The model is assumed to be in a state of plane strain (with no change in elastic strain in the vertical direction) and static mechanical equilibrium.Table 4
The hydraulic conductivity coefficientsKij, damage tensor Dii and Young’s modulus Eij, and shear modulus Gij of the fracture sample along joint plane direction.
Subscriptsij
K
i
j  (×10-6 m/s)
D
i
i
E
i
j (GPa)
G
i
j (GPa)
11
2.77
0.50
18.83
—
22
0.87
0.17
51.91
—
12
—
—
9.34
10.97Plane strain roadway model: (a) boundary conditions; (b) the finite element mesh.
(a)
(b)


## 4.6. Results and Discussion

### 4.6.1. Stress Distribution

Adverse performance of the rock mass in the postexcavation stress field may be caused by either failure of the anisotropic medium or slip on the weakness planes [3]. The elastic stress distribution around the roadways directly influences the deformation of rock and thus determines the design process.Figure17 shows the contour of first principal stress coupled with the seepage process. The orientation and magnitude of maximum principal stress controlled the distribution of the stress concentration in the heterogeneous media. The simulation result shows that the principal stress concentration zones appear mainly in rock surrounding the roadways. There exist maximum stress concentration areas in the arch foot and floor. Measures should be taken to control the deformation and assure the construction safety.Figure 17
Contours of the first principal stress forθ=15° case.To characterize the response of the stress to the hydraulic mechanics, a comparison of two scenarios is also presented as shown in Figure18. The first principal stress in the stress-seepage coupled model along the horizontal section A-A′ where y=27 is compared with a decoupled models. The result shows that the first principal stress increases when the seepage process is considered. Figure 19 shows a plot of normal stress in joint plane direction and hydraulic conductivity along the horizontal section A-A′ where y=27 as is shown in Figure 16. The coupled and decoupled model could be analyzed using Comsol Multiphysics code. In the coupled case, the governing equations for solid and fluid phase are solved in weakly coupled sense. For the sake of convenient contrast, the normal stress distribution as well as hydraulic conductivity is plotted when no seepage-coupled process is considered. The result shows that the normal stress will be underestimated when no coupling of seepage process is included. Moreover, the compressive stress leads to the decrease of permeability, and the hydraulic conductivity increases when seepage process is included.Figure 18
Contrast of first principal stress between coupled and decoupled model.Figure 19
Observed changes of normal stress in joint plane direction and hydraulic conductivity between coupled and decoupled models.


### 4.6.2. Seepage Distribution

Flow velocity with the angle of joint plane being 15° is shown in Figure20. All the process is not considered time dependent. Therefore, the velocity change is induced by the permeability variation caused by the compressive stress upon the joint plane. In this case study, water flows along the principal direction of joint plane where the permeability is largest in the model. An asymmetric seepage field is observed, and maximum of flow velocity is distributed on right top of the roadways roof. Moreover, the seepage pressure is also asymmetrically distributed as shown in Figure 21.Figure 20
Darcy’s velocity magnitude (m/s).Figure 21
The seepage pressure distribution (Pa).


### 4.6.3. Damage Zone

For underground mining, two principal engineering properties of the joint planes should be considered. They are low tensile strength in the direction perpendicular to the joint plane and the relatively low shear strength of the surfaces. The anisotropic strength parameters as tensile and compressive strength parameters for jointed rock mass are discussed by Chen et al. [46], Claesson and Bohloli [47], Nasseri et al. [48], Gonzaga et al. [49], and Cho et al. [50]. As discussed previously, the fluid pressure and velocity are sensitive to the joint plane angles θ. In this section, the Hoffman anisotropic strength criterion is used to assess the damage zone in this numerical model as shown in
(20)σ12XtXc-σ1σ2XtXc+σ22YtYc+Xc-XtXtXcσ1+Yc-YtYtYcσ2+τ122S2=1.The propertiesXt and Yt represent the tensile strength along joint plane direction and perpendicular to joint plane direction, respectively. Xc and Yc represent the compressive strength along joint plane direction and perpendicular to joint plane direction, respectively. S is the shear strength of the material along the joint direction. σ1 represents the normal stress along the principle direction of elasticity and σ2 represents the normal stress perpendicular to the principle direction of elasticity. τ12 represents the shear stress.The shear strength of the joints can be described by the simple Coulomb law(21)S=c+σ2tanϕ,
where c is cohesive strength and ϕ is the effective angle of friction of the joint surfaces.It should be noted that the mechanical parameters be acquired from laboratory or in situ tests. However, it is difficult to directly employ the strength parameters for jointed rock mass due to the inaccessibility of the tests for huge rock mass. Based on the in situ and laboratory tests of Heishan Metal Mine, the tensile, compressive, and shear strength parameters are listed in Table5.Table 5
Mechanical properties of rock mass in the direction of joint plane.
Tensile strength (MPa)
X
t
16.0
Y
t
4.0
Compressive strength (MPa)
X
c
120
Y
c
96
Shear strength (MPa)
c
1.5
ϕ
50The direction of joint plane is 15°, and thus the tensile strength in the direction perpendicular to joint plane is relatively low, compared with that in other directions. Figure22 shows the damage zone in this case study and the failure area mainly distributes in the direction perpendicular to the weakness plane, which gives an illustration of the roadway failure mode in tabular orebodies. Moreover, the failure of covered rock mass and the rock pillar in and between the two roadways does not influence each other in this case study, and thus the choice of roadway’s interval is proper from the aspect of the mechanics analysis.Figure 22
Damage zone of the simulating model (the direction of damage zone is approximately perpendicular to that of the joint planes whose angle is 15°).


## 4.6.1. Stress Distribution

Adverse performance of the rock mass in the postexcavation stress field may be caused by either failure of the anisotropic medium or slip on the weakness planes [3]. The elastic stress distribution around the roadways directly influences the deformation of rock and thus determines the design process.Figure17 shows the contour of first principal stress coupled with the seepage process. The orientation and magnitude of maximum principal stress controlled the distribution of the stress concentration in the heterogeneous media. The simulation result shows that the principal stress concentration zones appear mainly in rock surrounding the roadways. There exist maximum stress concentration areas in the arch foot and floor. Measures should be taken to control the deformation and assure the construction safety.Figure 17
Contours of the first principal stress forθ=15° case.To characterize the response of the stress to the hydraulic mechanics, a comparison of two scenarios is also presented as shown in Figure18. The first principal stress in the stress-seepage coupled model along the horizontal section A-A′ where y=27 is compared with a decoupled models. The result shows that the first principal stress increases when the seepage process is considered. Figure 19 shows a plot of normal stress in joint plane direction and hydraulic conductivity along the horizontal section A-A′ where y=27 as is shown in Figure 16. The coupled and decoupled model could be analyzed using Comsol Multiphysics code. In the coupled case, the governing equations for solid and fluid phase are solved in weakly coupled sense. For the sake of convenient contrast, the normal stress distribution as well as hydraulic conductivity is plotted when no seepage-coupled process is considered. The result shows that the normal stress will be underestimated when no coupling of seepage process is included. Moreover, the compressive stress leads to the decrease of permeability, and the hydraulic conductivity increases when seepage process is included.Figure 18
Contrast of first principal stress between coupled and decoupled model.Figure 19
Observed changes of normal stress in joint plane direction and hydraulic conductivity between coupled and decoupled models.


## 4.6.2. Seepage Distribution

Flow velocity with the angle of joint plane being 15° is shown in Figure20. All the process is not considered time dependent. Therefore, the velocity change is induced by the permeability variation caused by the compressive stress upon the joint plane. In this case study, water flows along the principal direction of joint plane where the permeability is largest in the model. An asymmetric seepage field is observed, and maximum of flow velocity is distributed on right top of the roadways roof. Moreover, the seepage pressure is also asymmetrically distributed as shown in Figure 21.Figure 20
Darcy’s velocity magnitude (m/s).Figure 21
The seepage pressure distribution (Pa).


## 4.6.3. Damage Zone

For underground mining, two principal engineering properties of the joint planes should be considered. They are low tensile strength in the direction perpendicular to the joint plane and the relatively low shear strength of the surfaces. The anisotropic strength parameters as tensile and compressive strength parameters for jointed rock mass are discussed by Chen et al. [46], Claesson and Bohloli [47], Nasseri et al. [48], Gonzaga et al. [49], and Cho et al. [50]. As discussed previously, the fluid pressure and velocity are sensitive to the joint plane angles θ. In this section, the Hoffman anisotropic strength criterion is used to assess the damage zone in this numerical model as shown in
(20)σ12XtXc-σ1σ2XtXc+σ22YtYc+Xc-XtXtXcσ1+Yc-YtYtYcσ2+τ122S2=1.The propertiesXt and Yt represent the tensile strength along joint plane direction and perpendicular to joint plane direction, respectively. Xc and Yc represent the compressive strength along joint plane direction and perpendicular to joint plane direction, respectively. S is the shear strength of the material along the joint direction. σ1 represents the normal stress along the principle direction of elasticity and σ2 represents the normal stress perpendicular to the principle direction of elasticity. τ12 represents the shear stress.The shear strength of the joints can be described by the simple Coulomb law(21)S=c+σ2tanϕ,
where c is cohesive strength and ϕ is the effective angle of friction of the joint surfaces.It should be noted that the mechanical parameters be acquired from laboratory or in situ tests. However, it is difficult to directly employ the strength parameters for jointed rock mass due to the inaccessibility of the tests for huge rock mass. Based on the in situ and laboratory tests of Heishan Metal Mine, the tensile, compressive, and shear strength parameters are listed in Table5.Table 5
Mechanical properties of rock mass in the direction of joint plane.
Tensile strength (MPa)
X
t
16.0
Y
t
4.0
Compressive strength (MPa)
X
c
120
Y
c
96
Shear strength (MPa)
c
1.5
ϕ
50The direction of joint plane is 15°, and thus the tensile strength in the direction perpendicular to joint plane is relatively low, compared with that in other directions. Figure22 shows the damage zone in this case study and the failure area mainly distributes in the direction perpendicular to the weakness plane, which gives an illustration of the roadway failure mode in tabular orebodies. Moreover, the failure of covered rock mass and the rock pillar in and between the two roadways does not influence each other in this case study, and thus the choice of roadway’s interval is proper from the aspect of the mechanics analysis.Figure 22
Damage zone of the simulating model (the direction of damage zone is approximately perpendicular to that of the joint planes whose angle is 15°).


## 4.7. Further Discussion

These simulations were performed to develop an understanding of the mechanics of joints and influence on stress and seepage fields and to gauge the ability of the proposed transversely anisotropic model to capture the response of jointed rock mass. For this purpose, a total of six scenarios with the joint plane anglesθ ranging from 0° to 150° with an interval of 30° are simulated in order to examine the effect of joint plane directions; see Figure 7 for the definition of θ. And the anisotropic properties of seepage field and damage zones are examined in these simulations.


### 4.7.1. Seepage Field

Figure23 presents the fluid pressure distribution and flow field arrows with different joint plane directions. It can be seen that the maximum pressure, which is located on the top boundary, equals the initial fluid pressure (2.21 MPa). The fluid pressure distribution on top of the roadways differs with the increase of directions.Contour of the fluid pressure and flow velocity vectors in different joint plane directions.
(a)
(b)
(c)
(d)
(e)
(f)When the joint plane is parallel or perpendicular to the floors of roadways, the fluid pressure distributions and flow arrows are all axially symmetric, which agrees with expectations. The fluid pressure in the 30° or 60° case is distributed unsymmetrically as shown in Figures23(b) and 23(c).The flow velocity iny direction along the horizontal section A-A′ where y=27 m as shown in Figure 16(a) curves are plotted in Figure 24 for θ=0°, 30°, 60°, and 90°, respectively. In all cases the absolute value of velocity increases approximately exponentially above the roadways. As shown in Figure 24, the flow velocity in scenario where θ=0° is in the lowest level. The reason is that the joint plane in this scenario is horizontally distributed and the compressive stress leads to the decrease of permeability. When angle of joint plane increases to 30°, flow velocities above both roadways increase and Roadway I has a more higher velocity than Roadway II. The scenario where θ=60° has similar performance. However, when the joint plane is vertically distributed, flow velocity above Roadway I is lower than that where θ=60°. The flow velocity in the case where joints are vertically or horizontally distributed is symmetric.Figure 24
The flow velocity iny direction along the horizontal section A-A′ for different numerical models, at θ = 0°, 30°, 60°, and 90°.


### 4.7.2. Damage Zone

The effect of joint plane angle on the damage zone is illustrated by using the Hoffman anisotropic strength criterion as shown in (23). The damage zones in different angles of joint planes are shown in Figure25 and could be an index to visualize the potential failure mode of the roadway.The damaged zone under different principal elastic directions (red arrows represent the flow velocity vector).
(a)
(b)
(c)
(d)
(e)
(f)It is clear from Figures25(a) to 25(f) that an increase of the joint plane angle has significantly influenced the shape and size of damage zone. When the joints are horizontally distributed (θ=0°), the tensile strength in the direction perpendicular to joint planes is much lower and thus the damage zone mainly concentrates within roof and bottom of roadways. This failure mainly manifests as roof falling and floor heave. Similarly, when the angle increases to 30°, the damage zones mainly concentrate in the rock mass of left roof and right floor within the roadways which is also in the direction perpendicular to joint planes. Compared with the result in Figure 25(b), the direction of damage zone with θ=60° shown in Figure 25(c) rotated significantly and the area also increases. When the joints are vertically distributed (θ=90°), the principal direction of damage zone is nearly horizontal. Lateral rock mass surrounding the roadways stabilizes with the increase of joint plane angle in this study, and the pillar between two roadways is also stable. The scenarios where θ=120° or θ=150° have the similar response to the joint plane direction with that of 60° or 30°. Qualitatively, the simulation results are in good agreement with the results by Jia et al. [51] and Huang et al. [52]. The model in this study could to some extent be effectively used to analyze the anisotropic property for jointed rock mass.


## 4.7.1. Seepage Field

Figure23 presents the fluid pressure distribution and flow field arrows with different joint plane directions. It can be seen that the maximum pressure, which is located on the top boundary, equals the initial fluid pressure (2.21 MPa). The fluid pressure distribution on top of the roadways differs with the increase of directions.Contour of the fluid pressure and flow velocity vectors in different joint plane directions.
(a)
(b)
(c)
(d)
(e)
(f)When the joint plane is parallel or perpendicular to the floors of roadways, the fluid pressure distributions and flow arrows are all axially symmetric, which agrees with expectations. The fluid pressure in the 30° or 60° case is distributed unsymmetrically as shown in Figures23(b) and 23(c).The flow velocity iny direction along the horizontal section A-A′ where y=27 m as shown in Figure 16(a) curves are plotted in Figure 24 for θ=0°, 30°, 60°, and 90°, respectively. In all cases the absolute value of velocity increases approximately exponentially above the roadways. As shown in Figure 24, the flow velocity in scenario where θ=0° is in the lowest level. The reason is that the joint plane in this scenario is horizontally distributed and the compressive stress leads to the decrease of permeability. When angle of joint plane increases to 30°, flow velocities above both roadways increase and Roadway I has a more higher velocity than Roadway II. The scenario where θ=60° has similar performance. However, when the joint plane is vertically distributed, flow velocity above Roadway I is lower than that where θ=60°. The flow velocity in the case where joints are vertically or horizontally distributed is symmetric.Figure 24
The flow velocity iny direction along the horizontal section A-A′ for different numerical models, at θ = 0°, 30°, 60°, and 90°.


## 4.7.2. Damage Zone

The effect of joint plane angle on the damage zone is illustrated by using the Hoffman anisotropic strength criterion as shown in (23). The damage zones in different angles of joint planes are shown in Figure25 and could be an index to visualize the potential failure mode of the roadway.The damaged zone under different principal elastic directions (red arrows represent the flow velocity vector).
(a)
(b)
(c)
(d)
(e)
(f)It is clear from Figures25(a) to 25(f) that an increase of the joint plane angle has significantly influenced the shape and size of damage zone. When the joints are horizontally distributed (θ=0°), the tensile strength in the direction perpendicular to joint planes is much lower and thus the damage zone mainly concentrates within roof and bottom of roadways. This failure mainly manifests as roof falling and floor heave. Similarly, when the angle increases to 30°, the damage zones mainly concentrate in the rock mass of left roof and right floor within the roadways which is also in the direction perpendicular to joint planes. Compared with the result in Figure 25(b), the direction of damage zone with θ=60° shown in Figure 25(c) rotated significantly and the area also increases. When the joints are vertically distributed (θ=90°), the principal direction of damage zone is nearly horizontal. Lateral rock mass surrounding the roadways stabilizes with the increase of joint plane angle in this study, and the pillar between two roadways is also stable. The scenarios where θ=120° or θ=150° have the similar response to the joint plane direction with that of 60° or 30°. Qualitatively, the simulation results are in good agreement with the results by Jia et al. [51] and Huang et al. [52]. The model in this study could to some extent be effectively used to analyze the anisotropic property for jointed rock mass.


## 5. Conclusion

The main purpose is to introduce the anisotropic model of seepage and stress and apply the model to the jointed rock mass. The coupled finite element analysis was performed, and the effects of joint planes direction on stress field, flow velocity, and Darcy’s velocity were verified. A more reliable model for the stability analysis of rock engineering and risk evaluation of water inrush was provided. Based on the results of a series of numerical simulations under different scenarios, the following conclusions are drawn.(1)
A linkage of digital information of fractures and mechanical analysis is realized. Relations between digital images involving detailed geometrical properties of rock mass and the quantitative determination of hydraulic parameters as well as elastic properties were realized. The results show that the scale for both damage tensor and permeability tensor of the rock mass’s REV in northern slope of Heishan Metal Mine is 7 m× 7 m.(2)
We examine the model for seepage-stress coupled analysis on anisotropic properties. The numerical simulations in this study have indicated that the existence of joint planes greatly affects the seepage properties and stress field. In anisotropic rock, water flows mainly along the joint planes, and water pressure is asymmetrically distributed where the angle of joint plane is 15° in the northern slope of Heishan Metal Mine.(3)
The influence of fractures cannot be neglected in stability analysis of rock mass. The numerical results visualized the damage zones in different directions of joint planes. The direction of damage zones was found perpendicular to that of joint planes which agrees well with the field observations and theoretical analysis.The work reported in this paper is an initial effort on the influence of joint orientation on the anisotropic property of seepage and stress in jointed rock masses. A model for more complex jointed rock mass remains to be quantified.



---

*Source: 420536-2013-07-10.xml*