## Abstract

A new method to study representative volume element (RVE) size is developed based on three-dimensional (3D) fracture numerical network modelling and a stochastic mathematical model. A 3D fracture numerical network is established, and the network factors that affect RVE size, such as fracture density, diameter, strike, and dip angle, are collected. This paper utilizes these four factors to obtain the RVE size based on the contingency table method. It is concluded that the side length of the RVE is 2–8 times larger than the average trace length. Fracture frequency is chosen to validate the RVE size through 3D fracture numerical network modelling. RVE size calculated through fracture frequency is consistent with the contingency table method. It is proven that the method for calculating RVE size through the contingency table method is feasible and credible.

Rock mass parameters, such as deformation modulus and permeability coefficient, change as a function of their sizes (Drugan & Willis 1996; Schultz 1996; Castelli *et al*. 2003; Chen *et al*. 2008). An increase in sample size results in a gradual change of parameters. When the size reaches a critical value, these parameter values no longer change with increase in sample size. The critical value is called the representative volume element (RVE), which can be used to analyse this complicated size effect. The results for a specimen can represent the entire rock mass only when its size is no less than the RVE. At the same time, the size also determines if the continuum model can be used in numerical simulations (Neumann 1998).

Mathematically, the RVE is a finite length size limit, where samples with the RVE size appear uniform and can represent the features of the whole rock mass (Min & Jing 2003). The first formal definition of the RVE states that it must be structurally typical of the whole mixture on average, and the sample should contain a sufficient number of inclusions for overall moduli to be effectively independent of the surface values of traction and displacement, so long as these values are macroscopically uniform (Hill 1963). Researchers have studied the RVE in different ways. Kulatilake (1985) studied the elastic modulus and RVE using the finite element method (FEM), and proposed that it is feasible to analyse RVE size using the discrete element method (DEM). Subsequently, the DEM, distinct element code (DEC), discontinuous deformation method (DDM), and particle flow code (PFC) were successfully employed to investigate the influence of the scale effect of fracture rock masses (Kulatilake *et al*. 1993; Grenon *et al*. 2001; Castelli *et al*. 2003; Min & Jing 2003; Esmaieli *et al*. 2010). Zhou & Yu (1999) analysed this issue statistically, and constructed a model to ensure a given accuracy of the overall estimated properties obtained by spatial averaging of the distribution characteristics and equivalent flexibility matrix; however, only the density of joints was considered. Lu *et al*. (2005) used a three-dimensional (3D) fracture numerical network model to create different samples with different sizes and confirmed the RVE size based on fractal theory. Pariseau *et al*. (2008) and Esmaieli *et al*. (2010) studied the RVE by adding joints until the samples do not change significantly with further joint addition, and solved the rock slope stability problem based on RVE analysis.

The anisotropic, inhomogeneous, and discontinuous features of rock masses are strongly affected by stochastically distributed discontinuities (Chen 2001; Hack *et al*. 2006; Brideau *et al*. 2009); thus, it is of great importance to analyse the RVE size considering the discontinuities in the rock mass. Density, strike, dip angle, and length are component elements of discontinuities. Therefore, it is important to use these factors to confirm the RVE size.

## Data acquisition

### General engineering situation

The studied rock mass is located within the dam area of the Baihetan hydropower station, which lies in the lower reaches of Jinsha River between Yunnan and Sichuan Provinces in China (Fig. 1). The station construction project is currently at the feasibility study stage. A concrete double curvature arch dam will be adopted for the Baihetan station. The dam is to be approximately 277 m high, with water surface elevation 590 m, dam crest 825 m, and storage water level 820 m. The hydroelectric generating capacity of the station is to be 4450 MW.

The valley in the dam area, which exhibits a mountain canyon geomorphology, has an asymmetric V shape. The cliffs on the left bank are at an elevation of 900–1100 m, whereas those on the right bank are at 1000–1400 m on the right bank. Mesas or gentle slopes occupy the top of the elevation.

No fault zone or large-scale fault lies along the river. The inner layer of the rock mass is unweathered (or only slightly weathered). The unloading degree of the rock mass is weak, and the weathering degree and depth gradually decrease from higher to lower elevations.

Basalt from the Permian period (P_{2}β) is the main component of the rock mass, and sandy shale from the Triassic period (T_{1}f) constitutes the overlying rock mass. The riverbed and mesa consist of Quaternary alluvium, diluvium, and eluvium. The site at which the dam will be built is mainly composed of basalt; therefore, the characteristics of basalt should first be considered in determining the stability of the dam.

The basalt is mainly composed of devonite, cryptocrystalline basalt, and amygdaloidal basalt. The basalt can be divided into 11 layers according to the intermittent and outbreak numbers. Each layer is composed of lava, breccia lava, and tuffaceous basalt from bottom to top. The rock formations strike NE–SW and dip 15–25° SE, and the columnar jointing in basalt develops inhomogeneously. The fractures in the basalt rock mass consist of high-angle fractures (70–85°) approximately parallel to the cylinders, and low-angle fractures (10–30°) normal to the cylinders. The trace lengths of low-angle fractures are 0.3–1 m, and those of the low-angle fractures are shorter than 0.2 m. The fracture traces are mainly flexuous, and the fracture surfaces are relatively rough. The fracture interspaces are filled with very thin calcite and quartz veins, usually well cemented, and a small fraction of the interspaces is filled with debris. In addition to the strong unloading zones and partially weak unloading areas, the fractures with steep dip angle are usually closed.

The river bed and the two sides are characterized by hard basalt. Locally weak rock burst phenomena exist, exfoliating the rock. Discal borehole cores 1–3 cm thick are distributed at an elevation of 385–470 m. The ground stress is low to medium. The maximum principal stress of the left bank is 5.0–12.4 MPa, with an average direction of 153° and an average plunge of 23°. The maximum principal stress of the right bank is 4.4–11.4 MPa, with an average direction of 7° and an average plunge of 12°. Earthquakes are mostly smaller than M6.0 according to historical records.

The annual average flow of the Jinsha River is 4110 m^{3} s^{−1}. The groundwater types are Quaternary pore water and fracture water, which are supplied by rainfall. Groundwater is discharged as springs, which converge with the gully and flow into the Jinsha River.

The discontinuities of exploration tunnel 1, which is on the left bank of the Jinsha River and is at an altitude of 675 m, are investigated. The strike of the tunnel is NE–SW and its dimensions are 100 m × 2 m × 2 m. Only structure fractures longer than 0.5 m were collected from the right lateral surface of the tunnel. The sampling window method was used in collecting the parameters of each fracture, such as the location, occurrence, trace length, filling, and shape. A total of 159 fractures were collected (Fig. 2).

## 3D fracture numerical network modelling

Numerical fracture network modelling can be used to represent a rock mass as an entity (Dershowitz & Einstein 1988). In this work, the discontinuities are considered as discs. The simulation procedures (Chen 2001; Starzec & Andersson 2002) are as follows.

Collection of fracture information, such as the location, occurrence, and trace length.

Preferential fracture sets demarcation; the following procedures aim at each fracture set respectively. The fracture sets were determined after the Kulatilake & Wu (1984

*b*) correction was performed (Fig. 3). Table 1 shows the occurrences of the four fracture sets. The Fisher constant was calculated to determine the degree of clustering (preferred orientation) within each fracture set (Fisher 1953; Priest 1993).Correction of the sampling bias on the fracture orientations. The sampling window method was used to collect the parameters of the fractures in the field. In practice, the size of the sampling window is finite and the fractures within the window were collected. The fractures of the sampling window display the following characteristics: (a) only one end of the fracture is surveyed; (b) both ends of the fracture are surveyed; (c) no fracture end is surveyed. These features cause the surveyed trace lengths to deviate; however, this deviation can be corrected using the method introduced by Kulatilake & Wu (1984

*a*). Fracture sets 1 and 3 show normal distributions, whereas fracture sets 2 and 4 show gamma distributions. The corrected trace lengths are listed in Table 1.Correction of the fracture trace length in the sampling window and determination of the probability density function of the corrected length. The method developed by Kulatilake & Wu (1986) was adopted to confirm the diameter of each fracture based on the corrected trace length. The diameters of the fracture sets show gamma or normal distributions, which are the same as those of the trace length (Table 1).

Simulation of the fracture density in 3D space. The spacing of each fracture set is highly important in generating the 3D fracture numerical model. The 3D spacing (total spacing) can be deduced from 2D spacing (normal spacing). The normal spacing of each fracture set can be obtained using the generated scan lines. On the other hand, the total spacing can be derived by dividing the normal spacing by cos(θ), where θ is the acute angle between the scan line and the line vertical to the average orientation of the fracture group. In 2D spacing, the frequency of each group can be calculated according to the method of Karzulovic & Goodman (1985). The number of fractures in the network can then be determined on the basis of the fracture density (Oda 1982).

Monte Carlo simulation of the fractures, which can be repeated until statistical agreement between the field data and simulated data is reached. A 140 m × 50 m × 50 m box, which is larger than the sampling window, was generated to accommodate the simulated fractures. Based on the fracture network parameters in the aforementioned steps, the Monte Carlo simulation can be used to create the fracture networks. RVE analysis can be conducted using a network that can statistically represent the actual rock mass. The parameters (i.e. orientation, trace length, spacing, and

*P*_{21}) in the field were confronted with those of the plane for the structural field mapping using the Kolmogorov–Smirnov (KS) test. A 5% significance level was used to determine whether the sampled data were consistent with those of the simulated model. The KS test for the simulated model should be continuously conducted before the significance level is reached. Subsequently, the model was used to represent the field rock mass. Generally, the features of the field fractures can be represented by their orientation, trace length, spacing, and intensity. Based on the KS test, the fracture network was stochastically similar to that in the field. Thus, the network can be used to simulate the actual fractures.

By using these procedures, fracture generations are achieved. As an example, the discontinuities of exploration tunnel 1 are expressed numerically in Table 2. Each discontinuity can be represented by a set of parameters, including coordinates (*x*, *y*, *z*), disc diameter, strike, and dip angle. The 3D fractures are visualized based on OpenGL (Fig. 4).

## Contingency table method and RVE calculation

The discontinuity information in the RVE is stochastically consistent with that of the overall rock mass (Zhou & Yu 1999). Therefore, based on the discontinuity information (i.e. density, diameter, strike, and dip angle), samples with RVE size are the same as each other as well as a sample of the overall rock mass. If the samples of the overall rock mass are stochastically the same, the sample size is equal to the RVE value. In this study, the overall rock mass was divided into samples with the same size, and the contingency table method was applied to check whether samples from the overall rock mass are the same. When the samples are the same considering discontinuity information, each sample includes sufficient information to represent the overall rock mass (i.e. the size of the sample is the RVE).

### Review of contingency table method

The contingency table method is evolved from the χ^{2} test, and is usually used to check whether samples are the same (Ge & Hu 2010; Kuroda *et al*. 2010). The χ^{2} value of tested samples is used to check the degree of difference of tested samples. If the difference between tested samples is significant, the χ^{2} value is large. When the χ^{2} value is larger than the critical value in the χ^{2} table, the samples are considered different; otherwise, the samples are considered the same.

Let (*S*_{1}, *S*_{2}, . . ., *S _{n}*) be samples from a population of unknown mean μ and of unknown variance σ

^{2}. The hypotheses of the contingence table method are described as follows.

*H*_{0} (null hypothesis): there is no significant difference between the tested samples.

*H*_{1} (alternative hypothesis): there is a significant difference between the tested samples.

The χ^{2} value can be calculated as

(1)

where *n* is the total number of analysed objects, *n _{ij}* represents the actual frequency of sample

*S*(

_{i}*i*= 1, 2, . . .,

*R*) in attribute

*j*(

*j*= 1, 2, . . .,

*C*),

*np*represents the theoretical frequency,

_{ij}*R*is the number of rows,

*C*is the number of columns,

*M*is the total number of column

_{j}*j*, and

*N*is the total number of row

_{i}*i*(Table 3).

The contingency table method may be unreliable, especially when some actual or theoretical frequencies are small (<5). This can be resolved by grouping the categories together. For example, if the disc number of a cube (1, 1, 1) is three, which is lower than five, this disc number can be added to that of another cube. The disc numbers of the two cubes are then considered the entity for the χ^{2} test. The value of χ^{2}_{(}* _{R}*−

_{1)×(}

*−*

_{C}_{1)},

_{α}in the χ

^{2}table can then be determined, where α is a significance level, which is always specified to be 0.05 or 0.1. The value should be compared with the calculated value χ

^{2}. If χ

^{2}≥ χ

^{2}

_{(}

*−*

_{R}_{1)×(}

*−*

_{C}_{1)},

_{α}, the null hypothesis

*H*

_{0}is rejected for the alternative hypothesis

*H*

_{1}. Otherwise the null hypothesis

*H*

_{0}is accepted.

## RVE calculation algorithm

The χ^{2} test was conducted by comparing the attributes of the objects (Table 3). For example, the overall rock mass can be divided into *n* samples, which represent the objects. Each sample contains fractures that are characterized by density, length, strike, and dip angle. Each characteristic is divided into several parts. For example, the strike range of the discontinuities (80°, 200°) can be divided into several parts, namely, (80°, 90°), (90°, 100°), …, (190°, 200°). The intervals of the characteristics are considered the attributes for the χ^{2} test. If the attributes of the *n* samples are the same considering the χ^{2} test, the samples are the same as the overall rock mass. The fracture characteristics in each interval can be counted (Table 4). If the size of the tested samples is the RVE, the number of each sample should be statistically the same in each interval.

Using the above-mentioned method, the RVE size can be calculated from the point of view of a contingency table. The computing steps are shown in Figure 5 and introduced as follows.

Parameter

*f*should be defined. Parameter*f*, which is not related to the RVE result, controls the program cycle. In the cycle,*f*can be equal to 0, 1, or −1.A cube, whose side length length is

*V*, should be determined. The rectangular box (rock mass) with size*l*×*w*×*h*should be divided by the cube. Therefore, the box will be divided into*n*parts as shown in the equation(2)

where ⌊ ⌋ rounds a number to the next smaller integer. The size of the box is 140 m × 50 m × 50 m. As shown in Figure 6, if the side length of the cube is 10 m, the box can be divided into 14 × 5 × 5 parts. Each part (cube) is represented by a code. For example, (1, 1, 3) indicates that this cube is the first one in the

*x*direction, the first in the*y*direction, and the third in the*z*direction.Assume that if a fracture disc centre is located in the cube (

*x*,*y*,*z*), the fracture is analysed in the cube no matter which one it touches. For every cube, the inner fractures should be extracted. All the fractures can be assigned to the corresponding cubes according the disc centre coordinates.For every cube, the range of the extracted fracture factors diameter, strike, and dip angle should be divided into no less than five equal intervals for each factor. The numbers of discs, diameters, strikes, and dip angles in different intervals can be calculated by a self-written program through Visual C++ (Table 4).

Calculate the χ

^{2}value. For a given cube size, the χ^{2}value can be calculated given by(3)

where

*n*is the number of cubes;*m*,_{d}*m*,_{dm}*m*and_{s}*m*represent the interval numbers into which disc numbers, diameters, strikes, and dip angles should be divided. For the numbers of discs,_{a}*m*is equal to unity;_{d}*m*,_{dm}*m*, and_{s}*m*can be arbitrary;_{a}*n*,_{dij}*n*,_{dmij}*n*, and_{sij}*n*represent the numbers of discs, diameters, strikes, and dip angles in the interval of (_{aij}*x*,_{j}*x*_{j}_{+1}) of cube*j*, respectively; and*n*, and_{dja}, n_{dmja}, n_{sja}*n*represent the average numbers of discs, diameters, strikes, and dip angles in the interval of (_{aja}*x*_{j},x_{j}_{+1}) of all cubes, respectively. The χ^{2}value is the sum of every divided interval of disc number, diameters, strikes, and dip angles in every cube.The program runs once. Because a value of zero is assigned to

*f*at the first step, the selection ‘*f*= 0’ is executed. If the calculated value is larger than the critical one in the χ^{2}table after the first cycle, a value of unity should be assigned to*f*(Fig. 5). A larger cube side length*V*+ 1 should be chosen to repeat Steps (3)–(5). Selection ‘*f*= 1’ runs automatically until the cube side length is*V′*, and the calculated value is smaller than the critical value; then the side length of the RVE is equal to*V*′. When the program has run once, if the calculated value is smaller than the critical value, a value of −1 should be assigned to*f*(Fig. 5), and*V*− 1 should be chosen to repeat Steps 3–5. Selection ‘*f*= −1’ runs automatically until the size is*V*′, and the calculated value is larger than the critical one; then the side length of RVE is equal to*V*′ + 1.

### RVE calculation based on practical engineering

The contingency table method may be unreliable, especially when some actual or theoretical frequencies are smaller than five. Therefore, if the number value of an arbitrary parameter for a cube in Table 4 is smaller than five, it should be added to the disc number value of another cube for the corresponding parameter until the sum is larger than five. The two combined cubes are then considered as entities for the χ^{2} test.

The χ^{2} values calculated by equation (3) are sensitive to the number of intervals *m _{dm}*,

*m*and

_{s}*m*. The larger the numbers, the larger the χ

_{a}^{2}value. However, only the difference between the calculated χ

^{2}value and the critical one can indicate the relationship between the tested samples. As an example, the RVE calculation results of exploration tunnel 1 are shown in Figure 7. The relative difference between the calculated value and critical value gradually decreases before the two lines in Figure 7 cross. The difference then increases after the crossing. According to the distributions of the discontinuity factors in different cubes, the first calculation cube size

*V*is set to 5 m.

*f*is set to zero, and the procedure runs once. The calculation result shows that the calculated value is larger than the critical one. Therefore,

*f*is altered to a value of unity, and the next larger size, 6 m, should be chosen to be computed, the cycle where

*f*= 1 runs. Until the side length

*V*is 8 m, the calculated value becomes smaller than the critical one. Then the distribution of discontinuities in each cube is the same, and the side length of RVE is 8 m. Therefore, the parameters of the sample with a side length equal to or larger than 8 m can represent those of the exploration tunnel 1. The equivalent continuum model can be used when the mesh size is equal to or larger than 8 m. Further analysis can be conducted based on the tested parameters and the selected numerical simulation.

Because the average trace lengths of different fracture sets may vary extremely, different fracture sets should be analysed separately to determine the relationship between the mean trace length and the RVE size. Table 5 shows that the side length of the RVE is approximately 2–8 times larger than the average trace length. If there is a high degree of similarity between the discontinuity factors (trace length, strike, and dip angle) in the analysed set, the multiple will be smaller; for instance, in set 4 in the exploration tunnel 1, the side length of the RVE is 1.8 times larger than the average trace length. Otherwise, the multiple will be larger; for instance, in set 2 in exploration tunnel 3, the side length of the RVE is 7.6 times larger than the average trace length. To be consistent with the traditional research method, the population of the discontinuities is also analysed. Table 5 shows that the RVE size of all discontinuities is closest to the RVE size of the set with the largest number of fractures.

The exploration tunnels in the Baihetan dam are characterized by inhomogeneous features, whereas each of the exploration tunnels is structurally homogeneous. Therefore, the test sample size or mesh size (for the continuum analysis) should be equal to or larger than the corresponding RVE values of the exploration tunnels.

## Fracture frequency calculation

To check the accuracy of the calculation, a parameter of the rock masses should be chosen to validate the RVE size. If the parameter value changes gradually with increasing size, and remains fixed when the size is larger than a critical value, then the critical size should be equal to the RVE. Based on the 3D fracture network, the parameter fracture frequency is chosen to calculate and validate the RVE size.

To build a relationship between the fracture frequency value and the RVE size, scan lines with different lengths are set through the computer program. Because the fractures are distributed stochastically, it makes sense to set numerous scan lines to acquire the mean fracture frequency value. In this study, the procedures for calculation of the fracture frequency in the *x* direction are as follows.

Determine the scan line group numbers

*n*and_{y}*n*in the_{z}*y*and*z*directions, respectively.*n*×_{y}*n*scan lines will be set in the_{z}*x*direction. The length of the rock mass in the*y*direction is 50 m. If*n*is equal to six, the_{y}*y*coordinate can be divided into six parts, namely, 0, 10, 20, 30, 40, and 50. Similarly, if the*z*coordinate can be divided into 0, 10, 20, 30, 40, and 50, 36 scan lines can be set along the*x*direction traversing (*x*_{0}, 0, 0), (*x*_{0}, 0, 10), … , (*x*_{0}, 50, 50), where*x*_{0}is an arbitrary*x*coordinate that the scan lines cross.Choose a scan line length

*l*. Each scan line in step (1) can then be divided into 140/*l*parts. If 140/*l*= 2, the fracture frequency value in (0, 70) and (70, 140) can be calculated. Through this step,*X*/_{t}*l*×*n*×_{y}*n*scan lines can be set._{z}Determine the intersections between a scan line and a fracture according to the equation

(4)

where (

*x*_{0},*y*_{0},*z*_{0}) represents the disc centre of a fracture,*y*and_{c}*z*represent the_{c}*y*and*z*coordinates as mentioned in Step (2), α represents the dip angle of the fracture, β represents the strike of the fracture, and (*A*,*B*,*C*) represents the normal vector of the fracture plane. Equation (4) has a solution if the fracture is crossed by the scan line.Calculate fracture frequency. For the entire fractures, different intersections along a scan line can be confirmed by equation (4). Fracture frequency can be confirmed by the intersection number divided by the scan line length.

Repeat Steps (3) and (4). All fracture frequency values of different scan lines can be calculated. Because the rock mass is inhomogeneous, the values could vary extremely. The statistical fracture frequency should be represented by the average value.

Change the scan line length

*l*, and repeat Steps (3)–(5). The fracture frequency value will change when the scan line length is increased, and size effect could be analysed.

Take exploration tunnel 1 as an example; the size of the studied rock mass is 140 m × 50 m × 50 m. *n _{y}* and

*n*are set to 10 and 10, respectively. The change of fracture frequency value with the increased scan line length is determined (Fig. 8).

_{z}If the fracture value changes with the increase of the scan line length and remains the same when the length is larger than a critical value, the critical length can be taken as the RVE size. As shown in Figure 8, the fracture frequency value decreases when the scan line length is increased. The longer the length, the smaller the change rate. When the scan line length is smaller than 8 m, fracture frequency decreases obviously. When the length is larger than 8 m, fracture frequency remains unchanged by and large. According to the definition, RVE size is equal to 8 m; this is consistent with the contingency table method to calculate RVE size (Table 5). Therefore, it is credible to validate RVE size through the contingency table method.

## Conclusion

This study is conducted based on 3D fracture numerical network modelling and stochastic mathematics (contingency table method). Through a 2D field investigation and procedures of the 3D simulation described above, all fractures are located, and the occurrences and sizes of all discontinuities are confirmed. This is useful in analysing issues in 3D space. The discontinuity factors (density, length, strike, and dip angle) are taken into account to determine the RVE size quantitatively. The contingency table method is applied in RVE size analysis. By calculating the numbers of disc, diameter, strike, and dip angle, the discontinuity characteristics are considered comprehensively. The test model is established, and RVE size is confirmed. The calculation method and conditions are described. It is concluded that the RVE size is 2–8 times larger than the average trace length. If there is a high degree of similarity between the discontinuity characteristics for the analysed discontinuities, the multiple will be smaller. A program is created, the fracture frequency is calculated, and the RVE size is validated based on the size effect of the fracture frequency. The RVE size determined by fracture frequency calculation is consistent with that found by the contingency table method. Therefore, the calculation of RVE size using the contingency table method is proved feasible and credible.

## Acknowledgments

This work was supported by the Natural Science Foundation of China (grant numbers 40872170 and 40902077), Doctoral Program Foundation of Higher Education of China (grant number 20090061110054), Jilin University’s 985 Project (grant number 450070021107) and Basic Research (grant number 200903202).

- © 2013 The Geological Society of London