A new estimation method and an anisotropy index for the deformation modulus of jointed rock masses

Bohu Zhang ,Junyan Mu ,Jun Zhng ,Qing Lv ,Jianhui Dng

a State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation,Southwest Petroleum University,Chengdu,610500,China

b School of Geoscience and Technology,Southwest Petroleum University,Chengdu,610500,China

c Department of Civil Engineering,Zhejiang University,Hangzhou,310058,China

d Center for Balance Architecture,Zhejiang University,Hangzhou,310007,China

e State Key Laboratory of Hydraulics and Mountain River Engineering,College of Water Resource and Hydropower,Sichuan University,Chengdu,610065,China

Keywords:Deformation modulus Analytical method Anisotropy index Jointed rock masses Mechanical behavior Discrete fracture network(DFN)

ABSTRACT The deformation modulus of a rock mass is an important parameter to describe its mechanical behavior.In this study,an analytical method is developed to determine the deformation modulus of jointed rock masses,which considers the mechanical properties of intact rocks and joints based on the superposition principle.Due to incorporating the variations in the orientations and sizes of joint sets,the proposed method is applicable to the rock mass with persistent and parallel joints as well as that with nonpersistent and nonparallel joints.In addition,an anisotropy index AI dm for the deformation modulus is de fined to quantitatively describe the anisotropy of rock masses.The range of AI dm is from 0 to 1,and the more anisotropic the rock mass is,the larger the value of AI dm will be.To evaluate the proposed method,20 groups of numerical experiments are conducted with the universal distinct element code(UDEC).For each experimental group,the deformation modulus in 24 directions are obtained by UDEC(numerical value)and the proposed method(predicted value),and then the mean error rates are calculated.Note that the mean error rate is the mean value of the error rates of the deformation modulus in 24 directions,where for each direction,the error rate is equal to the ratio of numerical value minus predicted value to the numerical value.The results show that(i)for different experimental groups,the mean error rates vary between 5.06%and 22.03%;(ii)the error rates for the discrete fracture networks(DFNs)with two sets of joints are at the same level as those with one set of joints;and(iii)therefore,the proposed method for estimating the deformation modulus of jointed rock masses is valid.

1.Introduction

The mechanical behavior of jointed rock masses should be appropriately assessed when designing foundations,slopes,underground openings,and anchoring systems(Zheng et al.,2014a;Asadizadeh and Hossaini,2016;Bahaaddini and Hosseinpour,2019;Hoek and Brown,2019;Lai et al.,2020;Mo et al.,2020).The deformation modulus is an important parameter of rock masses.Numerous numerical methods for studying the stress and displacement of rock masses are performed by using the deformation modulus as an input parameter(Palmström and Singh,2001;Panthee et al.,2016;Alnedawi et al.,2019).In addition,the deformation modulus is a parameter that best represents the prefailure mechanical behavior of a rock mass(Kayabasi et al.,2003).

In nature,a rock mass is a discontinuous medium with natural fissures,fractures,joints,bedding planes,and faults(Zheng et al.,2014b,2017;Mohammad et al.,2015),which makes it much more dif ficult to determine the deformation modulus of a rock mass than that of an intact rock.Therefore,many scholars have conducted in-depth studies on how to estimate the deformation modulus of a rock mass.These estimation methods can be categorized into two main types:in situ test methods and indirect test methods.

The in situ test methods for determining the deformation modulus of a rock mass mainly include plate loading tests(PLTs),plate jacking tests(PJTs),dilatometer tests,and Goodman jack tests.Among these methods,the most common in situ tests are PLTs and PJTs.Both PLTs and PJTs are often referred to as plate bearing tests(Kavur et al.,2014).PLTs measure displacements only on the rock surface;conversely,PJTs,which use a central instrumentation hole and multipoint borehole extensometers(MPBXs),measure displacements at different depths within rock masses(Vrkljan et al.,1995).In dilatometer tests,the volume of the tested area is too small and the tensile stresses are involved in the borehole;hence,the deformation modulus values calculated through dilatometer tests are usually 2-3 times lower than those calculated through PLTs and PJTs(Rocha,1974;Bieniawski,1978).The deformation modulus values calculated through Goodman jack tests are also lower than those calculated through PLTs and PJTs due to the deformation of the loading platens in hard rock.Bieniawski(1989)noted that few previous projects featured a suf ficient number of different tests to allow a meaningful comparison of in situ test data types.In situ test results may be highly variable depending on different test methods.Hence,it is not reliable to accept the results from any in situ method,i.e.two or more methods are desirable to be used to crosscheck the test results.Therefore,the in situ tests for determining the deformation modulus of a rock mass can be considered to be expensive,time-consuming,and dif ficult to conduct.

The indirect test methods mainly include empirical methods and analytical methods.In contrast to in situ measurements,due to the advantages of simplicity,rapidity and low cost,empirical methods have gradually become highly popular.Empirical methods mainly use existing field data to establish a relationship between deformation modulus and rock classi fications,such as the rock mass rating(RMR),tunneling quality index(Q),rock quality designation(RQD),and geological strength index(GSI).Bieniawski(1978)first proposed an empirical equation of the deformation modulus that used RMR as the parameter.However,this empirical equation can be applied only if the value of RMR is greater than 50.Sera fim and Pereira(1983)suggested another equation that can be used when RMR is less than or equal to 50.Following these works,many scholars have studied the deformation modulus and developed some empirical formulae(Nicholson and Bieniawski,1990;Hoek and Brown,1997;Zhang and Einstein,2004;Hoek and Diederichs 2006;Sonmez et al.,2006;Chun et al.,2009;Martins and Miranda 2012;Khabbazi et al.,2013;Alemdag et al.,2015).In recent years,in addition to using traditional methods to obtain empirical formulae,fuzzy mathematics,arti ficial intelligence,multiple and polynomial regression analyses,and neural networks have also been used to predict the deformation modulus(Ali Bashari et al.,2011;Alemdag et al.,2015;Mohammad Rezaei et al.,2015;Asadizadeh and Hossaini,2016;Bahaaddini and Hosseinpour,2019;Adeyemi and Wang,2019).For the aforementioned indirect estimation methods,a range of differences in deformation moduli of the rock massEmobtained from different equations were observed for the same rock mass(Hoek and Diederichs,2006;Ajalloeian and Mohammadi,2014;Panthee et al.,2016).Moreover,neither in situ tests nor empirical equations can re flect the anisotropy of the deformation modulus of a rock mass.It is of great signi ficance to evaluate the anisotropy of the deformation modulus of a rock mass for engineering design and construction.

The deformation modulus of a jointed rock mass can also be estimated by analytical methods.Gerrard(1982)studied the deformation modulus of a rock mass with several sets of joints based on orthorhombic layer theory.Unfortunately,the joints in the same set were assumed to be planar and approximately equally spaced.Fossum(1985)deduced a formula for estimating the deformation modulus of a rock mass with random persistent joints based on a simple geometric averaging method.However,the stiffness of the joints in the study was set to a constant value.Amadei and Savage(1993)obtained a general solution of the deformability of rock masses with regular joints and noted that the deformation modulus of jointed rock masses is anisotropic.Zhang(2010)regarded a heavily jointed rock mass as isotropic and developed a method to estimate the deformation modulus of heavily jointed rock masses.However,Zhang(2010)assumed that all normal stiffness and shear stiffness values were equal,which may deviate from the actual conditions.Li(2001)developed a method to estimate the deformation modulus of jointed rock masses in three-dimensional(3D)space and presented it in a hemispherical projection.This method re flected the anisotropy of the deformation modulus of jointed rock masses,and different values of stiffness were considered for different sets of joints.However,the joints in the study were assumed to be in finite,which limited the applicability of the method.Generally,previous research on analytical methods of estimating the deformation modulus of jointed rock masses does not fully re flect the randomness of joints,and the assumptions of joints are idealized.

In summary,although determining the deformation modulus of rock masses is a classic issue,there are still some defects in existing solutions.This study aims to propose a new way to estimate the deformation modulus of jointed rock masses.This new method has the following advantages:(i)the deformation modulus of rock masses with non-persistent and nonparallel joint networks can be quickly estimated;(ii)the approach can be applied to jointed rock masses at the engineering scale,which is different from the relevant research achievements obtained by laboratory tests and numerical methods for indoor test-scale rock samples with regular joints(Singh et al.,2002;Wong and Einstein,2009;Zhou et al.,2018);and(iii)the deformation modulus can be obtained in any direction.Therefore,this method can be used to characterize the anisotropy of the deformation modulus of jointed rock masses.Hence,an anisotropy index for the deformation modulus is de fined to describe the anisotropy of the deformation modulus of jointed rock masses.

2.Analytical solution of the deformation modulus of a rock mass with persistent joints

2.1.The case with a single persistent joint set

First,a square rock mass is taken as the research object(Fig.1).Morland(1974)proposed that the deformation of a rock mass is the superposition of deformation of the joints and rock matrix,which can be expressed as

whereε,εrandεaare the strains of the rock mass,the rock matrix and the joints,respectively.Note that in this study,“εa”,not“εj”,is used to represent the strain of the joints in order to avoid the confusion with the joint set numberjwhich is used in the later sections.For simplicity,the rock matrix is considered to be a linearly elastic material,thusεrcan be expressed as

whereEris the deformation modulus of the rock matrix andσis the normal stress acting on the upper and lower boundaries of the rockmass.Thus,the key to calculating the strain of the rock mass is to calculate the strain of the joints.As shown in Fig.2,the deformation of joints can be decomposed into shear deformation and normal deformation(Morland,1974).For a single joint,the shear deformation and normal deformation can be represented by the constitutive model proposed by Goodman(1968),wherein these components are expressed as

whereδnandδsare the normal and shear deformation of the joint,respectively;σnandτsare the normal and shear stresses acting on the joint,respectively;andknandksare the normal and shear stiffnesses of the joint,respectively.

Note that the final deformation modulus of a rock mass is not affected by the whole deformation of the joints;only the deformation of the joints in the loading direction contributes to the final deformation modulus.Hence,the deformation of joints that can contribute to the deformation modulus can be expressed as

whereuais the deformation of the joints that can contribute to the deformation modulus.The strain corresponding to this component of joint deformation can be expressed as

When the normal stress acts on the upper and lower boundaries of the rock mass,based on a force equilibrium analysis,the normal and shear stresses acting on the joint can be expressed as(Goodman,1989):

Therefore,the deformation of a persistent joint under normal stress is calculated by combining Eqs.(3)-(5)and(7)and(8),and this deformation can be expressed as(Li,2001):

Note that Eq.(9)is only applicable to persistent joints(Fig.1)that do not intersect the upper and lower boundaries.However,some persistent joints intersect the upper or lower boundaries,as shown in Fig.3.For those joints,the deformation can be obtained by modifying Eq.(9),and this modi fied form can be expressed as(Yang et al.,2011):

Therefore,as shown in Fig.4,the deformation of the joint represented by the black line in the loading direction can be calculated by Eq.(9),while the joint represented by the red line can be calculated by Eq.(10).However,for persistent joints that do not intersect the upper and lower boundaries(represented by black line),lcosθ/L=1,and Eq.(10)reduces to Eq.(9).Hence,Eq.(10)can be applied to all types of persistent joints,regardless of whether they intersect the upper and lower boundaries.

Fig.3.Persistent joints intersecting adjacent boundaries or upper and lower boundaries.

Fig.4.Models with(a)persistent parallel and(b)nonparallel joints.Note that the black lines represent the joints that do not intersect the upper and lower boundaries,whereas the red lines represent the joints that intersect the upper or lower boundaries.

2.2.The case with multiple persistent joint sets

For simplicity,first,a set of persistent joints is considered(Fig.4),which is labeledj.The linear frequency is denoted as λmjin the mean normal direction of joints,and the acute angle between the loading direction and the mean normal direction of the joints is denotedθmj.According to Fig.4,the number of joints in the model is equal to the number of joints intersecting left and lower boundaries or right and upper boundaries(Fig.5).Therefore,the number of joints in the model can be expressed as

whereNjis the number of joints from setj,λmjcosθmjis the linear frequency in the direction of the left or right boundary,and λmjsinθmjis the linear frequency in the direction of the lower or top boundary.

The deformationuajof all persistent joints from setjin the loading direction can be expressed as

whereXjrepresents theXvalue of a joint from setj,whereinX=kn,ks,θ,orl;andE(·)orE{·}is the expected value of the function.We assume that the distributions ofθjand the trace lengthljof the joints are independent from each other,and the values ofknjandksjof the joints from the same set are constant.Therefore,Eq.(12)can be rewritten as

For simplicity,Eq.(13)can be rewritten as

Moreover,according to the de finition of the expectation,Eq.(14)can be rewritten as

wheregθj(θj)is the probability density function ofθj,glj(lj)is the probability density function of the trace lengthlj,andwjis the range of values forθj.Then,substituting Eq.(11)into Eq.(16)can give

Furthermore,if there areJsets of persistent joints in the area,according to Eq.(17),the total deformationuaofJsets of joints can be written as

The strain of all joints in the loading direction can be given by

Then,the deformation modulus of the rock mass can be expressed as

Therefore,the deformation modulus of the rock mass can be calculated by combining Eqs.(1),(2)and(18)-(20)as

According to the de finition of the expectation,Eq.(21)can be rewritten as

In practice,we can easily estimate the deformation modulus of a rock mass with multiple persistent joint sets according to Eq.(22)based on the geometric and mechanical parameters of joints and intact rocks.

3.Analytical solution of the deformation modulus of a rock mass with non-persistent joints

3.1.The case with a single non-persistent joint

It is well known that there are a very limited number of actual cases where the joints are persistent,i.e.joints are usually nonpersistent.Therefore,this section will discuss a method for estimating the deformation modulus of rock masses that contain nonpersistent joints.

There is no doubt that if the method for calculating the deformation modulus of rock masses with persistent joints is applied to rock masses with non-persistent joints,the obtained deformation modulus will be smaller than the actual value.This underestimation occurs because,for a rock mass with persistent joints,all deformation along the plane of a given joint is due to the contribution from the persistent joints,whereas for a rock mass with non-persistent joints,the deformation is a combination of contributions from the non-persistent joints and the rock bridges.Furthermore,rock bridges are stiffer than joints under the same load.Hence,the existence of a rock bridge can constrain the deformation of non-persistent joints and result in a negative contribution to the deformation of non-persistent joints.

It is obvious that the normal and shear stiffnesses can describe the deformation capacity of joints in the normal and shear directions,respectively.Similarly,the deformation and shear moduli can also describe the deformation capacity of a rock bridge in the normal and shear directions,respectively.However,there is an essential difference between the two mechanical parameters.When the stress is divided by the stiffness,the deformation is obtained,whereas when the stress is divided by the modulus,the result is strain.Therefore,further processing is needed in order to use the deformation and shear moduli to describe the deformation capacity of the rock bridge.Kulatilake et al.(1992)proposed the concept of fictitious joints to represent the rock bridges for modeling a rock mass containing non-persistent joints using discrete element method.As stated by Kulatilake et al.(1992),it is necessary to generate fictitious joints when they are combined with the non-persistent joints and discretize the problem domain into polygons in 2D space or into polyhedra in 3D space.The proposed fictitious joint(Fig.6)is a kind of joint,whose mechanical properties are identical to those of intact rock.Kulatilake et al.(1992)also suggested a method for estimating the normal and shear stiffnesses for the fictitious joint,which is expressed as

Fig.5.Intersection situations of joints and boundaries:(a)All joints;(b)Joints intersecting the left and lower boundaries;and(c)Joints intersecting the right and upper boundaries.The numbers of joints in(a),(b)and(c)are the same.

Fig.6.Adding fictitious joints to represent a non-persistent joint.l b is the length of the rock bridge.

wherekfsandkfnare the shear and normal stiffnesses of the fictitious joint,respectively;andGrandErare the shear and deformation moduli of the rock matrix,respectively.Fortunately,because the rock bridge is a part of the rock,the normal and shear stiffnesses of the fictitious joint can be applied to quantitatively describing the deformation capacity of the rock bridge.After this process,the mechanical parameters that quantitatively describe the deformation capacity of joint and rock bridge are uni fied.

Similar to Eq.(10),the deformation of the rock bridge can be given by

whereubis the rock bridge deformation that can contribute to the deformation modulus.

Jennings(1970)computed the weighted mean for the shear strength parameters of joint and rock bridge by regarding the connectivity rate as the weight,and then the weighted mean was taken as the shear strength of the rock masses.The detailed processes in this computation can be expressed as

wherecaandcbare the cohesions of the joint and rock bridge,respectively;φaandφbare the friction angles of the joint and rock bridge,respectively;andkis the connectivity rate of a nonpersistent joint.Eqs.(26)-(28)are called the Jennings strength criteria,which are widely used due to their simplicity and accuracy.Similarly,the weighted mean can also be applied to calculating the deformation of non-persistent joints.More speci fically,this process can be described in three steps.For simplicity,a single nonpersistent joint is considered.First,the connectivity rate should be calculated.Subsequently,the deformation of the non-persistent joint without constraints and the deformation of the rock bridge corresponding to the non-persistent joint in the loading direction should be represented by Eqs.(10)and(25),respectively.Finally,the weighted mean of the two deformation components is calculated using the connectivity rate as the weight,which is taken as the deformation of the non-persistent joint in the loading direction and can be expressed as

whereunpais the deformation in the loading direction of the nonpersistent joint.Substituting Eqs.(10)and(25)into Eq.(29)results in an equation of the deformation of the non-persistent joint in the loading direction,which is expressed as

For persistent joints,k=1 and Eq.(30)reduces to Eq.(10).Therefore,Eq.(30)is applicable to any type of joint,including persistent and non-persistent joints.

3.2.The case with multiple non-persistent joint sets

For simplicity,first,a set of non-persistent joints is considered,which is labeledj.Then,a scanline along the mean normal direction of joint setjis set,and the linear frequency along the scanline and the joint plane density are denotedλmjandλpj,respectively.Moreover,the relationship betweenλmjandλpjcan be expressed as(Zheng et al.,2020):

where njis the unit normal vector of the joints and lmjis the unit vector of the scanline,which is in the mean normal direction of joint setj.Therefore,the number of joints from setjcan be expressed as

Similar to Eq.(12),based on Eq.(30),the deformation of all nonpersistent joints from setjin the loading direction can be expressed as

wherekjis the connectivity rate of one non-persistent joint from setj.We assume that the distributions ofθjand the trace lengthljof the joints are independent from each other,and the values ofknjandksjof the joints from the same set are constant.Then,Eqs.(33)and(34)can be written as

Moreover,according to the de finition of expectation,substituting Eq.(32)into Eq.(35)can give

whereglb(lb)is the probability density function of the rock bridge.Furthermore,if there areJsets of non-persistent joints in the area,Eq.(37)can be written as

In summary,the deformation modulus of rock masses with joints can be estimated by combining Eqs.(1),(2),(19),(20)and(38),the result of which can be expressed as

According to the de finition of the expectation,Eq.(39)can be rewritten as

In practice,we can easily estimate the deformation modulus of a rock mass with multiple non-persistent joint sets according to Eq.(40)based on the geometric and mechanical parameters of joints and intact rocks.

4.Anisotropy index for the deformation modulus of rock masses

Generally,it is not easy to quantitatively represent the anisotropy of the deformation modulus of a rock mass,since it is dif ficult to obtain the values of the deformation modulus of a rock mass in different directions.As described above,this study provides a method to estimate the deformation modulus of a rock mass in any direction,which makes it possible to quantitatively represent the aforementioned anisotropy.Consequently,an anisotropy index for the deformation modulus is de fined as

whereEmmaxandEmminare the maximum and minimum values of the deformation modulus of the rock mass in different directions,respectively;andAIdmis the anisotropy index of the deformation modulus.The following conclusions can be drawn from Eq.(41):(i)ifEmmax=Emmin,all the values of the deformation modulus in the different directions are equal,thusAIdm=0;(ii)ifEmmax≫EmminorEmmin=0,thenAIdm=1;and(iii)the range ofAIdmis from 0 to 1,and the more anisotropic the rock mass is,the larger the value ofAIdmwill be.Note that a similar index was de fined by Zheng et al.(2020)to quantitatively represent the anisotropy of the permeability of rock masses.

5.Numerical validation of proposed method

5.1.A rock mass with persistent joint sets

5.1.1.Procedureforthenumericalexperiments

To evaluate the effectiveness of the proposed method for rock masses with persistent joint sets,some experiments were conducted with the universal distinct element code(UDEC)(Itasca,2015).UDEC was originally formulated for evaluating fluid flow paths through fractured rock masses and is currently used in numerical experiments for determining the mechanical properties of fractured rocks(Min et al.,2004).Several material behavior models have been built in UDEC for both the intact blocks and joints,which permit the simulation of response representative of jointed geologic(or similar)materials.Based on a“Lagrangian”calculation scheme,UDEC is well-suited to model the large movements and deformations of a rock mass system.The joints in this study are assigned to joint area contact model,which provides a linear representation of joint stiffness and yield limit.

Eight 10 m×10 m discrete fracture network(DFN)models(i.e.eight groups in Table 1)were generated by Monte Carlo simulations in MATLAB and the coordinate data of fractures were saved.In UDEC,it is not easy to mesh complex networks including many nonparallel joints.Thus,in the simulations,the trace angles of the joints from the same set were assigned constant values,i.e.parallel joint sets were used in the simulations.Since the purpose of this study is to investigate the deformation modulus in all directions,a 5 m×5 m block in the 10 m×10 m generated model rotates 15°around the center of the zone in every rotation,from 0°to 345°,for each DFNmodel(Fig.8).Therefore,there are 24 blocks of 5 m×5 m from each 10 m×10 m model(i.e.a group in Table 1)to be imported to UDEC.Since the coordinate data of original joints were saved before,the joint coordinate data within the range of 5 m×5 m after a certain angle of rotation can be obtained through coordinate transformation.The related code for clipping and coordinate transformation was developed in MATLAB.By applying a unidirectional load on the upper and lower boundaries and recording the corresponding stress and strain data,the stress-strain curve of each 5 m×5 m DFN model can be obtained.It should be noted that the value of the deformation modulus obtained by UDEC is the slope(secant modulus)from the origin point to the point on the stress-strain curve where the stress is equal to half of the uniaxial compressive strength,σc(Fig.7).Some mechanical parameters needed in this study are listed in Table 2.

Table 1Geometric parameters of the persistent joints for generating the DFN models.

Table 2Mechanical parameters of joints and intact rock.

Fig.7.Calculation of deformation modulus based on the stress-strain curve from UDEC(a schematic diagram).

5.1.2.Resultsandanalysis

For Group 1,which has a set of persistent joints,the predicted and numerical values of directional deformation modulus are shown in Fig.9.Note that(i)the directional deformation modulus represents the deformation modulus in different directions;(ii)the predicted values are obtained according to Eq.(22);(iii)the values of the directional deformation modulus obtained from UDEC are denoted as the numerical values;and(iv)as described above,the 5 m×5 m block in the 10 m×10 m generated model rotates 15°around the center of the zone in every rotation,from 0 to 345°,for each DFN model(Fig.8),and so for each DFN model,there are 24 values of the deformation modulus in 24 directions.

Fig.8.The 5 m×5 m DFN models with two sets of persistent joints with an intersection angle of 30°(Group 4 in Table 1).Dimensions are in m.Note that there are twenty four 5 m×5 m blocks for each 10 m×10 m model,but eight models are shown in this figure due to limited space.

The results in Fig.9 show that the predicted values of the deformation modulus in different directions are close to the numerical values,and their change trends in different directions are similar.To quantitatively compare the predicted and numerical values of the directional deformation modulus,the mean error rates are calculated for the 24 predicted values for this DFN model and the mean value is 6.41%.Note that the mean error rate is the mean value of the error rates in 24 directions,where for each direction the error rate is equal to the ratio of numerical value minus predicted value to the numerical value.In addition,the predicted and numerical values ofAIdmfor Group 1 are calculated according to Eq.(41)and the values are both equal to 0.85.By the above analysis,it seems that the method for estimating the deformation modulus of a rock mass with one persistent joint set is valid.

Fig.9.Predicted and numerical values of the directional deformation modulus of the DFN models with one set of persistent joints(Group 1).The values of the directional deformation modulus from UDEC are regarded as the numerical values.The numbers around the circle are in degree.

Fig.10.Predicted and numerical values of directional deformation modulus of the DFN models containing two sets of persistent joints:(a)Group 2,(b)Group 3,(c)Group 4,(d)Group 5,(e)Group 6,(f)Group 7,and(g)Group 8.The numbers around the circle are in degree.

Fig.10 shows the results of the directional deformation modulus for the DFN models containing two sets of persistent joints with different intersection angles(Groups 2-8).The results also indicate that the predicted values are close to the numerical values,and their change trends in different directions are similar.The mean error rates are also calculated for 24 predicted values in each DFN model,and the results are shown in Fig.11.The following conclusions can be drawn:

(1)When the intersection angle of the two persistent joint sets varies between 0°and 90°,the mean error rate of the predicted values of the deformation modulus varies between 5.06%and 7.12%,all of which are at satisfactory levels.

(2)As the intersection angle of the two persistent joint sets increases,the mean error rate of the predicted values of the deformation modulus fluctuates by approximately 5.9%,i.e.the error rate seems to have no relation with the intersection angle of two joints.

(3)The error rates for the DFN models with two sets of persistent joints are at the same levels as those with one persistent set of joints;therefore,the error rate does not increase as the number of joint sets increases.

The predicted and numerical values ofAIdmfor Groups 2-8 are also calculated with Eq.(41),and the results are shown in Fig.12.The following conclusions can be drawn:

(1)For different groups,the predicted values ofAIdmare very close to the numerical values,and all the error rates are less than 12.71%.

(2)As the intersection angle of the two joint sets increases,the values ofAIdmdecrease,i.e.the larger the intersection angle of the two joint sets is,the less anisotropic the rock mass is,which is sensible for a rock mass containing two sets of joints.

(3)When the intersection angle of the two joint sets is 75°,AIdmreaches the smallest value,i.e.the rock mass has the weakest anisotropy.This is because the anisotropy of the deformation modulus of a rock mass is related to the normal and shear stiffnesses of joint sets in addition to the intersection angle of the two joint sets.Fig.12 further indicates that the proposed method has a good capacity to estimate the deformation modulus of rock masses with two persistent joint sets and can be the basis for studying the method of estimating rock masses with non-persistent joints.

Fig.11.Mean error rates of the predicted and numerical values for Groups 2-8.

5.2.A rock mass with non-persistent joint sets

5.2.1.Procedureforthenumericalexperiments

To evaluate the method developed for a rock mass with nonpersistent joint sets in this study,some experiments were also conducted with the UDEC(Table 3).Although it is easy to model rock masses with persistent joints in UDEC,it is dif ficult to model rock masses with non-persistent joints in UDEC.Fortunately,a fictitious joint is introduced to solve this problem.The fictitious joint as proposed by Kulatilake et al.(1992)is a kind of joint whose mechanical properties are identical to those of intact rock(Fig.6).Therefore,it is necessary to add fictitious joints for each nonpersistent joint.To improve ef ficiency,a MATLAB code is developed to add fictitious joints for each non-persistent joint in this study.The coordinates of the non-persistent joints are input into the code,and the resulting output can be input to UDEC.DFN models with non-persistent joints(Fig.13)were generated according to Table 3.Then,through the above code and the same process as described in Section 5.1.1,the deformation modulus of each 5 m×5 m DFN model can be obtained.

Table 3Geometric parameters of the non-persistent joints for generating the DFN models.

Note that the value ofGr/kfsin Eq.(23)is assigned as 0.01 while using Eq.(40)to calculate the value of deformation modulus,which is the middle of the range of values.Huang and Huang(1998)derived a formula to estimate the connectivity rate based on the measuring window method,and this formula can be expressed as

wheren1jis the number of joints that are visible at one end from setj,n2jis the number of joints that are visible at both ends from setj,andn0jis the number of joints that are not visible at both ends from setj.The connectivity rate of a rock mass can be estimated by Eq.(40)through the application of a measuring window.Hence,Eq.(42)was used to calculate the connectivity rate of a rock mass found in Eq.(40).

Fig.12.Comparison between the predicted and numerical values of AI dm for Groups 2-8.

Fig.13.The 5 m×5 m DFN models with two sets of joints with an intersection angle of 30°(Group 16).Dimensions are in m.Note that there are twenty four 5 m×5 m blocks for each 10 m×10 m model,but eight models are shown in this figure due to limited space.

5.2.2.Resultsandanalysis

The results in Fig.14 show that the predicted values of the deformation modulus in different directions are very close to the numerical values,and their change trends in different directions are similar.To quantitatively compare the predicted and numerical values of the directional deformation modulus,the mean error rates are also calculated for the 24 predicted values for each DFN model,and the results are shown in Fig.15.The following conclusions can be drawn:

(1)When the expected value of the trace length of the joint set varies from 1 m to 5 m,the mean error rate of the predicted deformation modulus varies between 12.16%and 20.75%.

(2)As the expected value of the trace length increases,the mean error rate of the predicted deformation modulus first increases and then remains relatively stable.This is likely caused by the boundary effect of the joints(Zheng et al.,2015),which increases as the size of the joints increases.This boundary effect indicates that some of the joints that intersect the 5 m×5 m block have center points outside of the 5 m×5 m block.

(3)In addition,the reason that the results of Groups 11-13 have relatively large mean error rates than those of Groups 9 and 10(Fig.15)is that most of the deformation moduli in different directions of Groups 11-13 are relatively smaller than those of Groups 9 and 10(Fig.14),and a same absolute error will bring a larger relative error rate.Actually,the absolute errors of Groups 11-13 seem to be smaller than those of Groups 9 and 10(Fig.14).

In order to evaluate the error caused by this method in estimating the deformation modulus of rock mass,the results obtained by the proposed method and some empirical and test methods are compared.As for in situ tests,Bieniawski(1989)stated that‘verydifferentin situresultsmaybeobtaineddependingonthetestmethod.Eveninanextensivein situtestprogramin fairlyuniformandgoodqualityrockmassconditions,deformability datamayfeatureadeviationof25%,orasmuchas10 GPa’.For laboratory tests,Farmer and Kemeny(1992)found that the results from laboratory tests are 5-20 times higher than those of the in situ test,which may due to the failure of the specimen size to reach the representative elementary volume(REV)of the rock mass.Even for empirical methods,the error is inevitable.With respect to a rock mass withGSI=70,disturbance factorD=0,and intact rock deformation modulusEr=50 GPa,Shen et al.(2012)applied three empirical equations proposed by Carvalho(2004),Sonmez et al.(2004),and Hoek and Diederichs(2006)to calculate the deformation modulus of the rock mass,and the results are equal to 21.7 GPa,25.6 GPa,and 36.6 GPa,respectively.The difference between the maximum and minimum values is up to(36.6-21.7)/21.7=69%.Besides,some scholars(Aksoy et al.,2012;Karaman et al.,2015)have also studied the differences between different empirical methods,which can even reach several times.Therefore,the mean error rate of approximately 16%for all cases with one set of joints indicates that the proposed method for estimating the deformation modulus is valid.

Fig.14.Predicted and numerical values of the directional deformation modulus of the DFN models with one set of joints:(a)Group 9,(b)Group 10,(c)Group 11,(d)Group 12,and(e)Group 13.The numbers around the circle are in degree.

In addition,the predicted and numerical values ofAIdmfor Groups 9-13 are calculated according to Eq.(41),and the results are shown in Fig.16.This figure shows that for different groups,the predicted values ofAIdmare very close to the numerical values,and all the error rates are less than 17%.Moreover,with an increase in the expected value of the trace length of the joints,the values ofAIdmincrease,i.e.the larger the joint sizes are,the more anisotropic the rock mass is,which is sensible for a rock mass containing one set of joints.Fig.16 presents a comparison between the predicted and real anisotropies,which further indicates that the developed method has a good capacity to estimate the deformation modulus of jointed rock masses with one set of joints.

The predicted and real values ofAIdmfor Groups 14-20 are also calculated using Eq.(41),and the results are shown in Fig.19.The following conclusions can be drawn:

Fig.15.Mean error rates of the predicted and numerical values for Groups 9-13.

Fig.16.Comparison between the predicted and numerical values of AI dm for Groups 9-13.

Fig.17 shows the results of the directional deformation modulus for the DFN models containing two sets of non-persistent joints with different intersection angles(Groups 14-20).The results also indicate that the predicted values are close to the numerical values,and their change trends in different directions are similar.The mean error rates are calculated for 24 predicted values in each DFN model,and the results are plotted in Fig.18.The following conclusions can be drawn:

(1)When the intersection angle of the two joint sets varies between 0°and 90°,the mean error rate of the predicted values of the deformation modulus varies between 12.98%and 22.03%,all of which are at satisfactory levels.

(2)As the intersection angle of the two joint sets increases,the mean error rate of the predicted values of the deformation modulus fluctuates by approximately 17%,i.e.the error rate seems to have no relation with the intersection angle of two joints.

(3)A comparison of Figs.15 and 18 shows that the error rates for the DFN models with two sets of joints are at the same level as those with one set of joints.Hence,the error rate does not increase as the number of joint sets increases.

(1)For different groups,the predicted values ofAIdmare very close to the numerical values,and the mean errors rates are less than 12.71%.

(2)Overall,as the intersection angle of the two joint sets increases,the values ofAIdmdecrease,i.e.the larger the intersection angle of the two joint sets is,the less anisotropic the rock mass is,which is sensible for a rock mass containing two sets of joints.

(3)When the intersection angle of the two joint sets is 75°,AIdmreaches the smallest value,i.e.the rock mass has the weakest anisotropy.This is because the anisotropy of the deformation modulus of a rock mass is related to the normal and shear stiffnesses of joint sets in addition to the intersection angle of the two joint sets.

(4)The anisotropy variation laws with the intersection angle for deformation modulus of rock mass with nonpersistent joints are similar to that of the rock mass with persistent joints.Fig.19 further indicates that the developed method has a good capacity to estimate the deformation modulus of jointed rock masses with two non-persistent joint sets.

In conclusion,the above error analysis shows that the analytical method for estimating deformation modulus proposed in this study is validated,but for parts of the cases,the error rates reach 20%.Therefore,to obtain reliable results of deformation modulus of rock masses,we recommend combining the multiple methods.

6.Discussion

The proposed method can only be applied to a rock mass having a REV.Because if a REV does not exist,the equivalent continuum model cannot be used to analyze the deformation behaviors of the rock mass.Therefore,the size of the analyzed rock mass should not be less than its REV when using the developed method.

The proposed method requires the values of the normal and shear stiffnesses,which can be obtained from uniaxial compression tests and shear tests.These test methods are relatively cheap and easy to perform,which increases the practicality of the method proposed in this study.In addition,the developed formula incorporates the variations in the orientations and sizes of joint sets.Note that in the numerical validation experiments,parallel joints are used for the same sets,because when a complicated nonparallel and nonpersistent joint network is input to UDEC,it requires fine meshing,and the computation process is likely to fail.Hence,UDEC does not have suf ficient capacity to deal with a complicated nonparallel and non-persistent joint network.However,it is very convenient to use the proposed method to deal with the aforementioned complicated networks.For example,when the trace angle in Group 9 is normally distributed with a standard deviation of 5°,the predicted values of the directional deformation modulus of this case are plotted in Fig.20.Moreover,the proposed method can re flect the directional deformation modulus and quantitatively describe the degree of anisotropy of rock masses.

Considering that it is not easy to perform a validation procedure for 3D complicated joint networks,the proposed method is developed in 2D space.However,actual rock masses exist in 3D space,and their deformation behaviors usually have obvious 3D effects.Therefore,the proposed method will be extended to 3D space in future research.

Fig.17.Predicted and numerical values of directional deformation modulus of the DFN models containing two sets of joints:(a)Group 14,(b)Group 15,(c)Group 16,(d)Group 17,(e)Group 18,(f)Group 19,and(g)Group 20.The numbers around the circle are in degree.

Fig.18.Mean error rates of the predicted and numerical values for Groups 14-20.

Fig.19.Comparison between the predicted and numerical values of AI dm for Groups 14-20.

Fig.20.Predicted values of the deformation modulus when the joint angle in Group 9 is normally distributed and the standard deviation is 5°.The numbers around the circle are in degree.

Finally,since it is extremely dif ficult to carry out a veri fication test of this method in practice,the proposed method was only evaluated by numerical modeling herein.Therefore,the proposed method must still be tested in practical applications.

7.Conclusions

A new method for estimating the deformation modulus of jointed rock masses is proposed.The developed formula(Eq.(39)or(40))incorporates the variations in the orientations and sizes of the joint sets and can accordingly be used to study the directional deformation modulus of jointed rock masses.Moreover,an anisotropy index for the deformation modulusAIdmis de fined to quantitatively describe the deformation modulus anisotropy.

Moreover,some numerical experiments were conducted with UDEC to evaluate the method proposed in this study.The results show that the predicted values of the deformation modulus are close to the numerical values,and the speci fic results are summarized as follows:

(1)For a rock mass with a set of persistent joints,the mean error rate of the deformation modulus in 24 directions is 6.41%.When the intersection angle of the two persistent joint sets varies between 0°and 90°,the error rate of the predicted values of the deformation modulus varies between 5.06%and 7.12%.

(2)For a rock mass with one non-persistent joint set,when the expected value of the trace length of the joint set varies from 1 m to 5 m,the mean error rate of the predicted values of the deformation modulus varies between 12.16%and 20.75%,all of which are at satisfactory levels.Moreover,the larger the joint sizes are,the more anisotropic the rock mass is.

(3)For a rock mass containing two non-persistent joint sets,when the intersection angle of the two joint sets varies between 0°and 90°,the mean error rate of the predicted values of the deformation modulus varies between 12.98%and 22.03%,all of which are at satisfactory levels.Moreover,the larger the intersection angle of the two joint sets is,the less anisotropic the rock mass is.

Therefore,the proposed method for estimating the deformation modulus is valid.In addition,this method lays a foundation for future studies on estimating the deformation modulus in 3D space.

Declarationofcompetinginterest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to in fluence the work reported in this paper.

Acknowledgments

This study was funded by the National Key R&D Program of China(Grant Nos.2017YFE0119500 and 2018YFC1505005),and the National Natural Science Foundation of China (Grant No.41972264).