Shear Creep Properties and Creep Model of Gravel Sliding Zone: A Case Study of the Zhoujia Landslide in China

Creep behavior of landslide sliding zones is closely related to the long-term stability and safety of landslides. In this paper, shear creep tests are carried out on undisturbed samples of the gravel sliding zone in the Zhoujia landslide. Creep properties, such as creep rate and long-term strength, of the sliding zone are studied. The result shows that the sliding zone has typical time-dependent behavior. The relationship between the steady strain rate and shear stress can be described by an exponential equation. The long-term strengths of the sliding zone under different normal stresses are determined by using the isochronous curve cluster method. A nonlinear viscoelastic-plastic creep model is developed based on the Nishihara model. The model is shown to be suitable for describing the accelerated creep deformation of the sliding zone. The results obtained are of practical significance for understanding the deformations of the Zhoujia landslide.


INTRODUCTION
Landslide is an important geological phenomenon, in which a soil or rock mass on a slope slips along the shear failure surface . As one of the most frequently occurring geological disasters, it often causes heavy casualties, economic losses, and even catastrophic consequences (Froude and Petley, 2018;Lin and Wang, 2018;Yan et al., 2019;Wang H. et al., 2020). On June 18, 1972, 67 people were killed and 20 injured in the landslide caused by heavy rain on Baoshan Road, Mid-levels, Hong Kong, China. On 22 March 2014 (Au, 1998), the Oso landslide in the United States occurred on a large scale with a volume of 8.3 million m 3 , resulting in 43 deaths (Stark et al., 2017). China is one of the most serious areas in the world with landslide hazard. Especially since the 1980s, occurrences of large landslides in southwest China have been increasing year by year (Huang, 2009;Xu et al., 2014).
For most landslides, the entire process from formation to instability is characterized by distinct creep deformation (Oberender and Puzrin, 2016;Sun et al., 2016), which is a critical factor to consider in the stability evaluation and prediction of landslides. The evolution of deformation and stability of creeping landslides are affected by combined effects of internal and external factors (Sun et al., 2017). On the one hand, the deformation and mechanical properties of landslide materials, especially soil in the landslide zones, have a time effect affecting the behavior of landslides. On the other hand, the seepage fields and stress fields of slopes are continuously changing affected by a number of factors, such as rainfall, rise and fall of reservoir water level, and ground load, which may induce landslide sliding. Therefore, studies on the deformation and failure of creeping landslides should be based on the creep properties of the landslide materials. In recent years, there have been increasing efforts aimed at studying creep properties and creep models of sliding zones. Xu (2012) summarized landslide deformation, failure behavior, and deformation-time curves of a large number of landslides, and found that macroscopic deformation and failure behavior of creeping landslides were mainly caused by the flow and rupture of mesoscale particles of rock and soil mass from the perspective of mesomechanics. Ren et al. (2021) analyzed the mechanism of shear failure of sliding zones from the mesostructures of shear planes, based on direct shear tests and high-precision 3 D laser scanning technology. From the results of direct shear tests of soil-rock mixture (SRM), Yu et al. (2021) found that the higher the rock content, the stronger the bond between soil and rock, and the greater the shear strength. Tang et al. (2020) conducted triaxial creep tests on loess under different water contents and pressure conditions; it was found that creep behavior of loess is significant at high moisture levels and that less time-dependent deformation occurs at high confining pressures. Wang L. et al. (2020) investigated the effect of shear rate on shear residual strength of slip zone soils. Although theoretical and experimental studies on homogeneous materials, such as clay and silty clay, have garnered a lot of attentions in recent years, there has been a lack of research on gravel soil. Significant influence of large-size gravel on characteristics of deformation and strength of soil was demonstrated in a study by Cheng et al. (Zhanlin et al., 2007). Therefore, large-scale creep tests on gravel sliding zones are necessary for better understanding of the deformation evolution and stability of creeping landslides.
Based on the results of large-scale creep shear tests of gravel sliding zones, the creep properties of the sliding zone of the Zhoujia landslide, which is of great importance for further understanding the deformation and failure mechanism of the landslide, are investigated in this paper. A nonlinear viscoelasticplastic creep (NVPC) model is proposed. Compared with the existing traditional models, this model is more suitable for describing accelerated creep deformation under shear load. The model can be applied in stability analysis of the Zhoujia landslide, providing an important reference for the risk control of the area and the safety of the reservoir.

BACKGROUND OF THE ZHOUJIA LANDSLIDE
The Zhoujia landslide is located on the right bank of Yalong River in Sichuan Province of southwest China, which is about 9.2-11.0 km away from the dam site of Kala Hydropower Station. The landslide is M-shaped and is divided into three zones: Zone A, Zone B1, and Zone B2, as shown in Figure 1 ). The inclinometer profile of borehole INzj1-1 shows that significant horizontal displacement has occurred at depths of 17 and 52 m. This implies a double shear-surface structure of the landslide.

Test Materials
The sampling location is in exploration adit TD37 (elevation 2,191.8 m) in Zone B1 of the Zhoujia landslide, 63.5-65.0 m away from the entrance of the exploration adit. The exposed sliding zone ( Figure 3A) is composed of mainly gray-yellow and graybrown gravel soil, and its thickness is about 0.3-1.0 m. The content of gravel is 50-60%, and the structure of sliding zone is dense. Experiments on basic physical properties were conducted on the sliding zone and the results indicate that the average density is 2.3 g/cm 3 and the water content is 10.3%; hence, the average dry density of sliding zone is 2.1 g/cm 3 . The distribution of particle size is shown in Figure 4. The undisturbed sample of the sliding zone were prepared in the field in Figure 3B. Considering the size limitation of the direct shear box (150 × 150 × 150 mm), the sampling location was carefully selected so that the maximum gravel size did not exceed 75 mm. If the size of gravel was larger than 75 mm, they were replaced with finegrained soil in the process of sample preparation. The samples were sealed with a thin plastic film to preserve the in situ moisture content.

Test Equipment and Procedure
Creep experiments were carried out on a CSS-3940YJ shear rheological testing machine developed by Changchun Testing Machine Research Institute, China. All tests are conducted in the laboratory condition with constant temperature (24 ± 0.5°C) and humidity. The samples are consolidated for at least 24 h until the vertical settlement is less than 0.05 mm/h. Direct shear tests are first carried out on the sliding zone. The normal stresses, which are the converted overburden earth pressures at each sampling location, are set at 851, 1,100, and 1,380 kPa. The normal stress is kept constant  until completion of the test, and the shear rate is 0.02 mm/min. Two or three groups of tests are repeated under each normal stress.
According to the results of direct shear tests, the shear creep test procedure is formulated. The samples are consolidated for at least 24 h until the vertical settlement tended to be constant. Shear creep tests are then carried out. Because the number of prepared samples is limited, the multi-loading method is adopted. Each shear loading increment is 1/11 to 1/5 of the peak strength of direct shear under the corresponding normal stress. At least five levels of shear load are applied, and the loading rate is 0.1 kN/ min. When the shear deformation is less than 5 × 10 −4 mm/day, the next level of the shear load is applied.

Direct Shear Test
A series of direct shear tests under different normal stresses are carried out first to study the shear mechanical properties of the sliding zone and to provide data support for estimating shear stress levels of shear creep tests. The stress-strain curves of the sliding zone under various normal stresses are presented in Figure 5B. It is obvious that the stress-strain curves have no distinct peaks and generally exhibit ideal elastic-plastic characteristics. With the increase of normal stress, the shear strength of the sliding zone increases gradually. Deformations of the samples can be divided into three stages. During the initial loading stage, the slope of the stressstrain curve is relatively small, which means that shear stress increases slowly. This can be attributed to the closure of pores and cracks between gravel and soil. With the increase of shear stress, the slope of the stress-strain curve increases rapidly. When the stress is close to the peak value, the soil in the sample reaches yield, and the slope of the stress-strain curve decreases gradually. Lateral expansion occurs in the sample, with soil or gravel extruding. When the shear stress exceeds the peak shear strength, the soilrock structure in the sample breaks down and recombines, and the stress-strain curve shows strain-hardening characteristics. In Figure 5A, the Mohr-Coulomb criterion is used to fit the relationship between peak shear strength and residual shear strength under different normal stresses. It can be seen that there are good linear relationships between the normal stress and the peak shear strength and between the normal stress and residual shear strength. The internal angles of friction are 33.42°and 33.02°and the corresponding cohesions are 91.49 and 80.64 kPa, respectively.

Shear Creep Test
Similar to most creeping landslides in reservoir areas, the Zhoujia landslide is moving slowly with obvious creep properties. Because creep properties of the sliding zone are critically important for the long-term safety and stability of Kala Hydropower Station during both construction and operation, multi-loading shear creep tests were carried out under different normal stresses. The stress conditions are listed in Table 1.

Shear Creep Properties
The classical shear creep properties of all samples during the loading processes are shown in Figure 6. The deformation curves are divided into three creep stages: transient creep, steady-state creep, and accelerated creep (Jia et al., 2018;Han et al., 2021). It is obvious that the samples exhibit characteristics of transient deformation at the beginning of each loading process, with large shear displacement occurring in a short period of time. Transient deformation is a large part of the total deformation. As time goes by, the shear strain rate gradually decays to a constant, and the sample enters a steady creep FIGURE 4 | Distribution of particle size of sliding zone. stage. Finally, the sample breaks down, and the accelerated creep stage takes place. It is found that shear strain increases step-by-step with step-wise increase of shear stress under the same normal stress. However, there are special cases, such as stage Ⅰ in Figure 6B and stage Ⅱ in Figure 6C. This is due to the friction of one gravel with another gravel on the failure surface in the loading process, resulting in occlusal effect, which prevents the deformation of samples. In Figure 7, white friction scratches appear on the failure surface. With the same shear stress increment, the creep deformation of the sample under normal stress of 1,380 kPa is greater than that under normal stress of 1,100 kPa, which is greater than that under normal stress of 851 kPa.

Shear Creep Rate
The variations of shear creep deformation and shear strain rate with time under shear stresses of 194.4 and 666.0 kPa are presented in Figure 8A,B, respectively, for sample HD-3. When the shear stress is less than the value of the last stage ( Figure 8A), the shear strain rate versus time curve exhibits an L-shaped characteristic: the shear creep rate is very large initially, then it decreases rapidly and approaches a small constant value, usually less than 0.06 × 10 −2 mm/h. At the final shear stress level ( Figure 8B), the shear strain rate versus time curve exhibits a U-shaped characteristic: the shear creep rate decreases rapidly to a small constant value, and then the creep rate increases abruptly until the test is stopped.
It is found that the steady-state creep rate increases with the shear stress level. The exponential function can be used to fit the experimental data ( Figure 9) as follows (Jia et al., 2018) v c α e βτc where α and β are constant parameters, τ c is the shear stress, and v c is the corresponding steady-state creep rate. The values of parameter α are 1.845 × 10 −4 , 9.207 × 10 −4 , and 2.501 × 10 −4 , the values of β are 0.0058, 0.0025, and 0.0042, and the corresponding R 2 values are 0.951, 0.966, and 0.998, respectively, under normal stresses of 851, 1,100, 1,380 kPa. The high values of R 2 indicate that the curves fit the experimental results well.

Long-Term Shear Strength
It is of great significance to obtain the long-term shear strengths of sliding zones for the long-term safety and stability of landslides (Liu et al., 2004;Shen et al., 2012). In this study, the isochronous curve cluster method (Li et al., 2010) is applied to determine the long-term strength of the sliding zone of the Zhoujia landslide. Shear strains corresponding to different shear stresses at a certain time interval are selected, and the shear stress-strain isochronous curves are plotted, as shown in Figure 10 for sample HD-3. It is obvious that these isochronous clusters of curves have great similarities. The shear stress corresponding to the turning point of clusters of curves from linear to nonlinear is defined as the long-term strength. Values of the long-term strength of the sliding zone are listed in Table 2. It can be seen that the long-term strengths of the three samples are lower than their instantaneous strength, and the reduction range is 13-18%.

SHEAR CREEP MODEL AND VERIFICATION 5.1 Describing Creep Behavior Using Fractional Calculus
The fractional calculus is a branch of calculus that studies differential and integral operators of any order. It is a mathematical tool for solving problems of physical and mechanical modeling effectively. Based on Riemann-Liouvelle's theory (Koeller, 1984;Adolfsson et al., 2005), Scott-Blair described the constitutive equation of rock or soil as (Blair, 1944) where τ(t) is the shear stress on rock or soil, γ (t) is the corresponding shear strain, η is the coefficient of viscosity, t is time, and n is between 0 and 1. When n is 0, the material is an ideal solid; whereas when n is 1, the material is an ideal fluid. Based on Eq. 2, a new sticky pot element, named Able (Zhou et al., 2012), was defined to describe the creep deformations of materials between ideal solids and ideal fluids. When the shear stress τ tends to a constant value, integrating both sides of Eq. 2 using fractional calculus yields where Γ(·) is the gamma function.

Nonlinear Viscoplastic Element
Rock or soil generally exhibits accelerated creep properties under high shear stress, and nonlinear creep elements are used to describe this stage. In this study, a nonlinear viscoplastic model is used to describe the creep behavior as follows (Xu et al., 2006)  where m is creep index, which reflects the rate of shear creep in the accelerated creep stage. τ s is the long-term shear strength of rock or soil, which is generally obtained by tests.

NVPC Model
The Nishihara model is a good and simple model to describe the creep properties fairly comprehensively. However, it can only be   Frontiers in Earth Science | www.frontiersin.org February 2022 | Volume 10 | Article 838183 7 used to describe the creep curve prior to the accelerated creep stage. As shown in Figure 11A, it is composed of a Hookean solid (H), a viscoelastic body (N/H), and a viscoplastic body (N/St.V) connected in series. G 1 represents the instantaneous shear modulus, G 2 represents the viscoelastic shear modulus, η 2 , and η 3 represent the viscoelastic coefficients, τ s is the long-term shear strength of the sliding zone, and γ e , γ ve , and γ vp are the shear displacements corresponding to their respective creep bodies.
Under the assumption that τ is the total shear stress and γ is the total shear displacement, the constitutive equation of the Nishihara model can be expressed as follows (Jiang et al., 2013;Wang et al., 2019) However, the Nishihara model cannot describe the nonlinear accelerated creep properties of sliding zones. In order to describe the creep properties of sliding zones more accurately, based on the theory of fractional calculus as expressed by Eq. 3, the Newtonian element in the Nishihara model is replaced by Able element, and the viscoplastic body model (N/St.V) is replaced by the nonlinear viscoplastic element model described by Eq. 4. Finally, the long-term strength of the sliding zone obtained by shear creep tests is measured as the threshold of nonlinear accelerated creep. Considering the strains of three parts of the model, NVPC model is established in Figure 11B, which can reflect the characteristics of three-stage creep of a sliding zone as where n is the value of fractional order, k is a non-negative integer, η NV is the nonlinear viscoplastic coefficients, τ s is the long-term shear strength threshold of the sliding zone that  depends on the long-term cohesion and the internal angle of friction, and m is creep index.

Identification of the Parameters
In order to verify the NVPC model proposed in this study, the creep parameters are identified based on the shear creep test data of sliding zone sample HD-3, as shown in Table 3. These parameters are substituted into the model to obtain theoretical results and compared with the experimental test values, as shown in Figure 12A. The result obtained using the Nishihara model is also plotted for comparison in Figure 12B. It is obvious that Nishihara model has errors in the accelerated creep stage. Whereas the NVPC model can describe the three shear creep stages, particularly the accelerated creep stage, indicating that the creep model proposed in this study is appropriate.

CONCLUSIONS AND FUTURE WORK
The following conclusions can be drawn from this study: (1) Shear creep tests of three undisturbed samples taken from the Zhoujia landslide are performed under normal stresses of 851, 1,100, and 1,380 kPa, respectively. The shear creep curves of the sliding zone exhibit three creep stages: transient, steady-state, and accelerated creep stages. The tested samples from the sliding zone exhibit nonlinear viscoplastic deformations under the shear stresses. The applied normal stress has a significant influence on these stages. (2) In the transient creep stage, the creep rate attenuates quickly from a large value to a small constant value with the increase of time and enters the steady creep stage, displaying an L-shaped curve of shear strain rate versus time. When entering the accelerated creep stage, the creep rate increases rapidly until the test is stopped, displaying a U-shaped curve of shear strain rate versus time. The relationship between the steady strain rate and shear stress can be described satisfactorily by an exponential equation. The empirical relation is useful for monitoring and forecasting creep deformations of the Zhoujia landslide.
(3) NVPC model is proposed to describe all three creep stages of the sliding zone of the Zhoujia landslide. Based on the longterm strength parameters obtained from the tests, the parameters are identified. The results from the theoretical model are in good agreement with the experimental values.
In this paper, the analytical nonlinear viscoelastic-plastic creep model is established, and numerical verification of the model is not performed, which will be done in future work.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
SC wrote the manuscript. SC and MS carried out relevant experiments. All authors discussed about the contents.

FUNDING
This study was supported by Research Grants (No. 51939004, 52109122) from the National Natural Science Foundation of China.