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 gray-brown 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³ and the water content is 10.3%; hence, the average dry density of sliding zone is 2.1 g/cm³. 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 fine-grained soil in the process of sample preparation. The samples were sealed with a thin plastic film to preserve the in situ moisture content.

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⁻⁴ mm/day, the next level of the shear load is applied.

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 stress-strain 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 soil-rock 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.

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.

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) τ ( t ) = η d n γ ( t ) d t n (2) 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 γ ( t ) = τ η t n Г ( 1 + n ) (3) where Г ( ⋅ ) is the gamma function.

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) γ ( t ) = { 0 ,      τ < τ s τ − τ s η t m ,       τ ≥ τ s (4) 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.

The Nishihara model is a good and simple model to describe the creep properties fairly comprehensively. However, it can only be 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 ₁ represents the instantaneous shear modulus, G ₂ 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 ,     γ v e , and γ v p 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) γ ( t ) = { τ G 1 + τ G 2 ( 1 − e − G 2 η 2 t ) ,      τ < τ s τ G 1 + τ G 2 ( 1 − e − G 2 η 2 t ) + ( τ − τ s η 3 ) t ,    τ ≥ τ s (5)

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 γ ( t ) = { τ G 1 + τ η 2 , β ∑ k = 0 ∞ ( − G 2 η 2 , β ) k t n ( 1 + k ) n ( 1 + k ) Г [ ( 1 + k ) n ] ,       τ < τ s τ G 1 + τ η 2 , β ∑ k = 0 ∞ ( − G 2 η 2 , β ) k t n ( 1 + k ) n ( 1 + k ) Г [ ( 1 + k ) n ] + ( τ − τ s η NV , n ) t m ,       τ ≥ τ s (6) where n is the value of fractional order, k is a non-negative integer, η N V 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.

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.