A nonlinear elasto‐viscoplastic model for clayed rock and its application to stability analysis of nuclear waste repository

The creep behavior is one of the most important topics on the stability analysis of rock mass in underground engineering. In order to characterize the creep deformation of clayey rock more accurately, this study proposes a nonlinear elasto‐viscoplastic (EVP) model based on a new yield surface and a modified Mohr‐Coulomb (MC) flow potential. First, the finite element (FE) formulation for the constitutive creep behavior is derived and a user subroutine to define the material's mechanical behavior (UMAT) is implemented in the commercial finite element package ABAQUS. The numerical simulation is validated against the triaxial creep tests of a cylinder, and the results show that the proposed model can effectively predict the creep responses of clayey rock and thus compensate for the deficiency of the creep model in ABAQUS. Then, the proposed model is used to study the long‐term stability of clayey rock formation surrounding an underground experimental site. Based on the results from triaxial creep tests and in situ monitoring of the lining deformation, the time‐dependent material creep parameters are obtained for rock by using the back analysis method. The predicted lining deformation is in good agreement with the field data with most of the relative errors less than 3.5%. The creep zone in the surrounding rock demonstrates that the horizontal scope is larger than the vertical one. The proposed model is proven to provide accurate prediction on the creep characteristics of clayey rock, showing a potential to assess the long‐term stability and reliability of large‐scale underground engineering structures.


| INTRODUCTION
In the framework of geological disposal of radioactive waste, clayey rocks are considered to be a potential geological barrier. Indeed, the confining capacity of this geological barrier is directly related to the time-dependent properties of clayey rock, in particular the creep deformation and damage evolution of the surrounding rock. The creep of soft rock or loose rock mass arises from the evolution of microstructures such as mud fillers and interlayer fragment zones. [1][2][3][4] A number of evidences show that stressed rocks undergo creep deformation, a slow movement of mass, and that the rock mechanical properties become time-dependent. 1,5,6 The hard rock has similar creep properties under high geostress conditions. 7 With the rapid development of expressway construction and the extraction of resources from deep mines, more and more | 151 JIA et Al. tunnels are being constructed under high geostress. The plastic deformation and nonlinear creep of deep rocks are so remarkable that they are the main problems encountered and concerned in deep resource mining. 6,8,9 Therefore, a good understanding of the creep characteristics of deep rock or soft rock and their impact on the long-term stability of a subsurface structure is very important.
One aim of the study on rheological properties of material is to develop a proper creep constitutive model based on the experimental results or in situ monitoring data. There are generally two ways to achieve this objective for rocks. The first one involves microscopic studies by using modern devices such as X-Ray CT and scanning electron microscope (SEM) to account for the microstructural changes in rock. 10 The other one uses macroscopic observation from the triaxial creep tests and in situ creep measurement to build the elastoviscoplastic (EVP) models. 5,9,11,12 As it is difficult to control the modern devices for long-time testing (say several years), the macroscopic approach is still widely used in studying the creep behavior of rock material.
Time-dependent properties and related applications in rock material have received a wide range of attention during the past several decades. Generally, the time-dependent mechanical models for analyzing the creep behavior of rock include mainly the following four types, 1,5,[13][14][15][16] that is, the component rheological (CR) model, yield surface rheological (YSR) model, empirical rheological (ER) model, and endochronic rheological constitutive (ERC) model. Because of simple concepts such as elastic, viscous, and plastic strains, the CR model is the most popular one in describing the rheological properties of rock. However, its deficiencies involve fast creep accelerating and static yield surface. 5,6 The YSR model uses dynamic yield surface which changes with time. The ER model can characterize the actual creep behaviors under different types of loading paths, but it requires many creep tests to calibrate. The ERC model is able to deal with both static and dynamic rheological behaviors under cyclic loading or vibration. However, due to complex programing implementation, it is less used than other models. In general, a good creep constitutive model should precisely describe the time-dependent behaviors such as failure of rock material. Among the above models, the ER model is the most widely used in the underground operations because of easy implementation and accurate prediction of the creep curves of rocks. 2,5,7,10,13 In general, the time-dependent inelastic deformation of rocks is described by the phenomenological viscoplastic theory. [17][18][19] The time-dependent deformation is attributed to the inherent viscous effect of material. On the modeling of creep deformation in engineering, the viscoplastic theory is largely used. Different viscoplastic models, based on either thermodynamic potential or overstress approach, have been developed, including overstress theory, 20 equivalent timeline model, 21 viscoplastic degradation model, 22 and internal state model. 23 However, although these models can provide powerful mathematical description of creep deformation, the time-dependent mechanical responses are not reproduced completely and identified clearly. Further, the viscoplastic deformation is described independently with the plastic deformation. Two distinct formulations will be needed, respectively, for the shortterm behavior (elastoplastic model) and for the long-term behavior (viscoplastic model). The Mohr-Coulomb (MC) failure criterion has been widely used for geotechnical applications. A large number of routine calculations are still performed for design using the MC criterion. To represent the viscoplastic characteristic of rock material, the MC criterion is combined with time-dependent models to simulate the creeping behavior. 4,6,24,25 Creep and plasticity of rock materials can be treated simultaneously, and the relevant governing equations can be solved in a coupled manner. 12,19,26,27 Because the creep in such materials is intimately linked to the plasticity behavior, it is necessary to define the creep flow potential through experimental or field data and plastic parameters such as cohesion, friction angle, and dilation angle. Classical creep flow potential for rock materials is defined mainly based on von Mises, 2,28 Drucker-Prager, and Cap models. 10,29 However, fewer numerical studies were carried out to describe the creep deformation based on the MC potential flow, although a number of laboratory tests and in situ measurements have proved that the MC criterion is more reasonable for geotechnical applications.
Clayey rock mechanics for engineering applications requires powerful constitutive modeling, and creep is one of the main deformation mechanisms. The creep of clayey rock is a highly nonlinear mechanical process closely related to complex stress conditions and plastic deformation. Based on the concept of equivalent stress surface, we propose a new EVP model using the MC flow potential. To enlarge of the scope of the creep models used in ABAQUS, a UMAT subroutine is developed and incorporated. This paper is organized as follows. Section 2 presents several mechanical models for equivalent stress surface used in failure analysis; the finite element formulation and solution method will be detailed in Section 3, which is followed by computational procedure in Section 4; the proposed model is validated in Section 5 and applied to the stability analysis of underground experimental site in Section 6; and in the end, brief conclusions are given in Section 7.

EQUIVALENT STRESS SURFACE
The Mises potential function is used by most FE models to describe the creep strain rate of rock material 2,13,28,30,31 that meet the nonassociated flow rule condition without consideration of dilation angle, but inappropriate for rocks which obey the associated flow rule or when the dilation angle is not equal to 0°.

| Drucker-Prager criterion
The Drucker-Prager (DP) criterion is often used for predicting the yield surfaces of rocks that exhibit pressure-dependent behaviors and has the form 27 where c and ϕ denote the cohesion and friction angle, respectively; p is the hydrostatic stress; and q is the modified von Mises stress defined as where q is the von Mises stress, J 3 is the third invariant of the deviatoric stress, and the parameter K (0.778 ≤ K ≤ 1) is used to match different stress values in tension and compression in the deviatoric plane.
The plastic/creep flow potential, which is continuous and smooth, is used to ensure that the strain rate direction in the strain spaces is uniquely defined. The flow potential function, G, in the DP criterion is defined as 29 where φ is the dilation angle, σ 0 is the initial equivalent yield stress, and denotes the eccentricity that defines the rate at which the potential function approaches the asymptote. It should be noted that the flow potential curve becomes a straight line when the eccentricity is equal to zero.
The creep iso-surfaces (or equivalent creep surfaces) are those where the stresses have the same creep intensity. When the rock is in plastic deformation, the creep iso-surface should coincide with the yield surface. Therefore, we can define the creep iso-surfaces for a specified stress state by homogeneously scaling down the yield surface in the p-q plane ( Figure 1). The equivalent creep stress, , in the DP criteria is expressed as The DP model becomes the Mises model when ϕ = φ = 0 and K = 1 in Equations (1)-(4).

| Modified hyperbolic MC criterion
The creep models based on the DP and Cap criteria for many rock materials have been provided by the nonlinear FE software ABAQUS. 27 However, this software does not include the case based on the MC criterion.
The modified MC criterion considers the actual tensile strength of rock, thus overcoming the limitation of the classical MC criterion which overestimates the tensile properties of rock. 32 Therefore, the proposed yield surface in this study is determined by two different criteria, that is, one shear criterion for the MC surface and one optional tension cutoff criterion for the Rankine surface. The yield function, F, has the form ( Figure 2).
The MC yield surface in the three-dimensional space is a hexagonal pyramid, which, due to gradient discontinuity at the six vertices, imposes difficulties for the numerical calculation. A suitable choice for M(θ), which rounds the vertices of MC surface at θ = ±30° in the octahedral plane, can make the numerical convergence easily ( Figure  3). In the vicinity of the vertices, where | | > T and T is a specified transition angle within the range of 0°-30°, an alternative form of M(θ) is defined, in which larger values of T giving better fit to the MC cross-section in the octahedral plane. In this study, a piecewise function, M � ( ), is chosen and adopted to make the modified MC yield surface close to the original one. The yield function is treated smoothly and continuously at the edges and corners. When 25° ≤ |θ| ≤ 30°, the function denoting the modified MC yield surface takes the form 6,33 where the coefficients A and B are defined as where θ T should not be too near 30° to avoid ill-conditioning of the approximation and a suitable value is 25° and the function sign(θ) Similarly, the equivalent creep stress in the MC criterion is written as follows ( Figure 4) where R mc = M′(θ)/cos ϕ. Figure 4 shows the creep iso-surface in the p-q space, without creep deformation because the equivalent creep stress is negative.
According to the modified MC yield function, the plastic/ creep potential function is written as Indeed, the classical MC function can be recovered by substituting t = c cot and θ T = 30°.

| Creep strains
The creep potential function is a smooth curve and guarantees the uniqueness of the creep flow direction. The creep strain increment, dε c , is given by 29 where dε c is the increment of creep strain vector, and d c is the increment of equivalent creep strain. For example, in the case of uniaxial creep, d c = | | d c,11 | | , here c,11 is the axial creep strain under the uniaxial loading condition. σ denotes the stress tensor, and the function f cr is defined as At present, many achievements have been made in the study of creep model for soft rock and the existing creep models are mainly divided into empirical model and component model. Among them, empirical model is widely used because of its pertinence, flexibility, and few model parameters. 34 The empirical creep equations, including time hardening, strain hardening, and Singh-Mitchell types, 5,6 can be chosen according to the creep deformation characteristics, as listed in Table 1.
The creep strain can be determined by using the functional relationship among equivalent creep strain rate ̇c , equivalent creep stress , and time t from the test data. 5 In particular, the time hardening equation can be defined for the equivalent creep strain of clayey rock and its rate where c is the total equivalent creep strain; ̇c is the rate of equivalent creep strain; t is time; and C 1 , C 2 , and C 3 are the creep parameters obtained by fitting the test data.

| Basic equations
The stress-strain relationship is expressed in incremental form 27 where D e ijkl denotes the elastic matrix, dσ ij and d e kl denote the increment of stress and elastic strain, respectively; , and d c kl denote the increments of the corresponding total strain, plastic strain, and creep strain, respectively. The plastic strain increments is defined as where d is the plastic multiplier, d = d pl g , g = 1 c σ: G σ , pl is equivalent plastic strain. According to the virtual displacement principle, the finite element formulation is obtained in increment form 6 where K e is the elastic stiffness matrix, Δu is the displacement increment array, ΔQ is the external load increment array, B is the geometric matrix, and D e is the elastic matrix.
Given the total strain increment, Δε, both the plastic strain increment, Δε p , and the creep strain increment, Δε c , are unknown and need to be determined by using iterative method. The iterative equation is written as where Δε (n) p and Δε (n) c are the increments of ε p and ε c , respectively, after n iterations at the current time step, and δu (n+1) is the value of Δu after (n + 1) iterations at the current time step.

| Calculation of the stress and strain increment at the Gauss integral points
Once the value of δu (n) for the displacement increment ∆u is obtained from Equation (17), the modified strain δε (n) and the total stain increment Δε (n+1) = Δε (n) + δε (n) can be obtained by using the geometric equation. Prior to the next step, the new values for , and Δσ (n+1) at every Gauss integral point can be determined through the following nonlinear equations where Δt is the time step.

Name Creep equation Parameter
Time hardeninġc = C 1 Creep parameters; T Temperature; Total strain Modified time hardeninġc = C 1

| Schemes of iterative solution
The creep analysis can be carried out by using explicit or implicit integration schemes. The explicit integration scheme is stable, but the time step must satisfy the stability condition.
The implicit scheme has the characteristic of unconditional stability, but requires a long computation time at each time step, which can be compensated by increasing the time step. The nonlinear equations in Equation (18) can be rewritten in iteration form where D ep is the elastoplastic matrix, k represents the iteration number and where ϑ is a time integral factor. When ϑ = 0, ϑ = 1, and ϑ = 1/2, the iteration scheme is an explicit, implicit, and implicit gradient method, respectively. The latter two methods have the characteristics of unconditional stability.

| Selection of time step
By comparing Δε c(k+1) with Δε c(k) between two iteration steps, the convergence criteria are defined as where err is a tolerance.
If the above convergence criterion is met, the creep and plastic strain increments can be updated at the (n + 1)th step If the calculated creep strain does not meet the above convergence criterion, Δε c(k+1) is substituted into Equation (19), and the next iteration continues until the convergence criteria are met. When all integration points meet the convergence criterion and the iteration calculation is completed at the (n + 1)th step, Δε (n+1) c and Δε (n+1) p are obtained by During the iteration process, the time step has a large effect on the convergence. The time step for stable convergence can be defined as 6 where el | |t and ̇c | |t denote the equivalent elastic strain and creep strain rate, respectively, at time t.
The initial creep stage is unstable with the creep rate changing dramatically and thus leading to significant stress redistribution. Therefore, the time step should be small enough for the early time and then increase gradually depending on the stress levels and the constants for creep properties.

| COMPUTATIONAL PROCEDURE
The communication between the UMAT creep subroutine and ABAQUS main program is described as follows. 27 First, the main program implements the equilibrium equation solver at time t n and obtains the new value δu (n) of displacement increment Δu. The strain increment Δε is obtained by using geometric equations. Then, the main program provides the strain increment Δε to the UMAT subroutine, which calculates the stress increment t n +Δt σ by solving the constitutive equations. The accuracy of Jacobian matrix affects the convergence rate of the UMAT subroutine. After that, the UMAT subroutine will pass the new stress increment t n +Δt σ on to ABAQUS main program to carry out the next equilibrium iteration at time t n+1 .
The UMAT creep subroutine works in the following way: 1. Based on the results at time t n , including the stress tensor t n σ, total strain t n ε, total strain increment t n Δε, and time increment Δt, provided by the ABAQUS main program, the UMAT subroutine calculates the new stress tensor, t n +Δt σ according to the constitutive equations. 2. Update the stress tensor in the UMAT subroutine. The creep strain increment is calculated by using the first equation of Equation (18). Then, the elastic-plastic strain increment and stress increment can be obtained by using Δε ep = Δε − Δε c and Δσ = JΔε ep , respectively. Here, J is the Jacobian matrix of the constitutive model. The stress tensor is updated for the next time step, that is, t n +Δt σ = t n σ +Δσ. 3. The strain tensor is updated by using t n +Δt ε = t n ε +Δε for the next time step. The UMAT subroutine provides the updated stress and strain tensors to ABAQUS main program. 4. The ABAQUS main program implements the iteration method for force equilibrium using the formula of Newton-Raphson at time t n+1 . If the current solution is convergent, the main program provides the stress tensor, total strain, total strain increment, and time increment to the UMAT subroutine and performs the iterative calculation. If the current solution is not convergent, the ABAQUS main program will decrease the time step Δt and continue the iteration until convergence is reached.
The UMAT creep subroutine includes the elastic-plastic analysis module, creep strain analysis module, and stress and strain increment modules. The flowchart of the UMAT subroutine for the viscoelasto-plastic model is shown in Figure 5.

| MODEL VALIDATION
In order to validate the present model, a comparison is made for the mechanical responses of rock during the conventional triaxial creep test between three models which are based on Drucker-Prager (DP), Mohr-Coulomb (MC), and Mises creep flow potentials. The elastic modulus, Poisson's ratio, friction angle, dilation angle, and cohesion of the rock are 300 MPa, 0.25, 18°, 18°, and 0.3 MPa, respectively. It should be mentioned that, when the Mises creep potential model is applied, both the dilation angle and friction angle are equal to 0. The time hardening formulation is adopted to describe the creep response of rock, with the material parameters in Equation (13) chosen to be C 1 = 5.02 × 10 −3 , C 2 = 3, and C 3 = −0.9.
In Figure 6, the axial compressive loading imposed on the cylinder is 0.5 MPa. Two cases, that is, one with only uniaxial loading and the other with axial loading and radial confining pressure of 0.1 MPa, are analyzed. Figures 7 and 8 show the creep curves under the two loading conditions. It can be seen that the creep results based on the Mises model are much larger than those obtained from the other two models. In the Mises model, the equivalent creep stress is equal to the Mises stress q and the creep strain exists when q > 0 according to Equation (13). The creep zone (q > 0) for Mises model in p-q space is larger than that of DP ( Figure 1) and MC models ( Figure 4). As for the creep iso-surface in the p-q space, there exists a noncreep zone for the DP and MC models when the equivalent creep stress is negative. However, in this noncreep zone, the Mises stress, q, is positive. Therefore, the Mises model overestimates the actual creep deformation of rock material and thus is not suitable for clayey rock. This model test shows that the creep strain evolution based on the MC flow potential is similar to that based on the DP creep flow potential under different stress conditions due to similarity of creep iso-surfaces and creep zones in both types, as shown in Figures 1 and 4. The creep strain obtained from the MC model is slightly larger than that from the DP model in uniaxial compression case, but the difference tends to increase when the confining pressure increases because of different potential functions and equivalent creep stress used. In the deviatoric plane (as shown in Figure 9), the DP model represents the circumcircle of MC model, and the yield surface and equivalent stress surface in the former are larger than those in the latter. This means the creep zone for MC model is slightly larger than that for the DP model.
Due to efficient determination of material parameters and easy implementation for geotechnical application, the MC model has been widely used for rock material. According to many experimental studies, 6,33,35,36 the MC model reflects properly the elasto-plastic mechanical behavior of clayey rock. Considering that the creep in such rock material is intimately associated with the plasticity behavior, it is necessary to apply a consistent potential function for viscoplastic behavior of clayey rock. Therefore, the creep UMAT subroutine describing the creep developed based on the Mohr-Coulomb criterion is feasible.

| Site overview
With the continuous development of nuclear plants over the past several decades, a large amount of radioactive waste needs to be buried underground. An intensive study on the time-dependent behaviors of rock formation is required for the safety and stability evaluation of the target repositories. Due to its homogeneous property and low permeability, clayey rock is regarded as one of the repository alternatives to other rock type such as granite and salt rock by many F I G U R E 6 Calculation model for the creep test countries. 9,30,[36][37][38] In this study, an underground experimental site in clayey rock is selected as example, which has a depth of 222.9 m, and the coefficient of the lateral earth pressure is 0.82. 6 A schematic view of the underground experimental site is provided in Figure 10. [36][37][38] It started by drilling an access shaft, and then, the test drift is tunneled by artificial excavation, with the excavation radius and length of the drift being 2.5 m and 67 m, respectively; the outer radius and thickness of the lining are 1.175 m and 0.6 m, respectively. To minimize the convergence, the connecting gallery was constructed by shield tunneling method. The deformation through diametrical convergence of lining 15, 29, 43, 71, 83, and 105 in test drift is monitored for about 10 years ( Figure  11). The creep characteristics of clayey rock will be studied by using the Mohr-Coulomb potential model.

| Determination of the creep parameters
The data used to determine the creep parameters are from two triaxial creep tests, that is, Test 1 and Test 2, which were performed under drained conditions. 6 The confining pressure σ 3 is loaded at 10 −5 MPa/s up to a set value, and then the deviatoric stress, σ 1 -σ 3 , is increased step by step and kept constant during the creep phase. In Test 1, the deviatoric stress is first increased to 1 MPa and then to 1.5 MPa. The confining pressure, deviatoric stress, and axial and radial strains are plotted in Figure 12. It can be seen that the creep strains are small when the deviatoric stress is 1 MPa, but cannot be neglected when the deviatoric stress becomes larger (say 1.5 MPa). In Test 2, the deviatoric stress is first increased to 1.5 MPa and then to 3.0 MPa. The variation of the axial strain with time under the confining pressure of 4.5 MPa is shown in Figure  13. The creep deformation is important within the testing time, and the deformation rate becomes larger when increasing the deviatoric stress. The creep parameters cannot be determined through the triaxial creep tests directly. However, they can be determined by back analysis method. A key to back analysis is to establish a reasonable error function, and it can be viewed as an optimization problem that the error function as the objective A combination of the mixed penalty function method and the Nelder-Mead algorithm is applied. Based on the error function above, the finite element program ABAQUS is used as the solver for the corresponding back analysis. 39 The triaxial creep test is simulated by using the Mohr-Coulomb criterion and time hardening creep model. The numerical model is equivalent to the real rock specimen, and the loading procedure is consistent with real test process in the triaxial creep test. By using the above test and back analysis method, the creep parameters in Equation (13) for clayey rock can be obtained, that is, C 1 = 1.1 × 10 −4 , C 2 = 3, and C 3 = −0.72. Figure 14 compares the creep strains of clayey rock between the measured values and the simulation results. It can be found that they are in good agreement. This means that the proposed Mohr-Coulomb creep model can provide a good prediction of the creep characteristics of clayey rock.

| Numerical model
The geometrical model for the finite element analysis is shown in Figure 15. The length and height of the studied domain are 100 m and 80 m, respectively. Initial vertical effective stress is 2.25 MPa, and pore pressure is 2.25 MPa. The laboratory and field tests demonstrate the permeability of clayey rock is 3 × 10 −19 m 2 and the porosity is 0.39. [36][37][38] The triaxial tests were performed with a confining pressure between 0.89 MPa and 5.42 MPa to estimate the mechanical properties of clayey rock. 6,40 The lining is assumed to be elastic, and its stiffness is reduced because of the joints of lining blocks and plates. 6,38 As its permeability is one order of magnitude larger than that of clayey rock, the drainage of lining must be considered. The mechanical parameters of clayey rock and lining are listed in Table 2.
The laboratory test showed that the concrete material of lining has obvious nonlinear viscoelasticity. 41,42 According to the results of creep tests on concrete material under a confining pressure close to the initial in situ stress, 6 the elastic modulus of concrete material is defined as a function of time where E 0 is the initial elastic modulus and t is time with a unit of day.
The main computational steps in this study are as follows. First, the equilibrium equation is solved under in situ stress field; second, tunnel excavation is realized by element death technique and the lining element is added by element birth technique; third, the extra-excavation element is added; and finally, the calculation is performed.

| Numerical results and discussion
In Figure 11, the diametrical convergence of L15 and L43 is the largest and smallest, respectively, among the linings. The diametrical convergence data of L29 are the most representative and thus selected as the monitoring data for this study. The time-dependent parameters, including C 1 , C 2 , and C 3 in Equation (13), can be determined by back analysis method by using finite element results and the measured convergence values. The accumulated relative error for the diametrical convergence of L29 between the numerical results and measured data is selected as the objective function, which is defined as where d t c and d n c are the tested results and numerical results of diametrical convergence of lining, respectively.
During the back analysis, the objective function gradually approaches a minimum with a convergence tolerance of 5%. The initial values of the parameters C 1 , C 2 , and C 3 for back analysis are set according to the preliminary estimate in Section 6.2, as listed in Table 3 as well as the range of these parameters.
The creep parameters obtained from the back analysis are shown in Table 3. By comparing them with those from laboratory tests, it can be found that the difference in the parameter C 2 is small while the difference in C 1 and C 3 is large. Determination of these material creep parameters from in situ back analysis is dependent on the excavation disturbance, lining deformation, seepage condition, and the surrounding rock properties. Figure 16 compares the diametrical convergence at L29 obtained from the back analysis and test results. They show an excellent agreement with a relative error less than 3.5%, indicating that the proposed model has the potential of capturing effectively the creep characteristics of clayey rock.
The creep strain evolution at three points, that is, A, B, and C, on the inner wall of L29 is shown in Figure 17. The creep strain at location B is much larger than that at locations A and C, and the creep strain at the bottom (point C) is slightly larger than that at the top (point A), which is caused by overexcavation and the lateral earth pressure. The maximum creep strain of the surrounding rock is 3% after 10 years. As the stress conditions in this case are close to those of Test 2 mentioned in Section 6.2, we can find that the creep strain evolution at the inner wall show similar trend as that in Figure 14, where the creep strain increases quickly in the early stage and gradually tends to be steady. Figure 18 shows the pore pressure variation with time in the surrounding rock. The early change of pore pressure is significant, especially near the inner wall of test draft, redistributing the stresses and thus affecting the deformation behavior of surrounding rock. With increasing time, the seepage field in the rock formation gradually tends to a steady state. The pore pressure at the inner wall of test draft declines slowly with time because of the drainage of lining. Figures 19 and 20 show the corresponding creep zone and maximum creep variation with time in the surrounding rock. It can be found that the creep strain and disturbance zone increase gradually with time in a butterfly shape and the effect of creep on the deformation of surrounding rock is important. The creep zone size in the horizontal direction is significantly larger than in the vertical direction, which is caused by overexcavation, the lateral earth pressure, and stress redistribution of surrounding rock during lining drainage and convergence deformation.

| CONCLUSIONS
In this work, a nonlinear visco-elastoplastic creep model, which considers the creep potential and equivalent creep stress of the modified Mohr-Coulomb (MC) criterion, is proposed for clayey rock. The finite element formulations are derived for creep deformation and incorporated into the ABAQUS software. A comparison is made for the conventional triaxial creep tests between the numerical models which are based on Drucker-Prager (DP), MC, and Mises creep flow potentials. In the end, the proposed model is applied to creep behavior of underground experimental site in clayey rock. Some conclusions are drawn: 1. The Mises model overestimates the actual creep behaivor of rock material, and the creep strain calculated using the MC model is slightly larger than that from the DP model. 2. For rock material that the creep behavior is coupled with plasticity behavior, it is necessary to apply the consistent potential function for viscoplastic behavior of rock. There is no conflict of interests exists in this paper.

DATA AVAILABILITY STATEMENT
The data used to support the findings of this study are included within the article.