A Nonlinear Creep Damage Model Considering the E ﬀ ect of Dry-Wet Cycles of Rocks on Reservoir Bank Slopes

: Water has a crucial e ﬀ ect on the time-dependent behavior of rocks. The long-term cyclical ﬂuctuations of reservoir water level lead to dry–wet (DW) cycles of rocks on reservoir bank slopes, making this inﬂuential factor more complex. To deeply understand the time-dependent behavior of rocks under DW cycles, argillite from the reservoir bank slope of Longtan Hydropower Station was used to perform a series of triaxial creep tests. Subsequently, based on analysis of creep test results after di ﬀ erent DW cycles, a damage nonlinear Burgers viscoelastic-plastic (DNBVP) model considering the e ﬀ ect of saturation–dehydration cycles was proposed by introducing a nonlinear viscoplastic body and a damage variable describing DW cycles. Then, the three-dimensional creep equations of the new model were derived and its creep parameters were identiﬁed. Comparison between the theoretical curves and the test results shows that the theoretical curves of the DNBVP model were able to describe rock creep tests results after di ﬀ erent DW cycles. Furthermore, by comparing classical creep models with the proposed model, it was found that the DNBVP model can accurately reﬂect the nonlinear characteristics of rocks at the accelerated creep stage. Finally, the sensitivity of the DNBVP model was analyzed and discussed, and three-dimensional central di ﬀ erence expressions necessary for secondary development of the new model were also derived in detail. The proposed new model with secondary development may provide a basis for improving the geotechnical design of reservoir bank slopes and the control of reservoir bank landslides.


Introduction
The time-dependent behavior of rocks is closely related to the safety of geotechnical engineering [1,2]. Numerical studies of engineering practices and other studies show that rock instability and failure have direct relationships with time [3,4], and variation in the mechanical properties of rocks over time can easily lead to natural disasters, especially landslides [5][6][7][8][9].
Water is considered a crucial factor in the time-dependent behavior of rocks [10][11][12][13]. With the construction and operation of numerous hydropower projects in China's western regions, a series of high and steep reservoir bank slopes have appeared [14,15]. With reservoir impoundment, the reservoir water level undergoes dramatic cyclic fluctuations every year, resulting in dry-wet (DW) cycles affecting rocks in the hydro-fluctuation belts of reservoir bank slopes. The effect is a cumulative process, with the damage and deterioration of rocks accelerating over time [16,17], leading to the development (iii) For creep models considering the effect of DW cycles, it is necessary to derive equations necessary for secondary development.
To address the problems above, argillite obtained from the reservoir bank slope of Longtan Hydropower Station was used to perform a series of triaxial creep tests after different DW cycles in this study. Firstly, the experimental data obtained in the creep tests were analyzed. On this basis, a nonlinear damage model accounting for the effects of saturation-dehydration cycles was proposed. Furthermore, the new model was verified and parameter sensitivity analyses of the proposed model were performed. Finally, the three-dimensional central difference expressions necessary for secondary development of the new model was derived in detail.

Preparation of Test Rock Specimens
Argillite was obtained from the reservoir bank slope of Longtan Hydropower Station in China (see Figure 1). In the dam site area, argillite makes up about 30.8% of the strata of the slope, and argillite in the "hydro-fluctuation belt" (Figure 1c) is subject to obvious creep under the influence of reservoir water level fluctuations. The sampling site is shown in Figure 1d. After the drilling of the boreholes, argillite samples (predominantly gray) were obtained in situ. Then, rock samples were packed and transported to the laboratory and processed into cylinders with a diameter of 50 mm and height of 100 mm, according to the GB/T50266-2013 standard for test methods of engineering rock mass [30], as shown in Figure 1e. Table 1 lists the basic mechanical parameters of the argillite samples.
Water 2020, 12, x FOR PEER REVIEW 3 of 17 properties of rocks after different DW cycles. (iii) For creep models considering the effect of DW cycles, it is necessary to derive equations necessary for secondary development.
To address the problems above, argillite obtained from the reservoir bank slope of Longtan Hydropower Station was used to perform a series of triaxial creep tests after different DW cycles in this study. Firstly, the experimental data obtained in the creep tests were analyzed. On this basis, a nonlinear damage model accounting for the effects of saturation-dehydration cycles was proposed. Furthermore, the new model was verified and parameter sensitivity analyses of the proposed model were performed. Finally, the three-dimensional central difference expressions necessary for secondary development of the new model was derived in detail.

Preparation of Test Rock Specimens
Argillite was obtained from the reservoir bank slope of Longtan Hydropower Station in China (see Figure 1). In the dam site area, argillite makes up about 30.8% of the strata of the slope, and argillite in the "hydro-fluctuation belt" (Figure 1c) is subject to obvious creep under the influence of reservoir water level fluctuations. The sampling site is shown in Figure 1d. After the drilling of the boreholes, argillite samples (predominantly gray) were obtained in situ. Then, rock samples were packed and transported to the laboratory and processed into cylinders with a diameter of 50 mm and height of 100 mm, according to the GB/T50266-2013 standard for test methods of engineering rock mass [30], as shown in Figure 1e. Table 1 lists the basic mechanical parameters of the argillite samples.  Notes: ρ = Natural density, ρd = Dry density, ρs = Saturated density, W = Natural water content, Ws = Saturated moisture content.

Experimental Procedure
Landslides frequently occur when the reservoir bank with a rock mass is in saturated conditions. Therefore, conventional triaxial tests and triaxial creep tests on rock samples after DW cycles were all performed in a saturated state in this study. In order to achieve saturation, a vacuum saturator was used to saturate the rock samples, according to the procedure mentioned in [31].
In the present study, specimens underwent 0 (saturated state), 1, 5, 10, 15 and 20 DW cycles, respectively. The DW cycles included a saturation process and a dehydration process. The saturation process of rock samples was as follows: firstly, a vacuum air pump was used to seal and saturate the rock samples, then a negative pressure of 0.1 MPa was maintained for 12 h to ensure that the samples were fully saturated. For the dehydration process, according to the standard testing method for engineering rock masses (GB/T50266-2013) [30], the samples were placed in an oven for 12 h at 104 • C to remove water completely and then were cooled naturally. A triaxial compressive strength test (with confining pressure of 4 Mpa) and a triaxial rheological test (with confining pressure of 4 Mpa) were performed after each saturation-dehydration cycle.
For the triaxial creep tests, a rock temperature-stress-seepage coupling triaxial creep meter produced by Top Industrie was used. The test procedure of this device is described in detail in [32]. During the experiment, the confining pressure was maintained at 4 MPa, and axial loads of the rock samples were applied level by level [33,34], which were approximately 30%, 45%, 60%, 75% and 85%, respectively, of the triaxial compressive strength. The loading scheme is shown in Table 2.  Figure 2 shows the results of creep tests after different saturation-dehydration cycles. Deviatoric stress in this figure is the difference in axial load and confining pressure.

Test Results
The strain-time and stress-time curves of argillite samples are presented in Figure 2. It can be seen that: (i) At the same level of deviatoric stress, the strain after stabilization of the creep curve increases nonlinearly as the number of DW cycle increases. For example, after the first level of deviatoric stress is applied, the values of stabilized creep strain after 0, 1, 5, 10, 15 and 20 saturation-dehydration cycles are about 0.0387%, 0.0404%, 0.0413%, 0.0421%, 0.0433% and 0.0440%, respectively. (ii) Under low axial loads, the rheological curves only exhibit two stages (decelerated stage and stable stage). However, three rheological stages-decelerated stage, stable stage and accelerated stage-may appear only when the axial load approaches or reaches the critical failure value of rocks [33,35]. (iii) As shown in Figure 2a-c, when the number of DW cycle is small, e.g., n = 0, 1 or 5, obvious shear failure zones are observed after sample failure, and conjugate shear fractures occur, with approximately X-shaped sample fracture surfaces. However, when n = 10, 15 or 20, the macromorphology of the samples after creep failure is very complex. As shown in Figure 2d-f, several shear fractures and obvious axial and horizontal tensile cracks and flaws were observed. A reasonable explanation of the above phenomenon is that the mechanical properties of the samples are degraded as the number of DW cycles increases, which also results in propagation of numerous micro-fissures inside the samples. In short, the macromorphology of the samples after creep failure becomes more complicated as the number of DW cycles increases. Water 2020, 12, x FOR PEER REVIEW 5 of 17  Figure 3 shows the axial strain/strain rate/time curves after one saturation-dehydration cycle when the fifth level of the axial load is applied. In Figure 3, S t is the moment when the rock samples enter the stable stage from the decelerated stage, and S ε is the strain value at that moment. P t is the moment when the rock samples enter the accelerated stage from the stable stage, and P ε is the strain value at the dividing point between the stable and accelerated stages. F t is the moment when failure of the rock samples occurs, and F ε is the strain value at that moment. At the decelerated stage, the strain continuously accumulates and the accumulation gradually increases, but the strain increase rate reduces continuously, i.e.,   Figure 3 shows the axial strain/strain rate/time curves after one saturation-dehydration cycle when the fifth level of the axial load is applied. In Figure 3, t S is the moment when the rock samples enter the stable stage from the decelerated stage, and ε S is the strain value at that moment. t P is the moment when the rock samples enter the accelerated stage from the stable stage, and ε P is the strain value at the dividing point between the stable and accelerated stages. t F is the moment when failure of the rock samples occurs, and ε F is the strain value at that moment. At the decelerated stage, the strain continuously accumulates and the accumulation gradually increases, but the strain increase rate reduces continuously, i.e., ε > 0, . ε > 0, .. ε < 0. At the stable stage, the strain continuously accumulates over time, but the strain increase rate is constant, i.e., ε > 0, . ε > 0, .. ε = 0. At the accelerated stage, the strain develops rapidly, and the strain increase rate shows a trend of increase, i.e., ε > 0, From the moment of t P , micro-fissures and pores in the rock samples propagate, converge and become connected. And at the moment t F , complete failure of the rock samples occurs.
Water 2020, 12, x FOR PEER REVIEW 6 of 17 accumulates over time, but the strain increase rate is constant, i.e., At the accelerated stage, the strain develops rapidly, and the strain increase rate shows a trend of increase, i.e.,

Axial strain
Strain rate

Establishment of Creep Model
The Burgers model, one of the most effective constitutive models, is widely adopted to describe the creep behavior of rocks [29,36]. This model is capable of accurately describing the creep law of rocks above long-term strength [37,38]. However, the elements in the Burgers model are linear, so it is not capable to reproduce the nonlinear viscoelastic-plastic creep properties of rocks reasonably [39,40]. It is noteworthy that the irreversible viscoplasticity is considered a critical factor resulting in creep failure of rocks [41,42]. The Burgers model does not have a yield threshold; thus, it cannot describe viscoplastic deformation when the applied stress exceeds the long-term strength [33]. Based on this understanding, we proposed a nonlinear rheological model that can reflect the viscoelastic-plastic properties of rocks was proposed based on the Burgers model, as shown in  As can be seen in Figure 3, the rheological curves of the rock samples at the accelerated stage increase nonlinearly. At this stage, considering the change of parameters of the creep model over time would be an ideal method for describing the creep properties of the rock samples [43,44]. Therefore, the viscosity coefficient 3 η in Figure 4 should decrease over time [34,45]. The viscosity coefficients in Figure 4 are expressed in section. At stages other than the accelerated stage, the

Establishment of Creep Model
The Burgers model, one of the most effective constitutive models, is widely adopted to describe the creep behavior of rocks [29,36]. This model is capable of accurately describing the creep law of rocks above long-term strength [37,38]. However, the elements in the Burgers model are linear, so it is not capable to reproduce the nonlinear viscoelastic-plastic creep properties of rocks reasonably [39,40]. It is noteworthy that the irreversible viscoplasticity is considered a critical factor resulting in creep failure of rocks [41,42]. The Burgers model does not have a yield threshold; thus, it cannot describe viscoplastic deformation when the applied stress exceeds the long-term strength [33]. Based on this understanding, we proposed a nonlinear rheological model that can reflect the viscoelastic-plastic properties of rocks was proposed based on the Burgers model, as shown in

Axial strain
Strain rate

Establishment of Creep Model
The Burgers model, one of the most effective constitutive models, is widely adopted to describe the creep behavior of rocks [29,36]. This model is capable of accurately describing the creep law of rocks above long-term strength [37,38]. However, the elements in the Burgers model are linear, so it is not capable to reproduce the nonlinear viscoelastic-plastic creep properties of rocks reasonably [39,40]. It is noteworthy that the irreversible viscoplasticity is considered a critical factor resulting in creep failure of rocks [41,42]. The Burgers model does not have a yield threshold; thus, it cannot describe viscoplastic deformation when the applied stress exceeds the long-term strength [33]. Based on this understanding, we proposed a nonlinear rheological model that can reflect the viscoelastic-plastic properties of rocks was proposed based on the Burgers model, as shown in  As can be seen in Figure 3, the rheological curves of the rock samples at the accelerated stage increase nonlinearly. At this stage, considering the change of parameters of the creep model over time would be an ideal method for describing the creep properties of the rock samples [43,44].
Therefore, the viscosity coefficient 3 η in Figure 4 should decrease over time [34,45]. The viscosity coefficients in Figure 4 are expressed in section. At stages other than the accelerated stage, the . Nonlinear Burgers viscoelastic-plastic (NBVP) rheological model. E 1 is the instantaneous modulus of elasticity; E 2 is the viscoelastic modulus; η 1 is the viscoelastic coefficient; η 2 and are the viscosity coefficients; σ S is the rock yield strength obtained from laboratory triaxial tests [43]; ε 1 , ε 2 , ε 3 and ε 4 are the strains corresponding to the rheological bodies shown in the Figure. As can be seen in Figure 3, the rheological curves of the rock samples at the accelerated stage increase nonlinearly. At this stage, considering the change of parameters of the creep model over time would be an ideal method for describing the creep properties of the rock samples [43,44]. Therefore, the viscosity coefficient η 3 in Figure 4 should decrease over time [34,45]. The viscosity coefficients in Figure 4 are expressed in section. At stages other than the accelerated stage, the viscosity coefficient η 3 remains constant. At the accelerated stage, the viscosity coefficient η 3 is subject to nonlinear processing, and the constitutive equation of the fourth part in Figure 4 is where t 0 is the unit reference time (the value of which is one), t 0 has the effect of maintaining dimensional uniformity of the time unit [46]; m is a parameter to be determined through rheological experiment; and <> is the conditional judgment switching function, in the equation:

One-Dimensional Creep Equations of the NBVP Model
One-dimensional creep equations of the NBVP model, based on the Laplace transform and the inverse Laplace transform [33], are as follows: (1) When σ 0 < σ S , (2) When σ 0 ≥ σ S , t ≤ t P ,

Three-Dimensional Creep Equations of the NBVP Model
The triaxial creep stress state is σ 2 = σ 3 which remains constant, and the confining pressure remains constant after each level of load σ 1 is applied. The rheological equation of the NBVP model under triaxial creep test is shown as follows: where K is the bulk modulus of argillite; G 1 and G 2 are the shear modulus corresponding to E 1 and E 2 in the three-dimensional stress state [47].

Analysis of the Damage Law under the Effect of Saturation-Dehydration Cycles
Defects such as holes, micro-cracks and micro-fissures in rocks propagate inevitably under the effect of external factors such as saturation-dehydration cycles [25]. The damage of rocks also means the development of micro-defects inside rocks [17].
Based on the theory of damage mechanics [48,49], the damage factor D was introduced. It is assumed that (i) the damage factor of the rock samples gradually increases as the number of saturation-dehydration cycles increases, and its value is less than one, and (ii) the change in the damage factor with the number of saturation-dehydration cycle is continuous. Then, where D 1 (n) is the damage factor of the instantaneous shear modulus of elasticity; G 1 (0) is the instantaneous shear modulus of elasticity when the number of saturation-dehydration cycles is zero. G 1 (n) is the instantaneous shear modulus of elasticity when the number of saturation-dehydration cycles is n; when n is zero, the rock samples have no damage, D 1 (n) = 0, and as the number of saturation-dehydration cycles increases, D 1 (n) increases but is less than one. Since the level-by-level loading method is used for triaxial creep tests in this study and the creep mechanical properties of the rock samples vary with the deviatoric stress [6], creep parameters corresponding to each deviatoric stress are also different. To determine the final creep parameters of the rock samples, the mean value method is widely adopted in previous research [26,50]. Table 3 lists the mean values of creep parameters of the NBVP model for argillite under different levels of deviatoric stress after different DW cycles, which were obtained before the rock samples reach the yield strength. Table 3. Values of creep parameters of argillite after different saturation-dehydration cycles. Note: G 1 , G 2 , η 1 and η 2 are respectively the mean values of the creep parameters G 2 , G 2 , η 1 and η 2 of the NBVP model; their damage factors are D 1 (n), D 2 (n), D 3 (n) and D 4 (n), respectively.
According to Table 3, the mean values of G 1 , G 2 , η 1 and η 2 of the creep parameters of argillite decrease as the number of saturation-dehydration cycles n increases. Referring to Equation (7), the damage factor of each creep parameter can be determined by fitting the data in Table 3 using the least square method. The equation expressing the relationship between the damage factor of each creep parameter and the number of saturation-dehydration cycles n is written as

A Nonlinear Creep Damage Model for Rocks Considering the Effect of Saturation-Dehydration Cycles
Creep damage of rocks is mainly due to the change of creep parameters, such as the shear modulus of elasticity and the viscosity coefficient [34,51]. Therefore, the introduction of a variable parameter (damage factor) into the creep constitutive model to express the damage evolution process provides a new method for studying the creep damage properties of rocks [50,52].
Based on the theory of damage mechanics, Equation (8) is introduced into the NBVP rheological model to establish a damage nonlinear Burgers viscoelastic-plastic (DNBVP) model, which is capable of considering the effect of saturation-dehydration cycles (see Figure 5). modulus of elasticity and the viscosity coefficient [34,51]. Therefore, the introduction of a variable parameter (damage factor) into the creep constitutive model to express the damage evolution process provides a new method for studying the creep damage properties of rocks [50,52].
Based on the theory of damage mechanics, Equation (8) is introduced into the NBVP rheological model to establish a damage nonlinear Burgers viscoelastic-plastic (DNBVP) model, which is capable of considering the effect of saturation-dehydration cycles (see Figure 5). (1) When (2) When Combining Equations (6) and (8), the rheological equations of the DNBVP model under triaxial creep tests are given as follows: (2) When σ 1 − σ 3 ≥ σ S and t ≤ t P , (3) When σ 1 − σ 3 ≥ σ S and t > t P , According to the stress state in rheological tests, relevant parameters of the DNBVP model can be obtained by combining Equations (9)-(11) and using the Origin fitting tool.

Discussion
In this section, the DNBVP model verification and parameter sensitivity analysis were conducted. Moreover, the equations necessary for secondary development of this model in FLAC3D were deduced, and the limitations of the DNBVP model were discussed.

Verification of the DNBVP Model
The rheological curves in Figure 2 were processed using the Boltzmann superposition principle [53,54], and parameters of the DNBVP rheological model were identified according to the creep curves. The specific steps are shown as follows: (1) When the deviatoric stress is less than the rock yield strength, the model has four rheological parameters, namely, G 1 , G 2 , η 1 and η 2 . First, according to the rheological test data, let t = 0, and the instantaneous shear modulus G 1 can be determined according to Equation (9). Then, calculate the other three rheological parameters, G 2 , η 1 and η 2 , by fitting the rheological test data. (2) When the deviatoric stress is greater than the rock yield strength, and t ≤ t P , there are five corresponding rheological parameters, namely, η 1 , η 2 , η 3 , G 1 and G 2 , which can be calculated by Equation (10) and the fitting of test curves. (3) When the deviatoric stress is greater than the rock yield strength, and t > t P , there are six corresponding rheological parameters, namely, η 1 , η 2 , η 3 , G 1 , G 2 and m. As can be seen in Equation (11), the six rheological parameters should be determined in order to solve the equation.
The rheological curves are fitted according to Step (2) when the t value is below t p , and then the corresponding parameters η 1 , η 2 , G 1 and G 2 can be calculated; parameters η 3 and m can be given by plugging the above parameters into Equation (11).
With the parameters obtained according to three steps above, the resulting data from the creep tests under the fifth level of load are shown in Figure 2-these were used as examples to verify the proposed DNBVP model. Figure 6 illustrate the comparison between the theoretical curves and the test data. The results show that the theoretical curves in Figure 6 can successfully describe the rheological test data of rocks after different saturation-dehydration cycles.
calculate the other three rheological parameters, G2, η1 and η2, by fitting the rheological test data.
(2) When the deviatoric stress is greater than the rock yield strength, and P t t ≤ , there are five corresponding rheological parameters, namely, η1, η2, η3, G1 and G2, which can be calculated by Equation (10) and the fitting of test curves.
(3) When the deviatoric stress is greater than the rock yield strength, and P t t＞ , there are six corresponding rheological parameters, namely, η1, η2, η3, G1, G2 and m. As can be seen in Equation (11), the six rheological parameters should be determined in order to solve the equation. The rheological curves are fitted according to Step (2) when the t value is below tp, and then the corresponding parameters η1, η2, G1 and G2 can be calculated; parameters η3 and m can be given by plugging the above parameters into Equation (11).
With the parameters obtained according to three steps above, the resulting data from the creep tests under the fifth level of load are shown in Figure 2-these were used as examples to verify the proposed DNBVP model. Figure 6 illustrate the comparison between the theoretical curves and the test data. The results show that the theoretical curves in Figure 6 can successfully describe the rheological test data of rocks after different saturation-dehydration cycles. The classic Kelvin-Voigt model [55,56] and Burgers model [29,[36][37][38][39][40] were compared with the DNBVP model by fitting the test data with n = 20 in Figure 6, respectively (see Figure 7). The results show that the Kelvin-Voigt model and the Burgers model can only describe the first two stages of creep curves, whereas the proposed DNBVP model can describe all creep stages of rocks, especially the accelerated creep stage. The correlation coefficients of the fitted curves obtained from the three models are 0.85, 0.90 and 0.99, respectively. Reasons for the above difference can be explained as follows: at the first two creep stages, the classic models had little difference from the proposed model in terms of fitted curves. However, after the rock samples enter the accelerated creep stage, the classic Kelvin-Voigt and Burgers models cannot describe the nonlinear properties of the rock samples because they lack corresponding nonlinear parameters, whereas the DNBVP model can accurately reflect the nonlinear properties of the rock samples at this stage because it introduced a nonlinear viscoplastic body. The classic Kelvin-Voigt model [55,56] and Burgers model [29,[36][37][38][39][40] were compared with the DNBVP model by fitting the test data with n = 20 in Figure 6, respectively (see Figure 7). The results show that the Kelvin-Voigt model and the Burgers model can only describe the first two stages of creep curves, whereas the proposed DNBVP model can describe all creep stages of rocks, especially the accelerated creep stage. The correlation coefficients of the fitted curves obtained from the three models are 0.85, 0.90 and 0.99, respectively. Reasons for the above difference can be explained as follows: at the first two creep stages, the classic models had little difference from the proposed model in terms of fitted curves. However, after the rock samples enter the accelerated creep stage, the classic Kelvin-Voigt and Burgers models cannot describe the nonlinear properties of the rock samples because they lack corresponding nonlinear parameters, whereas the DNBVP model can accurately reflect the nonlinear properties of the rock samples at this stage because it introduced a nonlinear viscoplastic body.

Analysis of Sensitivity of the DNBVP Model
With the same number of saturation-dehydration cycles, the most critical parameter in the DNBVP model that reflects the characteristics of the rock creep curves is the nonlinear parameter m. Curves at the accelerated creep stage of rocks are complex [27,33]. To better understand the effect of

Analysis of Sensitivity of the DNBVP Model
With the same number of saturation-dehydration cycles, the most critical parameter in the DNBVP model that reflects the characteristics of the rock creep curves is the nonlinear parameter m. Curves at the accelerated creep stage of rocks are complex [27,33]. To better understand the effect of the change of parameter m on the accelerated creep stage of rocks, we analyzed its sensitivity as follows.
Taking the test data in Figure 7 as an example, the obtained parameters of the DNBVP model are listed in Table 4 according to the model parameter identification method mentioned in Section 4.1. By keeping the other creep parameters constant in Table 4 and only assigning different values to the nonlinear parameter m, the tendency of change of the creep curve with the variation in m is obtained, as shown in Figure 8.  As can be seen in Figure 8, the value of creep strain and the creep rate strongly depend on the value of m. When m is greater than 0.1, the creep curve exhibits the accelerated creep stage, and the slope of the creep curve increases as the value of m increases, which means the creep deformation rate changes more rapidly.

Implementation of the DNBVP Model in FLAC3D
FLAC3D is an important tool for analyzing geotechnical engineering issues. To facilitate the As can be seen in Figure 8, the value of creep strain and the creep rate strongly depend on the value of m. When m is greater than 0.1, the creep curve exhibits the accelerated creep stage, and the slope of the creep curve increases as the value of m increases, which means the creep deformation rate changes more rapidly. From the above analysis, it can be found that, with the same number of saturation-dehydration cycles, the creep rate and the creep curve pattern of the DNBVP model are largely restricted by the parameter m, and the change of m fully reflects the nonlinear accelerated creep property of rocks. This suggests the flexibility of the nonlinear parameter m when describing the DNBVP model.

Implementation of the DNBVP Model in FLAC3D
FLAC3D is an important tool for analyzing geotechnical engineering issues. To facilitate the practical application of the DNBVP model, FLAC3D was used as the platform in this study to carry out secondary development of the DNBVP model according to the implementation process described in the user manual [57]. The secondary development of the constitutive model using FLAC3D is actually a process in which new stress is obtained according to the stress and strain increment within the last time-step by using the constitutive model equation established [52]. Therefore, it is crucial to derive three-dimensional central difference expressions necessary for secondary development of the constitutive model equation [47]. The three-dimensional central difference expressions of the DNBVP model are derived as follows: It is specified that S ij is the deviatoric stress; e ij is the total deviatoric stress; S N ij and S O ij are the new and old stress deviators within a certain time increment, respectively; and e N ij and e O ij are respectively the new and old strain deviators within a certain time increment. S S is the three-dimensional yield strength.
(1) When (S ij ) 0 < S S , the three-dimensional central difference expression for the new stress deviator is where e 2,O ij is the new strain deviator of the rheological element of the second part in Figure 5. (2) When (S ij ) 0 ≥ S S and t ≤ t P , the three-dimensional central difference expression for the new stress deviator is At this moment, the central difference scheme of spherical stress of the model is where σ N 0 and σ O 0 are the new and old spherical stresses within a certain time increment, respectively. ∆ε vol is the spherical stress increment within the time step ∆t.
(3) When (S ij ) 0 ≥ S s and t > t P , the three-dimensional central difference expression for the new stress deviator is At this moment, the central difference scheme of spherical stress of the model is In conclusion, the stress-strain relationship of the DNBVP model can be expressed by the three-dimensional central difference expressions for the new stress deviator in Equations (12), (15) and (16). A callable dynamic link library (.dll) is generated through VC++ programming [33], and secondary development of the DNBVP rheological model in FLAC3D can be realized.

Limitations of the DNBVP Model
In this study, a series of creep tests were performed on argillite after different saturation-dehydration cycles, and the DNBVP model, describing the creep properties after DW cycles, was proposed as well. However, the new model proposed still has some limitations. For example: (i) When the number of saturation-dehydration cycles is greater than 20, the rationality of the model cannot be verified, since no triaxial creep test is performed under such conditions. (ii) The expression of the model is very complex, in which six parameters need to be determined, and some parameters, including the value m, which are independent of the physical property of the sample, do not have definite physical significance. Thus, the variation of these parameters currently cannot be explained by linking the mineralogy, the clay content or other aspects of rock characteristics. (iii) The rock samples used in laboratory tests undergo saturation-dehydration cycles, whereas for rocks that are not in an unsaturated state and whose lithology is not argillite, the rationality of the model needs further verification. Therefore, more samples will be used in our future studies to discuss the above issues to verify and optimize our model. Moreover, the performance of creep tests on samples of other lithology in fully saturated and unsaturated states is also essential, in order to understand the creep behavior of rocks under reservoir water level fluctuations. Finally, it is also a challenge for us in the future to analyze practical engineering cases according to field monitoring data using FLAC3D and the outcomes of secondary development of the DNBVP model, in order to further verify the proposed model.

Conclusions
Triaxial creep tests were conducted on argillite after saturation-dehydration cycles and the law in creep tests was analyzed. Furthermore, a nonlinear creep damage model for rocks (the DNBVP model) considering the effect of saturation-dehydration cycles was proposed. The following conclusions can be drawn: