A Meso/Macroscale Theoretical Model for Investigating the Large Deformation of Soft Rock Tunnels Considering Creep and Anisotropic Effects

The rheological deformation of soft rock resulting from tunnel excavation can lead to significant construction and safety challenges. In this study, a multiphase numerical model was developed to simulate the rheological deformation of soft rock surrounding a tunnel after excavation. The developed model considers the coupled meso/macroscale creep and damage processes of the rock using the coupled discrete element method–finite element method (DEM–FEM). In particular, the damage and deformation accumulation at the mesoscale (i.e., initial phase before excavation, loading phase due to the disturbance of the excavation and creep-induced damage phase leading to large deformation) were incorporated into the model. The model predictions were validated using field monitoring data. By incorporating the coupled meso/macroscale deformation process of the rock into the model, the predicted time-dependent displacements of the tunnel face agree reasonably well with the monitoring data. In addition, the results demonstrate that tunnel brittle damage accumulated in mineral clusters severely leads to instantaneous deformation, which becomes less important in the creep evolution stage. Furthermore, the results indicate that the final deformation is characterized by a high sensitivity to the value of mesoscale modeling parameters. A multiphase numerical model was developed to simulate the rheological deformation of soft rock surrounding a tunnel after excavation. The developed model considers the coupled meso/macroscale creep and damage process of the rock using the discrete element–finite element coupling method (DEM–FEM). By incorporating the coupled meso/macroscale deformation process of the rock into the model, the predicted time-dependent displacements of the tunnel face agree reasonably well with the monitoring data. Bedding joints were considered in the model, and the influences of joints on the local movement of grains as well as the final converged deformation were assessed. The influences of mesoscale tensile strength, mesoscale cohesive strength and the fraction of grains experiencing creep on macroscale deformation were discussed. A multiphase numerical model was developed to simulate the rheological deformation of soft rock surrounding a tunnel after excavation. The developed model considers the coupled meso/macroscale creep and damage process of the rock using the discrete element–finite element coupling method (DEM–FEM). By incorporating the coupled meso/macroscale deformation process of the rock into the model, the predicted time-dependent displacements of the tunnel face agree reasonably well with the monitoring data. Bedding joints were considered in the model, and the influences of joints on the local movement of grains as well as the final converged deformation were assessed. The influences of mesoscale tensile strength, mesoscale cohesive strength and the fraction of grains experiencing creep on macroscale deformation were discussed.

• A multiphase numerical model was developed to simulate the rheological deformation of soft rock surrounding a tunnel after excavation. The developed model considers the coupled meso/macroscale creep and damage process of the rock using the discrete element-finite element coupling method (DEM-FEM). • By incorporating the coupled meso/macroscale deformation process of the rock into the model, the predicted timedependent displacements of the tunnel face agree reasonably well with the monitoring data. • Bedding joints were considered in the model, and the influences of joints on the local movement of grains as well as the final converged deformation were assessed. • The influences of mesoscale tensile strength, mesoscale cohesive strength and the fraction of grains experiencing creep on macroscale deformation were discussed.

Introduction
The rheological deformation of surrounding rock caused by excavating soft rock tunnels can lead to significant construction and safety challenges due to creep-induced damage (Xu et al. 2012;Kabwe et al. 2020a, b). Therefore, it is necessary to study the time-dependent deformation mechanisms of the soft rock surrounding rock caused by tunnel excavation (Manchao 2014;Wang and Cai 2022a, b). The ultimate large deformation of soft rock is a creep-damage coupled problem (Li and Tang 2015;Wang and Cai 2020;Xia et al. 2021a, b). Inspired by the fundamental understanding of rock creep and damage mechanics, previous studies have mainly focused on the macroscale level using creep-damage constitutive models to reveal the rheological process of soft rock tunnels (Ping et al. 2016;Liu et al. 2017;Huang et al. 2020a, b;Li et al. 2022aLi et al. , b, 2023. In addition, the different stages of long-term deformation of soft rocks were also studied from clay mineral slippage and disintegration (Fahimifar et al. 2010;Jiang et al. 2016;Liu et al. 2020a, b;Kovacevic et al. 2021;Maheshwari 2021;Cui et al. 2022). These results further demonstrate that the large deformation of soft rock tunnels is caused by the combined effects of damage and creep (Song et al. 2016;Yang et al. 2019;Kabwe et al. 2020a, b;Mu et al. 2020;Yu et al. 2020).
However, challenges remain in correlating the dissolution and slide of clay minerals at the mesoscale with crack propagation and coalescence at the macroscale. Long-term failure propagation in soft rock engineering is actually a multiscale issue (Li and Tang 2015;Xue et al. 2021). The mechanical properties and spatial orientations of rock joints could significantly influence tunnel excavation (Jia and Tang 2008;Deng et al. 2014). From a macroscopic perspective, joints may affect the local mineral grain flow and cause anisotropic deformation of the rock mass Cai 2021a, b, 2022a). A multiscale analysis needs to be incorporated into the large deformation analysis of soft rock. Using a mesoscopic representative volume element (RVE), Tang et al. developed a theoretical model for simulating the trans-scale progressive failure of rock based on statistical and continuum damage mechanics (Li and Tang 2015). In addition, renormalization group theory (Xia et al. 2021a, b) and fractal geometry (Wei et al. 2019;Zhou et al. 2019) are commonly used to analyze meso-to-macro cascading in the phase-changing process, but they have limited capabilities of accurately predicting the deformation of the rocks surrounding the tunnel. Several attempts have been made to model the creep and time-dependent damage of soft rock from trans-scale modeling (Potyondy 2007;Li et al. 2017Li et al. , 2018Cui et al. 2019;Song et al. 2019;Guo et al. 2020;Wang and Cai 2021a, b;Xia et al. 2021a, b). However, the challenges of the computational efficiency and anisotropic nature of soft rock and rock joints need to be addressed before theoretical models can be implemented in engineering practice.
In the present study, a multiphase numerical model of soft rock was developed to investigate the creep-induced damage behavior of the soft rocks surrounding a tunnel after excavation through coupled meso/macroscale simulation. The anisotropic mechanical behavior of the rock due to joints was also considered. The numerical results of the developed model were validated by using the field monitoring data of a tunnel case study.

Methods
The discrete element method (DEM) can reproduce and characterize rock mechanical behavior. As one of the most widely used DEM software programs, Particle Flow Code (PFC) shows excellent performance in simulating the fabric features of rocks using multiple-grain models (Potyondy and Cundall 2004;Bahaaddini et al. 2013).
Rock deterioration is considered a micro-to-macro trans-scale process, where the accumulated internal mesoscale damage could lead to a macroscale strength loss. In this paper, we want to correlate mesoscopic damage to macroscopic large deformation. For rocks, microcracks initiate at the grain scale. The parallel bond model (PBM) implemented in PFC can simulate the cementation of spherical grains at contacts. Burger's model adopted in PFC can mimic the sliding creep of grains. The bond damage and grain creep are termed mesoscale damage herein. The mesoscale damage initiated from the grain contact area is limited at the mesoscale (less than 1 mm). The grain assembly can represent rock and rock mass, which are referred to as the macroscale (1-10 cm) and the engineering scale (> 1 m), respectively.

3
As shown in Fig. 1, the damage and deformation of the soft rock surrounding rock resulting from tunnel excavation can be studied through coupled meso/macroscale modeling. The mechanical properties of soft rock at the macroscale are characterized by bedrock strength, joint distribution, and tectonic stress, which are anisotropic in nature. The mechanical behavior of the rock at the mesoscale can be modeled by the skeleton and mineral bonding contents (Pinyol et al. 2007). In addition, the damage of the soft rock at the mesoscale can be divided into three phases: (1) initial phase before the excavation, (2) loading

Modeling the Macroscale Damage and Deformation of the Surrounding Soft Rock
Resulting from Tunnel Excavation Figure 2 shows the macroscale modeling of the surrounding soft rock with consideration of the anisotropic state of the rock. The stress state of the soft rock surrounding the tunnel can be defined as: For hydrostatic conditions, the principal stresses of the surrounding rock can be described as: where m = 1 3 I 1 = 1 3 1 + 2 + 3 is the average principal stress. The stress tensor comprises the tectonic and gravity stresses of the rock surrounding the tunnel, which is composed of hydrostatic and deviatoric stresses. That is, where s kk is the deviatoric stress tensor and ij is the Kronecker delta (if i = j , ij = 1; if i ≠ j, ij = 0). The damage process of the soft rock surrounding the tunnel can be treated as an irreversible dissipation process. The thermodynamic potential function (Helmholtz energy) of the surrounding rock can be expressed as a function of elastic strain ( e ) and internal variables, v k (k = 1, 2, 3 … , n) , under isothermal conditions. That is, The anisotropic damage characteristics of the rock can be defined using the isotropic hardening variable v p for plastic deformation, the damage variable D for anisotropic damage, and the scalar variable for the current state of damage. That is, In an engineering case, the rock mass has initial damage from rock mass joints (Fig. 1). Based on the strain equivalence principle, the strain caused by stress acting on damaged material in the stressed state is equivalent to the strain caused by effective stress acting on undamaged material, and we obtain: where E 0 and E 1 are the elastic moduli of intact and damaged rock, respectively. The damage variable D can be introduced as the decay of E 0 (Xue et al. 2022), Based on the irreversible thermodynamic theory (Zhou and Zhu 2010), it is assumed that the plastic deformation and damage evolution of soft rock are mutually independent. Thus, the Helmholtz free energy function of the surrounding soft rock can be described as: The term d ( ) can be expressed in a simple form using the internal variable , which represents damage accumulation and mesoscale damage evolution. Therefore, the damage potential function at the macroscale can be expressed as: where d is the parameter denoting the damage consistency parameter. If F d < 0 , the soft rock is in the elastic stage, and the stress tensor can be defined as Eq. (4b). In the damage stage, according to the non-associated rule of orthogonal flow, we obtain: where Q denotes the damage yield function. is determined by the damage potential function F d and the hardening parameter.

Modeling the Mesoscale Damage and Deformation of the Surrounding Soft Rock Resulting from Tunnel Excavation
To analyze the damage accumulation of soft rocks during tunnel excavation at the mesoscale, the following assumptions are made in this study: (1) The tunnel rock is composed of skeleton particles and cementations (Pinyol et al. 2007;Li and Wong 2016).
(2) There are two types of bonding in soft rock, i.e., (1) strong bonding for tension, shear and torsion resistance and (2) weak bonding for creep resistance. (3) Strong bonding fails when the bearing capacity is exceeded. Soft rocks experience creep when the slippage and dissipation of weak bonding occurs (Liu et al. 2020a, b).

Initial Phase Before the Excavation
Before excavation, the rock surrounding the tunnel is subject to tectonic in situ stress ij , which is a stable and inelastic state. It is assumed that the total stress ( ij ) can be expressed by the sum of the strong bonding stress ( F ij ) and weak bonding stress ( f ij ). That is, where kl * F ij is the anisotropic stress, which is expressed as kl * f ij . kl and kl are tensors revealing anisotropy determined by the actual physical state of the rock surrounding the tunnel.

Loading Phase Due to Excavation-Induced Stress Redistribution
Due to tunnel excavation, the mechanical equilibrium of the surrounding rock is disturbed. As shown in Fig. 1, the cementations between grains play a role in resisting deformation caused by stress, preventing the relative shearing, deviation, and rotation of skeleton grains. The mesoscale damage can be described as the bonding failure of rock gains due to tension or shear, as shown in Fig. 3.
The parallel bonded model (PBM) (Potyondy 2007;Xia et al. 2021a, b) is used for modeling the bonding failure of rock gains at the mesoscale mainly based on the binding characteristics of the rock gains, i.e., shear strength c , tensile strength c , normal stiffness k n , tangential stiffness k s , and bonding radius R . That is, where F n and F s are the normal and tangential forces, respectively. M n and M s are the normal and tangential moments, respectively. U n , U s , n , and s are the relative displacements and rotation angles in the normal and tangential directions, respectively. A , I , and J are the area, moment of inertia, and polar moment of the bonds section, respectively, which can be described as The maximum tensile stress and shear stress on the grain surface can be calculated using the following equations: Mesoscopic damage to soft rock tunnels occurs when ≥ c or ≥ c , as shown in Fig. 2. The values of c and c determine the mechanical capacity of cementation, which can be obtained through a calibration process based on the peak strength of rock samples.

Creep-Induced Damage Phase Leading to Large Deformation
Clay minerals in soft rock may expand and disintegrate when exposed to water or constant loading. This long-term process can finally lead to the stress redistribution and damage accumulation of the rock surrounding the tunnel (Kontogianni et al. 2006;Tang et al. 2018;Xiong et al. 2020;Zhang et al. 2021). Because weak bonding is unable to carry tension, shear, and torsion, only the normal pressure is considered in the creep model. In this study, the creep behavior of soft rock is modeled using Burger's model, which combines the  Figure 4 shows the major components of Burger's model. The details are shown in Table 1. Figure 5 shows the creep-damage model for simulating the large deformation and viscoelasticity damage of the soft rock surrounding the tunnel. It is assumed that the strong bonding is linear elastic before brittle failure, where E 3 is the linear elastic modulus.
The governing equations of Burger's model can be described as where (20) For a parallel connection, as shown in Fig. 5, the strain ( c ) and stress ( c ) of the coupled creep-damage model can be described as: where 3 and 3 are the strain and stress of the damage model. By combining Eqs. (20)-(24), the coupled creep-damage model can be derived as: At the initial phase ( t = 0, = 0 ), Eq. (24) becomes In addition, the initial conditions of strain, stress, and creep rate can be described as: Equation (26) can be rewritten as where A, B and C are parameters related to the properties of the soft rock. By solving the second-order nonhomogeneous linear Eq. (30), the roots of the eigen-equation (i.e., 1 and 2 ) can be obtained as: Then, the large deformation of the soft rock ( c ) can be obtained as: Fig. 4 Burger's creep model The coefficients of C 1 and C 2 in Eq. (31) can be calculated when substituting Eqs. (27)-(29) into (30), and we obtain:

Coupled Meso/Macroscale Modeling of Damage and Deformation of the Surrounding Rock Resulting from Tunnel Excavation
The coupled meso/macroscale process resulting in the large deformation of soft rock can be described as: and d cd ij =̇c . By considering the anisotropic mechanical behavior of the soft rock, Eq. (32) can be written as: where ij represents the anisotropic mechanical behavior of the rock due to the tectonic stress condition at the macroscale. That is, Substituting Eq. (12) into Eq. (36), the multiscale coupled damage-creep model of soft rock surrounding a tunnel can be obtained, As shown in Fig. 6, Eq. (38) can be solved using the DEM based on a previously developed cell model that correlates the mesoscale structure and mechanical properties of soft rock (Li and Wong 2016;Xia et al. 2021a, b). The blue and green particles represent the normal and creep grains, respectively. It is assumed that there is a relatively strong bonding between the blue particles, whereas the bonding strength between green particles is relatively weak. The particles in the representative elementary volume (REV) are randomly arranged. In the numerical calculation, grain Between two grains, the most common stress states are compressive or tensile forces applied to the normal direction of the bonding surface. The deformation between grains involves instantaneous and rheological deformation. The PBM mimicking strong bonding has five components ( k n , k s , k n , k s and c ). PBM deformation is determined by the normal ( k n , k n ) and shear ( k s , k s ) stiffnesses of springs. The damage is determined by the strength component c , which complies with the threshold shown in Fig. 3. The Burger's model simulates weak bonding by four components ( E 1 , E 2 , 1 and 2 ). Weak bonding creep appears in the dashpot component ( 1 , 2 ) when a normal stress is applied. E 1 and E 2 are elastic springs herein.

Project Description
The soft rock tunnel selected in this study is located in northwest China with a total length of approximately 15 km and a MPa. In addition, the tectonic stress measured by the original hydraulic fracturing shows that the maximum and minimum horizontal principal stresses are 20 and 14 MPa, respectively (Zhang et al. 2020). The deformation of the surrounding soft rock after excavation is relatively large, with the convergent deformation exceeding 1 m. A total station and 3D laser scanner were used to monitor the section of the No. 2 inclined shaft of the tunnel. The deformation of the surrounding rock mass was measured by a BLSS 3D laser scanner. The 3D spatial deformation was calculated by distance, angle, and inclination correction, which was recorded by a laser distance sensor, position feedback encoder and inclination sensor, respectively. The measuring accuracy of the BLSS 3D laser scanner is 0.7 mm in a 10-m range. The final deformation is shown in Fig. 7a, b.
The monitored displacement of points "A", "B", and "C" was measured by a Leica TCA2003 total station. When measuring, the reference point is set up at a stable place, and reflective sheet targets are attached to the top of the vault. Three measurements were taken at each point for the average, which was continuously monitored for 25 days during the excavation. The monitoring results of the tunnel section investigated are presented in Fig. 8. The measuring range and accuracy of the Leica TCA2003 total station are 1000 and 0.1 mm, respectively.

Model Setup and Calibration
The numerical modeling process involves the following five steps: (1) Define the simulation domain of the grain model As shown in Fig. 9a, the size of the research domain is 60 m (width) × 60 m (height) × 10 m (length), considering the size of the tunnel. "A", "B", and "C" are the three monitoring points in the tunnel face.

(B) Distribute the grain particles into the simulation domain
As shown in Fig. 9b, the grains are uniformly distributed between the boundaries of the domain, and the mechanical system is iterated until equilibrium is reached, in which fine particles are distributed in the tunnel zone (Gutierrez-Ch et al. 2022). Considering the modeling convenience, the fine particle zone is selected in a square region. The density and modulus of the grain assembly are 2500 kg/m 3 and 600 MPa, respectively, based on the rock samples obtained from the tunnel.
(C) Tectonic stresses applied to the domain As shown in Fig. 9c, the tectonic stresses according to the in situ stress condition are applied to the simulation domain based on the servo method (Lin et al. 2022). The stress is calculated based on the overall contact force between grains and walls, while the velocity of walls is controlled by FISh (a computer coding language commonly used for DEM) until the dynamic stress balance is reached. The in situ stress is determined based on the field measurement results, i.e., the maximum, intermediate and minimum principal stresses are 20, 14 and 14 MPa, respectively.

(D) Combined finite and discrete element approach
In the numerical simulation, the model length is usually set to nearly eight times the tunnel diameter to avoid boundary effects. However, limited by computing efficiency, the DEM grain-based model cannot meet all the requirements.  Thus, as shown in Fig. 9d, a combined DEM-FEM approach (Huang et al. 2020a, b;) was used in this paper for modeling the continuous deformation in the inner zone and expanding the simulation range. The inner FEM zone is used for the support material. In the DEM-FEM coupled zone, the calculated unbalanced force of grains will be attributed to the adjacent FEM nodes. This will cause a velocity alteration of FEM nodes, which will be assigned to the adjacent grains. The above coupling process is completed in each timestep. Before the model implementation, bonding models were added, and it is assumed that the cemented grains form a stable system. The bonding strength between the grain particles was calibrated based on the average strength of the intact rock mass (i.e., approximately 30 MPa). It is important to note that the already existing rock support could influence the initial equilibrium. After excavation, the support element was applied with a null modulus. The real lining structure is installed instantly after Day 1 when creep deformation begins.

(E) Application of initial joints and boundary conditions
As shown in Fig. 9e, the deformation due to excavation can be modeled using the commercial DEM software Particle Flow Code (PFC). The modeling process of the joints is realized by using the smooth joint model in PFC, which can provide relative sliding friction (Hu et al. 2018). After the application of joints, the initial equilibrium before excavation was obtained under the boundary conditions shown in Fig. 9e. After excavation, the displacements at monitoring points "A", "B" and "C" were numerically predicted.
When a mesoscale modeling approach is used to model the mechanical response of rocks, the mechanical parameters of grain elements as well as bonding elements need to be well calibrated (see Fig. 10). In this paper, the calibration process is complex for considering joints and creep conditions. Typically, the strength and elastic modulus of rock are determined by a 'trial and error' process (Coetzee 2017). In this paper, we calibrate the creep, strength, and   Fig. 10, which is presented in Table 2. The trials involved were repeated three times under different random seeds. Based on field investigation, the basic rock properties are obtained from (Zhang et al. 2020), and the time-dependent creep test results are collected from . The mesoscale parameters presented in Table 2 are determined according to the following steps: (1) General parameters adopted in the DEM modeling In the grain-based model, we need to consider calculation accuracy and efficiency. According to a previous study, when L/R (where L is the model length and R is the mean grain radius) is greater than 68, the properties of the material are nearly constant (Li et al. 2014;Cheng and Wong 2020). Previous studies have shown that in a DEM system, the number of mesoscale elements can meet the requirements of accuracy and calculation efficiency when greater than 15,000   (Lin et al. 2022). In the rock sample and rock mass calibration, we choose an L/R value greater than 100 in this paper. A total of 31,516 grains were used in the simulation based on convergence analysis. The damping coefficient is set to a default value of 0.5 to dissipate kinetic energy.

(B) Loading test calibration
Static mechanical properties were calibrated prior to the creep deformation parameters. The loading test is the most commonly adopted method to characterize the elastic modulus and strength of rock samples. We conducted numerical tests to obtain reasonable bonding strength and stiffness. The calibration results are shown in Fig. 11a, where we find that the mechanical properties are well-calibrated compared to the real rock parameters. We compared the detailed mechanical properties of the rock samples in Table 3.

(C) Creep test calibration
Because of the rheological deforming scenario discussed in this paper, we need to consider two kinds of contacts as bonding. The ratio of PBM to Burger's model and their parameters will affect the final creep behavior of the rock. In this stage, the bonding strength is applied to a high value (tensile strength = 200 MPa, cohesive strength = 400 MPa), which will not damage the creep tests. The calibration process is finished by changing the parameter values after the 'trial and error' process.
The creep tests were carried out at 12, 24 and 36 MPa stresses (see Fig. 11). During the creep simulation process, we found that the results were reasonable when equating 5 × 10 5 steps to a real 24 h in the laboratory. We keep the same treatment in this paper in regard to the conversion of timesteps to real time. In the calibrated sample, creep grains make up 18.7% of the total grains, which are randomly distributed in the particle assemblage. The final calibration results are presented in Fig. 11b, c, where we can compare the creep curves. The detailed calibration error is shown in Table 4, and the mean error of creep strain in different loading forces is less than 5%.
(D) Verification on the engineering scale According to the mentioned conclusion, when the model size to grain size ratio is larger than 68, the mechanical behavior of the DEM model is constant. It is necessary to verify this conclusion to promote mesoscale parameters to an engineering scale. Therefore, we conducted an engineering scale (15 m × 15 m × 30 m) uniaxial loading test via the DEM. We find that when the model size to grain size ratio is fixed, the strength and elastic modulus of the rock are also nearly constant (see Fig. 11d).

(E) Bedding joint parameter calibration for rock mass
After rock sample calibration, we need to further determine the rock mass properties, where the joints are considered. Numerical back-analysis is a powerful method that can be used as a complementary technique to in situ or laboratory experiments to determine mechanical properties (Fakhimi et al. 2004;Li et al. 2022a, b). In situ failure, movement or convergence can be used to back-calculate the physical properties of rock joints. In this paper, the following steps were used for back-analysis: A. Strength of 29.6 MPa and elastic modulus of 7.6 GPa were used as determined from the laboratory loading test. B. A joint angle of 45 • and a spacing of 15 m were assumed. The value of angle and space was justified based on field mapping. C. Reasonable values for joint roughness, joint cohesive strength, and friction coefficients based on the results of in situ tests and engineering judgment were assumed. We modified these parameters until there was good agreement between the numerical and measured in situ failures. After the back-analysis process, we obtained a reasonable rock mass numerical mechanical property compared with the rock mass property (E = 0.6 GPa, ci =7.4 MPa) obtained from field investigation (Zhang et al. 2020).  Figure 12 shows the macroscopic deformation of the surrounding rock during the excavation. It shows that the predicted time-dependent displacements at "A", "B" and "C" agree with the monitoring data reasonably well.

Comparison Between the Predicted Displacement and Monitoring Data
In addition, the displacement at point "B" is larger than that at "A" and "C". For example, on Day 25 after the excavation, the displacement at point "B" is more than twofold larger than that at point "C". Furthermore, the maximum displacement surrounding the surface of the tunnel is approximately 526 mm, which could lead to tunnel damage.

The Asymmetrical Deformation of the Tunnel Face at Monitoring Points "B" and "C"
The deformation of the rocks surrounding the tunnel face is generally asymmetric due to the application of joints with different inclination angles (i.e., 30°, 45°, and 60°), and the asymmetric nature of the deformation becomes obvious at a joint inclination angle of 45° (Zhang et al. 2020). As shown in Fig. 13, the deformations at points "B" and "C" are symmetrical without initial joints, whereas the application of initial joints could lead to a significant difference in displacement between "B" and "C". Under the joint inclination angle of 45°, the rock mass is prone to slide along the joint surface, which ultimately could lead to a large displacement on the left side of the tunnel. Figure 14 demonstrates the movement of rock grains, which is captured from the initial equilibrium stage at cycle time step 10,000. The movement of rock grains demonstrates a symmetrical "X" shape, whereas the application of initial joints could result in the location of grains along the joints and two shearing zones. As monitoring point "B" is in one of the shearing zones, it has the largest displacement among the three measuring points. Ultimate displacements/ mm Predicted ultimate displacements at "A", "B" and "C" Model prediction with initial joints Model prediction without initial joints Fig. 13 The ultimate displacements at points "A", "B" and "C" with or without initial joints

Mesoscopic Scale Damage and Deformation of the Surrounding Rock at Different Phases
As shown in Fig. 15, the large deformation of soft rock surrounding the tunnel is a gradual process. In the initial phase, the whole surrounding rock is stable. Due to the excavation-induced stress redistribution, the surrounding rock undergoes the loading phase, in which the weak bonds in the rock experience creep under excessive local stresses. As the deformation accumulates, creepinduced damage develops around the tunnel. The coupled creep-damage process ultimately causes large deformation of the surrounding rock of the tunnel. From Days 1-4, the excavation-caused stress concentration around the tunnel causes severe mesocrack accumulation; the creep-induced stable increase in crack number is captured thereafter.

Model Sensitivity Analysis
Given the complexity of the scenario considered in this paper, we conducted a sensitivity analysis of the proposed model.
(1) Effect of the change in tensile/cohesive strength on the final displacement Figure 16 shows the predicted ultimate displacement at points "A", "B" and "C" under different bonding strengths of cementations quantified by the mesoscopic tensile and cohesive strength of the rock grains. The initial values of the bonding characteristics ( c =30 MPa, c = 61 MPa) are obtained through the calibration process based on the peak strength of the rock samples. It shows that the decrease of the bonding tensile strength by 50% could increase the displacement at the top of the tunnel (point "A") by 300% but have relatively little influence on the displacements at points "B" and "C". However, the same proportional change in cohesive strength has a smaller effect on the final deformation. Thus, the final deformation is more sensitive to the tensile strength.
(B) Effect of the fraction of grains experiencing creep on the model mechanical responses In this paper, we have coupled the damage and creep factors in a soft rock tunnel, and we also obtained a reasonable result. We refer to the grains that undergo creep as creep grains. We extended our discussion to some more scenarios with different creep grain ratios involving 18.7, 20.2, 21.7, 23.2 and 24.7% (Fig. 17). The rock strength, creep strain, and excavation deformation at three monitoring points "A" "B" and "C" are compared in Table 5 and Fig. 18 for sensitivity analysis. When the creep grain ratio rises to 24.7%, the error of the predicted deformation at the three monitoring points becomes large (greater than 10%). If this ratio continues to increase, the deformation prediction results obtained will differ significantly from the actual situation.
In a time-dependent comparison, we conducted creep tests for laboratory-scale rock samples under 24 MPa stress. Figure 17 shows that the ultimate creep strain and deformation at the left side of tunnel monitor point "B" increase with increasing creep grain ratio. When the incremental creep grain ratio increases from 4.5 to 6.0%, the incremental deformation at the three monitoring points nearly doubles.

Conclusions
We implemented a damage-based, time-dependent modeling approach to reproduce the entire evolution of a large-deformation soft rock tunnel. By incorporating the coupled meso/macroscale deformation process of the rock into the model, the predicted time-dependent displacements of the tunnel face agree reasonably well with the monitoring data.    The model is able to reproduce both instantaneous damage patterns of brittle failure and time-delayed creep deforming patterns. Tunnel excavation strongly affects the initial stability of the surrounding rock mass, and brittle damage accumulated in mineral clusters leads to severe instantaneous deformation, which becomes less important in the creep evolution stage.
The application of the theoretical model to engineering cases provided the opportunity to investigate the evolution mechanisms of large deformation. As an emerging property, the model simulates a wider range of scenarios with different rock static and creep parameters. The final deformation is characterized by a high sensitivity to the value of mesoscale modeling parameters.