The Specific Length of an Underground Tunnel and the Effects of Rock Block Characteristics on the Inflow Rate

Abstract

The specific length of a tunnel (STL) and a new analytical model for calculating the block surface area of the rock mass are introduced. First, a method for determining the appropriate length of a tunnel for a numerical simulation is described. The length is then used to examine the correlation between the inflow rate to the tunnel and the block volume, the block surface area, and the fracture intensity (P₃₂) through analytical and numerical modeling. The results indicate that the length of the tunnel should at least be equal to the least common multiple (LCM) of the apparent spacings of the joint sets at the wall of the tunnel to obtain the more reliable and immediate results for the inflow rate to a tunnel that is excavated in a fractured rock mass. A new analytical model was developed to calculate the block surface area and determine the essential joint set parameters, which include the dip, dip direction, and spacing. The determination of the rock block characteristics through numerical modeling requires considering the intact block for calculations. The results indicated that the inflow rate to the tunnel increased with an increase in fracture intensity and a decrease in block volume and surface area. The STL and the analytical model used for calculating the block surface area are validated through numerical simulations with 3DEC software version 7.0.

1. Introduction

Quality information on the inflow rate to an underground excavation is useful in a wide range of civil works. Taking into account the project requirements on the control of groundwater inflow during tunneling [1], assessment of the stability of the rock blocks at the wall of the tunnel [2], tunnel lining design [3], probable damage to the tunnel support structure and lowering the tunnel excavation rate [4], and settlement of the aboveground buildings [5,6] are some of the risks that relate directly to the rate of inflow to a tunnel and emphasize the significance of estimating the inflow rate. The inflow rate is also effectively a determining parameter for the cost of a civil and mining engineering project. Previous studies that aimed to estimate the inflow rate to a tunnel were developed through empirical, analytical, and numerical methods.

Many of the empirical equations that have been developed for estimating the inflow rate to a tunnel consider increasing depth, which results in the permeability of the rock mass decreasing and the hydraulic gradient at the wall of the tunnel increasing [7,8,9]. The developed empirical equations are based on data obtained from the double packer [10] or Lugeon tests [11] for the hydraulic conductivity of the rock mass. Several investigations have also been conducted for adjusting the relationships between the permeability of the rock mass and its geological indices [12,13], e.g., RQD, RMR, and GSI.

The proposed analytical models for estimating the inflow rate to tunnels were developed with the assumptions of substituting the fractured or porous media through equivalent homogeneous and isotropic formations. The application of various boundary conditions and solving the Laplace equation then allowed the development of a set of equations for the calculation of the inflow rate to a tunnel [14,15]. These equations could be applied to any tunnel excavated in a fractured rock mass using pre-determined values of its equivalent hydraulic conductivity using Packer or Lugeon tests. However, a minority of analytical methods that have been developed are not necessary for equivalent permeability [16], and, thus, conducting a field test is essential before the application of the analytical models.

A wide range of studies has also focused on the determination of the inflow rate to a tunnel using numerical simulations. In this regard, the effects of geometrical characteristics of the rock mass, i.e., joint spacing, aperture, and orientation [17], as well as the overburden load/stress [18] on the inflow rate to the tunnel have been considered using a two-dimensional [19,20] or three-dimensional [21,22] model. Although the inflow rate increases by increasing the fracture aperture and decreasing the spacing, the effects of joint set orientation on inflow rate remains unclear. Furthermore, by increasing normal stress on a fracture plane, the fracture aperture decreases and any shear displacement results in a lower increase in fracture transmissivity and inflow rate. Among the above-mentioned parameters, the effects of the block characteristics on the inflow rate to a tunnel have not been well studied. In addition, no criterion has been determined for the length of the tunnel to be considered in the numerical simulations of the inflow rate.

In this study, the relationship between the inflow rate to a tunnel and rock block characteristics including block volume, block surface area, and volumetric fracture intensity (P₃₂) have been investigated using Itasca 3DEC version 7.00 software. 3DEC is a three-dimensional command-driven numerical program based on the distinct element method for discontinuum modeling. It is used to evaluate the response of the fractured media to the static or dynamic forces, and to carry out hydromechanical coupling simulations. Because the actual geometry of a rock mass is usually too complex for being simulated, similar to all other numerical models for this purpose, several simplifying assumptions are considered in this study. For a rock mass that includes less than three joint sets, randomly generated joints define the block volume; if more than three joint sets are present, only the three prominent sets are considered to determine the block volume [23]. Thus, it is assumed that the rock mass includes three joint sets with different orientations and spacings, but each set has a fixed value of these parameters. The length of the tunnel is defined by a newly proposed analytical method, called the Specific Tunnel Length (STL), which is adapted to a three-dimensional model as used in the current study, and yields more reliable results. The STL represents the minimal length of a tunnel that is representative of its entire length regarding its hydraulic characteristics. The STL is introduced as an important parameter because in all of the 3D numerical simulations, the length of the tunnel is a critical input parameter that strongly affects the value of the inflow rate. The length of the tunnel smaller than the STL value has been proven to mostly produce an error in the value of the inflow rate because it is not representative of the complete hydraulic conditions observed along the tunnel. In contrast, the tunnel length larger than the STL will effectively increase the processing time of the numerical modeling without bringing more information. The optimum length of the tunnel for the numerical simulations in relation to hydraulic characteristics has been demonstrated to be equal to STL. The existence of the STL is supported by numerical simulations using 3DEC version 7.00 software. The last section of this paper focuses on the evaluation of the relationship between the tunnel inflow rate and three other parameters, which include block volume, block surface area, and volumetric fracture intensity. For this purpose, a series of numerical simulations have been conducted based on the already introduced specific length of the tunnel. In addition, an analytical model is developed for the calculation of the block surface area that is created by the intersection of three joint sets. Having the values of the inflow rate on the one hand and the numerically and analytically calculated block volume, block surface area, and P₃₂ on the other hand, allows the investigation the effects of each parameter on the inflow rate.

2. Specific Tunnel Length

In all numerical models, the first and the most important concern is to determine the minimum reliable size of the model that is representative of its unlimited sizes. For the rock mass, this size was addressed in the literature as Representative Elementary Volume (REV), which could be defined for the mechanical properties [24] or hydraulic characteristics [25]. As a general rule, to gain reliable results from the numerical simulations, the model has to be equal or greater than the REV of the rock mass. However, in the case of excavating a tunnel in a fractured rock mass, the predefined REV is not appropriate for use in the calculation of the inflow rate to the tunnel. Thus, for a better understanding of the problem, Figure 1 shows a rock mass that includes three joint sets. It is assumed that the intact rock is impermeable and the fluid flows only through the fractures. Figure 1A shows the predefined dimension of the REV and the horizontal tunnels excavated in two different directions in Figure 1B,C.

Based on Figure 1, the tunnel can be excavated in any direction in a REV of a rock mass. Despite the requirement that the REV of a specific formation should yield a unique hydraulic behavior, Figure 2 shows that the number and total length of the fractures that cross the tunnel vary by changing their direction. The intact rock is presumed to be impermeable and if a fixed level of the water table is set for both models of Figure 2c,d, by variation of the trace length of the fractures at the wall of the tunnel (Figure 2a,b), therefore, the inflow rate to the tunnel will differ for the same length of the tunnel.

The differences observed in Figure 1 and Figure 2 illustrate the necessity of revising the concept of the hydraulic REV for the cases that a tunnel excavated in the rock mass. For this purpose, the hydraulic REV of a tunnel can be defined as the shortest length of the tunnel for which the inflow rate is representative of a unit inflow rate obtained from all hydraulic expressions of the fracture sets. This length is called the specific tunnel length or STL in this paper. In other words, the STL is the shortest length of the tunnel that integrates the flow rate contribution of the entire fracture sets (inflow rate per each meter of the tunnel length). Multiplying the inflow rate provided by the fractures along with the STL value by the total length of the tunnel divided by the STL value yields the total inflow rate through the tunnel. For a specific formation, different STL values exist for each direction of the tunnel and, accordingly, the length of the tunnel in the numerical simulations should be equal to an integer multiplication of the STL to attain reliable results.

2.1. Determination of the STL

To determine the STL, the joint sets are assumed to be persistent, the intact rock is impermeable, and fluid flows only through the fractures. Furthermore, the values of the spacing, aperture, and orientation of the joint sets are fixed. Figure 3 shows the cross-section of the tunnel excavated in a rock mass having three joint sets. The apparent spacing that can be seen at the wall of the tunnel will always be equal or greater than the true spacing and the difference between the apparent and true spacings depends on the angle between the joint set and tunnel directions. As a reminder, true spacing refers to the shortest distance between two adjacent fractures of a joint set and apparent spacing refers to the observed spacing in any direction that is not essentially perpendicular to the plane of the joint set.

The fractures are assumed to be the sole channels for water flowing into the tunnel, and, thus, to define the STL, the minimum length of the tunnel whose arrangement of the joint sets repeats at any multiple of STL (i.e., 1 STL, 2 STL, …, n STL) should be defined. In this regard, Figure 4 illustrates the locations of the traces of each joint set at the wall of the tunnel. For simplicity, it is assumed that at the first point (zero points), the traces of all joint sets are overlapped. Thus, the same results for the value of the STL might be observed if the traces do not overlap at the initial point. As the apparent spacing of set 1 is 1, then it will be repeated at each multiple of 1, and similarly, joint sets 2 and 3 are repeated at each multiple of 1.5 and 2, respectively. According to the arrangement of the traces, all traces overlapped again at points 6, 12, and 18. Furthermore, the sequences of the traces from point 1 to 6 are repeated in the ranges between 7 to 12 and 13 to 18.

By considering the repetition points (6, 12, and 18), the first repetition point is found as the least common multiple (LCM) of the apparent spacings of joint sets at the wall of the tunnel, as in Figure 4, 6 is the LCM of 1, 1.5, and 2. Furthermore, the second and third repetition points are the integer multiples of the LCM. Therefore, the STL is equal to the LCM of the apparent spacings of joint sets at the wall of the tunnel. The existence of the STL is demonstrated graphically in Figure 5. In this figure, the rock mass includes three joint sets with apparent spacings equal to 0.1, 1, and 6 m. It is also evident that in each case, the sequence of the fractures in all spans along the tunnel direction with length equal to the STL is repeated in all next spans with the length equal to STL.

Before specifying the STL, the apparent spacings of each joint set at the wall of the tunnel should be determined. Based on Figure 6, the cross-section of the intersection of the cylinder (tunnel) and fracture plane is a circle perpendicular and ellipse in an oblique cross. The apparent spacing will increase accordingly by deviating from the perpendicular cross.

The true and apparent spacings of set i are named as S_i and Sp_i, respectively, and the angle between normal to joint set i and tunnel direction is called Ɵ_i. Equation (1) expresses the relationship between apparent and true spacings.

$${\mathrm{Sp}}_{\mathrm{i}}=\frac{{\mathrm{S}}_{\mathrm{i}}}{\cos {\mathsf{Ɵ}}_{\mathrm{i}}}$$

(1)

The angle between normal to joint set and tunnel direction Ɵ_i could be specified by identifying the dip/dip direction of the joint set and direction of the tunnel. Components of the unit normal vector to the joint set i (direction cosines of V_i in Figure 6) can be specified according to Equation (2).

$$\begin{gathered} {\mathrm{l}}_{\mathrm{i}}=\sin \left({\mathrm{DD}}_{\mathrm{i}}\right)\sin \left({\mathrm{D}}_{\mathrm{i}}\right) \\ {\mathrm{m}}_{\mathrm{i}}=\cos \left({\mathrm{DD}}_{\mathrm{i}}\right)\sin \left({\mathrm{D}}_{\mathrm{i}}\right) \\ {\mathrm{n}}_{\mathrm{i}}=-\cos \left({\mathrm{D}}_{\mathrm{i}}\right) \end{gathered}$$

(2)

where D_i and DD_i are the dip and dip direction of set i, respectively. Therefore,

$${ \overrightarrow{\mathrm{V}}}_{\mathrm{i}}=\left({\mathrm{l}}_{\mathrm{i}}\overrightarrow{\mathrm{i}},{ \mathrm{m}}_{\mathrm{i}}\overrightarrow{\mathrm{j}},{ \mathrm{n}}_{\mathrm{i}}\overrightarrow{\mathrm{k}}\right)$$

(3)

Furthermore, by identifying trend (T_T) and plunge (P_T) of the tunnel centerline, the components of the unit vector of the tunnel direction (V_T) are shown as Equation (4):

$$\begin{gathered} {\mathrm{l}}_{\mathrm{T}}=\sin \left({\mathrm{T}}_{\mathrm{T}}\right)\cos \left({\mathrm{P}}_{\mathrm{T}}\right) \\ {\mathrm{m}}_{\mathrm{T}}=\cos \left({\mathrm{T}}_{\mathrm{T}}\right)\cos \left({\mathrm{P}}_{\mathrm{T}}\right) \\ {\mathrm{n}}_{\mathrm{T}}=-\sin \left({\mathrm{P}}_{\mathrm{T}}\right) \end{gathered}$$

(4)

The same could be defined for the unit vector of the tunnel direction:

$${ \overrightarrow{\mathrm{V}}}_{\mathrm{T}}=\left({\mathrm{l}}_{\mathrm{T}}\overrightarrow{\mathrm{i}},{ \mathrm{m}}_{\mathrm{T}}\overrightarrow{\mathrm{j}},{ \mathrm{n}}_{\mathrm{T}}\overrightarrow{\mathrm{k}}\right)$$

(5)

Meanwhile, the angle Ɵ_i will be specified using the inner product of vectors V_T and V_i as

$${ \overrightarrow{\mathrm{V}}}_{\mathrm{T}}.{ \overrightarrow{\mathrm{V}}}_{\mathrm{i}}=\left|{ \overrightarrow{\mathrm{V}}}_{\mathrm{T}}\right|\times \left|{ \overrightarrow{\mathrm{V}}}_{\mathrm{i}}\right|\times \cos {\mathsf{Ɵ}}_{\mathrm{i}} \stackrel{\left|{ \overrightarrow{\mathrm{V}}}_{\mathrm{T}}\right|=\left|{ \overrightarrow{\mathrm{V}}}_{\mathrm{i}}\right|=1}{\to } \cos {\mathsf{Ɵ}}_{\mathrm{i}}={ \overrightarrow{\mathrm{V}}}_{\mathrm{T}}.{ \overrightarrow{\mathrm{V}}}_{\mathrm{i}}$$

(6)

Combining Equations (2)–(5) with Equation (6) yields the following:

$$\begin{array}{cc}\hfill \cos {\mathsf{Ɵ}}_{\mathrm{i}}={ \overrightarrow{\mathrm{V}}}_{\mathrm{T}}.{ \overrightarrow{\mathrm{V}}}_{\mathrm{i}}=& \left({\mathrm{l}}_{\mathrm{T}}\overrightarrow{\mathrm{i}},{ \mathrm{m}}_{\mathrm{T}}\overrightarrow{\mathrm{j}},{ \mathrm{n}}_{\mathrm{T}}\overrightarrow{\mathrm{k}}\right). \left({\mathrm{l}}_{\mathrm{i}}\overrightarrow{\mathrm{i}},{ \mathrm{m}}_{\mathrm{i}}\overrightarrow{\mathrm{j}},{ \mathrm{n}}_{\mathrm{i}}\overrightarrow{\mathrm{k}}\right)={ \mathrm{l}}_{\mathrm{T}}{\mathrm{l}}_{\mathrm{i}}+{\mathrm{m}}_{\mathrm{T}}{\mathrm{m}}_{\mathrm{i}}+{\mathrm{n}}_{\mathrm{T}}{\mathrm{n}}_{\mathrm{i}}\hfill \\ & =\left[\sin \left({\mathrm{T}}_{\mathrm{T}}\right)\cos \left({\mathrm{P}}_{\mathrm{T}}\right)\times \sin \left({\mathrm{DD}}_{\mathrm{i}}\right)\sin \left({\mathrm{D}}_{\mathrm{i}}\right)\right]\hfill \\ & +\left[\cos \left({\mathrm{T}}_{\mathrm{T}}\right)\cos \left({\mathrm{P}}_{\mathrm{T}}\right)\times \cos \left({\mathrm{DD}}_{\mathrm{i}}\right)\sin \left({\mathrm{D}}_{\mathrm{i}}\right)\right]+\left[\sin \left({\mathrm{P}}_{\mathrm{T}}\right)\times \cos \left({\mathrm{D}}_{\mathrm{i}}\right)\right]\hfill \end{array}$$

(7)

Finally, for a rock mass that includes three joint sets, the STL could be specified using Equation (8):

$$\mathrm{STL}= \mathrm{LCM}\left({\mathrm{Sp}}_1,{\mathrm{Sp}}_2,{\mathrm{Sp}}_3\right)=\mathrm{LCM}\left(\frac{{\mathrm{S}}_1}{\cos {\mathsf{Ɵ}}_1},\frac{{\mathrm{S}}_2}{\cos {\mathsf{Ɵ}}_2},\frac{{\mathrm{S}}_3}{\cos {\mathsf{Ɵ}}_3}\right)$$

(8)

where Sp₁, Sp_2, and Sp₃ are the apparent spacings and S₁, S_2, and S₃ are the true spacings of joint sets 1, 2, and 3, respectively. Furthermore, Ɵ₁, Ɵ_2, and Ɵ₃ are the angles between tunnel direction and normal to each joint set.

2.2. Validation of the STL

The existence of the STL in a rock mass that includes three joint sets is investigated and validated using the numerical simulation by 3DEC version 7.00 software. It should be emphasized that the rock mass could contain more than three joint sets, and, in this case, the same method is applicable for defining the STL. Numerical models with various values of the apparent spacings, Sp_i in Equation (1), were prepared to evaluate the existence of the STL. Afterward, the inflow rates to the tunnels with a length of half STL, 1 STL, 1.5 STL, 2 STL, and 3 STL were calculated. Figure 7 shows the models used for the numerical simulations. The x and z dimensions of the models are kept constant and only the y dimension varies to create the tunnels with different lengths.

Figure 8 illustrates the boundary conditions and the applied method for calculations of the inflow rate to the tunnel. As shown in Figure 8a, a fixed level of the water table is adjusted at the top of the model, and the pore pressure at the sides of the model parallel to the y-z plane is set to be equal to its initial values during the numerical calculations. The pore pressure at the wall of the tunnel is set to be zero to simulate the inflow to the tunnel, as shown in Figure 8b. In addition to the hydraulic boundary conditions described above, a set of mechanical and displacement boundary conditions have been applied to the model. In this regard, the displacement at all boundaries of the model is considered equal to zero. In addition, no stress, including gravitational stress, is applied to the model, and as a result, the hydraulic aperture is fixed in all parts of the model and is equal to its initial value. Because the inflow rate decreases by moving from the tunnel wall toward the boundaries of the model, the area close to the wall of the tunnel is selected for the calculation of the inflow rate as shown in Figure 8c. Figure 8 shows that the spatial discretization must increase where hydraulic gradients are higher, as is the case when we get closer to the tunnel to decrease the time of the calculation and to maintain the accuracy of the results.

In Figure 8, the flow plane is the planar polygon corresponding to face-to-face contact between solid blocks, the flow plane zone is a triangular discretization element of the flow plane, and the flowknot is the vertices of a flow plane zone that generally correspond to a sub-contact between solid blocks. To ensure that the model is in steady-state condition, the variation of the pore pressure with cycling steps is recorded at four points around the tunnel circumference (the nearest flowknot to the wall of the tunnel at 3, 6, 9, and 12 o’clock positions). The model reaches the steady-state condition when the pore pressure at all points does not vary anymore.

Based on the results of the numerical simulations, it is observed that in all cases, the average inflow rate (inflow rate per meter of tunnel length) in 1STL, 2STL, and 3STL is equal and differs from what was calculated for 0.5STL and 1.5STL. In this study, the direction of the tunnel is always in S-N and the hydraulic apertures of all joint sets are set to be equal. However, it could be simply proved that the different hydraulic apertures of joint sets and tunnel direction do not affect the length and existence of the STL.

For this purpose, the STL is evaluated for five different cases according to

Table 1. Joint set characteristics considered in the numerical models for assessment of the STL (DIP, DD, S, and Ap are the dip, dip direction, true spacings, and hydraulic aperture of the joint sets, respectively).

	Joint Set 1				Joint Set 2				Joint Set 3				Tunnel Diameter (m)	Water Head (m)
Case No.	DIP	DD	S (cm)	Ap (cm)	DIP	DD	S (cm)	Ap (cm)	DIP	DD	S (cm)	Ap (cm)
1	22	25	205.2	1 × 10−4	21	342	205.2	1 × 10−4	80	0	9.8	1 × 10−4	1	5
2	81	5	39.4	1 × 10−4	73	25	519.6	1 × 10−4	82	354	98.5	1 × 10−4
3	86	9	590.9	1 × 10−4	67	340	8.7	1 × 10−4	80	0	590.9	1 × 10−4
4	84	8	590.9	1 × 10−4	86	351	590.9	1 × 10−4	80	0	590.9	1 × 10−4
5	61	8	8.7	1 × 10−4	21	342	13.7	1 × 10−4	80	2	9.8	1 × 10−4

.

The hydraulic aperture and tunnel diameter in all cases were set to be equal to 1 × 10⁻⁴ (cm) and 1 m, respectively, to simplify the models and speed up the time of calculations. The level of the water table is fixed to 5 m above the centerline of the tunnel in all cases. The apparent spacings relevant to Table 1 and unit inflow rate to the tunnel for each model are listed in

Table 2. The apparent spacings (Sp1, Sp2, and Sp3) at the wall of the tunnel for each model, STL, and the value of the inflow rate to the tunnel for the tunnel length of 0.5 STL, 1 STL, 1.5 STL, 2 STL, and 3 STL.

					Inflow Rate (m3/s) per Meter of Tunnel Length
Case No.	Sp1 (m)	Sp2 (m)	Sp3 (m)	STL (m)	0.5 STL	1.0 STL	1.5 STL	2.0 STL	3.0 STL
1	6	6	0.1	6	2.15 × 10−8	9.35 × 10−10	1.02 × 10−8	9.35 × 10−10	9.35 × 10−10
2	0.4	6	1	6	7.15 × 10−10	3.58 × 10−10	9.41 × 10−10	3.58 × 10−10	3.58 × 10−10
3	6	0.1	6	6	1.28 × 10−8	1.51 × 10−8	1.36 × 10−8	1.51 × 10−8	1.52 × 10−8
4	6	6	6	6	5.96 × 10−10	8.98 × 10−10	1.50 × 10−9	8.98 × 10−10	8.98 × 10−10
5	0.1	0.4	0.1	0.4	7.15 × 10−10	8.28 × 10−10	5.96 × 10−9	8.28 × 10−10	5.96 × 10−10

.

Figure 9 is drawn based on the data provided in Table 2 and illustrates the equal inflow rate in 1 STL, 2 STL, and 3 STL (and generally n STL where n is an integer) as well as the different values for 0.5 STL and 1.5 STL. Regarding this figure and the spacings listed in Table 2, it is worth considering that the case numbers 1 to 5 are created by different arrangements of three apparent joint spacings, including 0.1, 0.4, and 6 m. It shows that the various values of the inflow rates may be obtained in different orientations of the tunnel. Furthermore, Figure 9 confirms that the value of the average inflow rate to each integer multiplication of the STL is equal to the average inflow rate to 1 STL and, accordingly, it is equal to the infinite length of the tunnel. The deviation of the inflow rate for 0.5 STL and 1.5 STL with 1 STL is also decreased, indicating that the deviation will tend to zero by increasing the length of the tunnel when the length is not equal to the n × STL.

3. Evaluation of the Effects of the Block Characteristics on the Inflow Rate

Using the concept of STL, which is explained and validated in Section 2, a series of numerical simulations are designed to study the relationship between geometrical characteristics of the rock block, i.e., rock block volume, block surface area, and volumetric fracture intensity on the one hand and the inflow rate to the tunnel on the other hand. This section aims to define the parameter with the most important impact on the inflow rate to the tunnel. That parameter could be an efficient representation of the geometrical characteristics of the rock mass, e.g., spacing, dip, and dip direction of the joint sets [26,27].

In this regard, numerical models comprised three joint sets with various spacings, dip, and dip directions. The tunnel is always excavated in the N-S direction and a fixed level of the water table is applied to the model. Because the orientations and spacings of the joint sets vary, the trace length of the fractures at the wall of the tunnel change, and the inflow rate will vary. Using a FISH function in the 3DEC commands, the block characteristics, as well as the inflow rate to the tunnel, are calculated numerically. The block volume and surface are calculated analytically using a formerly developed analytical equation for block volume specification [23,26,28] and an analytical equation developed in this study for the block surface computation. Rather than a comparison of the results of the numerical and analytical block volume and block surfaces, a discussion on the reliability of the block volume and surface calculation is presented at the end of this section.

3.1. Analytical Calculation of the Block Volume and Surface Area

For the evaluation of the effects of one parameter on the output of a model, it is mostly assumed that all variables are in their ideal conditions. For the case of a fractured rock mass, it is assumed that each joint set has a constant value of the spacing, orientation, and aperture, and that the persistence of the fractures is always 1, i.e., all the fractures cross the boundaries of the model. In the current study, it is assumed that the rock mass includes three persistent joint sets with fixed values of the spacings, dip, and dip direction for each joint set. In addition, the model does not contain random joints and the blocks are formed by the joint sets only. Based on the previous studies and the above-mentioned assumptions [23,26,28], the volume of the block created by the cross of three joint sets could be calculated using Equation (9):

$${\mathrm{V}}_{\mathrm{b}}^{\mathrm{A}}=\frac{{\mathrm{S}}_1\times {\mathrm{S}}_2\times {\mathrm{S}}_3}{\sin {\mathsf{\gamma }}_1\times \sin {\mathsf{\gamma }}_2\times \sin {\mathsf{\gamma }}_3}$$

(9)

where ${\mathrm{V}}_{\mathrm{b}}^{\mathrm{A}}$ is the analytically calculated block volume; S₁, S₂, and S₃ are the true spacings of joint sets; and γ₁, γ₂, and γ₃ are the angles between joint sets. For calculation of the surface area of a rock block created by three joint sets, an equation is mathematically developed using the methodology described briefly in Figure 10.

Except for the blocks cut by the boundaries of the model, the other blocks’ geometries are identical. An intact block is selected for model development as illustrated in Figure 10b. Assuming that each edge of the block could be shown by a vector (A, B, and C in Figure 10c), the surface area of each face of the block can be determined by the cross product of each pair of vectors A, B, and C. The surface area of each frontside faces of the block are equal because all fractures in a joint set are parallel. Accordingly, the surface area of the block could be calculated by Equation (10):

$${\mathrm{S}}_{\mathrm{b}}^{\mathrm{A}}=2\times \left(\left|\stackrel{\rightharpoonup }{\mathrm{A}}\times \stackrel{\rightharpoonup }{\mathrm{B}}\right|+\left|\stackrel{\rightharpoonup }{\mathrm{A}}\times \stackrel{\rightharpoonup }{\mathrm{C}}\right|+\left|\stackrel{\rightharpoonup }{\mathrm{B}}\times \stackrel{\rightharpoonup }{\mathrm{C}}\right| \right)$$

(10)

where ${\mathrm{S}}_{\mathrm{b}}^{\mathrm{A}}$ is the analytically calculated block surface and A, B, and C are the edge vectors of the intact block. Each edge vector has a direction and a magnitude. The direction of a vector could be specified by a unit vector, and for this specific case, it could be defined by Equation (11):

$$\begin{gathered} \stackrel{\rightharpoonup }{{\mathrm{u}}_{\mathrm{A}}}=\stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}2}}\times \stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}3}} \\ \stackrel{\rightharpoonup }{{\mathrm{u}}_{\mathrm{B}}}=\stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}1}}\times \stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}3}} \\ \stackrel{\rightharpoonup }{{\mathrm{u}}_{\mathrm{C}}}=\stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}1}}\times \stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}2}} \end{gathered}$$

(11)

where N_J1, N_J2, and N_J3 are the normal to joint set J₁, J₂, and J₃, and u_A, u_B, and u_C are the unit vectors of A, B, and C, respectively. On the other hand, by considering Figure 11, the magnitude of edge vectors is determined by Equation (12).

$$\begin{gathered} \left|\stackrel{\rightharpoonup }{\mathrm{A}}\right|=\frac{{\mathrm{S}}_1}{\cos {\mathsf{Ɵ}}_1} \\ \left|\stackrel{\rightharpoonup }{\mathrm{B}}\right|=\frac{{\mathrm{S}}_2}{\cos {\mathsf{Ɵ}}_2} \\ \left|\stackrel{\rightharpoonup }{\mathrm{C}}\right|=\frac{{\mathrm{S}}_3}{\cos {\mathsf{Ɵ}}_3} \end{gathered}$$

(12)

where ǀAǀ, ǀBǀ, and ǀCǀ are the magnitude of each edge vectors, S₁, S₂, and S₃ are the spacings, and Ɵ₁, Ɵ₂, and Ɵ₃ are the angles between normal to joint sets and direction of the edge vectors. These parameters are illustrated in Figure 11. The directions of the spacings S₁, S₂, and S₃ are parallel to N_J1, N_J2, and N_J3, respectively.

Based on Equation (12), to specify the magnitude of the edge vectors, Ɵ_i should be determined. According to Figure 11, Ɵ_i can be calculated by Equation (13):

$$\begin{gathered} \stackrel{\rightharpoonup }{{\mathrm{u}}_{\mathrm{A}}}.\stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}1}}=\left|\stackrel{\rightharpoonup }{{\mathrm{u}}_{\mathrm{A}}}\right|\times \left|\stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}1}}\right|\times \cos {\mathsf{Ɵ}}_1 \\ \stackrel{\rightharpoonup }{{\mathrm{u}}_{\mathrm{B}}}.\stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}2}}=\left|\stackrel{\rightharpoonup }{{\mathrm{u}}_{\mathrm{B}}}\right|\times \left|\stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}2}}\right|\times \cos {\mathsf{Ɵ}}_2 \\ \stackrel{\rightharpoonup }{{\mathrm{u}}_{\mathrm{C}}}.\stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}3}}=\left|\stackrel{\rightharpoonup }{{\mathrm{u}}_{\mathrm{C}}}\right|\times \left|\stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}3}}\right|\times \cos {\mathsf{Ɵ}}_3 \end{gathered}$$

(13)

Given that u_A, u_B, u_C, and N_J1, N_J2, and N_J3 are unit vectors and their absolute values are 1, by combining Equations (11) and (13), we can obtain the following:

$$\begin{gathered} \cos {\mathsf{Ɵ}}_1=\left(\stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}2}}\times \stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}3}}\right). \stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}1}} \\ \cos {\mathsf{Ɵ}}_2=\left(\stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}1}}\times \stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}3}}\right). \stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}2}}. \\ \cos {\mathsf{Ɵ}}_3=\left(\stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}1}}\times \stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}2}}\right). \stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}3}} \end{gathered}$$

(14)

Considering Equations (11), (12), and (14), vectors A, B, and C could be specified by Equation (15).

$$\begin{gathered} \stackrel{\rightharpoonup }{\mathrm{A}}=\left|\stackrel{\rightharpoonup }{\mathrm{A}}\right|\times \stackrel{\rightharpoonup }{{\mathrm{u}}_{\mathrm{A}}}=\frac{{\mathrm{S}}_1}{\left(\stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}2}}\times \stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}3}}\right). \stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}1}}}\times \left(\stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}2}}\times \stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}3}}\right) \\ \stackrel{\rightharpoonup }{\mathrm{B}}=\left|\stackrel{\rightharpoonup }{\mathrm{B}}\right|\times \stackrel{\rightharpoonup }{{\mathrm{u}}_{\mathrm{B}}}=\frac{{\mathrm{S}}_2}{\left(\stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}1}}\times \stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}3}}\right). \stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}2}}}\times \left(\stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}1}}\times \stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}3}}\right) \\ \stackrel{\rightharpoonup }{\mathrm{C}}=\left|\stackrel{\rightharpoonup }{\mathrm{C}}\right|\times \stackrel{\rightharpoonup }{{\mathrm{u}}_{\mathrm{C}}}=\frac{{\mathrm{S}}_3}{\left(\stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}1}}\times \stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}2}}\right). \stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}3}}}\times \left(\stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}1}}\times \stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}2}}\right) \end{gathered}$$

(15)

It should be mentioned that the unit normal vector to joint sets (N_J1, N_J2, and N_J3) could be determined knowing the dip and dip direction of the joint sets, according to Equation (16)

$$\begin{gathered} \stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}1}}=\left(\sin {\mathrm{DD}}_1.\sin {\mathrm{D}}_1\right) \stackrel{\rightharpoonup }{\mathrm{i}}, \left(\cos {\mathrm{DD}}_1.\sin {\mathrm{D}}_1\right) \stackrel{\rightharpoonup }{\mathrm{j}}, \left(-\cos {\mathrm{D}}_1\right) \stackrel{\rightharpoonup }{\mathrm{k}} \\ \stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}2}}=\left(\sin {\mathrm{DD}}_2.\sin {\mathrm{D}}_2\right) \stackrel{\rightharpoonup }{\mathrm{i}}, \left(\cos {\mathrm{DD}}_2.\sin {\mathrm{D}}_2\right) \stackrel{\rightharpoonup }{\mathrm{j}}, \left(-\cos {\mathrm{D}}_2\right) \stackrel{\rightharpoonup }{\mathrm{k}} \\ \stackrel{\rightharpoonup }{{\mathrm{N}}_{\mathrm{J}3}}=\left(\sin {\mathrm{DD}}_3.\sin {\mathrm{D}}_3\right) \stackrel{\rightharpoonup }{\mathrm{i}}, \left(\cos {\mathrm{DD}}_3.\sin {\mathrm{D}}_3\right) \stackrel{\rightharpoonup }{\mathrm{j}}, \left(-\cos {\mathrm{D}}_3\right) \stackrel{\rightharpoonup }{\mathrm{k}} \end{gathered}$$

(16)

where DD₁, DD₂, and DD₃ are the dip directions and D₁, D₂, and D₃ are dip of the joint set 1, 2, and 3, respectively. Finally, the surface area of the blocks could be specified by incorporating Equations (16) and (15) into Equation (10). For calculation of the block surface using Equation (10), an excel sheet is prepared and attached to this paper in Supplementary Materials section.

3.2. Validation and Comparison of the Analytical Models

As it is explained in Section 3.1, the analytical method for calculation of the block volume was previously developed by Equation (9) and a new analytical model is developed in this article for calculation of the block surface according to Equation (10). In this section, the outputs of the analytical models for 13 cases that are listed in

Table 3. Characteristics of the discontinuities used for comparison of analytical and numerical calculation of the block volume and block surface and also assessing the effects of the block characteristics on the inflow rate to the tunnel.

	Joint Set 1				Joint Set 2				Joint Set 3
Case	DIP	DD	Spacing (m)	Ap (m)	DIP	DD	Spacing (m)	Ap (m)	DIP	DD	Spacing (m)	Ap (m)
1	23	30	0.34	2 × 10−5	20	10	2.05	5 × 10−7	27	320	0.14	2 × 10−5
2	90	350	0.39	2 × 10−5	20	352	0.14	2 × 10−5	90	10	0.39	5 × 10−6
3	73	25	0.35	5 × 10−6	61	352	0.35	2 × 10−5	84	8	5.91	5 × 10−6
4	22	25	0.34	5 × 10−7	90	350	0.98	2 × 10−5	60	0	0.35	2 × 10−5
5	90	350	0.39	5 × 10−6	90	10	0.39	5 × 10−7	32	50	2.05	2 × 10−5
6	90	350	5.91	5 × 10−7	90	70	0.14	5 × 10−6	90	30	5.2	5 × 10−6
7	20	350	0.34	2 × 10−5	22	25	1.37	2 × 10−5	61	354	5.2	5 × 10−6
8	54	65	0.34	2 × 10−5	90	350	0.98	5 × 10−6	90	30	0.35	5 × 10−6
9	80	0	5.91	5 × 10−6	90	70	0.14	5 × 10−7	66	292	0.34	5 × 10−6
10	60	0	3.46	2 × 10−5	67	340	5.2	5 × 10−6	90	30	0.87	5 × 10−6
11	90	30	3.46	5 × 10−6	54	295	1.37	5 × 10−6	80	0	3.94	5 × 10−7
12	90	70	0.34	2 × 10−5	43	300	1.37	5 × 10−7	73	25	5.2	5 × 10−7
13	43	60	0.14	5 × 10−7	73	335	5.2	5 × 10−7	90	10	0.39	2 × 10−5

are compared with the results of the numerical simulations using the 3DEC version 7.0.

Based on the data listed in Table 3, a numerical model was prepared for the calculation of the representative block volume and surface for each case. For this reason, the volume and surface of the intact block, as illustrated in Figure 10c, should be specified. The most important point is to be certain that at least one intact block exists in the model. On the one hand, this block is representative of all blocks of the model when the dimension of the model is unlimited and on the other hand, the intact block is created and its number increase by extending the size of the numerical model. From now on, the RBL, RBLV, and RBLS will be used in this article for referring to the representative block, representative block volume, and representative block surface, respectively. It should also be emphasized that the RBL is the intact block shown in Figure 10c. In

Table 4. Comparison of the analytical and numerical calculation of the block volume (RBLV) and block surface (RBLS) for the cases of Table 3. N/A means that the parameter is unmeasurable.

	Block Surface (RBLS) (m2)		Block Volume (RBLV) (m3)
Case	Analytical (This Study)	Numerical (3DEC)	Analytical [23]	Numerical (3DEC)
1	320.2	320.2	4.31	15.14
2	1.63	1.63	0.07	0.07
3	54.15	54.16	8.55	4.6
4	49.73	49.71	0.36	3.65
5	12.08	12.07	1.03	1.08
6	N/A	N/A	10.57	N/A
7	131.54	131.53	24.90	17.02
8	4.21	4.20	0.25	0.31
9	21.35	21.35	0.45	1.04
10	278.38	278.39	89.33	85.36
11	317.93	318.04	39.91	124.94
12	54.75	54.75	3.80	7.09
13	10.15	10.14	0.53	0.51

, the results of the analytical and numerical calculations of RBLS and RBLV for the 13 cases of Table 3 are listed.

Based on Table 4, the analytical model for RBLS is in good accordance with the results of the numerical models; however, a remarkable discrepancy exists between the numerical and analytical calculations of the RBLV. Furthermore, for case number 6, the analytical and numerical model for RBLS and the numerical model for RBLV could not specify any value for these parameters because the constructed blocks are columnar and stretched in the z-direction (all dips are 90° with various dip directions). Therefore, by increasing the size of the numerical model, the RBLS and RBLV increase, and hence, the numerical RBLS and RBLV depend fully on the size of the numerical model. However, for such a case, the analytical model defines 10.57 m³ for RBLV. Because of this discrepancy and the existence of noticeable differences between numerical and analytical calculations of the RBLV, this method needs to be revised.

3.3. Evaluation of the Effect of Block Characteristics on the Inflow Rate

The average inflow rate to the tunnel that is excavated in the y-direction (S-N) of a rock mass is calculated using 3DEC version 7 software for the 13 cases of Table 3. To simplify the model and place more focus on the effects of the block characteristics on the inflow rate, the hydraulic aperture is assumed to be constant for all joint sets, and a fixed level of the water table is applied to all of the models.

The variations on the unit inflow rate to the tunnel with the numerical and analytical RBLV for the cases of Table 3 and Table 4 are shown in Figure 12 using the logarithmic scale of the block volume. For this purpose, a tunnel is excavated in the model and a fixed level of water table is applied to the formation, as well. The geometrical characteristics of the tunnel and join set as well as the hydrogeological condition of the rock mass are given in

Table 5. Geometrical and hydrogeological characteristics of the tunnel and rock mass for evaluation of the impact of rock geometries on inflow rate.

Case No.	1	2	3	4	5	6	7	8	9	10	11	12	13
Water head (m)	40	10	10	10	40	100	10	40	40	40	40	100	10
Tunnel radius (m)	4	1	4	4	2	4	1	2	4	1	1	4	1

.

The inflow rate decreases with the increase in the volume of rock blocks. However, this dependency is a little more relevant for the numerically calculated RBLV than the analytical ones (based on the R² of the trendlines). The reason is probably because of the accuracy of the numerical models for calculation of the block volume compared with the analytical model. However, as a general rule, the inflow rate to the tunnel is expected to decrease by increasing the block volume.

The relation between the logarithm of RBLS and inflow rate to the tunnel is illustrated in Figure 13. As the results of numerical and analytical models are almost the same, the diagrams in Figure 13 are identical. However, similar to Figure 12, the inflow rate to the tunnel decreases by increasing the surface of the block although the relationship between RBLV and inflow rate seems to be more relevant than the RBLS and inflow rate.

Another block characteristic assessed in this study is the 3D volumetric fracture intensity or P₃₂, which is defined by Equation (17). The P₃₂ represents the fracture area per unit volume of the block.

$${\mathrm{P}}_{32}=\frac{\mathrm{RBLS}}{\mathrm{RBLV}}$$

(17)

where P₃₂ is the 3D volumetric fracture intensity. Based on the data presented in Figure 12 and Figure 13, the relationship between the inflow rate to the tunnel and P₃₂ is illustrated in Figure 14. From analytical calculation of the block volume and surface, no logical relationship exists between the P₃₂ and the inflow rate to the tunnel. However, the inflow rate to the tunnel increases by increasing the numerically calculated fracture intensity.

The relationship between the unit inflow rate to the tunnel and P₃₂, which are calculated by analytical methods (Figure 14a) and numerical (Figure 14b), shows a significant difference. The reason is in Table 4, which shows the accuracy of the current analytical model for calculating block volume, and, therefore, the P₃₂ calculated by numerical and analytical models are not the same for a particular fractured rock mass. Accordingly, as shown in Figure 14, the relationship between inflow rate and P₃₂, which is calculated numerically, is more relevant than that calculated by analytical methods due to the lower accuracy of the analytical method for calculating block volume.

4. Discussion

In all geomechanical and hydrogeological studies of the rock mass, the first step is to determine the reliable dimension of the rock mass that is representative of its unlimited size. In the case of excavation of a circular tunnel in a rock mass and for evaluation of the inflow rate, as the flow takes place through the fractures at the wall of the tunnel, the length of the fractures per surface area of the tunnel is the determining factor in this regard. Meanwhile, the surface area is defined by the tunnel radius and length, and since the radius is assumed to be constant all over the tunnel, the length of the tunnel is the parameter that its variation affects the inflow rate. With this approach, STL is the shortest length of the tunnel that integrates the flow rate contribution of the entire fracture sets, and, hence, the inflow rate per meter of the tunnel length to the STL is equal to the average inflow rate to an unlimited length of the tunnel. As illustrated in Figure 15, the inflow rate to 1STL is equal to 7STL and generally, it will be equal to n × STL. As a result, the STL is the minimum length of the tunnel that could represent the hydrogeologic characteristics of the tunnel with unlimited length.

To specify the STL and simplify the numerical simulations, it was assumed that the hydraulic apertures of all joint sets are equal. However, different values of the hydraulic aperture do not affect the validity of the STL because that arrangement that exists in each STL, repeats at the subsequent multiple of STL and, as a result, the inflow rate will be again equal for all next lengths equal to STL.

All the apparent spacings at the wall of the tunnel are assumed to be integers to simplify the validation process of STL in Section 2.2. However, the apparent spacings are mostly decimal values in nature. For example, if the apparent spacings of the three joint sets are 1.3, 0.4, and 2.7 m, the LCM of the spacings or STL is 140.4. However, through a small variation of spacings to 1.5, 0.5, and 2.5 m, the STL is 7.5 m. In cases such as this, if the 7.5 m considered as the STL, the inflow rate varies in 2STL, 3STL, etc., compared to the STL. However, the variations of the inflow rate for 2STL and 3STL compared to the STL is small enough to be neglected and the domain of variation decreases with the increase of n. Therefore, in such cases, the inflow rate to 1STL could provide an acceptable estimation of the inflow rate to n × STL and, hence, it is recommended that the STL (or LCW) should be calculated by rounding the apparent spacings to an acceptable value.

STL is an important parameter for the evaluation of the inflow rate to the tunnel and is the minimum length of the tunnel that should be considered in the studies. Its importance is much more in three-dimensional numerical simulations because it is possible to excavate a tunnel in various directions in a rock mass and according to the direction, the probability of intersecting discontinuities by the tunnel will be different. As a result, the STL and inflow rate vary according to the tunnel direction, and, hence, to simulate the inflow rate numerically, a specific STL should be defined for each direction of the tunnel. As the only joint sets that are parallel to the tunnel direction (when normal to joint set is perpendicular to tunnel direction) are considered in 2-dimensional numerical simulations; in such cases, any length of the tunnel could be considered as STL, and, as a result, the absence of the STL is one of the major weak points of 2-D modeling. In practical works and to design seeped water evacuating systems, consideration of the STL could be effectively useful for an acceptable estimation of the inflow rate to the tunnel.

An analytical method is developed for the calculation of the RBLS that is created by the 3 persistent joint sets, knowing their dip/dip directions and true spacings. This model accurately defines the unlimited surface area for the blocks that are created by three vertical joint sets having different values of dip directions. However, analytically calculated RBLV for such cases is not unlimited and, thus, Equation (9) is not sufficiently comprehensive at all. Furthermore, the results of the developed model for the calculation of the RBLS are validated by numerical simulations using 3DEC software. However, the results of the previously developed model for specifying the RBLV, Equation (9), have a debatable difference with the 3DEC software results (Table 4).

To numerically specify RBLV and RBLS, the size of the numerical model increases until at least one block among all blocks of the model reaches its maximum possible size. If the block characteristics are of interest, in such a state, the model is not in its representative size. In this regard, the size of the model should be increased until the average block size tends to the RBLV. This size maybe not equal to the representative size of the model when the other parameters of the rock mass such as spacing, aperture, etc. are of interest. Therefore, it is essential to define a separate REV for each interesting characteristic of the rock mass and discontinuities, the same as what was defined for STL in Section 2.

The block volume and surface that are defined by the numerical simulations should be evaluated carefully before being used in the studies, as the average and maximum block size and surface area that is defined might be different for one specific model with various sizes. Furthermore, these characteristics are highly affected by the arrangement of the discontinuities, i.e., the start point for the generation of the joint sets in numerical models. For avoiding this error, either the size of the numerical model should be large enough to include mostly the RBLs, or just the analytically calculated block volume/surface should be used.

Based on Figure 12, Figure 13 and Figure 14, it is specified that a relevant relationship exists between the inflow rate to the tunnel on the one hand and numerical RBLV, numerical and analytical RBLS, and numerical P₃₂ on the other hand. According to the coefficients of determination (R²) of each diagram, it could be concluded that the relationship between P₃₂ and inflow rate is more relevant than the other parameters and, hence, P₃₂ could be used for the prediction of the inflow rate to the tunnel that is excavated in a rock mass. The advantage of using block characteristics for estimation of the inflow rate is that the RBLV, RBLS, and P₃₂ could be representative of spacing, dip, and dip direction of joint sets and, in this regard, various interesting parameters in the studies could be summarized in one parameter.

Like all other analytical methods, the models that are developed in this article have limitations that are related to the considered assumptions. The presence of three persistent joint sets with fixed values of spacing and orientation for each set constitute a limitation of the analytical model for block surface calculation. Regarding the STL, in addition to the limitations of block surface model, the cross section of the tunnel has a single and fixed value all along its length and the aperture is fixed in all fractures of a joint set. Generally, the presence of persistent joint sets is likely the most significant limitation when using the analytical models that are developed in this article.

5. Conclusions

Before the calculation of the inflow rate, it is necessary to determine the minimum length of a tunnel that should be considered to obtain representative values of the inflow rate for the entire tunnel length. This STL is equal to the LCM of the apparent spacings of the joint sets at the wall of the tunnel. Using a numerical simulation with 3DEC version 7.0 software, the unit inflow rate (inflow rate per each meter of the tunnel length) has been shown to be constant for any subsequent integer multiple of the STL (n × STL). Accordingly, the STL is also verified to be the representative length of the tunnel for hydrologic purposes. In developing the STL, the discontinuities are assumed to be persistent and the joint sets are in their ideal states, i.e., with fixed values of spacing, aperture, and orientation. Hence, the STL is recommended to be considered as the representative length of a tunnel to gain immediate and reliable results in all numerical studies relating to the evaluation of the inflow rate to a tunnel, especially in three-dimensional numerical simulations. In addition, more reliable estimates of the inflow rate could be obtained by considering STL in underground works, such as excavation of an underground tunnel. The STL value can be estimated when the dip, dip direction, and spacing of the joint sets are known and is independent of the hydraulic aperture.

Block surface area is an important parameter in the geomechanical and hydrological investigations of a rock mass. An analytical model is developed for calculating the block surface area when the orientations and spacings of the joint sets are known. The accuracy of the model is supported using the results of the simulation using the 3DEC software. The block surface area measured by the numerical models are shown to be possibly misleading because of the boundary effects as several blocks are often cut along the model boundaries. This error can be reduced by considering a simulation domain large enough to include a larger number of RBLS.

An analytical model that was previously developed for the calculation of the block volume for the case of three joint sets is used in this article. This model needs to be revised because a comparison of the results of the model with the output of the numerical simulation shows a significant difference. Block size is a very important parameter and is widely used for predicting the geomechanical and hydrological behaviors of the rock mass. Therefore, the acceptable accuracy of this parameter could guarantee the reliability of all other dependent parameters.

The relationship between the inflow rate to the tunnel and the RBLV, RBLS, and P₃₂ are investigated. We observed that the inflow rate to the tunnel decreases by increasing the RBLV. Similar but weaker relationships can be observed between the RBLS and the inflow rate to the tunnel. Furthermore, the inflow rate increases by increasing the volumetric fracture intensity (P₃₂), and this relationship is more relevant than the relationship between the RBLV, RBLS, and inflow rate. As a result, RBLV, RBLS, and P₃₂ could be representatives of several geometrical characteristics of the rock mass, e.g., joint set orientations and spacings. Hence, new relationships could be developed for the calculation of the inflow rate to the tunnel as a function of block characteristics. From a practical point of view, the rock block characteristics could be used in the estimation of the inflow rate to a tunnel. As a general rule, a higher inflow rate should occur in a tunnel excavated in a rock mass with smaller blocks.

Based on the material presented in this article, possible topics for further research include examining and calibrating the proposed specific length of the tunnel (STL) with field data obtained from the real case of the tunnels and, also, evaluation of the impact of discontinuities characteristics on the STL.