Long-term assessment of creep and water effects on tunnel lining loads in weak rocks using displacement-based direct back analysis: an example from northwest of Iran

The time-dependent stability of tunnels is an important and challenging topic, mainly when the tunnel is excavated in incompetent and weak rocks. The creep property of rock is one of the crucial mechanical properties of weak rock and the main factor affecting the long-term stability of rock masses. Also, water as an important environmental factor influences both the short-term and long-term behavior of rocks and is one of the causes of geotechnical engineering disasters, such as tunnel collapse, slope sliding, surface subsidence, etc. In this research, the effects of rock’s creep behavior and underground water on the long-term stability of the Shibli tunnels were analyzed. Geological maps and reports of Shibli tunnels show a highly jointed condition in the surrounding rocks which have been crushed by two orogenic stages. The Burger-creep visco-plastic model was used to simulate the tunnel host rock creep behavior. The model's parameters were adopted based on the displacement-based direct back analysis technique using a univariate optimization algorithm. In addition, the influence of underground water is assessed under the condition of the varying water table. Support capability diagrams were used to evaluate the loading created on the tunnel’s permanent lining due to the creep behavior of rock mass and underground water. This study suggests that the weak rock's creep behavior and underground water significantly affect the time-dependent stability of tunnels. Results show that the induced stresses due to the rock's creep behavior and underground water are more considerable in the tunnel spring-line. Also, the increasing 20 m in the water table approximately decreases ten years of tunnel lining stability time at the fault zone. Rocks creep behavior and underground water significantly affect the time-dependent stability of tunnels in weak rocks. Displacement-based direct back analysis using a univariate optimization algorithm was used to determine the CVISC model’s properties. Increasing 20 m in the water table approximately decreases ten years of tunnel lining stability time at the fault zone. Rocks creep behavior and underground water significantly affect the time-dependent stability of tunnels in weak rocks. Displacement-based direct back analysis using a univariate optimization algorithm was used to determine the CVISC model’s properties. Increasing 20 m in the water table approximately decreases ten years of tunnel lining stability time at the fault zone.


Introduction
The long-term behavior and stability of rock engineering structures such as tunnels have received considerable attention (Jia et al. 2020;Do et al. 2020;Sharifzadeh and Tarifard 2014;Pellet et al. 2009). A significant number of engineering applications and laboratory studies indicated the mechanical response related to weak and soft rock's time-dependent behavior is relatively complicated Yan et al. 2020;Török et al. 2019;Singh et al. 2018;Patil et al. 2018;Vlastelica et al. 2018;Lyu et al. 2017;Günther et al. 2015;Miščević and Vlastelica 2014;Zhou et al. 2014). One of the significant parameters affecting the underground spaces' longterm stability is weak rock mass creep response (Kovačević et al. 2021;Panthi and Shrestha 2018). In addition, underground water as an important environmental component influences weak rock's short-term and long-term behavior, especially in fault and karst regions (Verbovšek and Kanduč 2016).
In previous studies, some vital researches have been conducted on underground construction's time-dependent behavior. For the long-term study of underground spaces, most researchers focused on the creep behavior of rocks and proposed several valuable creep constitutive models (Wang and Cai 2020;Bérest et al. 2019;Sharifzadeh et al. 2013;Weng et al. 2010;Schenk et al. 2006;Pellet et al. 2005;Grgic et al. 2003). Recently, some researchers have assessed the time-dependent stability of the tunnel's permanent lining. Xu et al (2019) investigated the weak rocks creep behavior on the failure process of the final lining of the Dujia tunnel. Xu and Gutierrez (2021) studied the damage evolution of the final tunnel lining under the combined effects of corrosion of the temporary support and the creep behavior of the surrounding rock. Before the last decade, most of the suggested constitutive models focused on the first and second phase of rock creep, but in the previous decade, some researches were conducted on the numerical modeling of the third phase of rock creep (Liu et al. 2021;Zhang et al. 2020a;Wang et al. 2017;Sainoki et al. 2017;Haifeng and Zhigen 2016;Yang et al. 2015;Xu et al. 2012). Some researchers also simultaneously examine the effects of water and creep on a laboratory scale (Jia et al. 2018;Schenk et al. 2006).
Some experimental studies have shown that the water considerably impacts the weak rocks' creep behavior (Yu et al. 2015;Schenk et al. 2006;Hoxha et al. 2006). Hoxha et al. (2006) found that the creep strain rate in triaxial and uniaxial compression creep tests was substantially dependent on the relative humidity.
Most of the practical rock mechanics projects in the humidity region, such as tunnels, are constructed in dry or relatively dry conditions due to tunnel drainage and ventilation, but in this kind of underground space, the underground water level may be gradually increased during the tunnel service life. Therefore, a detailed evaluation of the creep behavior of rock under various underground water tables is essential to understand the time-dependent stability of the underground constructions.
Our research analyses the effect of rock mass creep behavior and underground water on the time-dependent stability of the Shibli tunnels. The surrounding rocks of the tunnels are mainly marl, calcareous shale, and gray to black shale. Our research concentrates on the part of the northern tunnel near faults 9 and 10 ( Fig. 2) in which the host rock is highly fractured, and there was a water rush during tunnel construction ( Fig. 4a) (Aradan Construction Co. 2009). In this paper, the CVISC model was selected to analyze the rock mass's creep behavior, and the tunnel behavior was assessed in terms of increasing underground water tables. The direct back analysis was used to adopt model parameters. After analyzing creep and underground water effects, the time-dependent stability of tunnel lining was examined.

Characteristics of Shibli tunnels
The Shibli tunnels were excavated in the westnorthern area of the central Iran geological formation. These tunnels as parts of the Tabriz-Zanjan highway connect Tabriz and Bostanabad. (Fig. 1). For the elimination of heavy traffic, the northern tunnel with a length of 2244 m and the southern tunnel with a length of 2289 m have been excavated. The Shibli tunnels' surrounding rocks mainly consist of marl, limestone, calcareous shale, and gray to black shale, that significantly crushed by two orogenic stages (Fig. 2). Based on the engineering geology studies, the tunnels' surrounding rocks were categorized in A, B, and C blocks (Aradan Construction Co. 2009). Block A, consisting of calcareous and sandy shale, is the most competent region of the tunnel's surrounding rocks. Block B is made up of gray to black shale and jointed limestone, and its RMR index is associated with the fair class. Furthermore, block C's host rocks include crushed black shale and marl (Sharifzadeh et al. 2013). Block C contains the poorest rocks that initiate the most challenging circumstances, particularly where high-water inflows and heavily tectonized zones coincide. Twelve small fault zones with a Northwest-Southeast strike deviate from the Tabriz fault pass through the Shibli tunnels. The northern and southern tunnels have maximum rock cover thicknesses of 180 and 188 m, respectively.
Rock mass classification methods like RMR, Q, and GSI were used to identify the rock mass's quality. Table 1 shows detailed information on rock mass rating. (Sharifzadeh et al. 2013(Sharifzadeh et al. , 2012Aradan Construction Co. 2009).
In Table 1, the UCS tests were conducted on ovendry conditions. As the most incompetent rock, the host rock of block C is mainly composed of crushed black shale. Based on previous researches, the softening coefficient for shale could be in various ranges. For   Guo et al (2012) calculated it as 0.78. A top-heading and benching approach was used to excavate the Shibli tunnels. Both tunnel crosssections are horseshoe-shaped, the width of the tunnels is 13.26 m, and their height is 10.5 m. In block C, the initial support system involves six meters of grouted rock bolts, shotcrete with a thickness of 25 cm, and steel rib with 50 cm spacing, which were placed after each tunneling sequence. The permanent lining with a thickness of 40 cm was installed one year after topheading excavation. Figure 3 shows the tunnel's cross-section, and Table 2 lists the proprieties of the temporary and permanent lining system.
Ground movement monitoring is necessary during tunnel building to ensure construction safety and obtain crucial data for back analysis. Instrumentation was used to monitor the deformations of the Shibli tunnels. Three convergence measurement stations were placed following an excavation of tunnel's upper part to monitor tunnel deformations. The array was extended to offer five-point convergence measurements with bench tunneling. Furthermore, three-point borehole extensometers were placed and recorded. (Sharifzadeh 2012). The research analyzes the long-term stability of the northern tunnel at chainage 27 ? 210. In this station, the lithology is composed of strongly jointed and crushed shale; four joint sets were measured, and the average spacing between joints was measured 60 mm (Aradan Construction Co. 2009). During tunnel construction in this location, due to faults 9 and 10 ( Fig. 2), water inrushed into the tunnel and created some difficulties in tunnel construction ( Fig. 4a) (Aradan Construction Co. 2009). In addition, high precipitation levels in areas near the tunnel (Fig. 4b) increase the possibility of underground water rise, especially in fault and karst regions. Figure 5 depicts the finite-difference grid used for this cross-section. The model's dimensions were selected greater than ten times the tunnel's span to reduce boundary impacts. Element sizes and aspect ratios are refined and gradually rise outwards near the tunnel wall.

Creep modeling
The CVISC model was adopted to assess the creep response of the tunnel's surrounding weak rock (Itasca Consulting Group 2002). Figure 6a represents the elasto-plastic volumetric behavior, and Fig. 6b describes the visco-elasto-plastic deviatoric behavior of this model. The Burger model refers to the viscoelastic constitutive law. The Mohr-Coulomb (MC) model refers to the plastic constitutive law that pseudo-simulates the creep's third phase. Because the plastic slider is not connected to a viscous dashpot, plastic yielding is not time-dependent and only depends on the stress. If the model is stressed above the slider's yielding stress, it behaves like an elastoplastic material. In contrast, if the model is stressed under the yielding threshold, it acts like the Burger body (Kabwe et al. 2020a).  For the CVISC model, Eq. (1) describes the deviatoric strain rate ( _ e ij ). Equations (2) -(4) express these three units' constitutive laws of deviatoric behavior, respectively. In addition, Eq. (5) represents the constitutive laws of volumetric behavior.
where _ e K ij is the Kelvin strain rate, _ e M ij is Maxwell strain rate, _ e P ij is Plastic strain rates, g K is Kelvin viscosity, G K is Kelvin shear modulus, g M is Maxwell viscosity, G M is Kelvin shear modulus, r 0 and e vol are volumetric components of the stress and strain tensors, k Ã is a multiplier that can be removed from the calculation later. For the MC unit, the plastic potential g and failure criterion f can be calculated as Eqs. (6) and (7), respectively: where r 1 and r 3 refer to the major and minor principal stresses,u is friction angle, c is cohesion, and w is dilation angle of the MC model. In this model, the Maxwell viscosity (g M ) and Kelvin viscosity (g K ) could be illustrated by purely viscous dampers. For the Maxwell part, an elastic element occurs instantly, similar to a spring, and relaxes immediately when the pressure is released. The viscous component of Maxwell expands when the stress is applied over time. A solid undergoing reversible viscoelastic strain is represented by the Kelvin component. When the material is subjected to constant stress, it deforms at a decreasing rate and reaches the steady-state strain asymptotically. The material relaxes and returns to its original shape by removing the stress.
For this model, the frictional slider predicts the instantaneous plastic deformation, yielding, and permanent strain rates in the accelerated creep phase. Furthermore, the anticipated plastic strains in this stage are strongly dependent on stress rather and not on time (Kabwe et al. 2020b). In fact, for the CVISC model, the plastic slider is independent of time because it is not attached to a viscous dashpot. As a result, plastic yielding in the tertiary creep stage captured using this model depends only on stress.

The time-step for creep modeling
In FLAC software, the time is the main discrepancy between time-independent and creep constitutive models. In creep modeling, the time-step indicates real time, but in static calculations, the time-step is employed as a virtual amount for stepping steady-state condition (Sharifzadeh et al. 2013;Itasca Consulting Group 2002). The default value for the time-step in FLAC software is zero, where the software considers the material as the linearly elastic (viscoelastic models) or the elasto-plastic (viscoplastic models) based on the circumstances. The zero value of timestep could be utilized to achieve equilibrium prior to commencing a creep modeling. Because the constitutive rules use the time-step in their calculations, it affects the modeling response. Suppose you want a system to stay in mechanical equilibrium all of the time (as in a creep analysis). In that case, the timedependent stress variations created by the constitutive law need to be small compared to the strain-dependent stress variations. In other respects, out-of-balance forces will be significant, and inertial effects may influence the solution.
A particular study of the time-step is required to ensure the creep modeling's stability. For creep analysis, it is essential to select an appropriate timestep because the constitutive laws employ it in the formulas, and it can influence the modeling's response. The maximum creep time-step for the CVISC model is determined as Eq. (8): 3.3 Sequential excavation of the tunnel and underground water modeling In this research, the effect of underground water is examined under the varying groundwater level. The underground water tables are assumed to fluctuate between 10 and 60 m above the tunnel's crown (Fig. 5). In this model, acceptable hydraulic boundary conditions are also specified. The top of the model was considered the seepage boundary, so the water pressure at this boundary is zero. The model's bottom, right, and left sides are fixed as impermeable boundaries. It is also assumed that the head of water is constant, and the water flow is stable. For the sequential excavation modeling, constructing the tunnel's top-heading and the bench was carried out by relaxing 70 percent of initial stress. The shotcrete and rock bolts are installed as a temporary support system, the tunnel is advanced to provide complete relaxation and pressure increases on the primary support. It is significant to mention that the shotcrete is adequately permeable to create minimal resistance to water movement when the studied water level is returned. After that, a final concrete lining is placed. A waterproof plastic membrane protects the concrete lining; the combination of concrete and membrane is not permeable.
After installing the final lining, the water level is increased to a different table above the tunnel crown. This simulation stage is done in two phases. First, to determine the pore-pressure distribution, a flow-only analysis was conducted. Then, because of the water pressure acting on the lining, a mechanical-only measurement was carried out to assess the loading and stress change in the lining.

Back analysis
This investigation used a displacement-based direct back analysis with a univariate optimization method. Back analysis techniques are now widely used as a helpful tool for detecting the rock mass's uncertain parameters, geometric systems, and initial-boundary conditions in geotechnical engineering problems with the assistance of deformation, stress, and strain monitoring through the construction of underground spaces. In general, there are two types of back analysis solutions: inverse and direct (Sharifzadeh et al. 2012). The direct method utilizes the trial quantities of uncertain variables as input data for the stress analysis algorithm, minimizing the discrepancy between measured and computed results. For complicated constitutive models, the direct back analysis is adaptable and practical. Following remarks are essential to choosing to perform a back analysis (Oreste 2005): (a) A computation model capable of evaluating the stress and strain in the rocks and analyzing the construction procedures. (b) A practical algorithm for decreasing the error between the numerical modeling results and the observed in situ measurements. (c) The error function.
In this research, the error function defined by Eq. (9) was employed to minimize the discrepancy between measured and computed deformations.
where u i and u m i , i = 1,2,…,n are the measurement data and corresponding computed values, respectively. In order to reduce the influence of measurement errors, a normalized error function is employed. It should be mentioned that the lost displacements were considered and updated in the modeling results because of the delay in the installation of the instruments. After that, the modeling results were evaluated to correct the measured data. Table 3 illustrates the final back analysis results for the parameters of the CVISC model and initial stresses ratio (K) at station 27 ? 210. Figure 7 compares the modeled and measured findings at various times. In fact, numerical modeling results agree well with measured displacements using back analysis results.

Long-term stability assessment of tunnel lining
The creep behavior of rock mass and supposed underground water at different tables have been modeled during the tunnel service life. Support capability diagrams were used to evaluate the stress generated on the final lining due to rock mass creep behavior and underground water (Chen et al. 2019). In the capacity diagram, the axial force-bending moment interaction (N-M interaction) diagram graphically represents the induced axial force and bending moment at the tunnel lining with the failure envelope. An axial-shear force interaction diagram shows the presentation of axial and shear forces on a tunnel lining, as well as the corresponding failure zone. Table 4 illustrates the limiting formula for creating an axial force-bending moment and axial-shear force interaction diagram (Carranza-Torres and Diederichs, 2009).
In Table 4, FS is a factor of safety, r c is maximum compressive stress and r t is maximum tensile stress. In our research, it is supposed that the permanent lining behavior is elastic. By long-term analysis of tunnel behavior, it is recognized that the axial force-bending moment interaction is more critical in this case study. The axial force and bending moment values are obtained by passing the time and fluctuating water. Tables 5 and 6 show numerical analysis results for different water tables and times corresponding to the tunnel crown and spring-line, respectively. In addition, Fig. 8 illustrates the N-M interaction diagram at the tunnel's crown where the water table is 40 m above the tunnel crown. Figure 9   Table 4 The limiting formula for the creation of capability diagrams Thrust-bending moment interaction

FS
The axial force limitation value (N) for compression results The axial force limitation value (N) for tension results The critical values of bending moment (M cr ) for both failure in compression and tension Thrust-shear force interaction N ¼ rcA The axial force limitation value (N) for compression results The axial force limitation value (N) for tension results The critical values of shear force (Q cr ) for both failure in compression and tension also shows the N-M interaction diagram at the tunnel's spring-line, where the water table is 10 m above the tunnel crown. Simulation results show that induced stresses can rise significantly over time for the crown and spring-line of the permanent lining. As shown in Fig. 9 and Table 6, for the tunnel spring-line, after 70 years and in the case of water table 10 m, concrete's compressive strengths would not sustain the pressures generated by the axial force and the bending moment. Table 6 represents that by increasing the underground water table, the stability time of tunnel lining decreases. In the tunnel crown, the lining is stable until the supposed underground water table becomes 40 m above the tunnel crown. In this case (Fig. 8), after 86 years, concrete's compressive strengths are insufficient to sustain the pressures created by the creep behavior of rock and underground water. Figure 10 represents the effect of fluctuating water tables on the stability time at the tunnel spring-line. As shown in Fig. 10, the increasing water table decreases the stability time of tunnel lining. The results show that when the water table is 10 m, the tunnel's springline stability time is 70 years and when the water  Figure 11 shows that the increasing water table significantly affects thrust force and bending moment at the tunnel's spring-line. For instance, for the specific tunnel life, for example, 60 years, numerical results show that tunnel lining is stable until the water table increases to 30 m above the tunnel crown. After this level, the concrete's compressive strengths will not withstand the pressures caused by the rising water table.

Discussion
The time-dependent deformation of tunnels and rock masses has been studied in detail for many years (Sulem et al. 1987). The behavior analysis is usually based on numerical modeling (Zhang et al. 2020a(Zhang et al. , 2020b, in situ measurements (Kovačević et al. 2021;Selen et al. 2019), or laboratory testing of deformation (Liu et al. 2021). The deformation of weak rocks is often linked to water and is especially critical when clay minerals are present (Jia et al. 2018;Lyu et al. 2017). In our back analysis approach, the monitored deformation patterns were compared with  (Fig. 7). Kovačević et al. (2021) also demonstrated that the in-situ monitoring data could be a valuable technique to determine rock mass time-dependent properties. Our findings indicate that tunnels constructed in the weak rock mass, especially for the faults zone, the tunnel lining's stresses rise considerably with time and fluctuating underground water. It is in good agreement with previous works, where time-dependent changes in viscoelastic rock masses and the interaction of liners were studied (Do et al. 2020). The introduction of underground water strongly modifies the system's stability, especially when weak rocks are present (Zhou et al. 2014). In our study, the tunnel crown is stable in the case of rock mass creep behavior, but by introducing underground water and in the case of water table equals 40 m, the pressures created by axial force and bending moment are too considerable for the tunnel crown to sustain (Table 5 and Fig. 8). In this condition, the creep behavior of rocks and underground water need to be investigated for the tunnels' long-term stability analysis (Xu et al. 2020). Tables 5  and 6 indicate that the tunnel's spring-line stresses are larger than the tunnel's crown. Furthermore, displacement measurements and numerical modeling (Fig. 7) show that displacement at the spring-line is more  table 10 m   water table 20 m   water table 30 m   water table 40 m   water table 50 m   water table 60 m   FS=1 FS=1.5 remarkable than the crown. These differences are mainly because of the anisotropic in-situ initial stresses and the tunnel shape. Xu et al. (2019) also showed that due to the creep behavior of the Dujia tunnel's host rocks, the cracks appeared earlier at the tunnel spring-line than the tunnel crown. The longterm analysis of the Shibli tunnel's final lining indicates that by passing time, the bending moment in the tunnel spring-line is more increased than the bending moment in the tunnel crown. It illustrates that besides the time-dependent behavior of rocks, the initial in-situ stress and the tunnel shape also influence the long-term induced stress in tunnel lining. Figures 10 and 11 indicate that the presence of underground water affects the long-term stability of the tunnel. For this case study, the increasing 10 m in the water table approximately decreases five years of tunnel lining stability time.

Conclusions
The Shibli tunnel's time-dependent behavior was examined by the finite-difference method. The creep behavior of rock mass and underground water effects were analyzed to study the time-dependent stability of the tunnel lining. This research shows that the creep behavior of weak rocks and underground water significantly influence the long-term behavior of tunnel lining. Direct back analysis with a univariate optimization algorithm was used, and the CVISC model's properties and initial stresses ratio were obtained. The comparison of measured and calculated results demonstrates that the back analysis techniques could adopt rock mass creep properties successfully. The effect of underground water was analyzed for the different groundwater tables, and long-term stability analyses of tunnel lining were conducted. The main findings of this study are as follows: • The stress distribution in the tunnel lining varies over time in response to the creep behavior of rocks and fluctuating subsurface water table. For the tunnel spring-line, after 70 years and, in the case of water table 10 m, concrete's compressive strengths will not sustain the pressures caused by induced axial force and bending moment and increasing water table. In addition, the critical stability time and water table for tunnel crown are 86 years and 40 m, respectively. • Simulation results indicate that the tunnel's springline stresses are larger than the tunnel crown. The average bending moment for the tunnel's springline is 2.5 times of the bending moment at the tunnel's crown during the tunnel service life. • Underground water decreases the time of stability of tunnel lining. The numerical modeling results reveal that the increasing 20 m in the water table approximately decreases ten years of tunnel lining stability time. • Increases in the water table have a more significant impact on the induced axial forces than on the bending moment. The axial force rises approximately 13 percent for the tunnel spring-line by increasing 10 m in the water table. • Induced tunnel lining stresses in the tunnel springline are higher than the tunnel crown. The reason is probably related to the initial stresses and the shape of the tunnel. • According to the current research, it is highly recommended that the time-dependent properties of rocks are considered in the long-term maintenance of the underground spaces constructed in the weak rock mass.