## Abstract

This paper presents an application of a 3D fracture network model and 3D persistence to confirm the representative volume element (RVE) in a fractured rock mass. Based on field fracture data collected from the Songta dam site, a large 3D fracture network is generated to analyse the RVE. The 3D persistence, which is governed by the comprehensive characteristics of fracture location, orientation, size and density, strongly affects the shear strength of fractured rock masses. Hence, it is valuable to confirm the RVE size considering the persistence, which combines the geometric and mechanical properties of a rock mass. The persistence values of the entire model and of cubes of different sizes, which are located in the centre of the model, are determined. A relative error of 10% is selected as the benchmark to evaluate the similarity between the cubes and the model. Finally, 15 m is confirmed as the size of the RVE in the study area.

The presence of discontinuities induces a significant size effect on the rock mechanical and hydraulic response (Song *et al.* 2015). Rock mass behavioural parameters, such as compressive strength, longitudinal elastic modulus and permeability coefficient, vary with the size of sample under study. When the sample size is relatively small, the rock mass parameters vary strongly with increasing size of the sample. The parameters no longer vary when the sample size has increased to a critical value, which is the size of the representative volume element (RVE) (Zhang *et al.* 2012). Knowledge of the RVE size contributes to the analysis of the behaviour of rock masses. Once the RVE size is determined, modelling of fractured rock mass problems at a scale equivalent to or larger than that of the RVE is possible using equivalent homogeneous properties (Esmaieli *et al.* 2010). Hence, the determination of the RVE size is an essential step before the generation of the continuum model for a fractured rock mass.

The first formal definition of the RVE was proposed by Hill (1963). According to this definition, the RVE should be structurally entirely representative of the whole mass, and the sample should contain a sufficient number of discontinuities, which can be used to represent the mean and variance of the entire discontinuities. Subsequently, various methods to confirm the RVE size were proposed.

Kulatilake (1985) adopted the finite element method (FEM) to study the size effect for rock masses and determine the RVE size. Castelli *et al.* (2003) studied the scale effects on the stiffness modulus based on the displacement discontinuity method (DDM) and the fictitious stress method (FSM). Min & Jing (2003) applied the distinct element method (DEM) to study the size effect on the equivalent elastic compliance tensor and to confirm the RVE size of a fractured rock mass. Pierce *et al.* (2007) generated a 3D fracture network using field data and applied 3D particle flow code (PFC3D) to estimate the RVE size. These methods assess the mechanical parameters of the RVE of fractured rock masses based upon numerical analyses.

The anisotropic, inhomogeneous and discontinuous features of rock masses are often characterized by stochastic variability. Hence, it is important to analyse the RVE size considering the variability of the discontinuities in a rock mass. Zhang *et al.* (2013*a*,*b*) generated a 3D fracture network model and determined the RVE size of the rock mass at the Baihetan dam site. In the 3D model, the volumetric fracture density (*P*_{32}), fracture dip direction, dip angle and size were considered in determining the RVE.

In this paper, a new method that considers the discontinuity persistence is presented for determining the RVE size. The discontinuity persistence is the degree of connectivity of the discontinuities along a given plane (line) (Fig. 1). In a 2D model, the discontinuity persistence is expressed as the ratio of the total length of the discontinuity traces to the length of a given line lying in a discontinuity plane (Jennings 1970; Zhang 1989). In a 3D model, the persistence is defined as the ratio of the total area of the discontinuity planes to the area of a given plane containing a combination of discontinuity planes and intact rock regions (Einstein *et al.* 1983). The persistence along a certain plane, which is decided jointly by the location, orientation, size and density of discontinuity around the plane, is one of the most significant discontinuity parameters and affects the shear strength of fractured rock masses and the stability of rock slopes (Park 2005). Hence, it is valuable to confirm the RVE size based on the persistence.

## Generation of a 3D fracture network

### Study area

This paper takes the Songta hydropower station as the study area to carry out the study of the RVE size for a granite body affected by discontinuities. This hydropower station, which is now being constructed, lies in the upper reaches of the Nu River in SW China. The flow direction of the Nu River at the dam site is *c.* 188° SW (Fig. 2).

The valley at the dam site exhibits an asymmetric ‘V’ shape, which is a typical mountain-canyon geomorphology (Fig. 3a). A concrete double-curved arch dam, with a maximum height of 318 m, is planned to be built in the valley. The elevations of the normal water level and storage level are 1700 and 1925 m above sea level (a.s.l.), respectively. The hydroelectric generating capacity of the station will be 3600 MW.

Biotite monzonitic granite and plagioclase amphibolite of the late Yanshanian (Cretaceous) period are exposed at the dam site. Biotite monzonitic granite is the predominant lithology and is mainly composed of quartz, plagioclase, potassium feldspar and biotite. Plagioclase amphibolite is the main component of the dykes, which have widths varying from 0.05 to 5 m. Although dozens of dykes are developed randomly at the dam site, these dykes are almost parallel and the general orientation is N20−35°W/SW∠78−87°. In the contact region between the two types of lithology, the plagioclase amphibolite has experienced slight argillization (Fig. 3b).

To investigate the engineering geological properties of the rock mass, adits were excavated in the dam abutment. The PDC3 adit, with an elevation of 1780 m a.s.l. on the right bank of the Nu River, has been used to determine the RVE (Fig. 4a). The strike of the adit is approximately perpendicular to the flow direction of the Nu River. The adit is 100 m long, with a cross-section 2 m wide and 2 m high. The window sampling method was utilized to collect the discontinuity data in the adit (Song & Lee 2001). All the discontinuities with trace lengths larger than 0.5 m that were found on the downstream rock wall of the PDC3 adit were measured. In the process of measurement, discontinuity information, including orientation, trace length, aperture, filling, roughness and weathering, was recorded. These data are used for the 3D fracture network simulation. The map of the 2D fracture traces is shown in Figure 4b. In the PDC3 adit, the fracture distribution is relatively homogeneous and the trace lengths are mainly concentrated between 0.7 and 1.5 m. In addition, the majority of the fractures dip outside the river valley.

### Fracture network modelling

Numerical fracture network modelling can be utilized to represent the geometric characterization of rock masses (Dershowitz & Einstein 1988). In this study, the discontinuities are considered to be thin circular discs with a finite diameter. The simulation procedures are as follows (Baecher *et al.* 1977).

(1) Demarcation of preferential fracture sets.

The fracture sets in the PDC3 adit were derived according to the method proposed by Chen *et al*. (2005*b*). The fractures were divided into three groups, as shown in Table 1 and Figure 5.

The sampling bias on orientation of each fracture set was then corrected in the light of the method introduced by Kulatilake & Wu (1984*a*). In this method, the weighting function *W _{i}* was utilized to correct the observed frequency of orientations collected using the sampling window:
(1)where

*w*and

*h*are the width and height of the sampling window,

*θ*is the dip angle of the

_{i}*i*th fracture,

*d*is the diameter of the

_{i}*i*th fracture and

*β*is the angle between the strike of the sampling plane and the dip direction of the

_{i}*i*th fracture.

The Fisher constant *K*, which is used to reflect the clustering degree of each fracture set, was obtained according to the method proposed by Priest (1993):
(2)where *M* is the number of fractures and *r _{i}* is the unit normal vector of the

*i*th fracture.

(2) Determination of trace length, diameter and density.

Owing to sampling errors on the finite window, the observed trace lengths are different from the real trace lengths, and should be corrected in accordance with the method proposed by Kulatilake & Wu (1984*b*):
(3)where *R*_{0} is the proportion of fractures with both ends censored within the sampling window, *R*_{2} is the proportion of fractures with both ends observable within the sampling window, *θ* is the dip angle of fractures and *g*(*θ*) is the probability density function of *θ*. The Kolmogorov–Smirnov goodness-of-fit test was used to identify the best probability distribution to represent the distribution of trace lengths. The corrected trace lengths, which showed a Gamma distribution, are listed in Table 1.

Usually, the fracture is deemed to be a disc. Based on the corrected trace lengths, the diameter of each fracture set can be easily determined using the method of Kulatilake & Wu (1986). The goodness-of-fit test was utilized to determine the diameter probability distribution according to the distribution of the corrected trace lengths (Table 1).

The 2D fracture network can be determined according to the fracture information of the starting point coordinate, terminal point coordinate and orientation. Scanlines are generated to calculate the spacing of each fracture set in the 2D fracture network. Based on the relationship between 2D and 3D fracture spacing, the volume density of each fracture set can then be calculated by (Oda 1982)
(4)where and are the volume density and the line density of the *j*th set and *E*(*D*^{2}) is the second moment of the fracture diameter distribution.

(3) Monte Carlo simulation.

The fracture number can be determined according to the density and network size. As required by the engineering, the size of the 3D fracture network was designed to be 100 × 50 × 50 m. Based on the parameters of each fracture set, including orientation, diameter and density, the 3D fracture networks can be generated using the Monte Carlo simulation (Chen 2001).

(4) Model verification.

Using the Monte Carlo simulation, a number of fracture network models can be generated. A relatively suitable model should be chosen for RVE analysis. The generated 3D fracture network models can be cut along a plane, the direction of which should be the same as that of the surveyed window. The parameters of fractures in the cut plane should then be compared with those in the surveyed window. A model with a high degree of similarity to the investigated data can be selected as the optimal one, as shown in Table 2.

Based on the above steps, an optimal model was chosen to provide the basic data for RVE analysis. In the model, the fractures can be expressed using the following parameters: centre coordinates (*x*, *y*, *z*), diameter, dip direction and dip angle (Table 3). To provide a visual display for the fracture network model, a cube of size 20 m is extracted from the centre of the model and is illustrated in Figure 6.

Through the above statistics and analysis, three groups of fractures collected from the PDC3 adit have the following characteristics. The development law of the fracture orientation, trace length and spacing, is mainly affected by the stress state and tectonic history of the fractured region. Set 1 and Set 2 are fractures with steep dip angles, which are rarely and stochastically developed in the process of rock formation and tectonic movement. Set 3 consists of fractures dipping 120° SE at gentle angles, which are heavily and intensively developed because of the unloading effect of the river valley. Hence, the fracture spacings of Set 1 and Set 2 are similar and are larger than that of Set 3, and the Fisher constants of Set 1 and Set 2 are similar and are larger than that of Set 3. In addition, the average trace lengths of the three sets are similar because the trace lengths are mainly concentrated between 0.7 and 1.5 m.

## RVE determination

### Persistence calculation

The 3D persistence of the rock mass, *P*, is expressed as the ratio of the total area of the discontinuities along a given plane. *P* is governed by 3D fracture location, orientation, size and density, and is a very important parameter that influences the mechanical behaviour of fractured rock masses. In this study, *P* is calculated using the method proposed by Chen *et al*. (2005*a*). The calculation procedure is as follows.

(1) Determination of ideal cross-section.

The persistence, which possesses not only magnitude but also direction, is a vector. Hence, a cross-section should be determined in the 3D fracture network before calculating the persistence. Normally, the persistence within the ideal cross-section is calculated and confirmed according to the geological conditions and the engineering requirements.

(2) Selection of fractures around the cross-section.

For a cross-section, not all of the fractures have influence on the persistence. Therefore, the fractures that lie within a certain distance from the cross-section should be selected. In nature, the failure surface of a rock mass is not smooth and flat but results from fractures of different orientations. Hence, the fractures selected above are reselected to include only those for which the angle between the fracture and the cross-section is smaller than 15°.

(3) Projection of selected fracture onto the cross-section.

There is a certain distance and angle between the selected fracture and the cross-section. To calculate accurately the persistence, the selected fractures should be projected onto the cross-section. The location relationships between fractures and cross-section can be divided into three types (Fig. 7): (a) fracture surfaces with intersection; (b) fracture surfaces without intersection, the projections of which on the cross-section have an overlapping area; (c) fracture surfaces without intersection, the projections of which on the cross-section do not have an overlapping area (Lu *et al.* 2004).

(4) Elimination of overlapping fracture.

Overlapping regions between projection planes of fractures should be common (Fig. 8). In addition, these regions do not play a multiple role in the failure of rock masses. Hence, the overlapping area should be calculated only once and other overlapping areas should be eliminated.

(5) Persistence calculation.

The 3D persistence can be calculated according to the following:
(5)where *n* is the total number of the selected fractures, *a _{i}* is the projected area of the

*i*th selected fractures,

*a*

_{0}is the eliminated area of the overlapping regions and

*A*is the area of the ideal cross-section.

### RVE calculation

Compared with the strength and deformation of intact rock, fractures within the rock mass are weaker (Zhang *et al*. 2013*a*). Hence, fractures play a key role in determining the properties of rock masses. It is important that the RVE size is determined considering the fracture characteristics. The 3D persistence, which is governed by fracture location, orientation, size and density, seriously affects the shear strength of fractured rock masses. Therefore, a new method based on the 3D persistence, which combines geometric and mechanical properties of rock masses, is proposed for determining the RVE size. This method involves the following stages (Fig. 9).

First, to eliminate the spatial effect, a series of cross-section groups with different directions are arranged in the entire model. The cross-sections are parallel and have equal spacing in each group. When the mean of the 3D persistence values within a certain cross-section group reaches a maximum, the cross-section group will be selected as the optimal one. In this study, the dip direction of each cross-section is approximately equal to the preferential dip of each fracture set and the dip angle of each cross-section changes around the preferential angle of each fracture set (Table 1). The 3D persistence values of the cross-sections with different dip angles are calculated and are illustrated in Figure 10. The direction and persistence of the ideal cross-section are listed in Table 4.

Second, some cubes with side lengths varying from 1 to 25 m are extracted from the centre of the entire model (Fig. 11). The determined cross-sections of each fracture set are arranged in each cube. The 3D persistence values of the ideal cross-sections within different cube side lengths are counted. The size effect can be shown by changing the sample size (Fig. 12). When the sample size is relatively small, the persistence varies strongly with increase in sample size. However, the value tends to stabilize when the sample size has increased to a certain critical value.

Third, the size effect of each fracture set can be determined according to the error rate (Kanit *et al.* 2003; Pelissou *et al.* 2009):
(6)where *ε*_{rel} is the relative error of the 3D persistence within different cube side lengths, *P _{j}* is the persistence of the ideal cross-section in a cube of size

*j*metres and

*P*is the persistence of the ideal cross-section in the entire model. The relative error of each fracture set is shown in Figure 13. In this study, an error rate of 10% is selected as the benchmark to confirm the RVE size. Figure 13 shows that when the cube side lengths are greater than or equal to 14, 15 and 11 m, the error rates of Sets 1, 2 and 3, respectively, are less than 10%.

By comparison, the cube side size determined based on the persistence is related to the geometric size of the fracture set. When the fracture set developed in a rock mass has a similar size, the cube side size is mainly affected by the orientation and density of the fracture set. The smaller the Fisher constant and spacing of the fracture set, the smaller is the determined side size.

When determining the RVE of a fractured rock mass, the size effect on all the preferential orientations should be comprehensively considered. Hence, the largest value, 15 m, is selected as the RVE size at the Songta dam site.

## Conclusions

This paper presents an application of the 3D fracture network model and the 3D persistence to confirm the representative volume element (RVE) at the Songta dam site in SW China. The window sampling method was applied to collect field fractures. Through the procedures of the 3D simulation described above, the locations, orientations and sizes of all fractures are determined. This is helpful to analyse the size effect of rock masses in 3D space.

The 3D persistence, which is governed by the comprehensive characteristics of fracture location, orientation, size and density, seriously affects the shear strength of fractured rock masses; that is, the persistence can comprehensively reflect the geometrical and mechanical behaviour of fractured rock masses. Hence, it is valuable to confirm the RVE size considering the 3D persistence. Based on the 3D fracture network, the persistence of the entire model is calculated using the projection method.

Cubes with side lengths varying from 1 to 25 m are extracted from the centre of the entire model. The persistence of each cube is then determined. An index, relative error, is used to evaluate the similarity between the cubes and the entire model. In this study, an error rate of 10% is selected as the benchmark to confirm the RVE size. When the cube side lengths were greater than or equal to 14, 15 and 11 m, the error rates of Sets 1, 2 and 3, respectively, were less than 10%. When determining the RVE of a fractured rock mass, the size effect on all the preferential orientations should be comprehensively considered. Finally, the largest value, 15 m, is selected as the size of the RVE at the Songta dam site.

## Acknowledgements and Funding

This study was supported by National Natural Science of China (Grants 41330636, 41402242 and 41402243), 2010 non-profit scientific special research funds of Ministry of Water Resources (Grant 201001008), Jilin University's 985 project (Grant 450070021107) and Graduate Innovation Fund of Jilin University (Grant 2015013).

- © 2017 The Author(s)