Study on Unloading Relaxation Characteristics of Columnar Jointed Rock Masses Based on Displacement Back Analysis

We aim to understand the relaxation of columnar joint rock masses during the excavation process of the diversion tunnel of the Baihetan hydropower station. This paper inverts the deformation parameters of the relaxed columnar joint rock based on the displacement monitoring data, and introduces a relaxation factor to describe the deterioration degree of anisotropic parameters of the relaxed columnar jointed rock. The equivalent strain is proposed as the criterion of unloading relaxation and the threshold is also given. Based on the software Flac3d, a program for calculating anisotropic elastoplastic model is developed. The distribution of the relaxation zone of the diversion tunnel after excavation is simulated, and compared with the results of the acoustic detection to verify the correctness and rationality of the program, which can provide a necessary reference for the design and construction of hydropower projects.


INTRODUCTION
The unloading relaxation of rock masses is mainly caused by the redistribution of the surrounding stress. After excavation, the stress level in the area near the excavation face decreases sharply. In lowstress and tensile stress areas, the connectivity rate of the joints, fractures, and other geological structural surfaces in the rock mass increases, resulting in the deterioration of mechanical properties of the rock mass, which has a very negative impact on the stability and safety of rock engineering. Therefore, the study of excavation relaxation has become an important topic in the field of rock engineering (Nguyen et al., 2001;Martino et al., 2004;Ma et al., 2019), especially the determination of the unloading relaxation range caused by excavation as this can provide reference for the support treatment of the relaxed rock in engineering, which has important practical significance in engineering construction.
During the last several decades, researchers have studied the formation mechanism of unloading relaxation (Cai et al., 2004;Sayers, 1990), and found that the relaxation formation is microscopically manifested as the initiation, expansion, and penetration of joint fissures in rock. In terms of engineering practice, Lu et al. studied the spatial distribution of the unloading relaxation of the high rock slope of Xiluodu Hydropower station based on an acoustic emission test (Lu et al., 2013). Sheng et al. estimated the unloading relaxation zone in the slope of the Three Gorges ship lock and believed that the weakening degree of rock mass in the relaxation zone is between 23 and 45% (Sheng et al., 2002). Li et al. studied the evolution characteristics of rockburst in the diversion tunnels of the Jinping II hydropower station through in-situ experiments (Li et al., 2012). Qian et al. determined the scale and distribution of the unloading relaxation zone of surrounding rock in the diversion tunnel of the Jinping II hydropower station through numerical simulation (Qian et al., 2009).
In the above studies, the rock mass is considered as an isotropic material. However, during the excavation of the diversion tunnel of the Baihetan Hydropower Station, a large section of columnar jointed rock mass was encountered. This rock mass exhibits obvious anisotropic mechanical properties and had significant unloading relaxation characteristics after excavation. Therefore, the traditional isotropic constitutive models are no longer applicable to the rock mass. It is necessary to study the anisotropic constitutive model reflecting the relaxation characteristics of the columnar jointed rock mass, and then determine the scale of the relaxation area caused by excavation. Pietruszczak and Morz extended the classical isotropy criterion to anisotropy criterion through combining the spatial distribution of the strength parameters (Pietruszczak and Mroz, 2000;Pietruszczak et al., 2002), and established an anisotropic yield criterion based on the microstructure tensor, which is a satisfactory solution to the anisotropy problem. On this basis, the paper takes account of the excavation relaxation characteristics of the columnar jointed rock mass, introduces a relaxation factor to describe the deterioration degree of the anisotropy parameters of the relaxed rock mass, and establishes a relaxation criterion based on the equivalent strain. An anisotropic elastoplastic calculation model embedded in software Flac3d is also developed. The distribution of the unloading relaxation zone of the diversion tunnel of Baihetan Hydropower Station after excavation is simulated, and compared with the in-situ acoustic test results to verify the rationality of the model, which can provide necessary reference for the design and construction of hydropower engineering.

Engineering Geology
The Baihetan Hydropower Station is located in Qiaojia County, Yunnan Province, China. It is a hydropower station of over ten million kilowatts similar to the Three Gorges Hydropower station or Xiluodu Hydropower Station. The bedrock of the dam site is dominated by cryptocrystalline basalt, almond basalt, and breccia lava, with hard lithology; the columnar joints are developed in some lithologic sections. The rock masses of the abutment and foundation of the dam are poorly developed, and the joints, interlayer, and intralayer dislocation zones are widely distributed. The faults with steep inclination and structural planes cut each other to form multiple wedges, which makes the rock masses exhibit high anisotropy characteristics (as shown in Figure 1),   and affects the stability of the arch dam and safety of the underground cavern group. As shown in Figure 2, the Baihetan hydropower station has five diversion tunnels; three diversion tunnels are set on the left bank and two diversion tunnels are set on the right bank, with elevations ranging from 574 to 605 m. The red part in Figure 2 represents the section of the diversion tunnel where the columnar jointed basalt is exposed, and its range is relatively large. Figure 3 shows the typical columnar jointed basalt in the diversion tunnel. The irregular columns have obvious contours, and their diameters range from 13 to 25 m. The inclination angle of columns are 70°-80°. The cross section of the columns are mainly irregular pentagonal and quadrilateral. Parallel cylindrical longitudinal micro-cracks are developed in the columns.

Relaxation Characteristics
In the excavation process of the diversion tunnels, unloading phenomena such as structural plane opening, tensile fracture, and rock mass relaxation appear in the surrounding rock near the excavation face. After the adjustment of the stress field near the surrounding rock, the original structural planes open and new fractures form. The widths of unloading fractures gradually decrease from the entrance to the end of the tunnel. Figure 4 shows the comparison of relaxation characteristics of the columnar jointed basalt in different sections of the diversion tunnel. The unloading fractures in the entrance section are generally developed, and are filled with gravel, rock debris, and secondary mud ( Figure 4A). The relaxed columnar jointed basalt shown in Figure 4B

Analysis of Diversion Tunnel Displacement
Taking the No.4 diversion tunnel as an example, this paper studies the unloading relaxation of columnar jointed basalt. The diversion tunnel is excavated in three steps from top to  bottom. The upper layer is 9 m high, the middle layer is 9 m high, and the lower layer is 6.2 m high. After excavation, the three-point displacement gauges are used to monitor the deformation of multiple sections of the tunnel. The deformation displacement was recorded at the depth of the hole at 0, 2, and 9 m to obtain excavation response characteristics of the columnar jointed rock mass. The layout and deformation of the multi-point displacement gauges of the 1 + 075 section are shown in Figure 5. After excavation, the top arch has the largest deformation, the shallow deformation reaches 40.71 mm, and the deep deformation also reaches 29.89 mm, which is due to the stress concentration at the arch top. The deformations of the side walls are smaller than that of the arch top. The displacement gauge in Myd6-2 measuring hole is damaged and has no recording. The displacements recorded by gauges in the remaining three monitoring holes show that the deformations of the side walls gradually decrease from shallow to deep. And the deformation of the lower part of the side wall is slightly higher than that of the upper part of the side wall, and the maximum deformation appears in the lower part of the left side wall, which is 22.38 mm. Figure 6 shows the displacement sequence curves at the lower part of the left-side wall and the upper part of the right-side wall. At the initial stage after excavation, the deformation rate of the surrounding rock continues to increase. After a period of time, the deformation rate slows down and gradually becomes stable, and the deformation rate of the deep part is lower than that of the shallow part.

Back Analysis of Deformation Parameters
Considering the columnar jointed rock mass is a transverse isotropic material, the above rock mass has 10 mechanical parameters including deformation modulus E 1 , E 3 , Poisson's ratioc 12 , c 13 , shear modulusG 12 , G 13 , cohesionC, internal friction angle φ, and anisotropic azimuths dd, dip, where the anisotropic azimuths are constant values. Therefore, there are still eight unknown mechanical parameters. For the convenience of engineering application, the parameters in the model need to be selected before back analysis.
Through parameter sensitivity analysis, we know that deformation modulusE 1 , E 3 and shear modulusG 13 are more sensitive to displacement, so these three mechanical parameters are selected as the parameters of back analysis.
The displacements of monitoring points MYD6-1 (0 m), MYD6-1 (2 m), and MYD6-5 (2 m) were selected as the target values of the back analysis. The three parameter levels of 85, 60, and 35%, and the orthogonal sample scheme shown in Table 1 are considered. The anisotropic elastoplastic model is used for calculation, and the simulation results of the sample are shown in Table 2.
The mapping relationship between deformation modulusE 1 , E 3 , shear modulusG 13 and displacement is established by support vector machines (Zhao et al., 2003). On this basis, a differential evolution algorithm is used as an optimization tool to search and monitor the deformation parameters of columnar jointed rock mass corresponding to displacement (Storn and Price, 1997). Table 3 is the target value of back analysis, and Table 4 is the result of back analysis.
The results of back analysis are used as mechanical parameters to simulate the diversion tunnel deformation, and the calculated displacements of the monitor points are recorded and compared with the measured displacements (Figure 7). The comparison results show that the displacements calculated by the results of back analysis are very close to the measured displacements, therefore, the mechanical parameter obtained by back analysis can be used to determine the relaxation scale of the diversion tunnel.

Anisotropic Elastoplastic Model
Pietruszczak and Morz extended the classical isotropy criterion to anisotropy by combining the spatial distribution of strength parameters and established anisotropy yield criterion based on the microstructure tensor (Pietruszczak 1999). In the model, each scalar of anisotropy parameter is composed of mixed variables of  stress and microstructure tensor. This means the physical and geometric properties of materials and the loading direction can be used as a set of variables. The microstructure tensor a ij is used to mark the material structure, and its deviatoric tensor is: where a kk is the principal value of the microstructure tensor a ij . The loading direction can be calculated by the stress component: The generalized loading vector l i can be expressed as: The anisotropy parameter η a ij l i l j is introduced to represent the projection of the microstructure tensor on the loading direction. Ning (2008) proposed an anisotropic yield criterion considering tensile strength: where, N φ 1+sin φ 1−sin φ Anisotropic Tensile Stress Yields Function: The shear potential function of columnar jointed rock mass corresponds to the non-associated flow rule, which can be expressed as: where ψ is the dilatancy angle, N ψ 1+sin ψ 1−sin ψ The tensile potential function g t corresponds to the associated flow law, which can be expressed as follows: The strain increment can be expressed as the superposition of elastic part and plastic part: The incremental stress-strain relationship of the elastic part is:  (9) where K ijkl is the elastic stiffness matrix.
The strain increment of the plastic part is determined by the plastic flow rule: Therefore, the incremental stress-strain relationship can be expressed as:

Relaxation Factor
Relaxation factor is introduced to describe the deterioration of elastic modulus and shear modulus of rock after excavation.
When the rock has no relaxation, the relaxation factor is 0; when the rock mass is completely relaxed, the relaxation factor is 1. For the rock in engineering, the relaxation factor is a value in the range of 0-1. The expression of the degradation of elastic modulus used in this paper is E (1 − S)E 0 (Eberhardt et al., 1999); S is the relaxation factor, and E 0 is the initial elastic modulus. The relaxation factor describing the relaxation degree of rock mass can be defined as: S 1 − E E0 . According to the aforementioned results of back analysis, the deformation modulus E 1 , E 3 and shear modulus G 13 are affected by relaxation. Therefore, three relaxation factors are needed to describe the different deterioration rates of these three mechanical parameters: where S i (i 1, 2, 3) are relaxation factors, E 0 i and E i (i 1, 2) are the initial and relaxation values of the two deformation moduli of rock mass, and G 0 13 and G 13 are the initial and relaxation values of the shear moduli of rock mass, respectively.

Relaxation Evolution Equation
Considering the influence of the three principal stresses comprehensively, the relaxation threshold is defined from the respect of equivalent strain. When the equivalent strain of rock is greater than the relaxation threshold, the columnar jointed basalt enters the relaxation stage. The equivalent strain is ε ε 2 1 + ε 2 2 + ε 2 3 .Here, ε i (i 1, 2 and 3) are the values of the three principal strains of rock, and the relaxation criterion is: whereε p is the relaxation threshold. Referring to the research results of Zhou et al. ), a relaxation evolution equation directly related to strain is obtained, S 1 − e −m(ε−ε p ) , wheremis the parameters of columnar jointed basalt related to relaxation. Three relaxation evolution equations can be obtained as: There are four parameters (m 1 、m 2 、m 3 、ε p ) that need to be determined in the relaxation evolution equation.

Parameter Determination
The corresponding relaxation factor values can be calculated according to the results of back analysis in Table 4. Table 5 is the calculated relaxation factor, which can be used as the target value of back analysis of relaxation evolution equation.
Considering the limitation of the length of this paper, the specific method of back analysis is not given a more detailed description, although one can be found in the literature (Shen and Chun, 2014). The optimal results of back analysis arem 1 14, m 2 20, m 3 14, ε p 1E−05.
The average value of three relaxation factors S 1 3 (S 1 + S 2 + S 3 ) is used to determine the scale of unloading relaxation zone, and the value is 0.49.

MODEL VALIDATION
Based on the development environment of Visual Studio 2008, the customized anisotropic elastic-plastic model with relaxation is compiled into dynamic link library DLL files for software Flac3d to call and execution. Figure 8A shows the calculation model of the section of the diversion tunnel in software Flac3d; the model has a total of 14,012 elements. The model size is 60 m*60 m, and the ground stress was applied to the boundaries. The parameters of the model are shown in Table.6.

Result Analysis
Figures 8B-E shows the scale of the relaxation zone of the surrounding rock after each step of excavation. After the first step of excavation, the relaxation zone is mainly distributed on both sides of the top arch ( Figure 8B). After the second step of excavation, a small amount of relaxation zone appears on the upper part of the top arch, and the relaxation zone on the side walls expands rapidly ( Figure 8C). After the third step of excavation, the relaxation zone of the top arch expands into a through relaxation ring with a maximum thickness of 0.65 m. The relaxation zone on the side walls is thicker, and the thickest part of the relaxation zone appears in the middle of the side walls. The maximum thickness of the relaxation zone on the left-side wall is 6.17 m, and the maximum thickness of the relaxation zone on the right-side wall is 6.20 m ( Figure 8D). Figure 8E shows the expansion process of the relaxation zone after each step of excavation. It can be seen that the next step of excavation will cause the expansion of the relaxation zone to induce by the previous step of excavation. Table 7 shows the relaxation depth of the diversion tunnel obtained by acoustic detection. It can be found that the thickness and distribution of the relaxation zone calculated by the customized anisotropic elastic-plastic model are well consistent with the results of the in-situ acoustic detection, especially at the side walls. The simulation results are in good agreement with the measured values. At the top arch, the relaxation depth calculated by the model is slightly different from the measured values, but the error is still acceptable, and fully meets the engineering requirements.   1) After the excavation of the diversion tunnel, the relaxation zone is mainly distributed on the side walls of the tunnel, especially in the middle of the side wall. The maximum thickness of the relaxation zone is about 6.20 m; the side walls of the tunnel need to be paid attention to in the subsequent support treatment. 2) The thickness and distribution of the relaxation zone calculated by the customized anisotropic elastic-plastic model are well consistent with the results of the in-situ acoustic detection. At the top arch, the relaxation depth calculated by the model is slightly different from the measured values, but the error is still acceptable, and fully meets the engineering requirements. 3) According to the results of parameter sensitivity analysis, this paper believes that the relaxation of rock is mainly due to the deterioration of deformation parameters, and the deterioration of the strength parameters is ignored. In further research, the deterioration of the strength parameters can be verified according to in-situ experiments.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
JY was responsible for defining the research aims, collating and reviewing the literature and previous research for the discussion, conceptualizing the article structure, writing the manuscript draft, and figure conception and execution. QZ, WY, and RW jointly agreed upon the research objectives; provided supportive analyses and interpretations; and contributed, reviewed, and edited text. HZ contributed and edited text.