Study on Surrounding Rock Deformation Laws of an Argillaceous Soft Rock Roadway Based on the Creep Damage Model

The indexes of argillaceous soft rocks in Western China, such as rock strength and softening coefficient, are much smaller than those of soft rocks commonly seen in other areas. Such argillaceous soft rocks are susceptible to a more serious creep phenomenon if encountering water. Accordingly, the study theoretically constructed a constitutive creep model considering the weakening coefficient and damage of surrounding rocks. Next, the viability of this model was verified by combining numerical simulation and field monitoring. Accordingly, this study conducted a comparative analysis on the creep range and cross-sectional deformation of argillaceous surrounding rocks of a roadway in a titled stratum under dry and damp states. Results showed that the Cvisc model considering the weakening coefficient of surrounding rocks was applicable to the decay creep and uniform creep stage, and that considering the damage of surrounding rocks was applicable to the accelerated creep stage. During the creep process, the argillaceous surrounding rocks were prone to viscoelastic deformation. The damp state had an obvious time effect on the damage of surrounding rocks. Specifically, the creep range was expanded under a damp state at 18 days compared with that under a dry state, and the change in the left and upper displacements was especially significant. The roadway cross-section changed asymmetrically under both dry and damp states, where it was changed into an approximate deflected rhombic shape under the dry state and into an approximate deflected trapezoidal shape under the damp state. Additionally, the roadway cross-sectional change under the humid state was evidently greater than that under the dry state.


INTRODUCTION
Rocks are primarily divided into magmatic, sedimentary, and metamorphic rocks, where sedimentary rocks play a dominant role, including limestones, sandstones, shales, and mudstones, which are especially susceptible to creep when encountering water (Xu, 2007;Zhang et al., 2011;Li et al., 2014). With the passing of time, the surrounding rocks on the roadway surface were gradually eroded by the moisture in the atmosphere, thus losing self-supporting ability. Consequently, the roadway is subjected to serious deformation, which further causes stability problems in terms of numerous aboveground and underground rock engineering. Therefore, argillaceous soft rock roadways need repeated renovation and are usually stuck in a vicious circle of "roadway formationfailure-repair-secondary failure-secondary repair." Roadway support is a major problem confusing mine safety and efficient production and one of the technical problems requiring urgent solutions in engineering (Wang and Wang, 2000;Chen et al., 2007;Lu et al., 2007).
With regard to the deformation instability problem of soft rock roadways, rock creep and soft rock characteristics have been extensively investigated by domestic (Chinese) and foreign scholars. Wang H J et al. explored the influence of different horizontal stresses on roadway deformation via FLAC 3D and pointed out that confining pressure was positively correlated with the creep rate of rock masses (Wang et al., 2014). Ma G Y et al. studied the mechanical properties of double-layer initial support for deep high-stress tunnels by constructing a damage evolution model considering the elastoplastic deviatoric strain decrement; they subsequently concluded that the rock deformation could be well reflected by the constitutive elastoplastic creep damage evolution model . Xu W Y et al. stated that if introduced into nonlinear functions, the generalized Bingham creep model could reflect the decay, steady-state, and accelerated creep stages during the creep process well . Xie Y P et al. probed the total creep characteristics of carbonaceous slate specimens by establishing a hardening damage model and pointed out that the Cvisc model improved on the basis of damage theory which could reflect rock creep well (Xie et al., 2021). Yang D S et al. studied the relationships of gas permeability and creep properties of mudstones with the moisture content (Yang et al., 2010;Yang et al., 2011). Huang Y B et al. explored the failure mechanism of deep soft rock chambers via numerical simulation and proposed an optimized recovery scheme (Huang et al., 2021). Liu G et al. analyzed the changes in the crustal stress before and after the excavation of a soft rock roadway and the influence on the change and failure of surrounding rocks; then, they discussed about the deformation failure mechanism of high-stress soft surrounding rocks from the angles of the changes in the confining pressure and strength of rock masses and engineering rock masses (Liu et al., 2000). Wu G J et al. built a nonlinear creep damage equation and indicated the positive correlation between damage and deviatoric stress (Wu et al., 2010).
The above study results have a very good guiding effect on the deformation of soft rock roadways, but the constitutive surrounding rock creep model thoroughly considers the damage and weakening coefficient of surrounding rocks. Additionally, the literature has scarcely explored the deformation laws of argillaceous soft rocks under a damp environment. In this study, a constitutive soft rock model considering the damage and weakening coefficient of surrounding rocks was constructed on the basis of the classical Cvisc constitutive model. Next, the reasonability of this constitutive model was verified by combining numerical simulation and field monitoring. On this basis, this research comparatively analyzed the creep range of surrounding rocks in an argillaceous soft rock roadway and the roadway cross-sectional deformation laws under dry and damp states, expecting to provide a basis for the field support design and facilitate mine safety and efficient production.

CONSTITUTIVE DAMAGE MODELING FOR ARGILLACEOUS SOFT ROCKS
Mechanical studies show that rock deformation not only includes transient elastic and plastic deformation characteristics unrelated to time but also includes rheological deformation characteristics related to time (Zhang and Xu, 2006). During deep tunnel excavation, the deformation with time characteristic accounts for 70% of the total deformation, and rheological deformation plays a dominant role under deep conditions (Paraskevopoulou and Diederichs, 2018). To conveniently explore the rheological behaviors of rock masses, three basic units-elastomer (H), viscous body (N), and plastic body (S)are typically combined in different forms with the main rheological models and deformation behaviors as seen in Table 1 (Zhu et al., 2011;Zhao, 2017): The above classical rheological models can describe well the partial creep properties of rock masses, but transient plastic deformation and accelerated creep deformation are neglected. In deep environment, rock masses are generally characterized by typical whole creep stage characteristics, namely, decay creep, uniform creep, and accelerated creep. As pointed out in the literature , the Cvisc model is a combination of Burgers model and Mohr-Coulomb model, where the former can describe the constitutive viscoelastic relations in the creep process well, and the latter is capable of describing the constitutive elastoplastic relations. Figure 1 displays the Cvisc viscoelasticplastic model:  As shown in Figure 1, the Cvisc viscoelastic-plastic model can describe the transient deformation, decay creep, uniform creep, and plastic deformation in the creep process of soft rocks, with the following constitutive equation: where ε t is the total strain, ε ρ stands for the plastic strain, σ is the horizontal stress, and t represents time. EM, ηM, EK, and ηK represent the elasticity modulus and viscosity coefficient of Maxwell body and Kelvin body. During the creep process of soft rocks subjected to argillization when encountering water, their deformation is largely affected by the damage degree inside the material and the intensity of water environment. According to relevant literature (Zhao, 2017), rock damage occurs in the accelerated creep stage during the creep process, and the damage variable D t presents the Weibull distribution with time t, expressed by the following equation: The influence of water environment on creep was considered by introducing the weakening coefficient λ of surrounding rocks. Then, λ was substituted into Eq. 1 to obtain the constitutive equations of transient deformation, decay creep, uniform creep, and plastic deformation during the creep process as follows: The Cvisc model considering damage could describe the accelerated creep stage (Wang et al., 2014) (Cai and Lai, 2003;Chen et al., 2022). In addition, if damage was considered for the components subjected to decay creep, uniform creep, and plastic deformation in Eq. 1; Figure 2 presents the Cvisc viscoelasticplastic model considering damage:  The constitutive equation for the accelerated creep stage during the creep process could be obtained by combining Eqs 2 and 3 on the bases of effective stress model and strain equivalence principle: Eqs 3 and 4 are the constitutive equations for the whole rock creep process (transient deformation, decay creep, uniform creep, accelerated creep, and plastic deformation) and consider the effects of rock damage variable D t and surrounding rock weakening coefficient λ on creep. As pointed out in relevant literature (Paraskevopoulou and Diederichs, 2018), the displacement of surrounding rocks is in direct proportion to strain during the excavation of the underground tunnel. The aforementioned strain-time relations in Eqs 3 and 4 were used to obtain the displacement-time relations of the surrounding rocks during tunnel excavation as shown in Figure 3.
From Figure 3, the whole time-dependent displacement process around the surrounding rocks could be well reflected by different rheological components within different time periods. This phenomenon indicates that the displacementtime relations of surrounding rocks could be accurately described by the constitutive Eqs 3 and 4 of soft rock creep that considered the rock damage variable D t and surrounding rock weakening coefficient λ.

Modeling
In a copper mine, the exploitation depth already reaches 600 m, and the surrounding rocks in its roadway are mostly mudstones and sandstones, accompanied by bedding development, low strength, and proneness to argillization when encountering water (Shen, 2014;Yang et al., 2018;Wu et al., 2020). According to the geological information of this copper mine, a 2D plane strain calculation model ( Figure 4) of this deep roadway was constructed using the UDEC discrete element software. To eliminate the influence of boundary effect on the model calculation, both the length and height of the model were set as 5 times of the tunnel diameter (30 m), and the monitoring points were arranged at the vault and 0.3 m depth of the left and right walls, respectively.
The pressure applied above the model was P = γh, where γ is the average volume weight of rock strata above the model. The horizontally applied stress was μP (μ = 0.42), and the horizontal displacements at the left and right sides and the vertical displacement at the bottom were restricted. In this model, the random fractures conforming to the geological survey results of this copper mine and the horizontal and vertical fractures preventing the automatic deletion of isolated random fractures in the model calculation were generated through the Monte Carlo method (Yu, 2015). Table 2 shows the reduced physical and mechanical parameters of this roadway.

Inversion of Model Creep Parameters
A precondition for numerical model calculation lies in the accurate numerical calculation of parameters. The BP neural  network can realize any nonlinear relation between input parameters and target parameters and acquire target prediction parameters, meeting the accuracy requirement by adjusting the initial weight and threshold (Chen, 2020).
To restore the properties of surrounding rocks in the roadway as much as possible, the Cvisc constitutive model was adopted for the mudstones of this roadway during numerical model calculation and Mohr-Coulomb constitutive model for other rock strata (Zha et al., 2021). Finally, the mean square errors between the numerically calculated displacements at monitoring points and the actually monitored displacements were made into an input set and the creep parameters (EM, ηM, EK, and ηK) into an output set for the prediction through the BP neural network until they meet the accuracy requirement. Figure 5 presents the numerical calculation results and actual monitoring results.
By analyzing Figure 5, the numerical calculation results largely accorded with the field measurement results, proving that the inversion parameters could accurately reflect the deformation of surrounding rocks in the roadway. Table 3 presents the parameters finally obtained by inversion.

ANALYSIS OF CALCULATION RESULTS
In the literature (Zhang et al., 2009), the softening coefficient of surrounding rocks in the roadway of this copper mine was obtained through uniaxial compression experiments as 0.12-0.21, indicating that the strength of surrounding rocks was largely affected by water. Given that groundwater prevention measures were already taken in the roadway excavation process of this mine, the strength of mudstones was weakened using the softening coefficient of 0.21 considering the influence of water on mudstone strength. To analyze the influence of dampness-induced  argillization on the deformation of surrounding rocks, the creep range of surrounding rocks and roadway crosssectional deformation laws under dry and damp states were comparatively analyzed using the numerical model in Section 3.1.

Verification of Constitutive Model
The creep parameters obtained through inversion in Section 3.2 were applied to the numerical model calculation, where the weakening coefficient λ of surrounding rocks was taken as 0 under the dry state. After the roadway excavation was stabilized, the strain and time of the roadway vault, left wall, and right wall were monitored. Next, the monitored data were fitted using Eqs 3 and 4, in which the stress level was the average stress within each time period. Table 4 lists the main fitted parameters under the dry state, and Figure 6 displays the fitting results. As seen in Table 4, the creep parameters EM, ηM, EK, and ηK basically had consistent orders of magnitude with those in Table 3, accompanied by high correlation coefficients. Thus, Eqs 3 and 4 could accurately describe the surrounding rock deformation-time relations under the dry state. The plastic deformations ε p at roadway vault, left wall, and right wall were all small, indicating that surrounding rocks were subjected to viscoelastic deformation in creep stages under the dry state. Table 5 lists the main parameters fitted through Eqs 3 and 4 under the damp state, and Figure 7 illustrates the fitting results. As could be known from Table 5, the creep parameters EM, ηM, EK, and ηK basically had consistent orders of magnitude with those in Table 3, accompanied by high correlation coefficients. This finding manifests that Eqs 3 and 4 could accurately describe the surrounding rock deformation-time relations under the damp  state. Moreover, the plastic deformations ε p at roadway vault, left wall, and right wall were all smaller than those under the dry state, indicating that due to the increase in water content of the surrounding rocks under the damp state, the viscoelastic deformation of surrounding rocks was more obvious in the creep stages. The λ value was kept at approximately 0.6, indicating that the surrounding rocks substantially deteriorated by the dampness-induced argillization so that the final deformation of surrounding rocks was greater than that under the dry state. As shown in Figures 6, 7, the monitoring data under dry and damp states were fitted very well by Eqs 3 and 4. Under the damp state, the maximum strains at roadway vault, left wall, and right wall were 0.29, 0.41 and 0.23%, respectively, and those under the dry state were 0.3, 0.35, and 0.21%, respectively. This finding indicates that the dampness-induced argillization of surrounding rocks impacted the left wall of the roadway in a titled stratum but slightly influenced the roadway vault and right wall. Figure 8 displays the strain-time relations under dry and damp states. Figure 8 shows that under both dry and damp states, the creep of surrounding rocks at different positions presented typical decay creep (3-18 days), uniform creep (18-46 days), and accelerated creep (46-52 days) with the passing of time. During the creep process, the strains at the roadway vault and right wall under the damp state were slightly greater than those under the dry state, and the strain at the left wall under the damp state was significantly greater than that under the dry state.
As aforementioned, the creep rates of surrounding rocks at different positions in the strain decay and accelerated stages were both higher than that in the uniform stage, thereby using the creep rate in the uniform stage as reference for analysis. Table 6 displays the creep rates of the surrounding rocks in different stages under dry and damp states. Under the dry state, the average creep rate of surrounding rocks at different positions in the decay stage was 1.97-3.74 times of that in the uniform stage, and that in the accelerated stage was 12.69-18.46 times of that in the uniform stage. Under the damp state, the average creep rate of surrounding rocks at different positions in the decay stage was 2.14-6.78 times of that in the uniform stage, and that in the accelerated stage was 12.73-20.14 times of that in the uniform stage. By comparison, the strain of surrounding rocks under the damp state was increased compared with under the dry state. Thus, the dampness-induced argillization accelerated the damage inside surrounding rocks. Finally, the damp state was positively correlated with the strain of surrounding rocks.

Determining the Creep Range for Surrounding Rocks
After the roadway excavation was stabilized, surrounding rocks creeped, thus slowly reducing the roadway cross section. As could be known from Figure 5, the horizonal convergence deformation at two roadway walls was greater than the subsidence deformation at the roadway vault. Accordingly, the horizontal displacement was taken as an example to analyze the creep range of surrounding rocks. Figure 9 exhibits the horizontal displacement nephograms of surrounding rocks in the roadway under dry and damp states within different time periods.
As could be observed from Figure 9, the scope of influence of horizontal displacement was gradually enlarged with the lapse of creep time. At 3-18 days, the scope of influence of horizontal displacement of surrounding rocks was small, and that under the dry state was basically consistent with that under the damp state. This finding indicates that under both states, the surrounding rocks deformed slowly at 3-18 days, during which the deformation of surrounding rocks had little to do with the dampness-induced argillization state of surrounding rocks. At 18-46 days, the scope of influence of horizontal displacement was gradually enlarged, and that under the damp state was slightly larger than that under the dry state. This finding indicates that the damage of surrounding rocks under the damp state was slightly greater than that under the dry state. At 18-46 days, the deformation of surrounding rocks was positively correlated with the dampness-induced argillization state of surrounding rocks. At 46-52 days, the scope of influence of horizontal displacement of surrounding rocks substantially expanded, and that under the damp state was larger than that under the dry state. This finding manifests that the damage of surrounding rocks under the damp state was greater than that under the dry state. Consequently, it took less time for surrounding rocks to experience the same deformation. In this stage, the deformation of surrounding rocks showed a significant positive correlation with their damp state. The above analysis results showed that after the argillization of surrounding rocks when encountering water, their deformation was enlarged. However, the dampness-induced argillization state had an obvious time effect on the damage of surrounding rocks, and surrounding rocks were damaged more evidently under the damp state after a certain time point. The stability of surrounding rocks could be reflected by their deformation rate. When the deformation was smaller than 0.02 mm/d, the surrounding rocks could be considered to be basically stable (Hou, 2017). Under the dry state, Figure 10 displays the horizontal displacement-time relations at different depths on the left side of surrounding rocks in this roadway. As could be known from Figures 6, 7, after 46 days, the surrounding rocks were in an accelerated creep stage. Thus, the creep range of surrounding rocks was distinguished by keeping their deformation rate at 0.02 mm/d (cumulative horizontal displacement: 0.00092 m) within 46 days.
As observed from Figure 10, the horizontal displacement of surrounding rocks was in a direct proportion to time and inversely proportional to the depth, and the cumulative horizontal displacement of surrounding rocks within the depth range of 1.0-5.5 m of surrounding rocks on the left side was greater than 0.00092 m. Thus, the creep range of surrounding rocks on the left side of this roadway was 1.0-5.5 m. Similarly, under the dry state, the creep ranges of surrounding rocks above the roadway at its right side and beneath it were 1.0-5.0, 1.0-4.5, and 1.0-3.5 m, respectively. Hence, the creep range of surrounding rocks was asymmetric under the dry state, and the overall creep range presented a clockwise reduction from the surrounding rocks on the left side of this roadway.
Under the damp state, the creep ranges of surrounding rocks on the left side of this roadway, above it, on its right side, and beneath it were 1.0-5.5, 1.0-5.5, 1.0-4.5, and 1.0-4.0 m, respectively. This finding manifests that the creep range of surrounding rocks under the damp state was relatively consistent with that under the dry state, but that under the damp state was larger than that under the dry state.

Analysis of Roadway Cross-Sectional Deformation
Owing to the creep of surrounding rocks in this roadway, the roadway cross section was reduced. Additionally, its shape also deviated from the original design, thus impacting the normal basic functions of this roadway. Figure 11 presents the timedependent changes in relative horizontal and vertical displacements at main monitoring points on the rectangular roadway cross section under the dry state.
As could be known from Figure 11, under the dry state, the relative horizontal displacement of roadway cross section in a tilted stratum under the dry state was -1.5 to 1.5 cm, and the relative vertical displacement was -2.2 to 1.5 cm. In addition, the relative displacement at each monitoring point was in a direct proportion to time. At monitoring point d, relative vertical displacement played a dominant role, which led to floor heaving of this roadway. At monitoring points b and f, relative horizontal displacement dominated, and the horizontal displacement at point b was greater than that at point f. Thus, the left wall largely caved than the right wall did. At monitoring points a, c, e, g, and h, both relative horizontal and vertical displacements occurred, leading to the inward reduction of roadway cross section, and the displacement at point h indicated the vault settlement and deviation to the left. At monitoring points b, f, h, and d, the relative displacement substantially influenced the cross-sectional change in this roadway, the overall cross-sectional change was asymmetric (approximate to a deflecting rhombus), and finally, the crosssectional shape was basically identical with that in relevant literature . Figure 12 shows the time-dependent changes in relative horizontal and vertical displacements at main monitoring points on the rectangular roadway cross section under the damp state.
As could be known from Figure 12, the main factors influencing the relative displacement-time relations and crosssectional changes of this roadway in a tilted stratum under the damp state were basically consistent with those under the dry state, the relative horizontal displacement was -1.2-1.6 cm, and the relative vertical displacement was -2.5 to 1.6 cm. The relative displacement range in both horizontal and vertical directions was larger than that under the dry state. Thus, under the damp state, the damage in the surrounding rocks was aggravated, the overall roadway cross section changed asymmetrically (approximating a deflecting trapezoid), and the final roadway cross-sectional shape basically coincided with the relevant literature .

CONCLUSION
1) The Cvisc model considering the weakening coefficient of surrounding rocks is applicable to the decay and uniform creep stages, and that considering damage is applicable to the accelerated creep stage. During the creep process, surrounding rocks were subjected to viscoelastic deformation; 2) At 3-18 days, surrounding rocks basically creeped under the dry state identically with that under the damp state. After 18 days, the creep range under the damp state expanded than that under the dry state, and the damage in surrounding rocks under the damp state had an obvious time effect; 3) Under the dry state, the maximum creep ranges of surrounding rocks on the left side of the roadway, above it, on its right side, and beneath it were 5.5, 5.0, 4.5, and 3.5 m, respectively. Those under the damp state were 5.5, 5.5, 4.5, and 4.0 m, respectively. The overall creep range was asymmetric, being reduced clockwise from the left side of this roadway; 4) Under the dry state, the relative horizontal displacement of this roadway was -1.5 to 1.5 cm, and the relative vertical displacement was -2.2 to 1.5 cm. The cross section of this roadway changed asymmetrically, resembling a deflecting rhombus or trapezoid under the dry state. Additionally, the change laws under the damp state were basically consistent with those under the dry state, but the roadway cross section largely changed under the damp state in an approximately deflecting trapezoid.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
Conceptualization, FC; methodology, LY; software, RS; validation, XY; formal analysis, HJ; investigation, CZ; writing-review and editing, RS. All authors have read and agreed to the published version of the manuscript.