Tunnel stability in fault zone base on 3D numerical simulation

Tunnel construction in complex geological environments such as fault zones and discontinuities are some problems in tunnel stability. This research focused on rock mass characterization and 3D numerical modelling of the finite element method to evaluate the diversion tunnel stability at the Bener Dam in Purworejo Regency. Rock deformation and stress distribution around the tunnel was the main discussions in this study to comprehend the design of a support system that could overcome these geological conditions. The diversion tunnel was located on the right side of the dam, where the fault was found and extended in the area around the intake, passing through the tunnel path. Geological surveys, point load tests, and rock engineering properties tests were then carried out, showing the rock in the location was dominated by andesite breccias with GSI values 30-68 and moderate to high weathering levels. The fault plane obtained in the study area had an N305°E/87° orientation. The simulation result showed a redistribution of displacement within the fault zone. Displacement and maximum stress occurred in the tunnel that intersected the fault zone, resulting in a decrease in the stability of rock around the fault. Therefore, an adjustment of the support system is required.


Introduction
One of the most critical issues in tunnel design and construction is ensuring the stability of the tunnel opening because most tunnel failures can be predicted from the stability of the tunnel opening itself during the construction period [1]. Tunnel stability issues are pivotal when the construction is in a complex geological environment, such as in karst areas or across fault zones and weak layers.
Fault as geological structures that cut the continuity of the layers and are accompanied by displacements in the fault plane can cause significant issues in geotechnical engineering. The presence of slip planes and the separation of faults that pass through the tunnel can result in relaxation and stress redistribution. Under conditions of high in situ stress, rock-bursts and the risk of water seepage may occur. Thus, faults that pass through the tunnel can be the main cause of the instability of tunnel excavation during construction and operation. Therefore, it is essential to comprehensively and extensively understand the influence of faults on the stability of the rock surrounding the tunnel.
Several methods to analyze rock tunnel stability include the physical model experiment method, the analytical method, and the numerical simulation method. Each of these methods provides advantages and disadvantages. For instance, an experimental model requires much power and material resources. Analytical solutions usually employ several theoretical simplifications and also cannot be applied if there are complex geological conditions.
With the rapid advancement of computer technology, the numerical simulation method has become the most popular method because of its efficiency in solving complex geotechnical and rock engineering problems compared to other methods [2]. In the numerical methods, several conceptual models have been developed to address tunnel stability and rock problems, such as finite element models (FEMs) [3][4], boundary element models (BEMs) [5], and discrete element models (DEMs) [6]. Among the three methods, the finite element method is the most popular one. The advantage of this method is that it can In this study, a three-dimensional RS3 numerical model from Rocscience is used to simulate the effect of faults on tunnel stability. In this software, discontinuities due to faults are modeled as a liner interface which is used to determine the effect and interaction between faults in tunnel excavation. The Bener Dam diversion tunnel is selected as the topic of study in this research. The results of this study are expected to provide a reference and evaluation of the existing reinforcement system in the presence of complex geological conditions.

Engineering background and geotechnical conditions
This research took place in the Bener Dam diversion tunnel located in Central Java Province, precisely in Guntur Village, Bener District, Purworejo Regency. The geographical location of Bener Dam is 7 o 35'54.59" SL dan 110 o 1'12.84" EL or with coordinates X = 391927.12 mE and Y = 9159958.76 mS. Figure 1 shows the profile of the tunnel path and the geological cross-section of this diversion tunnel. The length of this tunnel is 851.63 m, with the shape of a horseshoe.  The results of geological and geotechnical investigations showed that the diversion tunnel path consists of fresh breccia and andesite rocks with moderate to high weathering levels. Based on the surface mapping on rock outcrops with various levels of weathering, several lithologies were obtained, including andesite breccias, polimic breccias, tuffs, and andesites [7]. Faults manifestations were also obtained at this location, precisely at 7°35'45.24" SL dan 110°1'2.29" EL or at coordinates 391598.49 mE and 9160252.58 mS with fault orientation N305°E/87° (Figures 1 and 2). The limitation of drilling information at this location, especially on geometry information and data on fault discontinuity areas, made the study in this research focused on a small area around the fault that intersected the tunnel line, precisely at the STA. +200 to STA. +300. The parameters used were observation data at the location and rocks around the surface to identify technical data on fault properties.

Estimation of strength parameter for numerical analysis
Core drilling data, laboratory tests, and mapping had been carried out to obtain information related to subsurface rock masses. The evaluation of rock mass in this study used the Geological Strength Index (GSI) developed by Hoek et al. [8]. The purpose of this rock mass development was to provide input data related to the technical properties of rock masses required as input in numerical analysis such as tunnels, slopes, or foundations in rocks. The GSI system is more practical to use than other classification systems such as the Rock Mass Rating (RMR) and Q-system, especially when dealing with jointed rock masses and very weak and sheared rock masses [9].
Based on the rock mass condition and laboratory test, rock mass strength could be estimated using the Generalized Hoek-Brown Criterion equation introduced by Hoek [ Where σ1 and σ3 are major and minor principal stress; σci is the uniaxial compressive strength of the intact rock; mi, mb, s, and a are rock mass material constants; D is a disturbance factor; Ei and Erm are the modulus of deformation of intact rock and rock mass (Mpa).
The fault orientation in this study was N305°E/87°. In this modelling, the fault properties used slip criterion parameter in the form of Mohr-Coulomb equivalents, namely c ' as the cohesion peak and ' as the peak friction angle obtained from equations (6) and (7) and stiffness parameter in the form of Normal Stiffness (Kn) and Shear Stiffness (Ks) which were calculated using equations (8) and (9) developed by Barton [11].
In equations (8) and (9), Ei and Gi are modulus of elasticity and shear of intact rock; Em and Gm are the modulus of elasticity and shear of the rock mass, and L is the average distance between the discontinuity planes.

Tunnel excavation and installation of supports
Based on the planning data, the construction of the Bener dam tunnel used the New Austrian Tunnelling Method (NATM) construction method. NATM was developed in 1957 and 1965 in Austria. The principle of this method is to utilize the stress of the rock mass around the excavation to stabilize the tunnel itself [12]. The geometry of the diversion tunnel is horseshoe-shaped. Meanwhile, the main support system used is shotcrete and rock-bolt in the first layer, and the secondary support for the layer below is a concrete lining, as illustrated in Figure 3. The material properties of the support system can be seen in Table 2.  Figure 3 The support system of the diversion tunnel

3D numerical modeling
In order to determine the deformation and stress characteristics of the rock around the tunnel, a numerical model simulation was performed using RS3 software version 4.017 from Rocscience. RS3 is a finite element method-based software that has been widely used in analyzing geotechnical structure problems for civil and mining engineering, especially in the analysis of underground excavations such as tunnels.
As illustrated in Figure 4, a numerical model with complex geometries with a size of 130 m in the xdirection, 90 m in the y-direction, and 157 m in the z-direction have been used in this analysis. The modeling included stratigraphic layers based on geotechnical characteristics (Table 1) (Table 2) is modeled with several construction stages that represent the construction steps. The results of the analysis that will be discussed are at the final stage of construction, where excavations and support systems have been installed. Several review points will be analyzed as shown in Figure 5, namely in the midplane area, fault plane, and four concerning points on the excavation side.

Analysis of total displacement
The results of the analysis of the total displacement of the tunnel at the time of construction are complete in the form of the contours of the total displacement along the tunnel, the maximum total displacement on the crown, invert, left, and right sides and the deformation vector can be seen in Figure 6-11.   Figure 8 The total displacement of the midplane before excavation Figure 9 The total displacement of the midplane after excavation Figure 10 The total displacement of the fault plane before excavation Figure 11 The total displacement of the fault plane after excavation In Figure 6, it can be seen the contours of the total displacement along the tunnel. The distribution of the maximum total displacement value is located around the faults and on the upper left side of the tunnel from the center to the fault plane. The maximum total displacement on the sides of the tunnel excavation can be seen in Table 3. Based on data from Table 3, the maximum total displacement value is most significant on the left sidewall, where the overburden stress due to the soil above it on this site is greater than on the other side. The distribution of the maximum value of total displacement is concentrated in the area of intersection with the fault. The deformation vector along the tunnel can be seen in Figure 7. On the inverted side of the tunnel, the direction of the vector increases, which indicates an uplift so that the total displacement value in this area is smaller. Overall, the movement of the relative deformation vector towards the right sidewall is in line with the contours of the model geometry in the form of slopes.
The distribution of the total displacement in the midplane and the fault plane before and after excavation can be seen in Figures 8, 9, 10, and 11. It can be known that the distribution of total displacement changes is located on the top and left sides of the tunnel.

Analysis of stress distribution
During tunnel excavation, the stress on the rock mass around the tunnel will be redistributed due to increased stress. Before excavation, the in-situ stresses in the rock around the tunnel, namely v, h1, and h2 are evenly distributed over the rock mass. When the tunnel excavation is carried out, the stress in the rock mass around the tunnel will increase and will be distributed into new stresses. The main stresses are 1, 2 dan 3. The following is a graph of the stress changes from the top side of the tunnel crown at the midplane and fault plane: It can be seen from Figures 12 and 13 that the stress changes occurred both at 1 and 2 before excavation and after construction is complete. The area of influence due to stress changes on the upper side of the tunnel crown occurs up to 2 times the tunnel span or about 18.2 m in the midplane, while the fault plane is 10 m from the top side of the tunnel crown. When compared to the midplane and the fault plane views, the main stress 1 in the fault plane is 4073 kPa while in the midplane is 2930 kPa.

Analysis of strength factor
It can be seen in Figures 14, 15, 16, and 17 the distribution of the strength factor in the midplane and fault plane before and after tunnel excavation. Analysis of elastic condition obtained that strength factor in the rock mass decreased around the tunnel with a strength factor value of < 1. This value indicates that the excavation will be less stable. In the fault plane, the distribution of the strength factor values has a smaller tendency than in the midplane. Figures 18 and 19 show the distribution of the strength factor values in the upper side of the tunnel crown.

Conclusions
In this paper, a 3D numerical simulation using the finite element method has been carried out to study the stability and the effect of the fault as a discontinuous plane in Bener dam diversion tunnels. A numerical model with complex geometries and geotechnical conditions has been used in this analysis. With the fault crossing the tunnel patch, the maximum total displacement distribution was located around the fault and in the upper left side of the tunnel, which was 0.01165 m. Overall, the movement of the deformation vector tended to the right sidewall in line with the contours of the model geometry in the form of slopes. An increase in pressure occurred in the rock mass around the tunnel during the tunnel excavation. The maximum stress that occurred in the fault plane that cut the tunnel was 4073 kPa. In comparison, the stress that occurred in the midplane was 2930 kPa. The area of influence of the increased stress on the top side of the tunnel crown was 18.2 m long or two times the tunnel span in the midplane and 10 m in the fault plane. Based on the analysis of elastic conditions, the strength factor decreased after excavation. The distribution of the smallest strength factor values lied in the fault plane with a strength factor value of 0.71. This value is smaller than one, which indicates that the tunnel is less stable. Therefore, it is required to study and adjust the existing reinforcement system.