Study on Mechanical Properties and Nonlinear Strength Model of Deep Water-Sensitive Rock in Xianglu Mountain Tunnel

The complex mechanical properties of deep water-sensitive rocks during construction bring great challenges to the construction and stability analysis of underground engineering in water-sensitive strata. Taking Xianglu Mountain Tunnel in the Central Yunnan Water Diversion Project as the engineering background, the mechanical properties of two types of water-sensitive rocks including limestone and silty mudstone were studied through laboratory tests, and a nonlinear strength model based on Hoek-Brown criterion was proposed. The results show that (1) the stress-strain curves of deep limestone and silty mudstone have signi ﬁ cant strain softening characteristics and brittle-ductile transformation trend. Compared with limestone, silty mudstone has stronger ductility and lower critical con ﬁ ning pressure for brittle-ductile transformation; (2) the dilation e ﬀ ect of rock is a ﬀ ected by both con ﬁ ning pressure and plastic deformation. Under the same con ﬁ ning pressure and plastic deformation, the dilation characteristic of limestone is more evident than that of silty mudstone; (3) the Hoek-Brown criterion can describe the nonlinear strength characteristics of deep rock well. The evolution law of strength parameters m b and s with modi ﬁ ed plastic shear strain conforms to the negative exponential function, while the relationship between strength parameter a and modi ﬁ ed plastic shear strain satis ﬁ es cubic polynomial function; (4) the simulation results based on the nonlinear strength model are consistent with the test results, and the established strength model has high reliability. The research results can provide a basis for the stability analysis of the Xianglu Mountain Tunnel and the optimal design of the supporting structure and can also provide a reference for the study of the mechanical characteristics of the deep rocks in water-sensitive strata.


Introduction
With the rapid development of global economy and the sharp expansion of human living space, many underground projects under construction and planning continue to enter the deep strata [1][2][3]. Karst landforms are widely distributed in China. The distribution area of carbonate rocks in China is about 1.3 million km 2 . Especially in Southwest China, the karst area accounts for about half of the national distribution area [4]. Generally, many limestone, shale, and mudstone are distributed in karst areas, and their mechanical properties will be significantly weakened after they react with water, that is, they have strong water sensitivity [5]. Engineering disasters such as large deformation of soft rock and water inrush occur frequently in water-sensitive strata, which brings great challenges to the construction and stability of deep underground engineering. For example, the Yuanliangshan tunnel of Chongqing-Huaihua Railway in China, which passes through the soluble limestone stratum, has experienced dozens of water inrush accidents during construction, causing heavy casualties and economic losses [6]. The Qiguding Tunnel of China Meida Expressway encountered carbonized mudstone that resembled the characteristics of "dark chocolate," which was easy to hydrate and weather when it was exposed to water, and the self-stabilization ability of the surrounding rock after excavation was poor, which brought great difficulties to the construction [7].
The Central Yunnan Water Diversion Project is a large trans-basin water diversion project approved by The State Council of China, with a total length of 663.23 km [8]. The Xianglu Mountain Tunnel with a total length of 63.4 km is located in the first section of the main line of the Central Yunnan Water Diversion Project. The formation lithology along the tunnel is complex, including limestone, sandstone, and mudstone. The tunnel passes through 5 karst development areas, and the soluble rock stratum accounts for about half of the total length of the tunnel. The elevation of the tunnel site is generally 2400~3400 m, and the maximum buried depth is 1512 m. Many tunnel sections are located in the high to extremely high stress water-sensitive strata, which is the key control project of the whole water diversion project. In order to accurately analyze the stability of Xianglu Mountain Tunnel and provide scientific guidance for the design and optimization of supporting structure, it is necessary to carry out the study of mechanical properties of deep water-sensitive rock around Xianglu Mountain Tunnel.
In recent years, the research on the mechanical properties of water-sensitive rock has become a hot topic in the field of underground engineering [9,10]. Martin and Chandler [11] carried out uniaxial and tri-axial cyclic loading and unloading tests on Lac du Bonnet granite and studied the influence of cumulative damage on rock strength based on the Mohr-Coulomb (hereinafter referred to as M-C) criterion. The relationship between cohesion, friction angle, and damage variables is established. Based on Martin's work, Hajiabdolmajid et al. [12] proposed a CW-FS (cohesion weakening and friction strengthening) mechanical model for deep rock failure under high in-situ stress and gave a method to determine the corresponding parameters. Zhu [13] employed RV-84 Rheological testing machine to study the creep characteristics of tuff under dry and saturated conditions, and found that the influence of water content on the ultimate creep deformation of rock is extremely significant. Yang et al. [14] investigated the influence of water content on the creep mechanical properties of shale. The results showed that with the increase of water content, the creep rate of shale increased significantly. Xu et al. [15] and Li et al. [16] studied the influence of water content on the acoustic emission characteristics of sandstone and analyzed the energy evolution law during the experiment. The results show that with the increase of water content, the energy storage capacity and strain energy release capacity of rocks decrease. Xu et al. [17] studied the mechanical properties of typical rocks in the Danjiangkou reservoir area under dry, natural, and saturated conditions. It was found that with the increase of water content of samples, the uniaxial compressive strength, elastic modulus, and Poisson's ratio of rocks decreased, and the softening coefficient of rocks was directly related to the void ratio of rocks. Through the laboratory water absorption and dehydration tests of argillaceous siltstone, Jia et al. [18] studied the characteristics of density, wave velocity, stiffness, and strength under different saturation condition and proposed the meso softening mechanism of pore water on argillaceous siltstone. In order to reveal the softening damage characteristics of deep carbonaceous phyllite, Chen et al. [19] carried out mechanical tests of phyllite under different water content. It was found that the mechanical properties and brittleness of phyllite were weakened and the macro failure angle increased after encountering water.
Although many experts have done work on mechanical properties of water-sensitive rock, most studies focus on M-C strength parameters, that is, the evolution law of cohesion and friction angle with water content is mainly investigated. However, a large number of experimental results show that the strength envelope of deep rock has significant nonlinear characteristics under high stress conditions, and the applicability of linear M-C strength criterion is poor. Therefore, taking Xianglu Mountain Tunnel as the research background, this paper selects two types of water-sensitive rocks including limestone and silty mudstone to carry out laboratory tests. The mechanical properties and failure characteristics of saturated rock under different stress conditions are analyzed and discussed, with emphasis on the nonlinear strength characteristics and evolution law. Finally, the nonlinear strength model suitable for deep water-sensitive rocks is established, which provides scientific basis for stability analysis and design and optimization of supporting structure for deep tunnel in water-sensitive strata.

Laboratory Test
2.1. Rock Sample Description. Figure 1 shows the rock samples drilled from the construction site of Xianglu Mountain Tunnel. It can be seen that the limestone is grayish white, while the silty mudstone is grayish black. Through X-ray analysis of rock mineral composition, it is found that the main mineral of limestone is calcite, while the main mineral of silty mudstone is quartz and chlorite, containing a small amount of illite and feldspar. The retrieved rock samples are cut, cored, and polished to meet the specimen size recommended by ISRM, i.e. cylindrical specimen with height of 100 mm and diameter of 50 mm. Through the laboratory mechanical tests of two kinds of rocks in natural state and saturated state, it is found that the water content has a strong weakening effect on their stiffness and strength [20]. For limestone, compared with the natural condition, the elastic modulus, deformation modulus, and tensile modulus decrease by about 30% on average, and uniaxial compressive strength and tensile strength decrease by about 15% and 20%, respectively. Similarly, for silty mudstone, the elastic modulus, deformation modulus, and tensile modulus under saturated state are reduced by about 38%, and the uniaxial compressive strength and tensile strength are reduced by about 30% and 40%, respectively. Therefore, we can conclude that these two types of rocks are typical watersensitive rocks.

Test Apparatus and Method.
Uniaxial compression and conventional tri-axial compression tests are carried out by MTS-815 testing machine shown in Figure 2. During uniaxial compression, the displacement-controlled loading mode is adopted, and the axial stress is applied at the rate of 0.1 mm/min. When approaching the peak stress, the loading mode is switched to circumferential strain-controlled mode, and the circumferential strain rate is controlled at 5e-6/min. For the tri-axial compression test, the confining pressure is first applied at the speed of 0.02 MPa/s until the confining pressure increases to the predetermined value and then 2 Geofluids remains unchanged. In order to ensure that the test is carried out under quasi-static conditions, the axial force is applied at the speed of 0.01 MPa/s. When the axial stress approaches the peak stress, the loading mode is also switched to circumferential strain-controlled mode until the end of the test. According to the in-situ stress measured on site, the maximum principal stress of the rock mass in the tunnel area is about 47.7 MPa. So the confining pressure values set in the tri-axial test are 10 MPa, 20 MPa, 30 MPa, 40 MPa, and 50 MPa, respectively. In order to reduce the dispersion of test results as much as possible, rock samples with similar wave velocity are selected to carry out the test. Each group of tests is repeated three times, and the test results are taken as the average value. Due to the large depth, both types of rocks are actually saturated. Therefore, the specimens in this study need to be put into a vacuum saturator and saturated at -0.1 MPa air pressure before the experiment.

Test Results Analysis
3.1. Stress-Strain Curve. According to the tri-axial compression test results of deep limestone and silty mudstone under different confining pressure conditions, the stress-strain curve of typical specimens is shown in Figure 3. Here, it is specified that the compressive stress and strain are positive. In Figure 3, the right side is the deviatoric stress-axial strain curve, and the left side is the deviatoric stress-lateral strain curve. The deviatoric stress here is the difference between the axial stress and the confining pressure, i.e., (σ 1 − σ 3 ).
When the confining pressure is 0 MPa, the development of deviatoric stress-axial strain curve includes five stages. Taking the silty mudstone in Figure 3(b) as an example (see Figure 3(c)), they are void compaction stage (OA section), linear elastic deformation stage (AB section), pre-peak nonlinear deformation stage (BC section), post-peak softening stage (CD section), and residual stage (DE section). Under the condition of tri-axial compression, the initial compaction stage of the deviatoric stress-axial strain curve gradually disappears. The inclination of the pre-peak deviatoric stress-axial strain curve of two types of rocks is significantly affected by the confining pressure. The greater the confining pressure, the steeper the curve. It indicates that the increasing confining pressure can significantly increase the elastic modulus and deformation modulus of rock. However, for the deviatoric stress-lateral strain curve, the inclination of the pre-peak curve is almost not affected by the confining pressure. They basically coincide in the pre-peak stage.

Geofluids
In addition, the stress-strain curves of two types of rocks contain obvious post-peak softening stage and residual stage. When the confining pressure is low, the deviatoric stress will fall rapidly after reaching the peak strength. With the increase of confining pressure, the inclination of deviatoric stress-axial strain curve and deviatoric stress-lateral strain curve in the post-peak stage gradually slows down, and the ductility characteristics are significantly enhanced. This implies that the two types of rocks are gradually transformed from brittleness to ductility with the increase of confining pressure. For hard limestone, the post-peak curve still shows an obvious downward trend at 50 MPa confining pressure. However, when the confining pressure is greater than 30 MPa, the post-peak curve of silty mudstone becomes very gentle or even close to the horizontal. Therefore, the critical confining pressure for brittle-ductile transformation of silty mudstone is obviously smaller than that of limestone. Table 1 lists the deformation parameters of limestone and silty mudstone under different confining pressures. It can be seen that the deformation modulus of two types of rock under different confining pressures is always less than the elastic modulus.

Deformation Characteristics.
At lower confining pressure, the difference between deformation modulus and elastic modulus is large. With the increase of confining pressure, the gap between the two becomes smaller and smaller. This is because the rock is very dense under high confining pressure, which leads to the very close deformation modulus and elastic modulus. Although the elastic modulus and deformation modulus increase gradually with the increasing confining pressure, the increment value decreases gradually. When the confining pressure is lower than 30 MPa, the elastic modulus and deformation modulus increase linearly with the increase of confining pressure. At high confining pressure, they show a slow nonlinear growth, and the whole trend roughly conforms to the negative exponential function model shown in Eq. (1) below. Compared with elastic modulus and deformation modulus, Poisson's ratio is less affected by confining pressure. Due to the constraint of confining pressure, Poisson's ratio decreases slightly with the increase of confining pressure.
where the unit of coefficients A and C is GPa, and B is dimensionless quantity. The unit of σ 3 is MPa. In order to  5 Geofluids eliminate the dimensional inconsistency before and after the equal sign, σ 0 is equal to 1 MPa.
By fitting the data in Table 1 with the least square method, the fitting equation of rock elastic modulus and deformation modulus can be obtained.
Elastic modulus formula of limestone: Deformation modulus formula of limestone: Elastic modulus formula of silty mudstone: Deformation modulus formula of silty mudstone: Compare the curves drawn by Eqs.
(2)-(5) with the test value, as shown in Figure 4. The scatter points in the figure represent the test values. It can be seen that the fitted curves are very close to the test results. It indicates that the negative exponential model shown in Eq. (1) can well describe the nonlinear relationship between rock modulus and confining pressure.

Dilatancy Characteristics.
In order to study the dilatancy characteristics of deep rock in the post-peak stage, the variation curves of volumetric strain of limestone and silty mudstone with axial strain are drawn in Figure 5. Figure 5 shows that the variation trend of volume strain of the two types of rocks during loading is basically identical. Both of them have two stages: volume compression stage (OA section) and volume expansion stage (AC section). Furthermore, the volume expansion stage is divided into rapid expansion stage (AB section) and slow expansion stage (BC section) according to the expansion rate. These three stages correspond to the pre-peak stage, post-peak softening stage, and residual stage of the deviatoric stress-strain curve in Figure 3, respectively. In the pre-peak stage, the volume strain of rock decreases linearly. In the post-peak softening stage, the volume strain increases rapidly in a nonlinear manner. After entering the residual stage, the volumetric strain increases slowly at a near linear rate. When the con-fining pressure is 0 MPa, the growth rate of volume strain is very fast. With the increase of confining pressure, the constraint of confining pressure makes the growth rate of volume strain slow down gradually, which is shown as the slope of the curve in Figure 5 gradually slows down. In particular, when the confining pressure is greater than 40 MPa, the volume expansion of silty mudstone does not occur, indicating that the confining pressure has an obvious inhibitory effect on rock dilatancy.
In addition, the slope of AB section decreases gradually with the increasing confining pressure. It shows that the volumetric strain of rock is more sensitive to low confining pressure than high confining pressure. At the same time, it can also be concluded that under the same confining pressure, the growth rate of volumetric strain of limestone is significantly faster than that of silty mudstone. It indicates that the dilatancy characteristics of limestone are more significant than that of silty mudstone.
In order to describe the volume expansion effect during rock compression, scholars put forward the concept of dilation angle. In the conventional tri-axial test [21], the formula for calculating the dilation angle ψ is where dε p 1 and dε p 3 are the plastic strain increments in the axial and lateral directions, respectively. Referring to previous studies [16], the elastic-plastic coupling is not considered. Thus, the plastic strain increment in the axial and lateral directions can be separated according to the stressstrain curve in Figure 3, and then the dilation angle of rock can be calculated.
The plastic shear strain is defined as where ε p 1 and ε p 3 are the plastic strains in the direction of the maximum principal stress and the minimum principal stress, respectively, corresponding to the axial and lateral strains in the tri-axial compression test. Figure 6 shows the variation curve of dilation angle of limestone and silty mudstone with plastic shear strain under different confining pressures.
It can be seen from Figure 6 that the dilation angle of both limestone and silty mudstone decreases gradually with the development of plastic shear strain. When the plastic shear strain is small, the dilation angle decreases faster. As the plastic shear strain increases to a certain value, the dilation angle gradually tends to be stable. Furthermore, the confining pressure has a significant influence on the dilatation effect of the deep rock. Under the same plastic shear strain, the larger the confining pressure, the smaller the dilatation angle. Under the condition of uniaxial compression, the dilatation angle of rock can be about twice the dilatation angle at confining pressure 50 MPa. In addition, with the increase of the confining pressure, the distance of the dilatation angle variation curve is becoming closer and closer. It shows that the reduction of the dilation angle is smaller under the high confining pressure, and the dilatation characteristic of the deep rock is more sensitive to the low confining pressure. By extracting the critical plastic shear strain under different confining pressures in Figure 6, it is found that with the increase of confining pressure, the critical plastic shear strain of deep rock from softening stage to residual stage gradually increases. It roughly conforms to the variation trend of quadratic function shown in Eq. (8) below.  where γ p c is the critical plastic shear strain of the rock from the softening stage to the residual stage, the coefficients A, B, and C are all dimensionless quantities, and the meaning of σ 0 is the same as the above. Using the least square method, the relationship between the critical plastic shear strain and the confining pressure of the two types of rock is fitted according to Eq. (8).
The critical plastic shear strain equation of limestone is: The critical plastic shear strain equation of silty mudstone is: γ p c = 0:006 σ 3 /σ 0 ð Þ 2 + 0:113σ 3 /σ 0 + 5:657: ð10Þ Figure 7 shows the comparison between the curves drawn by Eqs. (9)-(10) and the experimental curves. It can be seen that the fitted curves are very close to the experimental results. It shows that Eq. (8) can describe the nonlinear relationship between critical plastic shear strain and confining pressure well.
In order to quantitatively analyze the influence of plastic shear strain on the dilatancy characteristics of surrounding rock, the plastic shear strain under different confining pressures is normalized, and the modified plastic shear strain is defined as 8 Geofluids Figure 8 shows the variation curve of dilation angle of limestone and silty mudstone with modified plastic shear strain. Figure 9 depicts the variation curve of dilation angle with confining pressure. It can be seen that the dilatation angle of the rock is affected by the modified plastic shear strain and confining pressure. That is, the dilatation angle is a binary function of the modified plastic shear strain and confining pressure. Under the same confining pressure or the same modified plastic shear strain, the dilation angle decreases gradually with their increase and finally tends to be stable. On the whole, the relationship between the dilation angle and the modified plastic shear strain and confining pressure meets the negative exponential change trend.
Therefore, the dilation angle model shown in Eq. (12) is constructed.
The dilation angle formula of silty mudstone:

Strength Characteristics.
According to the deviatoric stress-strain curve of rock shown in Figure 3, the peak strength, increment of peak strength, residual strength, and increment of residual strength of limestone and silty mudstone under different confining pressures can be obtained as shown in Table 2. Figure 11 depicts envelopes of peak and residual strength of limestone and silty mudstone. 10 Geofluids According to the analysis of Table 2 and Figure 11, the peak strength and residual strength of rock show obvious nonlinear characteristics. The existence of confining pressure greatly improves the bearing capacity of deep rock. With the increase of confining pressure, the peak strength and residual strength of rock exhibit a nonlinear growth trend, and the growth rate gradually slows down, which is consistent with the strength envelope described by the generalized Hoek-Brown (hereinafter referred to as H-B) criterion [22].
where σ 1 and σ 3 are the maximum and minimum principal stresses, respectively; σ ci is the uniaxial compressive  Table 2. Figure 11 shows the comparison between the fitted strength envelope and the strength envelope obtained from the test. It can be seen that they show good consistency, indicating that the gener-alized H-B criterion can well describe the nonlinear strength characteristics of deep rock.
3.5. Failure Mode. Figure 12 shows the typical failure modes of the specimens under different stress conditions. It can be observed that both types of rocks fracture along the weak position of the specimen under uniaxial tension. Moreover, the fracture is relatively flat and there is no rock debris, which is pure brittle fracture failure. Under uniaxial compression, splitting failure occurs in both types of rocks, and  13 Geofluids the splitting crack is roughly parallel to the axial direction. Compared with limestone, silty mudstone has worse integrity and more obvious fragmentation after damage. For the tri-axial compression test, the rock exhibits a composite failure mode of overall shear failure and partial splitting failure under the condition of low confining pressure. In addition, the inclination angle of the fracture surface is steep, and some rock debris is generated due to the shear slip of the fracture surface. When the confining pressure is moderate, the rock mostly shows shear failure of a single main crack surface, and the inclination angle of the fracture surface is relatively slow. However, the rock shows a conjugate "X" mode of multi-crack plane shear failure under the condition of high confining pressure. Therefore, the confining pressure 14 Geofluids has a great influence on the failure mode of the specimen.
With the increase of confining pressure, the failure mode develops as axial split failure-shear and split composite failure-monoclinic fracture plane shear failure-conjugate shear failure. Strain softening refers to the continuous deterioration of strength parameters of rock with the accumulation of plastic strain in the post-peak softening stage. That is, the yield criterion should include not only the stress tensor but also the plastic softening parameter in the post-peak softening stage. Therefore, the yield function of the nonlinear strength model of deep rock based on generalized H-B criterion can be expressed as

Nonlinear Strength Model
where κ is the softening parameter, and the strength parameters m b , s, and a are all functions of the softening parameter κ. Obviously, the key to the establishment of nonlinear strength model lies in the determination of strength parameters in the post-peak softening stage. At present, there are two main methods to determine the strength parameters in post-peak softening stage based on generalized H-B criterion. One is to select the GSI (Geological Strength Index) as the softening parameter, and the post-peak strength parameters are determined according to the post-peak softening law of GSI and the relationship between GSI and strength parameters [23,24]. The other is to select internal variables such as plastic shear strain, equivalent plastic strain, or plastic volumetric strain as softening parameters, and the post-peak strength parameters are determined by    15 Geofluids artificially assuming the evolution law of strength parameters with these internal variables [25,26]. For simplicity, it is usually assumed that the strength parameters vary linearly with the plastic internal variable. The evolution law of strength parameters in post-peak stage determines the rationality and accuracy of the model. Although the two commonly used research methods are simple, intuitive, and easy to understand, they all have some defects. For example, the first method ignores the limitation of the physical value of GSI, which easily leads to a large deviation in the estimation of the H-B strength parameter [27]. The second method artificially assumes the evolution law of the strength parameters and lacks a real and reliable experimental basis, and the determined strength parameters are highly subjective. In fact, when the stability assessment of the supporting structure is carried out, the reliability of the results largely depends on the real evolution law of the rock strength parameters in the post-peak softening stage.  Tables 3 and 4, respectively. Figure 13 shows the strength envelopes of the two types of rocks at different post-peak stages. From the analysis of Tables 3 and 4 and Figure 13, it can be seen that with the increase of confining pressure, the strength of rock in different post-peak stages increases significantly and shows obvious nonlinear characteristics. With the continuous development of post-peak softening stage, the peak strength envelope of rock gradually decreases and transits to the residual strength envelope. The strength envelopes with different softening degrees have significant similarities. It indicates that generalized H-B criterion can not only well describe the peak strength and residual strength envelope but also has good applicability to the strength envelope of different post-peak softening stages. Compared with limestone, the area between the peak strength envelope and residual strength envelope of silty mudstone is much narrower, indicating that the stress drop of silty mudstone in the post-peak softening stage is smaller and the ductility characteristics are more significant.

Evolution Law of Strength Parameters in Post-Peak
Stage. Through the nonlinear regression method, the H-B strength parameters corresponding to different modified plastic shear strains in the post-peak softening stage of the two types of rocks are obtained, as shown in Table 5. Figure 14 shows the evolution of strength parameters m b , s, and a with the modified plastic shear strain.
From Table 5 and Figure 14, it can be seen that the strength parameters m b and s gradually decrease with the increase of the modified plastic shear strain, and the decline speed is first fast and then slow, which roughly conforms to the negative exponential function model. Compared with parameter m b , parameter s decreases relatively faster. The strength parameters m b and s of limestone are larger than those of silty mudstone, indicating that the rock mass quality of limestone is significantly higher than that of silty mudstone. In addition, the influ-ence of lithology on the parameter m b is greater than that of s. However, with the emergence of the residual stage, the influence of lithology on the strength parameters becomes weaker and weaker. With the development of plastic deformation, the strength parameter a of the two types of rocks gradually increased, and the growth rate showed a slow-fast-slow trend, which roughly increased according to the law of the cubic function. In the whole post-peak softening stage, the parameter a of limestone is smaller than that of silty mudstone, indicating that the degree of nonlinearity of the strength characteristics of limestone is higher than that of silty mudstone.
From the above analysis, it can be concluded that the evolution model of the H-B strength parameter in the post-peak softening stage of the deep rock is as follows: where the coefficients A1~A3, B1~B3, and C1~C4 are all dimensionless quantities, which can be obtained by the fitting method. Putting this formula into Eq. (16), the nonlinear strength model of deep rock based on H-B criterion can be obtained

Program Development of Nonlinear Strength Model.
Since the H-B criterion was proposed, a large number of scholars have tried to develop the calculation program of the elastic-plastic model based on the H-B criterion [28]. Cundall et al. [29] used radial, associative, and composite plastic flow laws for different states such as tensile stress, low confining pressure, and high confining pressure, respectively, and developed an elastic-plastic finite difference calculation program based on the H-B criterion. In order to overcome the numerical singularity of the yield surface of H-B criterion at curved edges and sharp points, Wan [30] uses elliptic function to smooth the H-B yield function and realizes the finite element program development of H-B criterion based on Euler implicit backward integration algorithm. Based on the generalized H-B criterion, Clausen and Damkilde [31] and Chen et al. [32] divided the principal stress space into different regions, and used the implicit backward Euler algorithm to return mapping the stress to the edges or apex of the yield surface, and developed the corresponding elastic-plastic finite element calculation program. Compared with the traditional method of smoothing the irregular yield surface to solve the numerical singularity, the method proposed by Clausen is completely carried out in the principal stress space. The form is simpler and convenient for programming. Moreover, the idea of regional reflection of principal stress space can accurately realize the correction of plastic stress. However, most scholars currently focus on the elastic perfectly plastic model in programming research, and there are few studies on the post-peak softening model. In order to reasonably and effectively simulate the mechanical behavior of deep rock in the post-peak softening stage, the principal stress return mapping algorithm proposed by Clausen is adopted in this paper. Based on the user material subroutine UMAT of ABAQUS platform, the nonlinear strength model of deep rock obtained from the above research is programmed. Figure 15 shows the flow chart of principal stress return mapping algorithm. When the principal stress is reflected, it is necessary to calculate the partial derivative of the potential function. For simplicity, the selection of potential function in this paper is consistent with the yield function, that is, the associated flow law is adopted. User material subroutine UMAT is a secondary development program interface provided by ABAQUS to users to define specific material properties [33]. With the help of UMAT, users can easily define the constitutive relationship of materials and develop material models that are not in ABAQUS model library. When UMAT is combined with user-defined field variable program USDFLD, the material parameters of the model can be updated gradually in the calculation process to realize the dynamic evolution of material behavior. Since there is neither a material model based on the H-B criterion nor a material model that can reasonably simulate the post-peak softening stage of the material in ABAQUS, the modified plastic shear strain defined above is taken as a field variable, and the deformation of the rock is calculated by combining USDFLD and UMAT. By embedding the formula of deformation parameters and H-B strength parameters, the program development of the aforementioned nonlinear strength model can be realized. Figure 16 shows the flow chart of program development.

Model Verification.
In order to verify the reliability of the proposed nonlinear strength model and development program, tri-axial compression numerical tests were carried out on the above limestone and silty mudstone, and compared with the partial stress-strain curves obtained from laboratory tests. A cylindrical finite element model with the same size as the rock specimen (50 mm ×100 mm) is established as shown in Figure 17.  (20). Figure 18 shows the comparison of the results of deviatoric stress-strain curves of limestone and silty mudstone obtained from numerical simulation and test. The scattered points in the figure are the test results and the curve is the simulation results. Obviously, the deviatoric stress-axial strain curve and deviatoric stress-lateral strain curve obtained by simulation have post-peak softening stage and residual stage, which are close to the test results. On the whole, the nonlinear strength

Conclusions
Taking deep saturated limestone and silty mudstone as test objects, this paper systematically studies the characteristics of stress-strain curve, deformation and strength characteristics, and failure mode of two types of water-sensitive rocks under different stress conditions and obtains the nonlinear evolution law of post-peak strength parameters. The main conclusions are as follows: (1) Confining pressure has a significant influence on the deformation characteristics of rock. The elastic modulus and deformation modulus increase in a negative exponential function with the increase of confining pressure. The dilation angle decreases with the increase of confining pressure and the development of plastic deformation. The relationship between dilation angle and confining pressure and modified plastic shear strain satisfies a negative exponential function (2) The peak strength envelope, residual strength envelope, and the strength envelopes at different postpeak stages of the deep rock have obvious nonlinear (3) The failure modes of deep rock are significantly different under different stress conditions. In direct tension, the rock exhibits brittle tensile fracture failure, and in uniaxial compression, it shows axial multicrack splitting failure. With the increase of confining pressure, the failure mode gradually develops from splitting-shear composite failure to single inclined fracture plane shear failure and conjugate "X" multi-fracture plane shear failure (4) Based on the generalized H-B criterion, a nonlinear strength model of deep water-sensitive rock was established, and the corresponding calculation program was developed. The model is verified by carrying out numerical tests and comparing with the results of laboratory mechanical tests. The research results can provide a scientific basis for the stability analysis of Xianglu Mountain Tunnel and the design and optimization of the supporting structure and can also provide a reference for the study of the mechanical properties and models of the deep water-sensitive rock in similar projects

Data Availability
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest
The authors declare no conflict of interest.