Deformation Prediction of Tunnel-Surrounding Rock Considering the Time Effect of the Viscosity Coefficient: A Case of an NATM-Excavated Tunnel

Considering engineering problems such as complicated stress and the difficulty in controlling large deformation while a tunnel passes through a soft rock stratum, a theoretical prediction model of convergence deformation of tunnel-surrounding rock is proposed. Based on the longitudinal displacement profile curve reflecting the “space effect” of the excavation surface, the Hoek formula with better applicability was introduced to analyze and theoretically deduct the “time–space effect comprehensively.” By taking the influence of the “time effect” coefficient into account, an improved Nishihara model was established to derive the analytical equation of the viscoelastic–viscoplastic convergence of surrounding rock. Taking the Dingxi Tunnel of Wujing Expressway in Hunan Province, the physical and mechanical parameters of surrounding rock in the tunnel were firstly determined then they were used to calculate and predict the vault down of three typical sections with the scoping equation of surrounding rock deformation. Based on the calculation results, the causes of the differences between the measured and theoretical values were analyzed; moreover, it is indicated that the minimum root-mean-square error is 1.68, the minimum average error is 1.3%, and the correlation coefficient is 0.99. The comparison shows that the theoretical prediction results agree well with the corresponding field test results. The improved Nishihara model can accurately predict the final deformation of the surrounding rock. Simultaneously, it is also proved that the relevant parameters and the hypothesis and correlation of the nonlinear viscosity coefficient equation are reasonable, with particular effectiveness and applicability in practical engineering. This study is significant for further studying the tunnel-surrounding rock’s stability and accumulating theoretical and practical experience to develop rheological theory.


INTRODUCTION
In China, the center of tunnel construction is gradually shifting to the southwest and northwest. The construction difficulties are also shifting to high-stress areas with complex geology. Among them, carbonaceous shale is widely distributed in the strata of Southwest China. Because carbonaceous shale has the characteristics of obvious bedding development, low strength, and soft and easy disintegration, it has prominent rheological characteristics under a high-stress environment. Rheology is one of the important mechanical properties of rock, which is closely related to the long-term stability of rock mass engineering. Soft rock has distinct rheological properties. When tunnels cross such a layer, engineering problems, such as the complicated force and difficulty to control large deformation, are usually accompanied, severely restricting the construction and affecting the support's long-term stability (Sun, 2007;Li et al., 2018). Therefore, scientific evaluation and prediction of surrounding rock deformation and stability according to the rheological characteristics of tunnel-surrounding rock are essential to guide on-site construction and feedback design and a necessary condition to ensure tunnel construction safety. At present, the rock rheological constitutive models used to describe the tunnel rock mass under excavation conditions mainly include  the Kelvin, Bingham, Maxwell, Burgers, and Nishihara models. Among them, the Nishihara model can consider the viscoelastic characteristics of surrounding rock and the important characteristic of viscoplastic deformation. Therefore, the Nishihara model can describe the deformation of rock masses in tunnel engineering, such as sandstone, limestone, sandy shale, clay shale, and carbonaceous shale. To sum up, it is important to conduct further research on the deformation the mechanism of tunnelsurrounding rock under complex geological conditions, explore new methods for surrounding rock stability analysis, and establish a scientific prediction model, which can facilitate predictions of tunnel-surrounding rock deformation.
Some scholars have established constitutive creep models considering various effects by introducing different variables based on the Nishihara model (Shen et al., 2014;Li et al., 2014;Zhang et al., 2015;Li Z. Q. et al., 2017;Zhang J. X. et al., 2019). However, these models only consider the influence of a single variable and ignore the impact of time. At the same time, an improved Nishihara model was established based on the Nishihara model by considering different influencing factors, and it was also verified through the application. In the above aspects, a lot of important research work has been done to put forward many modified models for specific rock strata Cao et al., 2017;Wang et al., 2018;Zhu et al., 2019;Jin et al., 2019;Wang et al., 2019;Li et al., 2020;. However, their general applicability is limited to a certain extent. In addition, based on the traditional Nishihara model, some scholars use a nonlinear viscoplastic body or clay tank to replace clay tank and viscoplastic body or use the series element method to establish the improved Nishihara model (Xu et al., 2015;Zhang et al., 2016;Ou and fang, 2017;Li Y. G. et al., 2017;Pu et al., 2017;Lin et al., 2020).
Nonlinear improvements have been made to the existing models, but no research has been conducted on the time effects under nonlinear conditions. Considering the influence of time and stress on creep, some scholars established a nonlinear creep damage model and the corresponding constitutive equation Liu W. et al., 2020;Yan et al., 2020;Liu Y. et al., 2020). However, they did not study the influence of this model on the damage and deformation of surrounding rocks.
Some scholars use the Nishihara model to study tunnelsurrounding rock deformation and deduce the analytical expression of tunnel-surrounding rock deformation based on the improved Nishihara model (Zhang, 2016;Xiao et al., 2017;Yu, X.Y. et al., 2018;Yu et al., 2019;Yu et al., 2020;Yu et al., 2019). In recent years, the prevalent artificial intelligence method Zhou et al., 2017;Chen et al., 2019;Huang et al., 2020;Zhang K. et al., 2020;Zhu et al., 2020) provides a more potential development direction for tunnel construction deformation prediction. Combined with the application of deep learning algorithms and nonlinear theory, research on the multi-factor tunnel excavation system in complex geological environments shows excellent advantages (Melis et al., 2002;Huang et al., 2021;Huang et al., 2022). However, it is necessary to further combine the constitutive rock mass model and the viscoelastic plastic mechanics theory to give full play to this advantage.
In summary, previous scientific research achievements have made significant contributions to the optimization and improvement of the Nishihara model and the development of the prediction theory of tunnel-surrounding rock deformation, but there are still some shortcomings. In tunnel excavation, the surrounding rock will carry out stress redistribution from two dimensions of time and space. The time-space effect is an essential factor that must be considered in tunnel-surrounding rock deformation prediction. Therefore, based on the advantages of previous prediction methods, combined with the LDP curve reflecting the "spatial effect" of the excavation face, this study introduces the Hoek formula (Hoek, 2001) with better applicability to comprehensively analyze and theoretically deduce the "time-space effect". Considering the time-space effect of the plastic viscosity coefficient, an improved Nishihara model is established to study the deformation process of surrounding rock in the process of tunnel excavation and support. The theoretical model is then established by rheological theory, and the convergence calculation formula of surrounding rock deformation is deduced. Finally, a comparative study is carried out combined with engineering cases to provide theoretical support and technical reference for the design and construction of tunnel engineering.
Rheological properties of the weak surrounding rock are more susceptible to external forces, exhibiting nonlinear characteristics. The influence of viscous coefficients is not negligible, and it is closely related to the stress-strain state of surrounding rock. This study derived the analytical equation of the viscoelastic-viscoplastic convergence of surrounding rock based on the improved Nishihara model with the Matlab software.

Improved Nishihara Model
The traditional Nishihara model (Nishihara, 1957;Zhou et al., 2010) is composed of a Hooke body, Kelvin body, and Bingham body, as shown in Figure 1.
The viscosity coefficient η 2 in the traditional Nishihara model is replaced by the nonlinear viscoplastic coefficient η(t) related to time, so the model becomes an unsteady Nishihara model. This study divides the model into two parts: the generalized Kelvin model and the nonlinear Bingham model. The viscoelastic and viscoplastic regions of the surrounding rock are analyzed. The improved Nishihara model is shown in Figure 2:where η 1 , η 2 , and η(t) represent the coefficients of viscosity. E 1 and E 2 indicate the elastic modulus. σ s stands for the ultimate frictional resistance of St. Venant's body.

Basic Assumption
In this study, the calculation model assumes that the surrounding rock is a continuous isotropic rock mass mechanics model and does not consider the action of groundwater, and when the lateral pressure coefficient is λ(t), the viscoelastic plastic problem of stress and deformation of surrounding rock of a circular tunnel is assumed, and the original rock stress field is axisymmetric along the vertical direction. The viscoelastic-plastic surrounding rock mechanical model is shown in Figure 3. The radius of excavated tunnels is R 1 , R 1 -R 2 , and R 2 -R 3 , representing ranges of the crack zone and plastic zone, while the elastic deformation zone should be from R 3 to infinity in theory.
In addition, many scholars have made theoretical assumptions and derivations based on experiments for the research on the expression of the nonlinear viscosity coefficient (Xiong et al., 2010;Yan et al., 2010). However, the results do not have universal adaptability. According to the nonlinear viscosity coefficient characteristics, this experiment is firstly expressed as Eq. 1: (whether it is a reasonable hypothesis, rationality will be verified in the sixth part of this study.) where A and B are constants related to the surrounding rock properties, which can be obtained from experimental data.

VISCOELASTIC DEFORMATION ANALYSIS 3.1 Generalized Kelvin Constitutive Equation
Based on the derivation method of the tunnel lining displacement in the Zhao et al. (2016), the theory and the relationship between "spatial effects" were improved by combining LDP curves, in which Hoek's formula (Nishihara, 1957) was also introduced to conduct a comprehensive analysis and theoretical derivation of "space-time effects." The generalized Kelvin body analyzes the viscoelastic deformation. Considering that the shape changes only when the geometry undergoes elastoplastic deformation, it can be calculated because the volume deformation is almost negligible. According to the elastic mechanics, in the cylindrical coordinate system (r, θ, z), due to the axial symmetry of the tunnel calculation model, all stress components are independent of the coordinate θ, and the stress-strain is only a function of r, u. Therefore, the relationship between strain and displacement along the (r, θ) plane is as follows: where A(t) is a function of time, and substituting Eq. 4 into Eq. 2 yields: The generalized Kelvin constitutive equation is transformed into a stress bias form, and the boundary condition, σ r σ θ p, σ p when r → ∞, is considered, and the axisymmetric condition can be expressed as follows: where σ r , σ θ , σ-radial, circumferential, average stress; ε r , ε θ -radial strain and circumferential strain; E 1 , E 2 , η 1 -instantaneous, viscous elastic modulus, and viscosity coefficient.
For the periphery of the tunnel, the "equivalent initial geostress" p(x) corresponding to the displacement was introduced and substituted into Eq. 6 in the form of p(x) p 0 λ(t).
Considering the boundary condition of r r 0 , σ 0, the following relationship can be obtained: where p 0 is geo-stress.

Space Effect
To facilitate the study of the "space effects" of the working surface, we first consider the effects of elastic media. The "space effect" of the excavation surface is directly reflected in the radial displacement of the tunnel along the tunneling direction. This distribution curve is the longitudinal deformation profile curve, which is regarded as the LDP curve. At present, many scholars have made in-depth studies on LDP curves within the scope of "space effect." Panet (1995) summarized the following approximate expressions by applying the finite-element analysis method of elasticity theory and combining the relationship between the measured radial displacement around the tunnel and the distance between the driving faces: where R is the radius of the tunnel cave and x is the distance between the calculation section and the excavation face.
Hoek (1999) analyzed and fitted the field monitoring data of the Mingtam underground cavern project and summarized the following empirical expressions by using the related theory analysis method: For the applicability of the above two formulas, Sun (2007) thinks that the Panet method overestimates the radial displacement of the tunnel, making it is easy to underestimate the supporting load in practical engineering and inducing the design unsafety. The result of the Hoek method agrees with the numerical analysis in the tunnel vault, waist, sidewall, and other parts, which shows that Hoek's formula has stronger applicability and can better describe the "space effect" of the tunnel excavation surface. The radial displacement described by Hoek's formula varies with the excavation surface, as shown in Figure 4.
Therefore, this study transformed Hoek's empirical formula: , substitution Eq. 9: The relationship between the radial displacement u m measured around the tunnel and the cross-section distance x measured at a distance from the working face is as follows:  There is a relationship between u m∞ and constants χ as shown in Figure 5 below: Since the radial free displacement u 0 of surrounding rock had been generated before the excavation entry surface reached the cross-section of the survey point, the total displacement should be calculated as follows: If the equivalent initial geo-stress is p(x), the relationship between displacement and equivalent initial geo-stress can be expressed as follows: Therefore, the above displacement should be proportional to the initial geo-stress, and Eq. 12 should be rewritten as follows: where T is a constant relating to the characteristics of surrounding rock and the radius of the tunnel; Y is dimensionless.
, after substituting Eq. 13 into Eq. 11: Let λ 0 u 0 u ∞ , as shown by Eq. 9, when x/R 0, then λ 0 0.3078. Then Eq. 14 can be simplified as: According to Hoek's empirical formula, the value of λ 0 is 0.3078. According to the symmetry and continuity of Figure 5, we can see that:

Time Effect
Considering the influence of rock rheology, which has been neglected before, the Hoek formula is used to comprehensively analyze the "space effect" caused by the working face's excavation and the "time effect" caused by the viscous effect of the surrounding rock. Assuming that the tunnel is driven forward at a continuous and uniform speed, and the driving speed is υ, the simultaneous Eqs 15, 16 show that: When t > 0, the above λ e (t) expression is introduced into Eq. 7 to solve differential equations.
From the above formula, it can be seen that the displacement around the tunnel is u e A(t) r ; by taking into account r r 0 , Eq. 18 can be transformed into: According to Eq. 19, when t → + ∞, the limit displacement around the tunnel is u e ∞ (E1+E2)p0r0 E 1 E 2 by calculating the upper limit. During the tunnel's excavation, the surrounding rock in front of the working face is affected by the construction disturbance, promoting the redistribution of the stress in the surrounding rock, thus causing the initial displacement u 0 of the surrounding rock ahead.
To study the initial displacement, we should first discuss the case when t ≤ 0.
Equation 17 shows that when Equation 21 shows that excavation speed directly affects the initial displacement u 0 . If the tunneling speed is plodding, the viscous stress and strain of the surrounding rock will have enough time to be fully released. So, it can be concluded that as υ → 0, T a1 χ υ → ∞. Then transform Eq. 21, then u 0 On the contrary, if the tunneling face is speedy, it is instantaneously annoying from infinity to the study's monitoring section. The viscous stress-strain of the surrounding rock does not have sufficient time to release completely, and only the elastic deformation is emancipated. At this point, let us say: υ → ∞, then T a1 χ υ → 0, so there is u 0 → p0λ0r0 E 1 . According to the above analysis, the initial displacement u 0 of the surrounding rock expressed by Eq. 20 is between the two extremes, namely: Let E ∞ E1E2 E1+E2 , the upper form is transformed into: Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 10 | Article 843545 6 Therefore, the viscoelastic deformation of the tunnel is S e 2(u e − u 0 ).

Viscoplastic Deformation Analysis
In the process of tunnel excavation, besides the instantaneous elastic deformation released from the working face, the secondary stress generated by the stress redistribution of surrounding rock exceeds the yielding stress of the rock mass at that point, which results in the viscoplastic state that the deformation of surrounding rock increases with time. The "time effect" of rock rheology is related to viscoplastic flow in surrounding rock mass (Sun, 2007). This kind of viscous flow usually occurs in the weak surrounding rock with high stress, especially for the surrounding rock with well-developed joints and fissures and dissolution. It is also found that the nonlinear Bingham model is suitable for the rheological study of this kind of surrounding rock.
Considering that only when geo-materials undergo plastic deformation and volume deformation, the change of shape can be neglected, it is assumed that the volume of the plastic zone is constant, that is, the volume strain is zero. Then there are: thus, So the Bingham viscoplastic constitutive equation can be rewritten into the form of stress deviation. Furthermore, the boundary conditions are considered when r → ∞, σ r σ θ p, σ p and the axisymmetric state can be expressed as follows: when σ < σ s , Because there is no creep and stress relaxation in this case, we only analyze the second case When σ > σ s , When r r 0 , the σ r 0, σ p, ε r − A p (t) r 2 0 and axisymmetric conditions are substituted into Eq. 27, so p is calculated in the form of p p(x) p 0 λ p (t): When t < 0, the tunnel has not been excavated to the monitoring section, and the viscoplastic deformation of the surrounding rock of the section is minimal, which can be almost ignored. Therefore, it is only necessary to study the deformation after excavation. when t > 0, Ta is substituted into Eq. 28 to get: The above Eq. 29 is solved to obtain the expression A p (t) as follows: is substituted into Eq. 30 can be obtained: When t → + ∞, u p ∞ A p (t)/r 0 r 0 (p 0 +σ s ) η(t) t Ar 0 (p 0 + σ s ). Therefore, the viscoplastic deformation of the tunnel is S p 2u p .

Viscoelastic-Viscoplastic Deformation Analysis
The viscoelastic theoretical solution of surrounding rock deformation mentioned above can be considered as the stress redistribution of surrounding rock after tunnel excavation is completed instantaneously with the elastic or elastic-plastic wave's propagation velocity. The viscoelastic convergence expression of surrounding rock is S e 2(u e + u 0 ), and the viscoplastic convergence expression is S p 2u p . Therefore, the analytical formula of viscoelastic-viscoplastic convergence of surrounding rock is developed as follows: where σ s is the yield stress; E ∞ E1E2 E 1 +E 2 ; η 1 is the viscosity coefficients; A and B are constants relating to the properties of the surrounding rock, and C is a constant. In Eq. 32, excavation speed υ cannot approach infinity in tunnel engineering; the elastic deformation cannot be completed instantaneously in actual engineering, and the initial deformation u 0 can hardly be measured. Therefore, in order to make the calculated value more in line with the field monitored data, the above expression is simplified as S s 2(u e + u p ) for calculation.

CASE ANALYSIS
According to the above analysis, if the initial stress state and mechanical parameters of surrounding rock before tunneling are given out, the deformation of the surrounding rock can be calculated according to the theoretical formula. This study verifies the rationality of theoretical derivation through the example of the Dingxi tunnel monitoring project of Wujing Expressway in Hunan Province and predicts the final deformation of surrounding rock.

Project Overview
The Dingxi tunnel of the Wujing Expressway in Hunan Province is located in Huaihua City. The tunnel area is in a low mountainous landform with large terrain fluctuations, where the midline elevation is 510-648 m and the maximum height difference is about 148 m, as shown in Figure 6. The starting pile number of the left tunnel is ZK66 + 940-ZK68 + 066, with a length of 1,126 m, while the right tunnel is K66 + 926-K67 + 999 with a length of 1,073 m.
In this case, three typical sections of ZK67 + 220 (Grade V), ZK67 + 500 (Grade VI), and ZK67 + 900 (Grade V) on the left line of the Dingxi tunnel are selected for research. The detailed geological evaluation is as follows: (1) ZK67 + 220 section: the surrounding rock is the intensely weathered sandy slate, with a small amount of moderately weathered sandy slate. The joints and fissures of the rock mass are well developed, and the rock mass is broken.
(2) ZK67 + 500 section: the surrounding rock is the intensely weathered sandy slate, with developed joints and fissures, hard rock, and fragmented rock mass. (3) ZK67 + 900 section: the wall rocks are primarily distributed in the fault fracture zone, and they were intensely weathered sandy slates.

Field Monitoring Test Scheme
The Layout of Measuring Points According to the type of surrounding rock, the tunnel's depth, and the excavation method, measuring points were arranged in the tunnel wall along the longitudinal direction. The arrangement of measuring points and baselines was adjusted according to the specific construction scheme, and the distribution of measuring points should be set in the section where surrounding rock grades generally change. The spacing distance between measuring points is 5-10 m. The layout of measuring points is shown in Figure 7.

Measurement Frequency
The monitoring frequency in the experiment is to measure the data once a day.

Instruments and Equipment
The measuring base point is buried in the area outside the tunnel excavation longitudinal and transverse (3-5) times of tunnel diameter. According to the standard benchmark embedding method, two base points are buried for the mutual check. Therefore, the base point should be connected with the nearby benchmark to obtain the original elevation. A high-precision level observes the arch crown's settlement, and the "four fixed principles" must be observed. In other words, it is necessary to refer to the following requirements: fixed construction personnel, fixed station position, fixed measurement duration, and fixed construction sequence, to ensure that the measurement accuracy meets the specification requirements. When the measuring points of the convergence measurement section are buried, the distance between the measuring points and the excavation surface should be less than 2 m. The first measurement should be carried out within 24 h after the last blasting. According to the relevant measurement frequency requirements, considering the influence of the ambient temperature, the data is read three times to ensure accuracy. The testing instruments and equipment are shown in Table 1.

Mechanical Parameters
Because the hydraulic fracturing method is simple and easy to operate in the testing process without needing to drill rock core samples or using exact electronic instruments, the testing depth can reach more than 5,000 m. Therefore, this study applies it to the initial in situ stress test of tunnel-surrounding rock.
The principle of the hydraulic fracturing method is to make the excellent wall break by increasing the water pressure based on the designed measurement depth of the hydraulic pump, then to determine the pressure and fracture direction of the characteristic point in the fracturing process, and calculate the initial stress state of the rock mass at the measurement point. Hydraulic fracturing is a two-dimensional testing method that can determine the maximum principal stress's magnitude and direction and the minimum principal stress perpendicular from the borehole plane. Figure 8 illustrates hydraulic fracturing test device.
According to the geological engineering report of the Dingxi tunnel, the in situ stress of three representative sections, ZK67 + 220, ZK67 + 500, and ZK67 + 900 on the left line, was measured with the hydraulic fracturing method. The self-weight stress value was calculated according to the designed thickness of overlying strata (σ v ρgh, rock bulk density ρ 2.7g/cm 3 ). The test results are shown in Table 2.
Because the surrounding rock in this study area is intensely weathered sandy slate, it is suitable for the viscoelastic-viscoplastic rheological model. The mechanical parameters of the surrounding rock of the Dingxi tunnel are shown in Table 3. The model parameters in Table 3 were derived from the creep data under the corresponding stress by the .
Because the design of the Dingxi tunnel is a three-center circular section, in order to ensure the accuracy of a theoretical calculation, it is necessary to convert the circular radius of the tunnel used in theoretical calculation into the equivalent circle radius using Eq. 33. Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 10 | Article 843545 9 FIGURE 8 | Schematic diagram of testing device for hydraulic fracturing method. 1-recorder, 2-high pressure pump, 3-flowmeter, 4-pressure gauge, 5-high pressure hose, 7-pressure gauge, 8-pump, 9-packer, 10-fracturing section.
The equivalent radius of the Dingxi tunnel section is r 0 7.03 m, which was calculated by substituting b = 12.68 m and H = 10.08 m for the upper formula.

Calculation and Analysis
The vault settlement of ZK67 + 220, ZK67 + 500, and ZK67 + 900 sections was calculated by the analytic equation group S s 2(u e + u p ) and compared with the measured results.
The simultaneous Eq. 32 was used to calculate all the parameters into the equation group S s 2(u e + u p ), and the displacement expressions S s and ultimate displacement expressions S ∞ of surrounding rock deformation can be obtained as follows: The Section of ZK67 + 220 The surrounding rock grade in this section is V, taking A 0.19, B 0.58, and the movement speed of the working face was υ 1.2m/d. Monitoring can be started only when the driving distance of the working face exceeds 2 m of the cross-section of the measuring point. Therefore, when T a 3.5d, it can be calculated using the formula T a χ/υ (when x 2, χ 4.2). The theoretical final settlement calculated by Eq. 35 is S ∞ 49.30 mm. However, due to the secondary lining of the section when t 20d, the final actual settlement  Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 10 | Article 843545 12 measured in the field was S 48.97 mm. The comparison between the measured and theoretical values of the ZK67 + 220 section is shown in Figure 9.
The Section of ZK67 + 500 The surrounding rock grade for this section is IV, taking A 0.14, B 0.28, and the movement speed of the working face is υ 2.4m/d. Monitoring can be started only when the driving distance of the working face exceeds 2 m of the crosssection of the measuring point. Therefore, when T a 1.75d, it can be calculated using the T a χ/υ (when x 2, χ 4.2). The theoretical final settlement value calculated by Eq. 35 was S ∞ 36.77 mm. However, due to the secondary lining of the section when t 20d, the final actual settlement measured in the field was S 36.69 mm. The comparison between the measured and theoretical values of the ZK67 + 220 section is shown in Figure 10.
The Section of ZK67 + 900 The surrounding rock grade for this section is V, taking A 0.29, B 0.49, and the movement speed of the working face is υ 1.2m/d. Monitoring can be started only when the driving distance of the working face exceeds 2 m of the crosssection of the measuring point. Therefore, when T a 3.5d, it can be calculated using the formula T a χ/υ (when x 2, χ 4.2). The theoretical final settlement value calculated using Eq. 35 was S ∞ 53.70 mm. However, due to the secondary lining of the section when t 20d, the final actual settlement measured in the field was S 53.31 mm. The comparison between the measured and theoretical values of the ZK67 + 220 section is shown in Figure 11.

DISCUSSION
In practical tunnel engineering, the instability and failure of surrounding rocks are closely related to the rheological properties of the rock mass, and the reasonable constitutive model is the key to studying the rheological properties. A new and improved Nishihara rheological model was presented, based on which the convergence formula of tunnel deformation was derived (Eq. 34) to be applied to the actual tunnel engineering. In this study, based on the residual calculation results, the R MSE and the Mean Error (R E ) were calculated using the residual formula to evaluate the consistency of theoretical calculation results with real results and its calculation accuracy. The smaller the test standard value is, the more accurate the prediction result is, and the closer the theoretical calculation value is to the measured value. If the average error tends to "0," the estimation is considered unbiased. The closer the correlation coefficient R 2 is to "1," the higher the linear correlation between the theoretical and measured curves.
The following results are remarkable according to the comparative analysis of relevant parameters in Table 4 and Table 5. As shown in Table 4, the three sections' final theoretical calculated values were all larger than the measured values, with a difference smaller than 0.73%, which can be considered very consistent. It can be seen from Table 5 that the root-mean-square error of the dislocation-time curves of the three sections was 7.26 mm, the maximum mean error was 5.01%, and the minimum correlation coefficient was 0.86. It was analyzed that firstly, in the field monitoring process, the secondary lining was applied on all three sections when T = 20 d, so the residual surrounding rock deformation cannot be effectively measured. The final measured settlement value was less than the final settlement value calculated by theory. Secondly, it was assumed that the surrounding rock's elastic deformation is completed immediately after excavating the surrounding rock in the theoretical derivation of the model. However, due to the "time effect," the viscoplastic deformation occurred slowly. "time effect" is also the fundamental reason for the massive difference in the settlement values at 12 d in Figures 10, 11. Thirdly, the faster driving speed of the section in Figure 9 and the shallower buried depth of the section in Figure 11 are also important reasons for the massive difference. Fourth, the existence of groundwater is significant to the stability of surrounding rock. The theoretical model in this paper does not consider the effect of groundwater, but there is gushing water in the tunnel during tunnel excavation. Fifth, the local variation of the excavation method and the uncertainty of construction disturbance are also factors that can not be ignored. In conclusion, the above accuracy evaluation indicates that the improved rheological model proposed in this study can provide a reference and theoretical basis for analyzing the stability of surrounding rock under similar conditions in other places. However, there are some limitations in this study to be further discussed. Firstly, due to the limited testing time and the number of test samples, this model's validation was performed only by testing the deformation of the Dingxi tunnel. Therefore, some creep tests should be carried out on the rock mass of different tunnels in subsequent studies to verify the proposed improved rheological model. Secondly, rock mass fractures, groundwater, blasting, and other related factors were not considered in calculating the selected tunnel section. Considering these factors to make the improved rheological theory possible, it is a subject worthy of further study. Thirdly, the theoretical calculation results can only have particular reference significance in the tunnel's surrounding rock with relatively single lithology, stable rock structure, and good integrity. However, in engineering practice, the deformation, instability, and failure of surrounding rock are complex processes controlled by several interrelated factors. To evaluate and predict the qualitative risk of the stable surrounding rock, the improved Nishihara model should further consider the coupling effect of many factors, such as joints and fractures in the rock mass, groundwater, lithologic differences, blasting, and lousy geology.

CONCLUSION
Based on the monitorization and measurement of the Dingxi Tunnel on Wujing Highway in Hunan Province, this study establishes an improved Nishihara model to forecast the surrounding rock deformation of the Dingxi Tunnel. The conclusions are as follows: (1) Model establishment and theoretical derivation. By simultaneously considering the "time effect" of plastic viscosity coefficient and the "time-space effect" of the Hoek formula, an improved Nishihara model was established, and the viscoelastic-viscoplastic convergence analytical equation of surrounding rocks was derived.
(2) Application of the model. The comparison and verification show that the improved Nishihara model can accurately predict the final deformation of surrounding rock and further validate the assumptions and correlations of nonlinear viscosity coefficient equations. The correctness of the parameters indicates that this theoretical model has absolute validity and applicability in practical engineering. (3) Limitations and prospecting of the model. The model is affected by excavation speed and buried depth under different geological conditions, and the forecast accuracy has certain limitations. Therefore, the critical point of research in the future is to deeply excavate and improve the model by combining in-depth learning and other methods.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
DH was responsible for determining the research purpose, organizing and reviewing the literature and previous research for discussion, determining the article structure, writing the draft manuscript, graphic conception, and implementation. XL, YL, and YW jointly agreed on the research objectives; provided supporting analysis and interpretation. LJ provided some modification suggestions and edited local text.