# Modified Rheological Model for Deep-Buried Silty Mudstone and Support Time Analysis Application

Authors: Ankui Hu; Shuai Zhou; Mengkun Wu; Xiaohui Liu; Rui Zhao
Journal: Mathematical Problems in Engineering
 (2023)

---

## Abstract

The rheological properties of soft rocks should be considered in the long-term design and maintenance of deep-buried tunnels using uniaxial single-stage loading and graded incremental cyclic-loading methods. In this article, creep tests were performed on deep-buried silty mudstone from a specific water conveyance tunnel in China, with a buried depth in the range of 1650–2320 m and subjected to high in situ stress. The creep curves of silty mudstone under different loading stresses were obtained, showing evident rheological mechanical behavior under complicated external environmental conditions. Based on the classic Burgers rheological model, a new nonlinear creep model was established based on the creep properties of deep-buried silty mudstone in the project area. Typically, the designated rheological models for certain projects are unsuitable or inadequate. A nonlinear dashpot was calculated using the Levenberg–Marquardt (L–M) method coupled with origin to account for the deterioration trend in the strength of the silty mudstone over time. With the determined parameters, the modified Burgers model exhibited good qualitative consistency with field monitoring data. The user-defined material mechanical behavior (UMAT) subroutine of the modified Burgers model was successfully achieved after it was implemented in the numerical code ABAQUS. Based on the full-rheological effect, the proper supporting time of a deep-buried tunnel was studied, and it was proposed that a second lining should be cast in situ approximately 150 days after the excavation of the tunnel. The outcomes of the proposed modified rheological model can accurately represent the creep behavior of deep-buried silty mudstones in a specific engineering instance. The research results can provide a basis for the rheological behavior and supporting time of deep-buried silty mudstone.


---

## Body

## 1. Introduction

It is well known that the creep properties of rock materials significantly influence on the stability of underground structures [1–3]. Particularly, in deep-buried underground caverns and roadways, large deformations and rock-mass failures occur along with rock-mass unloading because of complicated external conditions. Therefore, the creep mechanical behavior of such rocks is closely related to the long-term stability of the engineering structure, which should be carefully studied in the design, operation, and maintenance of tunnels [4–6].As an important rock mechanics characteristic, the rheological properties of some types of rocks have been comprehensively studied through numerous specified creep, relaxation, and quasistatic compression tests [7, 8]. Laboratory testing is one of the most important methods for studying the mechanical properties of rock [9]. Compared with on-site measurements, the advantages of laboratory tests are greatly apparent, including easier long-term observation, stricter control over experimental conditions, elimination of secondary factors, ability to repeat the number of tests, and lower cost [10]. Thus, the research and analysis of rock creep constitutive relationship and parameters are indispensable for evaluating the long-term stability of rock engineering [11, 12]. Moreover, the component element model (CEM) [13] has been widely applied to the study of rheological characteristics in numerous engineering practices because of its clear physical meaning and good applicability in numerical calculations [14].Many studies have been conducted on the development of rock creep constitutive models, and numerous research results have been obtained. For example, the Burgers and Nishihara models are widely used in geotechnical engineering [15, 16]. Hou et al. [17] developed a superposition model for the prediction and analysis of ground dynamic subsidence in mining areas with thick loose layers based on the theory of the Kelvin viscoelastic rheological model and probability-integral method. Fan and Gao [18] established a nonlinear creep equation which can better describe the creep deformation characteristics at each stage. In an analysis of the two models, Tao et al. [19] found that the Burgers model described the creep behavior of rock better than the elasto-viscoelastic model. He et al. [20] proposed a new creep model based on the fractional derivatives of granite analysis at various temperatures. Wang et al. [21], through triaxial unloading creep experiments on layered-rock specimens, proposed a nonlinear creep-damage model of rock under unloading conditions. Liu et al. [22] conducted triaxial creep tests on sandstone and established a time-dependent creep constitutive model based on an unsteady fractional order, which better described the accelerated creep deformation law of rocks. Therefore, the application of an appropriate constitutive relation is crucial in practical engineering [23–25]. Furthermore, the designated rheological model for a particular project is typically unsuitable and inadequate for other models directly owing to various types of rock and the difference in their stress states. However, the CDM has limitations. It can describe the creep attenuation and steady-creep properties of rock but cannot reflect the acceleration-creep properties. Although there are many CEMS to describe the creep mechanical behavior of a rock mass, it is still difficult to develop a new constitutive model by the user-defined material mechanical behavior (UMAT) based on the secondary development platform of finite-element software. Only a few successful cases have been reported to date. The establishment of rheological models requires further study and development, particularly, for the silty mudstone of deep-buried tunnels, whose buried depth is in the range of 1650–2320 m with high in situ stress. Based on the secondary development platform of ABAQUS, the UMAT subroutine of a modified rheological model was successfully realized, which was in good accordance with the creep properties of deep-buried silty mudstone.This study presents the details of a specifically modified rheological model for deep-buried silty mudstone, which has a good qualitative consistency with field monitoring data, and is applied to the problem of proper supporting time for a specific deep-buried water conveyance tunnel. Figure1 shows a brief framework of the entire study. The rheological properties of deep-buried silty mudstones were studied using creep tests in a laboratory. The rheological model is described in the following section. The establishment of the modified Burgers model and parameter estimation is provided in Section 3, including the framework of the technological process of the secondary development platform in ABAQUS coupled with the UMAT subroutine. In Section 4, the application for the research of proper supporting time of a tunnel is described, based on taking full account of the rock rheological effect in Section 2. Finally, a discussion is presented in Section 5, and Section 6 concludes the study. The rheological model proposed in this article can accurately describe the overall time-varying characteristics of deep-buried silty mudstone and provide a reliable basis for creep numerical simulation for predicting the long-term stability and support time of projects.Figure 1
Flowchart of establishment of the rheological model for deep-buried silty mudstone.


## 2. Study on Rheological Properties of Deep-Buried Silty Mudstone

As the CEM is widely used to establish rheological models for a certain rock with the corresponding rheological equation determined, it is imperative to adopt an appropriate rheological model to describe the rock creep characteristics. Furthermore, the parameters were determined, followed by the laboratory creep tests [25, 26]. Finally, the exact form and equation of the rheological constitutive model were created.


### 2.1. Unified Rheological Mechanical Model

The unified rheological model consists of four basic and eleven complex rheological models. According to the loading and unloading characteristics of rock creep to test curves under different stress levels, one or two models suitable for specific rocks can be effectively selected according to Figure2 and common rheological models to satisfy engineering practice [26, 27].Figure 2
Unified rheological mechanical model.


### 2.2. Selection Principle for the Rheological Model

The creep mechanical behavior of rocks must be studied when the rheological mechanical model is applied to a specific practical problem. Therefore, the laboratory creep tests or on-site measurements of rock samples are commonly applied to analyze creep characteristics and geological regularity. Thereafter, a practical rheological mechanical model is selected or established through typical uniaxial-and triaxial-compression experiments of rocks with different confining pressures [28].


### 2.3. Creep Properties of Silty Mudstone

#### 2.3.1. Situation of Laboratory Creep Tests and the Specific Project

Because test equipment should maintain a constant strain for a long time during the creep tests [29], it maintains strict requirements for the confining pressure and control systems. In particular, the TLW-2000 triaxial rheology testing system [30] was adopted to study the creep mechanical behavior of deep-buried silty mudstone. In addition, the laboratory temperature was maintained at 20 ± 2°C, and the humidity remained constant. The UPS power was installed on the rheometer to eliminate the impact of power-supply interruptions.The samples for creep tests were gathered from silty mudstone from a specific water-conveyance tunnel in China, located in northwestern Tianshan, a significantly deep, buried, and long tunnel. The total length of this tunnel is approximately 42 km, whereas the proportion of the length is approximately two-thirds of its buried depth between 500 and 2268 m. The Institute of Crustal Stress of the State Seismological Bureau conducted geo-stress measurements for four boreholes, JDZK3 (460 m), JDZK4 (600 m), JDZK9 (460 m), and JDZK8 (886 m), which were successively completed along the water-conveyance tunnel. See Table1 for the relevant measurement results. According to the inferred results of geological exploration of the site, the maximum horizontal stress was approximately 46–64 MPa.Table 1
Actual geo-stress measurement results.
NumberDepth of measuring section (m)Main stress value (MPa)σHσhσvJDZK3237.85–238.7518.4311.436.55431.55–432.4530.4817.8811.90JDZK4297.00–297.6020.5513.837.74315.86–316.4620.5312.768.24JDZK8437.80–438.6024.7915.7911.38593.82–594.6229.3218.3215.43711.92–712.7233.9820.7818.50JDZK9150.49–151.3411.327.514.07428.42–429.2719.1012.2911.58Note.σH is the maximum horizontal principal stress, σh is the minimum horizontal principal stress, and σv is the vertical principal stress.Finally, the actual measurement results were combined with the creep test conducted by Zhang et al. [11]. All five samples used in the creep tests were naturally air-dried, with an air-dried density of 27.5 kg/m3. In addition, the samples were standard cylinders with a diameter of 50 mm and a height of 100 mm. The laboratory creep tests were performed on five rock triaxial creep test machines according to the actual buried depth of the water conveyance tunnel. Only one test specimen was used for the uniaxial-compression creep test through uniaxial single-stage loading with a confining pressure of 0 MPa, whereas the other four specimens were used for the common confining pressure triaxial-compression creep tests using the graded incremental cyclic-loading method. After the above loading preparations, the creep phase for the silty mudstone was initiated and each step was maintained for 125 h, as shown in Table 2.Table 2
Creep testing procedure for silty mudstone (unit: MPa).
Confining pressure1st load2nd load3rd load4th load5th load6th load0369121518551015202510102030405020153045603018365472


#### 2.3.2. Results of Creep Tests

Typical uniaxial- and triaxial-compression creep curves of the deep-buried silty mudstone are shown in Figure3, where σ is the deviatoric stress and ε is the axial strain of the samples. Figure 3 shows that all the samples exhibited instantaneous elastic deformation, which means that Hooke’s body should exist in the rheological mechanical model for deep-buried silty mudstone. Moreover, at a low deviatoric stress level, the creep curve shows evident attenuation creep characteristics, whereas with an increase in deviatoric stress, the steady-state creep phenomenon becomes significant. Thus, both decay and steady creep characteristics were evident under the condition of high-loading stress.Figure 3
Relation between axial strain and time under different confining pressures: (a) uniaxial compression under no confining pressure; (b) triaxial compression under confining pressure of 5.0 MPa; (c) triaxial compression under confining pressure of 10.0 MPa; (d) triaxial compression under confining pressure of 20.0 MPa; (e) triaxial compression under confining pressure of 30.0 MPa.
(a)(b)(c)(d)(e)Figure3 shows when the deviating stress reached a certain value, the strain rate first increased before subsequently decreasing and tended to be constant over time. Based on the variation law of the creep strain rate, the creep curve could be divided into two stages: transient and steady. It can be inferred that the strain rate decayed in the first stage, stabilized in the second stage, and accelerated in the third stage. Thus, the strain rate increased with an increase in the deviatoric stress.To confirm the influence of the confining pressure on the creep tests, typical compression stress-strain curves of the deep-buried silty mudstone are shown in Figure4. When the strain exceeded the limit value, the stress-strain curve exhibited a plastic platform. It has been reported that the peak strength increases with increasing confining pressure. Quasi-linear and reversible stress-strain relationships were obtained during the initial loading. Based on the stress-strain curves, it can be concluded that the deformation modulus of the samples is closely related to the nonlinear deformation under the initial load.Figure 4
Typical compression stress-strain curves.Figure5 shows the creep failure patterns under confining pressures of 0 MPa and 10 MPa. The sample failure is typically caused by pore compression and crack coalescence. Under time and loading effects, the high axial plastic strain and high strain rate caused by the long-term accumulation of creep effects are the main characteristics of failure. Therefore, the brittle failure is not observed in this sample.Figure 5
Typical creep failure patterns of silty mudstone: (a) no confining pressure; (b) confining pressure of 10.0 MPa.
(a)(b)Because the loading duration at all levels was not sufficiently long, the relative information of the third-stage creep was missing from the creep curves, which is called tertiary creep. Moreover, no unloading tests could be performed. Thus, the Burgers and Nishihara models remain adaptive to be the basic rheological models for deep-buried silty mudstone, according to the existing information from creep tests. However, it is difficult to determine the rheological lower limit value of the Nishihara model because of the missing tertiary creep section in the available creep test curves. Therefore, the Burgers creep model, which combines the Maxwell and Kelvin models, was selected.


## 2.1. Unified Rheological Mechanical Model

The unified rheological model consists of four basic and eleven complex rheological models. According to the loading and unloading characteristics of rock creep to test curves under different stress levels, one or two models suitable for specific rocks can be effectively selected according to Figure2 and common rheological models to satisfy engineering practice [26, 27].Figure 2
Unified rheological mechanical model.


## 2.2. Selection Principle for the Rheological Model

The creep mechanical behavior of rocks must be studied when the rheological mechanical model is applied to a specific practical problem. Therefore, the laboratory creep tests or on-site measurements of rock samples are commonly applied to analyze creep characteristics and geological regularity. Thereafter, a practical rheological mechanical model is selected or established through typical uniaxial-and triaxial-compression experiments of rocks with different confining pressures [28].


## 2.3. Creep Properties of Silty Mudstone

### 2.3.1. Situation of Laboratory Creep Tests and the Specific Project

Because test equipment should maintain a constant strain for a long time during the creep tests [29], it maintains strict requirements for the confining pressure and control systems. In particular, the TLW-2000 triaxial rheology testing system [30] was adopted to study the creep mechanical behavior of deep-buried silty mudstone. In addition, the laboratory temperature was maintained at 20 ± 2°C, and the humidity remained constant. The UPS power was installed on the rheometer to eliminate the impact of power-supply interruptions.The samples for creep tests were gathered from silty mudstone from a specific water-conveyance tunnel in China, located in northwestern Tianshan, a significantly deep, buried, and long tunnel. The total length of this tunnel is approximately 42 km, whereas the proportion of the length is approximately two-thirds of its buried depth between 500 and 2268 m. The Institute of Crustal Stress of the State Seismological Bureau conducted geo-stress measurements for four boreholes, JDZK3 (460 m), JDZK4 (600 m), JDZK9 (460 m), and JDZK8 (886 m), which were successively completed along the water-conveyance tunnel. See Table1 for the relevant measurement results. According to the inferred results of geological exploration of the site, the maximum horizontal stress was approximately 46–64 MPa.Table 1
Actual geo-stress measurement results.
NumberDepth of measuring section (m)Main stress value (MPa)σHσhσvJDZK3237.85–238.7518.4311.436.55431.55–432.4530.4817.8811.90JDZK4297.00–297.6020.5513.837.74315.86–316.4620.5312.768.24JDZK8437.80–438.6024.7915.7911.38593.82–594.6229.3218.3215.43711.92–712.7233.9820.7818.50JDZK9150.49–151.3411.327.514.07428.42–429.2719.1012.2911.58Note.σH is the maximum horizontal principal stress, σh is the minimum horizontal principal stress, and σv is the vertical principal stress.Finally, the actual measurement results were combined with the creep test conducted by Zhang et al. [11]. All five samples used in the creep tests were naturally air-dried, with an air-dried density of 27.5 kg/m3. In addition, the samples were standard cylinders with a diameter of 50 mm and a height of 100 mm. The laboratory creep tests were performed on five rock triaxial creep test machines according to the actual buried depth of the water conveyance tunnel. Only one test specimen was used for the uniaxial-compression creep test through uniaxial single-stage loading with a confining pressure of 0 MPa, whereas the other four specimens were used for the common confining pressure triaxial-compression creep tests using the graded incremental cyclic-loading method. After the above loading preparations, the creep phase for the silty mudstone was initiated and each step was maintained for 125 h, as shown in Table 2.Table 2
Creep testing procedure for silty mudstone (unit: MPa).
Confining pressure1st load2nd load3rd load4th load5th load6th load0369121518551015202510102030405020153045603018365472


### 2.3.2. Results of Creep Tests

Typical uniaxial- and triaxial-compression creep curves of the deep-buried silty mudstone are shown in Figure3, where σ is the deviatoric stress and ε is the axial strain of the samples. Figure 3 shows that all the samples exhibited instantaneous elastic deformation, which means that Hooke’s body should exist in the rheological mechanical model for deep-buried silty mudstone. Moreover, at a low deviatoric stress level, the creep curve shows evident attenuation creep characteristics, whereas with an increase in deviatoric stress, the steady-state creep phenomenon becomes significant. Thus, both decay and steady creep characteristics were evident under the condition of high-loading stress.Figure 3
Relation between axial strain and time under different confining pressures: (a) uniaxial compression under no confining pressure; (b) triaxial compression under confining pressure of 5.0 MPa; (c) triaxial compression under confining pressure of 10.0 MPa; (d) triaxial compression under confining pressure of 20.0 MPa; (e) triaxial compression under confining pressure of 30.0 MPa.
(a)(b)(c)(d)(e)Figure3 shows when the deviating stress reached a certain value, the strain rate first increased before subsequently decreasing and tended to be constant over time. Based on the variation law of the creep strain rate, the creep curve could be divided into two stages: transient and steady. It can be inferred that the strain rate decayed in the first stage, stabilized in the second stage, and accelerated in the third stage. Thus, the strain rate increased with an increase in the deviatoric stress.To confirm the influence of the confining pressure on the creep tests, typical compression stress-strain curves of the deep-buried silty mudstone are shown in Figure4. When the strain exceeded the limit value, the stress-strain curve exhibited a plastic platform. It has been reported that the peak strength increases with increasing confining pressure. Quasi-linear and reversible stress-strain relationships were obtained during the initial loading. Based on the stress-strain curves, it can be concluded that the deformation modulus of the samples is closely related to the nonlinear deformation under the initial load.Figure 4
Typical compression stress-strain curves.Figure5 shows the creep failure patterns under confining pressures of 0 MPa and 10 MPa. The sample failure is typically caused by pore compression and crack coalescence. Under time and loading effects, the high axial plastic strain and high strain rate caused by the long-term accumulation of creep effects are the main characteristics of failure. Therefore, the brittle failure is not observed in this sample.Figure 5
Typical creep failure patterns of silty mudstone: (a) no confining pressure; (b) confining pressure of 10.0 MPa.
(a)(b)Because the loading duration at all levels was not sufficiently long, the relative information of the third-stage creep was missing from the creep curves, which is called tertiary creep. Moreover, no unloading tests could be performed. Thus, the Burgers and Nishihara models remain adaptive to be the basic rheological models for deep-buried silty mudstone, according to the existing information from creep tests. However, it is difficult to determine the rheological lower limit value of the Nishihara model because of the missing tertiary creep section in the available creep test curves. Therefore, the Burgers creep model, which combines the Maxwell and Kelvin models, was selected.


## 2.3.1. Situation of Laboratory Creep Tests and the Specific Project

Because test equipment should maintain a constant strain for a long time during the creep tests [29], it maintains strict requirements for the confining pressure and control systems. In particular, the TLW-2000 triaxial rheology testing system [30] was adopted to study the creep mechanical behavior of deep-buried silty mudstone. In addition, the laboratory temperature was maintained at 20 ± 2°C, and the humidity remained constant. The UPS power was installed on the rheometer to eliminate the impact of power-supply interruptions.The samples for creep tests were gathered from silty mudstone from a specific water-conveyance tunnel in China, located in northwestern Tianshan, a significantly deep, buried, and long tunnel. The total length of this tunnel is approximately 42 km, whereas the proportion of the length is approximately two-thirds of its buried depth between 500 and 2268 m. The Institute of Crustal Stress of the State Seismological Bureau conducted geo-stress measurements for four boreholes, JDZK3 (460 m), JDZK4 (600 m), JDZK9 (460 m), and JDZK8 (886 m), which were successively completed along the water-conveyance tunnel. See Table1 for the relevant measurement results. According to the inferred results of geological exploration of the site, the maximum horizontal stress was approximately 46–64 MPa.Table 1
Actual geo-stress measurement results.
NumberDepth of measuring section (m)Main stress value (MPa)σHσhσvJDZK3237.85–238.7518.4311.436.55431.55–432.4530.4817.8811.90JDZK4297.00–297.6020.5513.837.74315.86–316.4620.5312.768.24JDZK8437.80–438.6024.7915.7911.38593.82–594.6229.3218.3215.43711.92–712.7233.9820.7818.50JDZK9150.49–151.3411.327.514.07428.42–429.2719.1012.2911.58Note.σH is the maximum horizontal principal stress, σh is the minimum horizontal principal stress, and σv is the vertical principal stress.Finally, the actual measurement results were combined with the creep test conducted by Zhang et al. [11]. All five samples used in the creep tests were naturally air-dried, with an air-dried density of 27.5 kg/m3. In addition, the samples were standard cylinders with a diameter of 50 mm and a height of 100 mm. The laboratory creep tests were performed on five rock triaxial creep test machines according to the actual buried depth of the water conveyance tunnel. Only one test specimen was used for the uniaxial-compression creep test through uniaxial single-stage loading with a confining pressure of 0 MPa, whereas the other four specimens were used for the common confining pressure triaxial-compression creep tests using the graded incremental cyclic-loading method. After the above loading preparations, the creep phase for the silty mudstone was initiated and each step was maintained for 125 h, as shown in Table 2.Table 2
Creep testing procedure for silty mudstone (unit: MPa).
Confining pressure1st load2nd load3rd load4th load5th load6th load0369121518551015202510102030405020153045603018365472


## 2.3.2. Results of Creep Tests

Typical uniaxial- and triaxial-compression creep curves of the deep-buried silty mudstone are shown in Figure3, where σ is the deviatoric stress and ε is the axial strain of the samples. Figure 3 shows that all the samples exhibited instantaneous elastic deformation, which means that Hooke’s body should exist in the rheological mechanical model for deep-buried silty mudstone. Moreover, at a low deviatoric stress level, the creep curve shows evident attenuation creep characteristics, whereas with an increase in deviatoric stress, the steady-state creep phenomenon becomes significant. Thus, both decay and steady creep characteristics were evident under the condition of high-loading stress.Figure 3
Relation between axial strain and time under different confining pressures: (a) uniaxial compression under no confining pressure; (b) triaxial compression under confining pressure of 5.0 MPa; (c) triaxial compression under confining pressure of 10.0 MPa; (d) triaxial compression under confining pressure of 20.0 MPa; (e) triaxial compression under confining pressure of 30.0 MPa.
(a)(b)(c)(d)(e)Figure3 shows when the deviating stress reached a certain value, the strain rate first increased before subsequently decreasing and tended to be constant over time. Based on the variation law of the creep strain rate, the creep curve could be divided into two stages: transient and steady. It can be inferred that the strain rate decayed in the first stage, stabilized in the second stage, and accelerated in the third stage. Thus, the strain rate increased with an increase in the deviatoric stress.To confirm the influence of the confining pressure on the creep tests, typical compression stress-strain curves of the deep-buried silty mudstone are shown in Figure4. When the strain exceeded the limit value, the stress-strain curve exhibited a plastic platform. It has been reported that the peak strength increases with increasing confining pressure. Quasi-linear and reversible stress-strain relationships were obtained during the initial loading. Based on the stress-strain curves, it can be concluded that the deformation modulus of the samples is closely related to the nonlinear deformation under the initial load.Figure 4
Typical compression stress-strain curves.Figure5 shows the creep failure patterns under confining pressures of 0 MPa and 10 MPa. The sample failure is typically caused by pore compression and crack coalescence. Under time and loading effects, the high axial plastic strain and high strain rate caused by the long-term accumulation of creep effects are the main characteristics of failure. Therefore, the brittle failure is not observed in this sample.Figure 5
Typical creep failure patterns of silty mudstone: (a) no confining pressure; (b) confining pressure of 10.0 MPa.
(a)(b)Because the loading duration at all levels was not sufficiently long, the relative information of the third-stage creep was missing from the creep curves, which is called tertiary creep. Moreover, no unloading tests could be performed. Thus, the Burgers and Nishihara models remain adaptive to be the basic rheological models for deep-buried silty mudstone, according to the existing information from creep tests. However, it is difficult to determine the rheological lower limit value of the Nishihara model because of the missing tertiary creep section in the available creep test curves. Therefore, the Burgers creep model, which combines the Maxwell and Kelvin models, was selected.


## 3. The Modified Burgers Model for Deep-Buried Silty Mudstone

### 3.1. Methodology of Modified Burgers Model

The constitutive equation for the classical Burgers creep model is as follows [11, 29]:(1)ε=σE0+ση2t+σE11−e−E1/η1t,where σ, ε, and t represent the deviatoric stress, strain, and time, respectively, E0 and η2 are the elastic modulus and viscosity coefficient of the Maxwell body, respectively, and E1 and η1 are the elastic modulus and viscosity coefficient of the Kelvin body, respectively.As shown in equation (1), the permanent viscous deformation item τ=σ/η2⋅t is proportional to the loading time in the classic Burgers model. However, it does not always increase throughout the loading time, which means that the increment of viscous deformation (τ) gradually decreases with an increase in loading time and eventually tends to reach a steady state. Therefore, a nonlinear correction of pure viscous elements (η2) is made based on the classic Burgers model. In other words, η2 is converted into a nonlinear dashpot, and the corresponding equation can be expressed as(2)η2=AeBt,where A and B represent the parameters to be estimated for the nonlinear dashpot (η2).Figure6 shows the conceptual structure of the modified Burgers model.Figure 6
Conceptual structure of the modified Burgers model.Under loading conditions, the creep equation of the modified Burgers model is in the form:(3)ε=σE0+σAB1−e−Bt+σE11−e−E1/η1t.Assuming unloaded at the moment oft0, the corresponding equation under unloading conditions is of the form:(4)ε=σAB1−e−Bt0+σE11−e−E1/η1t0∗e−E1/η1t−t0.It should be mentioned that the nonlinear dashpot degenerates into a linear viscoelastic pot when the value of parameterB is close to zero, which means that the modified Burgers model is turned into the classic Burgers model again.


### 3.2. Parameter Estimation

According to equations (3) and (4), there are five parameters to be determined, E0, A, B, E1, and η1. First, the elastic modulus (E0) of Hooke’s body can be obtained directly from the creep curve, which can be represented by the following mathematical equation. In the uniaxial creep test,(5)E0=σε0,where ε0 represents the initial strain value.In the triaxial creep tests,(6)E0=σ1−μσ2+σ3ε0,where σ1 is the deviatoric stress, whereas σ2 and σ3 are the confining pressures, which follow σ2=σ3, and μ acts as Poisson’s ratio.Four parameters remain to be recognized in the modified Burgers model. Taking a static uniaxial creep test as an example, the corresponding creep equation at the moment ofti can be expressed as(7)εi=σE0+σAB1−e−Bti+σE11−e−E1/η1ti.The constrained optimization problem can be expressed as follows:(8)GE1,η1,A,B=∑i=1Nεi−εi˜2.Also,(9)∂G∂E1=0∂G∂η1=0∂G∂A=0∂G∂B=0,where εi˜ and εi represent the theoretical and measured values, respectively, at the moment ti; the function of G is the square of the difference between the measured value and the theoretical value of the strain. Thus, the goal of least-square method is to minimize the function min GE1,η1,A,B, which should satisfy equation (9). This is typically solved by using successive linearization. To obtain the real optimization parameters instead of the local minima point, it should be mentioned that the initial value of the iteration is determined based on the characteristics of the creep test curves, instead of a blind value. Moreover, the L-M algorithm is adopted in the iterative computation, which is an improved algorithm of the least-square method [1], as it has been widely used and not overly dependent on the initial value of the iteration. Additionally, it is not easy to rapidly converge to the local minima or attain a high iterative-convergence speed.Subsequently, the creep equation of the modified Burgers model was embedded into the function library of origin through its user-defined-function (UDF) to determine rheological parameters using the L-M method. Taking 0 and 10 MPa confining pressures as an example, the relevant parameters for modifying the Burgers creep model are determined through the iterative analysis of uniaxial and triaxial creep tests, as shown in Table3.Table 3
Creep parameters of deep-buried silty mudstone under different confining pressures.
Confining pressure (MPa)Deviatoric stress (MPa)E0 (GPa)η2 (GPa·h)E1 (GPa)A (GPa·h)B (/h)0378.9500.5100.126130.70.299690.9488.4139.625940.70.144990.0560.2175.618303.20.08191264.8688.3173.518171.30.05691561.5875.0179.314637.20.03691841.3914.6197.910743.80.0271101017.83818.1182.692372.80.4532040.81838.3171.042912.20.1993038.92554.9285.222444.30.08004051.61715.785.631830.30.03125046.92054.8109.617746.6−0.207


### 3.3. Secondary Development of Procedure

#### 3.3.1. Creep Solution Method for the Modified Burgers Model

During operation of a water-conveyance tunnel, the temperature of the rock mass is reduced to a certain extent owing to the water environment. However, the distribution of the seepage field and pressure of deep-buried silty mudstones can be significantly influenced and changed by drainage measures during the construction period. Thus, it can be inferred that the rheological parameters of the silty mudstone are not invariable with changes in the temperature and seepage pressure of the rock mass during the operation period. Therefore, the parameters of the materials can be changed by changing the temperatureT or seepage pressure P during the secondary development of the procedure. In other words, the developed program can be used directly for transient analysis. Specifically, when the rheological parameters of materials are set to a series of fixed values under the condition of a certain temperature or seepage pressure, it is typically referred to as the steady-state analysis in the numerical simulation. In the case of different rheological parameters under varying input temperatures or pressures, the developed program will use a linear interpolation to obtain the corresponding rheological parameters according to the temperature and pressure of the current time to conduct the transient analysis. Figure 7 shows an overview of the secondary development of the modified rheological model for the deep-buried silty mudstones.Figure 7
Flowchart of developed program for the modified Burgers model.It should be mentioned that the core objective of the UMAT program is to provide the Jacobian matrix of the specific constitutive model, which is the change rate of stress increment corresponding to variable increment, namely,∂Δσ/∂Δε. Subsequently, the available stress was updated using the ABAQUS program. For example, the results of σn and εn are both known beforehand, and the strain increment dεn+1 is given by the ABAQUS program. Eventually, a new stress state (σn+1) was calculated based on the Jacobian matrix DDSDDE using UMAT.The UMAT subroutine can indirectly support the proposed model. In view of the establishment of the modified Burgers rheological mechanical model, the process of cooperative work is further described in Figure8, based on a material unit-integral point combined with the ABAQUS program with the UMAT subroutine.Figure 8
Collaborative work process of the ABAQUS program and UMAT subroutine.


#### 3.3.2. Numerical Verification

To verify the accuracy of the UMAT procedure, a numerical simulation of laboratory creep tests of silty mudstone was conducted using the modified Burgers model. Figure9 shows that the computational model was highly consistent with standard cylindrical samples: 50 mm in diameter and 100 mm in height. Moreover, in this numerical calculation model, there were more than 3045 nodes and 2560 elements, all of which are hexahedrons. The bottom edge constraint was fixed, while the entire side of the cylindrical samples with normal constraints represented the confining pressure, and the top with normal constraints represented the deviatoric stress. Most of the physical and mechanical properties of silty mudstone used for numerical testing were derived from the previous laboratory creep tests and parameter estimation results, while other properties of silty mudstone are listed in Table 4.Figure 9
Numerical simulation model of creep test samples.Table 4
The properties of silty mudstone employed in numerical simulations.
PropertiesValueUnit weight,ρ (kg/m3)2525Unit elastic resistance coefficient,K0MPa/cm40Firmness coefficient.f4Poisson’s ratio,μ0.265Cohesion strength,c (MPa)1.10Fiction angle,φ (°)41Dilatancy angle,ψ (°)5.12Shear wave velocity,Vs (m/s)2450Longitudinal wave velocity,Vl (m/s)4200As the node at the center of the cylinder surface is a displacement-monitoring point, more emphasis was placed on the process of its axial displacement in the entire monitoring process. Subsequently, the contrast to the creep curve prediction using the modified Burgers model and the numerical simulation test results is shown in Figure10, with confining pressures of 0, 10, and 30 MPa. In Figure 10, σ is the constant axial pressure and R2 is the fitting correlation coefficient.Figure 10
Comparison of Burgers model predicted creep curves with test results: (a) confining pressure of 0 MPa; (b) confining pressure of 10 MPa; (c) confining pressure of 30 MPa.
(a)(b)(c)Figure10 shows the intuitive results that the modified Burgers creep model can predict the creep curve with higher accuracy than the measured curve under the confining pressure of 0, 10, and 30 MPa. The square of correlation coefficients (R2) are both above 0.98, which means that the model parameters could be identified effectively by the technological process of the secondary development platform of ABAQUS coupled with UMAT. Thus, a nonlinear creep model was established in accordance with the creep properties of deep-buried silty mudstone in the project area, which could describe the overall time-dependent behavior.


## 3.1. Methodology of Modified Burgers Model

The constitutive equation for the classical Burgers creep model is as follows [11, 29]:(1)ε=σE0+ση2t+σE11−e−E1/η1t,where σ, ε, and t represent the deviatoric stress, strain, and time, respectively, E0 and η2 are the elastic modulus and viscosity coefficient of the Maxwell body, respectively, and E1 and η1 are the elastic modulus and viscosity coefficient of the Kelvin body, respectively.As shown in equation (1), the permanent viscous deformation item τ=σ/η2⋅t is proportional to the loading time in the classic Burgers model. However, it does not always increase throughout the loading time, which means that the increment of viscous deformation (τ) gradually decreases with an increase in loading time and eventually tends to reach a steady state. Therefore, a nonlinear correction of pure viscous elements (η2) is made based on the classic Burgers model. In other words, η2 is converted into a nonlinear dashpot, and the corresponding equation can be expressed as(2)η2=AeBt,where A and B represent the parameters to be estimated for the nonlinear dashpot (η2).Figure6 shows the conceptual structure of the modified Burgers model.Figure 6
Conceptual structure of the modified Burgers model.Under loading conditions, the creep equation of the modified Burgers model is in the form:(3)ε=σE0+σAB1−e−Bt+σE11−e−E1/η1t.Assuming unloaded at the moment oft0, the corresponding equation under unloading conditions is of the form:(4)ε=σAB1−e−Bt0+σE11−e−E1/η1t0∗e−E1/η1t−t0.It should be mentioned that the nonlinear dashpot degenerates into a linear viscoelastic pot when the value of parameterB is close to zero, which means that the modified Burgers model is turned into the classic Burgers model again.


## 3.2. Parameter Estimation

According to equations (3) and (4), there are five parameters to be determined, E0, A, B, E1, and η1. First, the elastic modulus (E0) of Hooke’s body can be obtained directly from the creep curve, which can be represented by the following mathematical equation. In the uniaxial creep test,(5)E0=σε0,where ε0 represents the initial strain value.In the triaxial creep tests,(6)E0=σ1−μσ2+σ3ε0,where σ1 is the deviatoric stress, whereas σ2 and σ3 are the confining pressures, which follow σ2=σ3, and μ acts as Poisson’s ratio.Four parameters remain to be recognized in the modified Burgers model. Taking a static uniaxial creep test as an example, the corresponding creep equation at the moment ofti can be expressed as(7)εi=σE0+σAB1−e−Bti+σE11−e−E1/η1ti.The constrained optimization problem can be expressed as follows:(8)GE1,η1,A,B=∑i=1Nεi−εi˜2.Also,(9)∂G∂E1=0∂G∂η1=0∂G∂A=0∂G∂B=0,where εi˜ and εi represent the theoretical and measured values, respectively, at the moment ti; the function of G is the square of the difference between the measured value and the theoretical value of the strain. Thus, the goal of least-square method is to minimize the function min GE1,η1,A,B, which should satisfy equation (9). This is typically solved by using successive linearization. To obtain the real optimization parameters instead of the local minima point, it should be mentioned that the initial value of the iteration is determined based on the characteristics of the creep test curves, instead of a blind value. Moreover, the L-M algorithm is adopted in the iterative computation, which is an improved algorithm of the least-square method [1], as it has been widely used and not overly dependent on the initial value of the iteration. Additionally, it is not easy to rapidly converge to the local minima or attain a high iterative-convergence speed.Subsequently, the creep equation of the modified Burgers model was embedded into the function library of origin through its user-defined-function (UDF) to determine rheological parameters using the L-M method. Taking 0 and 10 MPa confining pressures as an example, the relevant parameters for modifying the Burgers creep model are determined through the iterative analysis of uniaxial and triaxial creep tests, as shown in Table3.Table 3
Creep parameters of deep-buried silty mudstone under different confining pressures.
Confining pressure (MPa)Deviatoric stress (MPa)E0 (GPa)η2 (GPa·h)E1 (GPa)A (GPa·h)B (/h)0378.9500.5100.126130.70.299690.9488.4139.625940.70.144990.0560.2175.618303.20.08191264.8688.3173.518171.30.05691561.5875.0179.314637.20.03691841.3914.6197.910743.80.0271101017.83818.1182.692372.80.4532040.81838.3171.042912.20.1993038.92554.9285.222444.30.08004051.61715.785.631830.30.03125046.92054.8109.617746.6−0.207


## 3.3. Secondary Development of Procedure

### 3.3.1. Creep Solution Method for the Modified Burgers Model

During operation of a water-conveyance tunnel, the temperature of the rock mass is reduced to a certain extent owing to the water environment. However, the distribution of the seepage field and pressure of deep-buried silty mudstones can be significantly influenced and changed by drainage measures during the construction period. Thus, it can be inferred that the rheological parameters of the silty mudstone are not invariable with changes in the temperature and seepage pressure of the rock mass during the operation period. Therefore, the parameters of the materials can be changed by changing the temperatureT or seepage pressure P during the secondary development of the procedure. In other words, the developed program can be used directly for transient analysis. Specifically, when the rheological parameters of materials are set to a series of fixed values under the condition of a certain temperature or seepage pressure, it is typically referred to as the steady-state analysis in the numerical simulation. In the case of different rheological parameters under varying input temperatures or pressures, the developed program will use a linear interpolation to obtain the corresponding rheological parameters according to the temperature and pressure of the current time to conduct the transient analysis. Figure 7 shows an overview of the secondary development of the modified rheological model for the deep-buried silty mudstones.Figure 7
Flowchart of developed program for the modified Burgers model.It should be mentioned that the core objective of the UMAT program is to provide the Jacobian matrix of the specific constitutive model, which is the change rate of stress increment corresponding to variable increment, namely,∂Δσ/∂Δε. Subsequently, the available stress was updated using the ABAQUS program. For example, the results of σn and εn are both known beforehand, and the strain increment dεn+1 is given by the ABAQUS program. Eventually, a new stress state (σn+1) was calculated based on the Jacobian matrix DDSDDE using UMAT.The UMAT subroutine can indirectly support the proposed model. In view of the establishment of the modified Burgers rheological mechanical model, the process of cooperative work is further described in Figure8, based on a material unit-integral point combined with the ABAQUS program with the UMAT subroutine.Figure 8
Collaborative work process of the ABAQUS program and UMAT subroutine.


### 3.3.2. Numerical Verification

To verify the accuracy of the UMAT procedure, a numerical simulation of laboratory creep tests of silty mudstone was conducted using the modified Burgers model. Figure9 shows that the computational model was highly consistent with standard cylindrical samples: 50 mm in diameter and 100 mm in height. Moreover, in this numerical calculation model, there were more than 3045 nodes and 2560 elements, all of which are hexahedrons. The bottom edge constraint was fixed, while the entire side of the cylindrical samples with normal constraints represented the confining pressure, and the top with normal constraints represented the deviatoric stress. Most of the physical and mechanical properties of silty mudstone used for numerical testing were derived from the previous laboratory creep tests and parameter estimation results, while other properties of silty mudstone are listed in Table 4.Figure 9
Numerical simulation model of creep test samples.Table 4
The properties of silty mudstone employed in numerical simulations.
PropertiesValueUnit weight,ρ (kg/m3)2525Unit elastic resistance coefficient,K0MPa/cm40Firmness coefficient.f4Poisson’s ratio,μ0.265Cohesion strength,c (MPa)1.10Fiction angle,φ (°)41Dilatancy angle,ψ (°)5.12Shear wave velocity,Vs (m/s)2450Longitudinal wave velocity,Vl (m/s)4200As the node at the center of the cylinder surface is a displacement-monitoring point, more emphasis was placed on the process of its axial displacement in the entire monitoring process. Subsequently, the contrast to the creep curve prediction using the modified Burgers model and the numerical simulation test results is shown in Figure10, with confining pressures of 0, 10, and 30 MPa. In Figure 10, σ is the constant axial pressure and R2 is the fitting correlation coefficient.Figure 10
Comparison of Burgers model predicted creep curves with test results: (a) confining pressure of 0 MPa; (b) confining pressure of 10 MPa; (c) confining pressure of 30 MPa.
(a)(b)(c)Figure10 shows the intuitive results that the modified Burgers creep model can predict the creep curve with higher accuracy than the measured curve under the confining pressure of 0, 10, and 30 MPa. The square of correlation coefficients (R2) are both above 0.98, which means that the model parameters could be identified effectively by the technological process of the secondary development platform of ABAQUS coupled with UMAT. Thus, a nonlinear creep model was established in accordance with the creep properties of deep-buried silty mudstone in the project area, which could describe the overall time-dependent behavior.


## 3.3.1. Creep Solution Method for the Modified Burgers Model

During operation of a water-conveyance tunnel, the temperature of the rock mass is reduced to a certain extent owing to the water environment. However, the distribution of the seepage field and pressure of deep-buried silty mudstones can be significantly influenced and changed by drainage measures during the construction period. Thus, it can be inferred that the rheological parameters of the silty mudstone are not invariable with changes in the temperature and seepage pressure of the rock mass during the operation period. Therefore, the parameters of the materials can be changed by changing the temperatureT or seepage pressure P during the secondary development of the procedure. In other words, the developed program can be used directly for transient analysis. Specifically, when the rheological parameters of materials are set to a series of fixed values under the condition of a certain temperature or seepage pressure, it is typically referred to as the steady-state analysis in the numerical simulation. In the case of different rheological parameters under varying input temperatures or pressures, the developed program will use a linear interpolation to obtain the corresponding rheological parameters according to the temperature and pressure of the current time to conduct the transient analysis. Figure 7 shows an overview of the secondary development of the modified rheological model for the deep-buried silty mudstones.Figure 7
Flowchart of developed program for the modified Burgers model.It should be mentioned that the core objective of the UMAT program is to provide the Jacobian matrix of the specific constitutive model, which is the change rate of stress increment corresponding to variable increment, namely,∂Δσ/∂Δε. Subsequently, the available stress was updated using the ABAQUS program. For example, the results of σn and εn are both known beforehand, and the strain increment dεn+1 is given by the ABAQUS program. Eventually, a new stress state (σn+1) was calculated based on the Jacobian matrix DDSDDE using UMAT.The UMAT subroutine can indirectly support the proposed model. In view of the establishment of the modified Burgers rheological mechanical model, the process of cooperative work is further described in Figure8, based on a material unit-integral point combined with the ABAQUS program with the UMAT subroutine.Figure 8
Collaborative work process of the ABAQUS program and UMAT subroutine.


## 3.3.2. Numerical Verification

To verify the accuracy of the UMAT procedure, a numerical simulation of laboratory creep tests of silty mudstone was conducted using the modified Burgers model. Figure9 shows that the computational model was highly consistent with standard cylindrical samples: 50 mm in diameter and 100 mm in height. Moreover, in this numerical calculation model, there were more than 3045 nodes and 2560 elements, all of which are hexahedrons. The bottom edge constraint was fixed, while the entire side of the cylindrical samples with normal constraints represented the confining pressure, and the top with normal constraints represented the deviatoric stress. Most of the physical and mechanical properties of silty mudstone used for numerical testing were derived from the previous laboratory creep tests and parameter estimation results, while other properties of silty mudstone are listed in Table 4.Figure 9
Numerical simulation model of creep test samples.Table 4
The properties of silty mudstone employed in numerical simulations.
PropertiesValueUnit weight,ρ (kg/m3)2525Unit elastic resistance coefficient,K0MPa/cm40Firmness coefficient.f4Poisson’s ratio,μ0.265Cohesion strength,c (MPa)1.10Fiction angle,φ (°)41Dilatancy angle,ψ (°)5.12Shear wave velocity,Vs (m/s)2450Longitudinal wave velocity,Vl (m/s)4200As the node at the center of the cylinder surface is a displacement-monitoring point, more emphasis was placed on the process of its axial displacement in the entire monitoring process. Subsequently, the contrast to the creep curve prediction using the modified Burgers model and the numerical simulation test results is shown in Figure10, with confining pressures of 0, 10, and 30 MPa. In Figure 10, σ is the constant axial pressure and R2 is the fitting correlation coefficient.Figure 10
Comparison of Burgers model predicted creep curves with test results: (a) confining pressure of 0 MPa; (b) confining pressure of 10 MPa; (c) confining pressure of 30 MPa.
(a)(b)(c)Figure10 shows the intuitive results that the modified Burgers creep model can predict the creep curve with higher accuracy than the measured curve under the confining pressure of 0, 10, and 30 MPa. The square of correlation coefficients (R2) are both above 0.98, which means that the model parameters could be identified effectively by the technological process of the secondary development platform of ABAQUS coupled with UMAT. Thus, a nonlinear creep model was established in accordance with the creep properties of deep-buried silty mudstone in the project area, which could describe the overall time-dependent behavior.


## 4. Application for Proper Supporting Time Analysis

A water-conveyance tunnel is difficult to design and construct because its buried depth is in the range of 1650–2320 m with high in situ stress. In addition, it has been proven by laboratory tests that rock mass can show momentous rheological characteristics with high stress [27, 31], which has a significant influence on the mechanical behavior of deep-buried rock masses, and can drastically change the engineering stability. Particularly, for the hole section with obvious rheological effects, it should be carefully studied to determine the supporting time of the lining. If the support is not suitable, the deformation of the rock mass will not be effectively limited. In contrast, the support system will be destroyed because its strength is exceeded when the support is provided too early.


### 4.1. Numerical Calculation Model

Figure11 shows that the computational model scope was 120 m in the x-direction and 120 m in the z-direction for proper supporting time analysis of the deep-buried tunnel, which is more than six times of the tunnel cross-section diameter. It was 20 m in the y-direction, that is, the axial direction of the tunnel. In this numerical calculation model, there were more than 53,285 nodes and 59,392 hexahedral elements. The cross-sectional dimensions are shown in Figure 12. In the longitudinal direction, as illustrated in Figure 13, the span of each excavation cycle was 2.0 m. The first lining and rock bolt were installed immediately after each excavation cycle, and proper supporting time of the second lining was determined by numerical simulation using the modified Burgers rheological model and estimated parameters.Figure 11
Numerical simulation model of deep-buried water-conveyance tunnel.Figure 12
Cross-sectional dimensions of the water-conveyance tunnel.Figure 13
Schematic diagram and numerical simulation of excavation cycle during construction.To facilitate our analysis, the tunnel convergence monitoring locations, including the displacements and plastic strain of the rock-mass at the crown and the maximum compressive stress of the crown segment, are schematically illustrated in Figure14.Figure 14
Convergence monitoring positions in numerical simulations.
(a)(b)In this article, ten excavation cycles (K20000–K20020) were simulated step-by-step using ABAQUS. The properties of the silty mudstones are listed in Table4, and the linings are listed in Table 5, which are measured in the laboratory. The postexcavation simulation provides an “initial state” for further rheological calculations, considering the time-varying changes in the geological parameters and excavation support in the construction process. In Table 5, E denotes the elastic modulus, σc denotes the compressive strength, σt denotes the tensile strength, and tc denotes the construction lining thickness.Table 5
Properties of the support system employed in numerical simulations.
PropertiesFirst liningSecond liningAnchorE (GPa)23.025.0200.0μ0.170.170.30σc (MPa)14.319.1300.0σt (MPa)1.431.71300.0tc (m)0.150.40N/ABased on the engineering best practice, the key analysis steps in the numerical simulation are schematically presented in Table6. In contrast to the numerical analysis of conventional static problems, the time step in the creep analysis should correspond to the actual time of construction and operational processes. In practice, the tunnel boring machine (TBM) method was adopted for full-face tunnelling. We assumed the driving speed of the TBM is 20 m/d, since the construction speed of the TBM method is rapid, which means that the excavation zones could be removed from the numerical model in one day, as shown in Step 2. In the process of shotcrete hardening, it was assumed that the shotcrete procedure was completed within one day, whose strength can reach half the standard values. In the subsequent 27 days, the strength increased linearly to standard values. Undoubtedly, the determination of the time interval is the key point of analysis between the first and second linings, which are achieved by numerical calculation to ensure the ultimate stability of the tunnel.Table 6
Key steps of proper supporting time analysis.
Key stepNumerical simulation contentStep 1Calculation of in situ stress field by applying lateral-pressure coefficientsStep 2Tunnel excavation through ten excavation cycles in one dayStep 3First lining conducted along with excavationStep 4Shotcrete hardening for 28 daysStep 5Time interval between first and second lining to be determined by numerical simulationStep 6Second lining subsequently completedStep 7Long-term analysis for three years in the operation period


### 4.2. Deformation of Tunnel without Any Support

First, the deformation characteristics of deep-buried silty mudstones in a tunnel should be studied; therefore, the supporting time of the lining was initially determined. Three buried depths were selected for comparison in the process of numerical analysis: 1700 m, 2000 m, and 2320 m. Generally, the lateral-pressure coefficient is equal to one at a significantly buried depth [26, 32]. The displacement contour of the rock-mass followed a circular distribution after excavation, and the convergence deformation law of each point around the tunnel was essentially the same. Figure 15 shows the curves of roof subsidence (ur) without any support at three different buried depths, and the monitoring points are shown in Figure 14(a).Figure 15
Curves of roof subsidence without any support.It can be seen from Figure15 that the deformation of the rock mass tends to be larger with an increase in the buried depth. Moreover, the instantaneous deformation were 1.41 cm, 2.09 cm, and 3.01 cm after tunnel excavation with a buried depth of 1700 m, 2000 m, and 2320 m, respectively. The final stable deformation was 3.49 cm, 4.54 cm, and 5.88 cm, respectively. Therefore, the initial deformation accounted for 40.4%, 46.0%, and 51.2% of the final stable deformation, respectively, indicating that the rheological deformation of the deep-buried silty mudstones was significant. Furthermore, it can be seen from Figure 15 that the rheological deformation of rock-mass occurred in the first six months. The displacement rate of the roof subsidence (ur) was much larger at the beginning and decreased approximately 90 days later. The deformation of the rock mass was kept in a proper stable state. Therefore, the suggested time of the second lining for the tunnel was selected as 90 to 150 days after excavation.


### 4.3. Proper Support Time Analysis

Based on the above analysis, the supporting time of the second lining of the tunnel was set as three solutions for comparison: 90 days, 120 days, and 150 days after excavation. To study the stress state of the supporting system and the deformation regularity of the rock mass, a maximum calculated depth of 2320 m was used as an example. Finally, the appropriate support time for deep-buried silty mudstones was determined. Figures16–18 show the curves of each monitoring variable with respect to time under three different solutions in the numerical analysis process, and the monitoring points are shown in Figure 14. It is worth mentioning that ur represents the roof subsidence of rock-mass, εmax represents the maximum plastic strain of rock mass, and σmax is the maximum compressive stress of the segment.Figure 16
Curves of roof subsidence with a designated support system.Figure 17
Curves of maximum plastic strain of roof with a designated support system.Figure 18
Curves of maximum compressive stress of segment.From Figures16 and 17, it can be inferred that the deformation of the displacement and plastic strain were well controlled in the tunnel when the second lining was completed, which meant that the deformation of the tunnel tended to remain stable. It should be mentioned that an instantaneous increment of the monitoring variables was generated because of the application of water pressure. In particular, the time axis started after the second lining was completed, and it could be clearly seen that the stress levels of the second lining had a significant increase in a relatively short time owing to the application of water pressure under the three different solutions. Consequently, the stress of the second lining continued to grow slowly owing to the rheological properties of the deep-buried silty mudstones, and the stress growth process occurred almost 150 days later. The maximum compressive stress of the second lining was 24.59 MPa, 20.39 MPa, and 16.34 MPa in the final steady state at 90 days, 120 days, and 150 days after excavation, respectively. In contrast, C40 concrete was used for the second lining, and the designated compressive strength was 19.1 MPa. The second lining should have a certain degree of safety to ensure the overall safety of the support system. Consequently, it was suggested that the second lining be cast in place approximately 150 days after excavation to ensure the safety of the deep-buried tunnel.


## 4.1. Numerical Calculation Model

Figure11 shows that the computational model scope was 120 m in the x-direction and 120 m in the z-direction for proper supporting time analysis of the deep-buried tunnel, which is more than six times of the tunnel cross-section diameter. It was 20 m in the y-direction, that is, the axial direction of the tunnel. In this numerical calculation model, there were more than 53,285 nodes and 59,392 hexahedral elements. The cross-sectional dimensions are shown in Figure 12. In the longitudinal direction, as illustrated in Figure 13, the span of each excavation cycle was 2.0 m. The first lining and rock bolt were installed immediately after each excavation cycle, and proper supporting time of the second lining was determined by numerical simulation using the modified Burgers rheological model and estimated parameters.Figure 11
Numerical simulation model of deep-buried water-conveyance tunnel.Figure 12
Cross-sectional dimensions of the water-conveyance tunnel.Figure 13
Schematic diagram and numerical simulation of excavation cycle during construction.To facilitate our analysis, the tunnel convergence monitoring locations, including the displacements and plastic strain of the rock-mass at the crown and the maximum compressive stress of the crown segment, are schematically illustrated in Figure14.Figure 14
Convergence monitoring positions in numerical simulations.
(a)(b)In this article, ten excavation cycles (K20000–K20020) were simulated step-by-step using ABAQUS. The properties of the silty mudstones are listed in Table4, and the linings are listed in Table 5, which are measured in the laboratory. The postexcavation simulation provides an “initial state” for further rheological calculations, considering the time-varying changes in the geological parameters and excavation support in the construction process. In Table 5, E denotes the elastic modulus, σc denotes the compressive strength, σt denotes the tensile strength, and tc denotes the construction lining thickness.Table 5
Properties of the support system employed in numerical simulations.
PropertiesFirst liningSecond liningAnchorE (GPa)23.025.0200.0μ0.170.170.30σc (MPa)14.319.1300.0σt (MPa)1.431.71300.0tc (m)0.150.40N/ABased on the engineering best practice, the key analysis steps in the numerical simulation are schematically presented in Table6. In contrast to the numerical analysis of conventional static problems, the time step in the creep analysis should correspond to the actual time of construction and operational processes. In practice, the tunnel boring machine (TBM) method was adopted for full-face tunnelling. We assumed the driving speed of the TBM is 20 m/d, since the construction speed of the TBM method is rapid, which means that the excavation zones could be removed from the numerical model in one day, as shown in Step 2. In the process of shotcrete hardening, it was assumed that the shotcrete procedure was completed within one day, whose strength can reach half the standard values. In the subsequent 27 days, the strength increased linearly to standard values. Undoubtedly, the determination of the time interval is the key point of analysis between the first and second linings, which are achieved by numerical calculation to ensure the ultimate stability of the tunnel.Table 6
Key steps of proper supporting time analysis.
Key stepNumerical simulation contentStep 1Calculation of in situ stress field by applying lateral-pressure coefficientsStep 2Tunnel excavation through ten excavation cycles in one dayStep 3First lining conducted along with excavationStep 4Shotcrete hardening for 28 daysStep 5Time interval between first and second lining to be determined by numerical simulationStep 6Second lining subsequently completedStep 7Long-term analysis for three years in the operation period


## 4.2. Deformation of Tunnel without Any Support

First, the deformation characteristics of deep-buried silty mudstones in a tunnel should be studied; therefore, the supporting time of the lining was initially determined. Three buried depths were selected for comparison in the process of numerical analysis: 1700 m, 2000 m, and 2320 m. Generally, the lateral-pressure coefficient is equal to one at a significantly buried depth [26, 32]. The displacement contour of the rock-mass followed a circular distribution after excavation, and the convergence deformation law of each point around the tunnel was essentially the same. Figure 15 shows the curves of roof subsidence (ur) without any support at three different buried depths, and the monitoring points are shown in Figure 14(a).Figure 15
Curves of roof subsidence without any support.It can be seen from Figure15 that the deformation of the rock mass tends to be larger with an increase in the buried depth. Moreover, the instantaneous deformation were 1.41 cm, 2.09 cm, and 3.01 cm after tunnel excavation with a buried depth of 1700 m, 2000 m, and 2320 m, respectively. The final stable deformation was 3.49 cm, 4.54 cm, and 5.88 cm, respectively. Therefore, the initial deformation accounted for 40.4%, 46.0%, and 51.2% of the final stable deformation, respectively, indicating that the rheological deformation of the deep-buried silty mudstones was significant. Furthermore, it can be seen from Figure 15 that the rheological deformation of rock-mass occurred in the first six months. The displacement rate of the roof subsidence (ur) was much larger at the beginning and decreased approximately 90 days later. The deformation of the rock mass was kept in a proper stable state. Therefore, the suggested time of the second lining for the tunnel was selected as 90 to 150 days after excavation.


## 4.3. Proper Support Time Analysis

Based on the above analysis, the supporting time of the second lining of the tunnel was set as three solutions for comparison: 90 days, 120 days, and 150 days after excavation. To study the stress state of the supporting system and the deformation regularity of the rock mass, a maximum calculated depth of 2320 m was used as an example. Finally, the appropriate support time for deep-buried silty mudstones was determined. Figures16–18 show the curves of each monitoring variable with respect to time under three different solutions in the numerical analysis process, and the monitoring points are shown in Figure 14. It is worth mentioning that ur represents the roof subsidence of rock-mass, εmax represents the maximum plastic strain of rock mass, and σmax is the maximum compressive stress of the segment.Figure 16
Curves of roof subsidence with a designated support system.Figure 17
Curves of maximum plastic strain of roof with a designated support system.Figure 18
Curves of maximum compressive stress of segment.From Figures16 and 17, it can be inferred that the deformation of the displacement and plastic strain were well controlled in the tunnel when the second lining was completed, which meant that the deformation of the tunnel tended to remain stable. It should be mentioned that an instantaneous increment of the monitoring variables was generated because of the application of water pressure. In particular, the time axis started after the second lining was completed, and it could be clearly seen that the stress levels of the second lining had a significant increase in a relatively short time owing to the application of water pressure under the three different solutions. Consequently, the stress of the second lining continued to grow slowly owing to the rheological properties of the deep-buried silty mudstones, and the stress growth process occurred almost 150 days later. The maximum compressive stress of the second lining was 24.59 MPa, 20.39 MPa, and 16.34 MPa in the final steady state at 90 days, 120 days, and 150 days after excavation, respectively. In contrast, C40 concrete was used for the second lining, and the designated compressive strength was 19.1 MPa. The second lining should have a certain degree of safety to ensure the overall safety of the support system. Consequently, it was suggested that the second lining be cast in place approximately 150 days after excavation to ensure the safety of the deep-buried tunnel.


## 5. Discussion

Figure3 shows that the creep behavior of deep-buried silty mudstones is extremely significant under high-deviatoric stresses. The time-dependent deformation decreased with an increase in confining pressure, indicating that less creep of the rock sample may occur at a high confining pressure. The creep strain rate of the rock samples varied with deviatoric stress. The strain rate tends to be close to zero over time under low-deviatoric stress. However, the strain rate increased with the deviatoric stress.Figure10 shows that the parameters determined using the modified Burgers model qualitatively agree well with the on-site monitoring data of the deep-buried silty mudstones. Subsequently, the UMAT subroutine of the modified Burgers model was compiled successfully based on the secondary development platform of the ABAQUS program, and the numerical analysis was validated using the experimental data. However, unless it can be verified experimentally, it cannot be concluded that the developed model is suitable for other load conditions.Moreover, a nonlinear creep model was established based on experimental data from loading creep tests instead of unloading creep tests. Further study is required to resolve the challenges in laboratory creep testing procedures, including the stress state of the rock, testing boundary conditions and loading patterns, and specimen size effect [33]. For instance, it is still difficult to maintain constant temperature and humidity conditions in the laboratory, and creep-test results are influenced by changes in the external environment during the lengthy period of the rheological test. Therefore, the impact of these factors on quantitative estimation needs to be studied further. In addition, the creep of soft rocks under high-stress levels are often highly nonlinear and does not satisfy the linear superposition principle, which is conducted using the graded loading method. Therefore, it is necessary to use appropriate methods.


## 6. Conclusions

In this article, uniaxial- and triaxial-compressive creep tests of deep-buried silty mudstones were conducted experimentally, and an improved Burgers constitutive model was proposed based on the classical Burgers model. The following conclusions were drawn as follows:(1)
The modified Burgers model was in good agreement with the experimental results when a nonlinear dashpot was introduced to consider the nonlinear time-dependent strength deterioration of the rock-mass. Based on the test results, the corresponding creep equations were embedded into origin through its UDF function to determine the rheological parameters using the L–M method. Subsequently, a new nonlinear creep model was established based on the creep properties of deep-buried silty mudstones for a specific water-conveyance tunnel.(2)
The UMAT subroutine of the modified Burgers model was successfully compiled, based on the secondary development platform of the ABAQUS program. The numerical results demonstrate the high precision of the modified Burgers creep model in predicting the creep curve compared with the measured curve. Thereafter, a subprogram was used to study the surrounding rock-deformation law and support timing of the deep-buried silty mudstone tunnel section. The results show that the proposed rheological model can accurately describe the overall time-dependent behavior of deep-buried silty mudstone in a specific engineering instance, which could provide a scientific and convincing basis for creep numerical simulations for predicting the long-term stability of a specific deep-buried project.



---

*Source: 3124359-2023-02-11.xml*