Study on Creep Constitutive Model of Stratified Siltstone and its Application to Environmental Hazards Treatment


 Mining causes environmental hazards which have an adverse effect on ecological system and human living environment. Rock is the main body of underground coal mining, and among them stratified siltstone is very common. Creep and instability of underground stratified siltstone is one of the main inducing factors of ground fissures and ground surface settlement. In this study, on the basis of fully considering its stratified structure features, we built and theoretically analyzed the viscoelastic-viscoplastic constitutive model of stratified siltstone. Creep tests and numerical simulation under multi-stress were conducted. Numerical simulation results under high, medium and low stress were consistent with the creep test results, it made up for the deficiency of traditional creep models in predicting plastic deformation at low stress level. Numerical simulation results showed that under different boundary conditions, the proposed model had advantage to more accurately simulate creep behavior of stratified siltstone compared with traditional models. Also, the numerical simulation results reveled creep characteristics of underground roadway surrounding rock, based on which a supporting scheme was put forward. In the stress adjustment stage, the strain rate dropped significantly from 14.6mm/d to 10.3mm/d. Then, it entered the deceleration creep stage. Creep of roadway surrounding rock and ground surface settlement stopped. The environmental hazards were solved to a great extent.


Introduction
Ground fissures and ground surface settlement caused by extraction of coal are man-made environmental hazards. They have a negative influence on ecological system and human living environment (Fig. 1), and threaten safe mining and energy sustainable development. In terms of siltstone with durative service life and native stratified structures, creep is a key factor inducing instability of underground engineering, and in turn aggravates environmental hazards [1][2][3][4][5]. Creep behavior of stratified siltstone, due to its special stratified structure, are different from those of other rock such as granite. Moreover, with further extension of underground engineering to a deeper level, stress state of stratified siltstone becomes worse. Stratified siltstone is in the state of coexistence of energy storage and energy release during creep. When energy storage is dominant, stratified siltstone becomes carrier of concentrated stress. Once the energy storage or concentrated stress reaches its threshold, failure of stratified siltstone will be transformed into dynamic disaster, which is easy to spread to ground surface, inducing ground surface disaster [6][7][8][9]. In particular, creep and synergistic dynamic instability of stratified siltstone are very likely to be encountered during underground mining. Therefore, the study of creep behavior of stratified siltstone is the key to forecast and control environmental hazards.  [10][11][12]. As the research moved along, different behavior of rock in various mechanical environments were gradually discovered, and the traditional creep constitutive model showed its limitations. In this case, a series of rheological models were proposed through in-depth research. These models include viscoelastic elements and elastoplastic elements, which are able to simulate aging and non-aging behavior of rock [13], and has a wider scope of application. However, as rock is a kind of natural material whose physical and mechanical properties has significant anisotropy, it is necessary to take this key factor into consideration when building constitutive model. In recent years, scholars have been devoted to studying constitutive model of various rocks, and obtaining beneficial conclusions.
Based on the theory of rock creep and fractional calculus operator theory, Pu [14] established integer and fractional percolation-creep model, and obtained the idea of underground surrounding rock stability assessment in a water-rich environment. Shan [15] constructed a creep constitutive model of layered red sandstone in frozen state under low temperature with the independently developed experimental equipment, which provided a method for studying safety evaluation of wellbore and related projects during freezing construction. Wang [16] used creep experiment and theoretical analysis methods to study creep behavior of coal-rock mass, built and verified a nonlinear creep constitutive model. Tomanovic J. Hickman [23] proposed a rate-dependent constitutive model that could predict the change of axial yield stress with the power function of applied axial strain rate, which was able to reproduce rate-dependent behavior of porous rocks under various loading conditions. D. Sterpi [24] discussed influence of viscoelasticity (primary) and viscoplasticity (secondary) on the rheological model and analyzed influence of creep on tunnel shrinkage. F. Pellet [25] connected viscoplastic constitutive relation with damage theory and established a rock constitutive model that considered viscoplastic behavior and damage evolution over time. Mei [26] constructed a nonlinear creep model with five elements by fractional calculus, and determined model parameters through rock creep test, fitting results showed that the creep model was reasonable and feasible. Based on the creep experimental data obtained for gneiss, Zhang [27] proposed a creep composite model, and obtained accurate creep behavior characteristics and mechanism of gneiss. B. Nedjar [28] proposed a local damage model based on the viscoelastic coupling continuous damage mechanical framework, and conducted in-depth research.
Heusermann. S [29] built a LUBBY2 constitutive model that could describe the nonlinear creep behavior of rock salt, and gave an example of the finite element analysis of the cavity, as well as the corresponding stability and availability criteria. S. Olivella [30] established a creep constitutive model based on two deformation mechanisms of fluid-assisted diffusion transfer creep and dislocation creep, this model was able to predict the creep strain rate, and the simulation results were highly consistent with those of experiments.
In this study, a viscoelastic-viscoplastic constitutive model was proposed and theoretically analyzed, creep tests on stratified siltstone were also conducted to explore creep behavior of stratified siltstone. In order to verify scientific nature and accuracy of the model, governing equation of the model was introduced into numerical calculation program to simulate a series of rock creep test. Numerical simulation results of the model was consistent with creep test data. We also compared the proposed model with the simple elastic-plastic model, the Kelvin model and the Nishihara model. Results are shown that under different boundary conditions, the former could more accurately simulate creep properties of stratified siltstone. Therefore, the proposed model has reference value for creep analysis and stability evaluation of stratified siltstone. In addition, the results of numerical simulation reveled creep characteristics of roadway, based on which a supporting scheme was put forward and environmental problems were solved to a great degree.

Theoretical analysis of the viscoelastic-viscoplastic constitutive model
Where, ε T is the total deformation, and ε e is the strain generated by an elastoplastic component, strain of this part is generated by the elastic element in the component, which is completely reversible and has no aging property.
Stress-strain relationship of spring element with elastic coefficient KM in the viscoelastic component is as The equation is related to external tension exerted by internal stress of each component, internal equilibrium equation of the viscoelastic-viscoplastic constitutive model can be obtained by equation For viscoelastic component, the internal stress is = . By combining strain decomposition component and internal constitutive equation, it can be obtained: Thus, the flow rule of ( ) is: Colonnetti minimum rule is used and the plastic strain is regarded as additional deformation. When both the plastic strain and the total strain are known, the differential equation of above equation can be expressed as: The convolution law of stress can be obtained: Variable mediation is conducted for the above equation, the following equation can be obtained: Where, EU is uniform elastic parameter, which is related to specific material type. ̇( ) is a known value available from creep tests. Viscoplastic strain rate ̇( ) is a part of mechanical response, the specific value of which cannot be determined by the above known conditions. Therefore, it is necessary to refer to classical plasticity theory to evaluate and obtain this parameter. In a mathematical sense, viscoplastic strain is the integral of viscoplastic strain rate on the time axis, and its mathematical expression is: Viscoplastic activation equation of the viscoelastic-viscoplastic constitutive model can be defined as: According to the standard plastic judgment theory, When A(σ,γ)＜0, viscoplastic component is inactive, and viscoplastic component is activated when stress increases to A(σ,γ)＞0. Stress acting on the viscoplastic element produces viscoplastic strain, and viscoplastic strain rate can be calculated by the flow law: The above equation defines viscoplastic strain rate of the viscoelastic-viscoplastic constitutive model.
From the above equation, it can be seen that the internal strain of the viscoplastic component increases with viscoplastic strain, which leads to increase of viscoplastic activation threshold.
In the numerical simulation, when calculation time step reaches TN+1, integral relation of the convolution law of stress is: Where, t0 is the start time of numerical simulation, viscoplastic strain of the viscoplastic component equals to 0, and integral term of the above equation is: Where, Equation (14) can update time integral function at each calculation step during numerical calculation. So, if the initial value ( ) is given, results of next step can be obtained by the next iteration increment Δ ( ) , which can be defined as: Where, △t is time step, △t=tn+1-tn. Δ , Δ and Δ represents total strain increment, elastic strain increment and viscoplastic strain increment respectively. The stress at the end of the calculation step can be calculated by the following equation: .
Equation (16) shows that updated stress at the end of calculation step can be calculated by three variables.
The first part is total stress increment +1 , this parameter is viscoelastic stress response to a given total strain increment. The second part is viscoplastic stress correction caused by the increase of viscoplastic strain. The third part is the elastic stress response caused by elastic strain. Bvp is equivalent time step viscoplastic modulus.
The viscoplastic activation equation at the n+1 step is as follows: In the calculation step n+1, the above equation can be expressed as: Viscoplastic strain increment Δ is needed to obtain stress update value at the end of each step.
Step increment of the viscoplastic flow equation is as follows: Step increment can be expressed as: Where, In Equation (20), Bvp, ω and △γ are positive, so when +1 is positive, step increment of viscoplastic strain are not identically zero. When increment of the viscoplastic flow equation is substituted into the equation, we can get: The increment △γ can be obtained by solving the above equation: When △γ remains positive, the increment between steps of viscoelastic-viscoplastic computation is solved. The viscoplastic strain increment can be obtained through Equation (19). And the viscoelastic-viscoplastic constitutive model established in this study can be determined through Equation (10)- (22). The calculation process follows traditional discrete integral method, while evolution relationship of change rate of viscoplasticity is solved by the Euler difference method.

Model parameter identification
When the viscoelastic-viscoplastic model is used for numerical calculation, it is necessary to identify parameters in the model. In the standard creep test, the ε-t curve is usually exponential with a constant e as base. In general, it is difficult to find time point of relaxation since the curve has good continuity and transition. To solve this problem, data difference magnification processing method can be adopted to find relative significant difference between data, and the processed data can be plotted to locate corresponding relaxation time point. As shown in and K3 can be obtained from Equation (23). In order to identify and obtain viscoplastic coefficient, it is necessary to conduct creep tests on stratified siltstone under various stress. In addition, conventional monotonic increasing stress experiments were carried out to obtain the slope mean of σ-ε curve in the plasticity domain, and the expression is: Plastic hardening coefficient ω can be obtained from Equation (24), and viscosity coefficient ηv can be calculated by the following equation: Where, μv is built-in viscoplastic time, and its solution method is the same as relaxation time μi.

Fig. 3
Strain-time curve after data amplification Fig. 4 Creep test results of stratified siltstone

Computation validation based on creep tests
Creep test of stratified siltstone was carried out to obtain creep parameters. The rock samples were shaped on a cutting machine, then, the samples were polished on both ends. Its positive-negative error should be controlled within 0.1mm. The sample diameter was 50.0mm, the height was 100.0mm, and the ratio of height to diameter was 2:1. 0.85σc(σc =35.8MPa) stress was applied to load the sample for 4320.0s. The experimental results were shown in Fig. 4. The curve corresponds to three relaxation time points, thus corresponding to three sets of viscoelastic parameters. In the model calculation part, three-dimensional finite difference program FLAC3D was used to establish numerical model with the same size as sample used in the standard creep test. Input parameters required for calculation were shown in Table 1 In the creep test, uniaxial pressure test machine was used. Samples were taken from Jiangou coal mine in Urumqi, Xinjiang province, China. The number of samples was 9. 3 samples in each group, which were numbered A, B and C respectively. The specific experimental scheme was as follows: 0.6σc (21.5MPa), 0.7σc (25.1MPa) and 0.8σc (28.6MPa) was applied to group A, B and C, respectively. Constant pressure was maintained until the samples were damaged (Fig. 5). By comparing and analyzing parameters obtained from samples, the average value of data within the acceptable range of discreteness was calculated according to the standard of applied stress. Then average values of all data were analyzed to reduce or avoid the possibility of large deviation data.
(a)σ=21.5MPa  The stratified siltstone samples showed combined failure modes of sliding failure along structural plane and shear failure through structural plane. When stress applied to stratified siltstone samples was high enough or there are defects in the lithospheric structural plane, failure path not only developed along soft materials between bedding planes, but also broke through structural plane containing defects. Thus, resulting in combined sliding failure along structural plane and shear failure through structural plane (Fig. 5). Fig. 6 illustrates variation between experimental data at the loading stage and numerical calculation results obtained from the viscoelastic-viscoplastic constitutive model. According to the experimental results, stress paths were numerically simulated, and creep curves under loading stress 0.6σc (21.5MPa), 0.7σc (25.1MPa) and 0.8σc (28.6MPa) were obtained. It can be seen that for instantaneous deformation and creep, data obtained from creep test on stratified siltstone samples are highly consistent with the data acquired by the model. The model has a good performance in terms of strain rate, which can be attributed to the fact that plastic zone is more likely to develop under action of high stress level in the creep process of stratified siltstone. The creep strain of samples under low stress were analyzed, results showed that hardening behavior of the model can compensate for the deficiency of viscoplastic model in predicting plastic deformation under low stress.

The influence of loading type on the viscoelastic-viscoplastic constitutive model
Loading time has key influence on behavior of the viscoelastic-viscoplastic constitutive model. Therefore, we adopted numerical simulation of short-time loading and long-time loading based on two-dimensional simulation. Samples were 100.0 cm in height and 50.0cm in diameter with a 2:1 ratio of height to diameter. In order to simulate mechanical behavior of stratified siltstone samples, the input mechanical parameters were obtained through creep tests (Table 2). By comparing numerical simulation results of the viscoelastic-viscoplastic constitutive model, the behavior evolution characteristics of the model were studied. Short-time loading and long-time loading were simulated to obtain influence of loading time on the viscoelastic-viscoplastic constitutive model. All numerical simulations were multi-stage simulations (18 stages). In the short-time loading numerical simulation, the viscoelastic-viscoplastic constitutive model was calculated by solving command in each stage until unbalanced force fell below specified threshold and operation state was automatically stopped. The short-time loading numerical simulation was used to determine the minimum time required for threshold to reach equilibrium at each load or strain stage. In the long time loading numerical simulation, number of the shortest solution steps was doubled to check mechanical behavior of the model over time under constant conditions. After the calculation of the loading stage was completed, the unloading stage was carried out. Calculation step (or time) of unloading stage was the same as that of loading stage. In addition, after each stress or strain retention period (i.e. load -hold -unload and repeat to the next loading and unloading process), the strain change pattern accumulated throughout the test period was checked. The increment of axial stress applied on the sample was 1.5MPa, with a total of 20 fixed incremental load values, until constant load threshold reached 30.0MPa. At each constant load threshold, axial (vertical) strain was recorded.
During the numerical simulation, it was observed that samples showed elasticity at the beginning and then started to yield. As shown in Fig. 7, strain continued to increase during unloading process. In short-term and long-term loading modes, with increment of load, strain also increased. However, with extension of load retention period in the loading step, strain at each load threshold in the long-term loading was larger than that of the corresponding short-term loading model. Strain of long-term loading model was 43% higher than that of short-term loading model. When the model ran unloading stage after the loading and holding stages, continuous strain of the models were less than the strain in the loading stage without repeated unloading stage. Unloading stage before transition to next load threshold allowed strain dissipation of the model, which verified the fact that some strains in the model were recoverable. However, total strain at the unloading stage was still larger than that at the stage of maximum loading, indicating that the model produced unrecoverable strain at the beginning of yield.
Loading method also has key influence on behavior of the viscoelastic-viscoplastic constitutive model. Therefore, we adopted numerical simulation of static load (stress control) and constant strain rate (strain control) based on two-dimensional simulation. In addition, in order to obtain the improved properties of the viscoelastic-viscoplastic constitutive model, Kelvin model, Nishihara model and viscoelastic-viscoplastic constitutive model were selected as the target models for comparative analysis.
(a) Short-term loading (b) Long-term loading Fig. 7 Stress-strain relationship of viscoelastic-viscoplastic constitutive model under different loading type Fig. 8 Deformation characteristics of models under various loading modes (K, N, and V represent the Kelvin model, the Nishihara model, and the viscoelastic-viscoplastic constitutive model, respectively. L and S represent long-term loading and short-term loading, respectively. "σ" and "ε" represent stress control and strain control, respectively.) Fig. 8 describes the variation characteristics of the axial strain with the calculation step in the long-term loading and short-term loading modes under different boundary conditions. Comparison objects include the Kelvin model, the Nishihara model and the viscoelastic-viscoplastic constitutive model. It can be seen from the figure that a longer loading time under the same axial stress or axial strain produced a larger strain. Since the Kelvin model only captured the first-stage creep, its strain was smaller than that of the viscoelastic-viscoplastic constitutive model and the Nishihara model, but the model did not yield. The strain of the Nishihara model was larger than that of the viscoelastic-viscoplastic constitutive model over a long period of time without instability. While the viscoelastic-viscoplastic constitutive model began to yield when the final load cycle was reached. The viscoelastic-viscoplastic constitutive model requires less stress to achieve the same strain under strain control. In the stress control experiments, once a constant target stress level was reached, stress level oscillated as a result of the model adjusting to a new boundary condition. And it was controlled by the model's damping (the stiffness of the grid, density and region). On the contrary, strain control experiments showed that when the target strain was reached, the stress would increase gradually without oscillation.

Performance analysis of the viscoelastic-viscoplastic constitutive model
In order to analyze the performance of each model in simulating underground surrounding rock, the circular roadway with diameter of 10.0m and depth of 400.0m in the stratified siltstone were simulated. The mechanical parameters assigned to the roadway model were the same as those studied in the above models. In order to avoid influence of stress reflection, 10 diameters of the roadway were designed from center of roadway to model boundary. Bottom displacement of the model was constrained. Stress boundary was set to reflect only the gravity load, and grid was arranged radially outward without influence of supporting equipment. Focus of the analysis is to compare mechanical behavior of the selected models to explore long-term behavior of an unsupported roadway. Fig. 9 shows displacement of surrounding rock in the roadway after excavation. Results show that displacement of the viscoelastic-viscoplastic constitutive model was similar to that of the Nishihara model, and the Nishihara model produced smaller strain. The Kelvin model did not produce more displacement with time, and surrounding rock of roadway produced a certain amount of displacement when the Kelvin stiffness component of the model played its role. Compared with the Nishihara model, the Kelvin model could reach a stable state with infinite strain potential. It should be noted that there is no obvious relationship between actual (creep) time and number of steps in the model, as it is directly dependent on the time required for software and model to reach equilibrium. The study also showed that deformation of the viscoelastic-viscoplastic constitutive model is similar to that of the Nishihara model. In terms of deformation, the viscoelastic-viscoplastic constitutive model had larger deformation than the Kelvin model, so it was considered that the viscoelastic-viscoplastic constitutive model generated additional displacement due to viscous potential of rock.
Under control of stress and strain, static load experiment scheme was numerically studied, and mechanical behavior characteristics of the simple elastic-plastic model, Kelvin model, Nishihara model and viscoelastic-viscoplastic constitutive model were analyzed in detail (Fig. 9). Compared with the other three models, the viscoelastic-viscoplastic constitutive model was viscoelastic before entering yield state and converted to be plastic after entering yield state. In general, the viscoelastic-viscoplastic constitutive model can simulate viscoplastic properties of rock (time dependence of the plastic shear) and relate plastic yield strength to the previous creep process. From the phenomenological point of view, plastic yield stage can be matched with the third creep stage observed in the actual situation, and evolution law of the plastic creep stage of rock can be accurately simulated in mechanics.

Application of the viscoelastic-viscoplastic constitutive model to field problems
Roadway surrounding rock is the key part to support underground stope. Creep and widespread instability of roadway surrounding rock lead to mechanical imbalance, resulting in the continuous progressive destruction. When failure range of the surrounding rock reaches a critical value, it will lead to surface settlement and ground fractures. Therefore, it is necessary to support the surrounding rock and control its creep before it is completely destroyed. Based on the discrete element calculation program 3DEC, a numerical model was established based on mining condition of the steeply inclined coal seams (the average dip angle of coal seams is 85.0°). The viscoelastic-viscoplastic model was imported into the calculation program, and the input parameters could be obtained through rock mechanical tests (Table 3). Results showed that the surrounding rock converged significantly to the interior, and the amount reached 1844.0mm (Fig. 10), which was consistent with the field measurement results. However, there were great differences in creep locations and fracture modes of the roadway's sides. Creep on the left side showed synthetic bucking beginning in the central position of excavation free face, and the maximum amount was about 708mm. Creep on the right side showed shearing slip in the corner position of free face, and the maximum amount was about 1136.0mm. The sliding rock mass were piled up on the floor, which seriously affected function of the roadway. Also, large deformation of roadway surrounding rock led to strong deformation and destruction of working site, which further caused ground fissures and ground surface settlement.

Application and evaluation of the optimization support scheme
In order to prevent creep of roadway surrounding rock from being too large, an optimization support scheme was put forward (Fig. 11). According to the regional creep characteristics of the surrounding rock on both sides of the roadway, equal strength anchor bolts were installed in the middle and bottom of the roadway on the left and right side to reinforce the surrounding rocks. Bolt size was φ 25.0mm× L3000.0mm and row spacing was 750.0mm. Pretension and anchoring length of bolt were 180.0N•m and 1500.0mm respectively. In order to prevent collapse and instability of rock at the top of roadway, steel strand (size 15.2 × L9000.0mm) was arranged in the middle line position at the top of the roadway. The azimuth was shifted 30.0° to the north, and the row spacing was set as 3000.0mm. Fig. 11 Optimization support scheme of roadway Under the condition of the original support scheme, a section 550.0-650.0m of No. 3-6-1 roadway was selected as monitoring section. Depth of detection hole of each station was designed to be 3.5m. A total of 11 stations were arranged in the area, with the interval of 10.0m, numbered 1~11. Monitoring period was from March 15 to April 15, 2019. It could be found that the surrounding rock began to produce relatively large deformation after the roadway was built. And the deformation generated in a short time accounted for about 37.4% of total deformation (Fig. 12a). Subsequently, from March 17 to March 21, the roadway surrounding rock entered the stress adjustment stage. The data from stress monitoring equipment showed that stress presented a rapid downward trend, and the corresponding strain rate continued to decrease rapidly, with the minimum value dropping to 14.6mm/d. During April 10 to April 15, the creep rate of surrounding rock at some stations such as 3#, 5#, 6#, 7#, 10# and 11# increased rapidly. The fractures in the surrounding rock were fully developed and interwoven fracture networks were formed, transforming the overall surrounding rock into loose structure. The surrounding rock at the two sides of roadway was extruded in a large range, and rock bolts at some positions were deformed or even fell off with surrounding rock. In order to evaluate support effect of the optimization scheme, the whole distance multi-point displacement monitoring was adopted for section 700.0-800.0m of No. 3-8-15 roadway after implementation of the optimization scheme. Depth of detection holes of each station was designed to be 3.5m. A total of 11 stations were arranged in the area, with the interval of 10.0m, numbered 1~11. By analyzing the monitoring data from May 20 to June 20, 2019, it could be found that large-scale deformation occurred in the initial stage, and the deformation generated at this stage accounted for about 43.8% of total deformation (Fig. 12b). From May 21 to May 24, the surrounding rock began to enter stress adjustment stage. At this stage, it was observed that the stress dropped in a leaping way and the strain rate continued to decrease correspondingly. On May 23 and May 24, it reached the minimum value of this stage, which was 10.3mm/d. Then, it entered the deceleration creep stage from May 25 to May 28. During June 16 to June 20, creep of roadway surrounding rock tended to stop. In addition to the rock bolt in the No.1 station fell off due to construction reasons, leading to instability of surrounding rock. Monitoring results of the other stations showed that the surrounding rock stopped producing deformation. The surrounding rock showed good anti-deformation ability after implementation of optimization support scheme, the creep of surrounding rock was alleviated, and the surface settlement and ground fissures of overlying strata were improved to a large extent. In order to acquire the specific information of the surface settlement and ground fissures of overlying strata, we built the exploratory trench. Length, width and depth of the exploratory trench are 20.0m, 7.0m and 10.0m respectively. The exploratory trench revealed five sets of main strata, which are plain fill, silt, silty silt, sand and clay. The main fracture has a dip angle of 73°. The number of dislocation generated by the main fracture to the formation gradually increases with the increase of the depth. The dislocation generated by the main fracture to the bottom of the plain fill and fine sand is 15.0cm and 170.0cm, respectively. The secondary fracture is located in the footwall and produces 20.0cm dislocation for the coarse sand layer and 36.0cm dislocation for the coarse sand interlayer in the clay layer. The strata near the fractures dip downward, showing a clear traction structure. Large-scale ground fractures were no longer generated in the overlying strata (Fig. 13), and ground surface settlement was significantly improved. On the whole, the optimization support scheme had a remarkable effect.

Conclusions
It is of practical significance to study the creep behavior of stratified siltstone for the underground safe mining and environmental protection. The establishment of scientific and accurate creep constitutive model of stratified siltstone is a key factor to solve the problem. In this study, the viscoelastic-viscoplastic constitutive model was established. Numerical simulation and rock creep test were adopted to verify scientific nature and accuracy of the model. Specifically, numerical simulation results of the viscoelastic-viscoplastic constitutive model at high, medium and low stress levels were consistent with the results of creep tests, which made up for deficiency of other similar creep models in predicting plastic deformation at low stress levels. Compared with the simple elastic-plastic model, Kelvin model and Nishihara model, the viscoelastic-viscoplastic constitutive model was able to simulate creep characteristics of stratified siltstone more accurately in various boundary conditions. Furthermore, creep characteristics of roadway surrounding rock was obtained according to numerical simulation. Accordingly, we proposed a supporting scheme and the evaluation results showed that environmental problems were solved. Declarations Funding: This research was funded by the General Project of Shaanxi Provincial Department of Education, grant number 21JK0952, the High-level Talents Special Fund of Xijing University, grant number XJ21B12.

Conflicts of Interest:
The authors declare no conflict of interest. The funders were responsible for a part of work in the design of the study.