Loess is extremely sensitive to water and its strength reduces rapidly upon wetting (Xing et al., 2016; Yan et al., 2018). Research has revealed that water sensitivity of loess is manifest by the collapse of the original soil structures and the decrease in the shear strength and stiffness after encountering water (Feda, 1988; Zhou et al., 2014b; Garakani et al., 2015); however, in terms of loess landslides, previous studies (Tu et al., 2009; Jin et al., 2012; Wang et al., 2014; Carey et al., 2017; Li et al., 2020; Wang et al., 2020; Pu et al., 2021) usually evaluated stability of a loess slope by using the general theory of unsaturated soil mechanics, such as the Bishop’s shear strength formula (Bishop et al., 1961) which considers pore-water and pore-air pressure in soils but does not consider the water sensitivity of mechanical parameters of the loess (Xu L. et al., 2018; Xu Y. et al., 2018). This analysis method usually overestimates the slope stability due to neglecting the water sensitivity of loess, therefore, loess slopes should consider not only the effects of rainfall infiltration on the increase of the pore water pressure but also the variations in the mechanical parameters with increasing water content.

To demonstrate such specific hydraulic-mechanical characteristics of loess in slope stability analysis, especially for loess slopes with seepage fissures under the action of rainfall infiltration and traffic load, a new simulation method is proposed, as shown in Figure 3.

First, the hydraulic conductivity curve and soil-water characteristic curve of the loess are experimentally derived, and the unsaturated steady seepage theory is used to inversely calculate the initial water content and initial pore water pressure of the slope based on the pore water pressure and the water table as measured on-site.

Second, the water content increment field of the slope under rainfall is calculated based on the fields of water content before and after infiltration, and then the mechanical parameter fields, such as the distributions of cohesion, angle of internal friction, deformation modulus, and unit weight, are fitted based on the measured hydraulic-mechanical properties of the loess and the distribution of the increment in water content.

Finally, a constitutive model suitable for the loess (presented in Consideration of Transverse Isotropy of Loess ) is adopted to conduct mechanical analysis of the slope based on the pore water pressure and fitted mechanical parameters of the slope. The fissure-deformation development mechanism of the slope is revealed by analysing the calculated results, and the slope stability is evaluated.

It is noted that this method is an equivalent method of mechanical parameters to consider the water sensitivity of loess. The changes in mechanical parameters of the loess after encountering water, such as the decrease in strength and stiffness and the increase in unit weight, are fully considered in the method, but the collapse of the soil structure and the change in volume are not involved.

Soil structures are the inherent determinants of the strength and stiffness of the soils, especially for loess (Jiang et al., 2014; Wen and Yan, 2014; Xu and Coop, 2016). Figure 4A depicts the Q 2 eol loess collected at the study site, which exhibits transverse isotropy due to the layered structures and micro-cracks in the soil. The loess is equivalent to the mechanical model due to this transverse isotropy as shown in Figure 4B. In the equivalent model, x ' , y ' , and z ' are the local coordinates, where the x ' and y ' -axes are in the plane of the layer and the z ' -axis is normal to it. The strength and stiffness of the soil are isotropic in the layer plane (in either the x ' or y ' -direction) (Jiang et al., 2014; Luo et al., 2018).

The Anisotropic-Elasticity Ubiquitous-Joint Model in Flac^3D (Itasca Consulting Group, 2018) was modified to simulate the transversely isotropic strength and stiffness of the loess. The modified model integrates the logic of the transversely elastic model (which is used to simulate the transverse isotropy of the deformation of the loess) with that of the ubiquitous-joint model (which is employed to simulate the yield of the loess in a specific direction, i.e., the anisotropic strength of the loess). In calculation using the model, the elastic stress and strain of the elements are first calculated by the transversely elastic model (i.e., elastically trial calculation). Then, the plastic elements are determined based on the elastic stress and the yield function, and accordingly, the stress is corrected based on the distribution of the plastic elements. Followed by the stress correction, the strain of the plastic elements is corrected using the plastic flow rules to determine the incremental plastic strains. The main modification of this model is the stress and strain correction in the x ′ , y ′ , and z ′ -directions.

The Mohr-Coulomb strength criterion is employed to judge the yield elements and correct the stress, and then the associated or non-associated plastic flow rule is accordingly used for the strain correction. Let τ i be defined as the magnitude of the tangential traction component in i ′ direction ( i ,   j , k ∈ { x ,   y ,   z } ), then we obtain τ i = σ i ′ j ′ 2 + σ i ′ k ′ 2 (1) where σ i ′ j ′ is stress component of the j ′ -direction in the plane normal to the i ′ -direction. Let γ i be the strain variable associated with τ i , then γ i = ε i ′ j ′ 2 + ε i ′ k ′ 2 (2) where ε i ′ j ′ is the elastic strain component of the j ′ -direction in the plane normal to the i ′ -direction. For shear failure in the i ′ -direction, the failure criterion f i s has the following form f i s = τ + σ i ′ i ′ t a n ϕ i − c i (3) where ϕ i and c i are the angle of internal friction and cohesion in the i ′ -direction, respectively. For tension failure in the i ′ -direction, the failure criterion f i t has the following form f i t = σ i ′ i ′ − σ i t (4) where σ i t denotes the tensile strength in the i ′ -direction. The potential functions, g s and g t , are used to define the shear plastic flow and the tensile plastic flow, respectively. The function g s corresponding to a non-associated law has the following form g i s = τ + σ i ′ i ′ t a n ψ i (5) where ψ i is the angle of dilation in the i ′ − direction. The function g t corresponding to the associated law is g i t = σ i ′ i ′ (6) It is noteworthy that the main modification of the Anisotropic-Elasticity Ubiquitous-Joint Model involves the yield judgment and stress correction in the x ′ , y ′ , and z ′ -directions based on the structural characteristics of the loess. The transversely elastic calculation method and the forms of the failure criteria and potential functions are still consistent with the existing model. Similar modification has been made for jointed rock mass and the rationality of such modification has been verified (Wang and Huang, 2014; Leng et al., 2021), so, the verification of the model modification is not presented here due to the existing similar verification having been published previously.

The undisturbed Q₂ loess collected from the Yangpoyao loess slope is brown yellow, hard, and dense, with layered structures and micro-cracks in the soil (Figure 4A). The physical properties were listed in Table 1. The main particle components of the loess are silt particles: the proportion of the loess particles with a diameter larger than 0.05 mm is 12.9%, the proportion of those from 0.005 to 0.05 mm is 71.5%, and the proportion of those smaller than 0.005 mm is 15.6%. The liquid limit and plastic limit of the soil are 31.4 and 16.7%, respectively. The unit weight of the loess varies with the water content. In the natural state, the lowest water content of the loess is about 14%, and the corresponding unit weight is 1.52 g/cm³; the saturated water content is about 30%, and the corresponding unit weight is 1.73 g/cm³. The relationship between the water content and the bulk unit weight is roughly linear: γ = 1.   33 + 1.   33 w (7) where γ and w are the bulk unit weight and water content of the loess, respectively.

To measure the hydraulic characteristics of the loess, the tests of soil-water characteristic curve were conducted using a volumetric pressure plate apparatus under pressures of 10, 50, 100, 200, 300, 500, 750, 850, and 1,000 kPa. The test data are shown in Figure 5A. Generally, the soil-water characteristic curve can be described by some models, such as van Genuchten function (van Genuchten, 1980), Fredlund-Xing function (Fredlund and Xing, 1994), and modified Kovacs function (Aubertin et al., 2003). In the present study, as the main particle components of the loess are silt particles, the Fredlund-Xing function (Eq. 8) is used to fit the soil-water characteristic curve, and the fitted function is displayed in Eq. 9. θ w = C ( φ ) θ S a t { l n [ ( e + φ a ) n ] } m (8) θ w = 0.   48 { l n [ ( e + φ 10.   65 ) 15.18 ] } 0.209 (9) where C ( φ ) is a correlation function, θ w is the volumetric water content, θ S a t denotes the saturated volumetric water content, φ is matrix suction (equal to the negative pore water pressure), and parameters a, n, and m represent the fitting parameters that control the shape of the volumetric water content function. θ w can be converted to water content w by equation w = θ w / ( 1.   5 + θ w ) × 100 % where 1.5 is the bulk unit weight of the dry loess.

The hydraulic conductivity curves can be fitted based on the existing models, such as the van Genuchten model (van Genuchten, 1980) and Fredlund-Xing-Huang model (Fredlund et al., 1994). In the study, Fredlund-Xing-Huang model (Eq. 10) is used to fit the hydraulic conductivity function based on the water-soil characteristic curve function and the soil properties, such as the saturated hydraulic conductivity, and saturated and residual water contents. The fitted hydraulic conductivity curve obtained by using the estimation techniques of Seep/W (part of the Geostudio suite) is demonstrated in Figure 5B. K w ( θ w ) = K s a t ∫ θ r e s θ w θ w − x φ 2 ( x ) d x / ∫ θ r e s θ s a t θ s a t − x φ 2 ( x ) d x (10) where K s a t is the saturated hydraulic conductivity, and x represents a dummy variable of integration representing the water content.

Hydraulic-mechanical characteristics and their mathematical expressions in the loess are necessary for the seepage-stress coupling simulation of the slope under rainfall and traffic load. Previous studies revealed some hydraulic-mechanical characteristics of loess, including the variation in the matric suction, shear strength and stiffness with the water content (van Genuchten, 1980; Fredlund et al., 1994; Aubertin et al., 2003; Haeri et al., 2014; Wang W. et al., 2018; Xu Y. et al., 2018; Liang et al., 2018). However, these studies rarely propose mathematical expressions of the variation in the mechanical parameters with the water content, so, it is difficult to apply the results to numerical simulation in any direct sense.

These mathematical expressions were further studied by using direct shear tests. The tests were conducted using the direct shear apparatus at a horizontal (shear) speed of 0.8 mm/min under the condition of the vertical stresses of 100, 200, 300, and 400 kPa, respectively. Six groups of the remoulded samples with the water contents of 14 (initial water content), 16, 18, 20, 22, and 24%, respectively, were prepared for the direct shear tests. The relationships between the shear stress and shear deformation of the samples with varying water contents tested under the vertical stress of 100 kPa (those at 200, 300, and 400 kPa had similar patterns and are not given) were revealed, as illustrated in Figure 6A.

The mechanical parameters, including the cohesion (Figure 6B), angle of internal friction (Figure 6C), and deformation modulus (Figure 6D), were obtained based on the results of the direct shear tests. The fitted functions of the variation in the mechanical parameters with the water content increment are shown in Eqs. 10–13, where c is the cohesion (kPa), Δ w is the change in water content (%), ϕ is the angle of internal friction (°), and E m denotes the deformation modulus (MPa). c = 80 − 689 Δ w (11) ϕ = 28 − 92 Δ w (12) E m = 19.   3 − 57.   4 Δ w (13) It is noted that due to the limitations imposed by the prevailing conditions at the study site, it is difficult to obtain undisturbed Q₂ loess samples with different water contents, such that remoulded samples with a 14–24% water content were formed in the laboratory. According to related studies (Wang Y. et al., 2019; Chang et al., 2021), the difference in mechanical parameters between the undisturbed loess and remoulded samples gradually decreases with an increasing water content and is acceptable when the water content is 14–24%. As the remoulded samples lose their structural characteristics compared with the undisturbed loess, it is necessary to convert the mechanical parameters of remoulded loess into those of undisturbed loess. Based on the transformation method proposed in previous studies (Jiang et al., 2014; Zhang L. et al., 2020; Liu et al., 2020), the strength and stiffness parameters of undisturbed loess is about 1.05 times that of remoulded loess in the isotropic plane and 0.9 times that of remoulded loess out of the isotropic plane. In general, this relatively small difference between the remoulded loess and undisturbed loess does not affect our investigation of the deformation development and destabilisation mechanism of fissured loess slopes.

The seepage boundary conditions are set as shown in Figure 7A based on the water table revealed in the site investigation and the rainfall conditions over the study area. In the seepage model, the boundary conditions at the bottom and on the highway are impervious boundaries with zero water flux (WF), those on the left and right sides below the water table are constant pore-water head (PWH) boundaries, those on the left and right sides above the water table are water unit gradient (WUG) boundaries and that on the upper surface is a water-infiltration boundary with a water flux (WF) of 50 mm/d (the seepage lasts for 1 day). To simulate the fissure development process, the depth of the fissure is set to 0, 0.5, 1.0, and 1.5 m, respectively. The boundary condition at the fissure is pore-water pressure head varying with the depth of the fissure. The saturated/unsaturated seepage theory is used in seepage calculations, and the material parameters, including the soil-water characteristic curve and hydraulic conductivity curve, are assigned to the model according to the fitted results (Unsaturated Hydraulic Characteristics of the Loess).

The mechanical boundary conditions are set as shown in Figure 7B. In the mechanical model, the horizontal displacements of the left and right boundaries and the vertical displacement of the bottom boundary are fixed. The boundary conditions of the highway (with six lanes) and temporary road (with two lanes) are stress boundaries with a normal stress thereon of 10.5 kN/m (per unit width out of the plane of analysis) per lane. Other boundary conditions of the model, including the slope surface and the fissure with varying depth, are free (stress boundary with zero stress). Moreover, the pore water pressure distribution is applied to the model based on the calculation results of the seepage. The modified constitutive model in Consideration of Transverse Isotropy of Loess is used for the mechanical calculation, and the material parameters, including the cohesion, friction, deformation modulus, and bulk unit weight of the loess, are fitted based on the seepage results and assigned to the model, as shown in Settings for the Seepage-Coupled Stability Simulation.

Figures 8A–D show the results of the seepage calculation on the slope with a 1.5-meter-deep fissure. It can be seen from the figures that the water content and pore-water pressure on the slope surface, especially at the fissure, have increased significantly after infiltration (Jin et al., 2012; Zhou et al., 2014a; Chang et al., 2021). The seepage results are used for the mechanical analysis of the slope by the method proposed in Consideration of Water Sensitivity of Loess and the steps are as described below.

Firstly, the pore-water pressure field after rainfall infiltration (Figure 8B) is directly applied to the mechanical calculation model of the slope.

Secondly, the water content increment can be obtained by subtracting the water content field before (Figure 8C) from that after (Figure 8D) water infiltration. The mechanical parameter fields (Figure 8E), including those pertaining to the cohesion, angle of internal friction, and deformation modulus, can be calculated via Eqs. 11–13 based on the water content increment. A similar method is used to fit the bulk unit weight field via Eq. 7 based on the water content distribution after water infiltration (Figure 8D).

Finally, the fields pertaining to the cohesion, angle of internal friction, deformation modulus, unit weight and pore-water pressure are assigned to the mechanical model for slope stability analysis.

Previous studies (Ma et al., 2016; Tan et al., 2017) divided the fissure-deformation development process of soil slopes under internal and external actions into the trailing-edge tensile fissure stage (initial deformation stage), side-edge shear fissure stage (uniform deformation stage), leading-edge tensile fissure stage (accelerated deformation stage), and fissure penetration stage (sharp deformation stage); however, in terms of loess slopes with seepage fissures, the fissure-deformation development process is different.

Figure 9 illustrates the schematic diagrams of the shear strength field, the distribution of the plastic zones and the displacement counter at different fissure development stages of the slope under rainfall and traffic load. The fissure-deformation development process can be back-analysed based on the data shown in Figure 9.

Initial stage (Figure 9A): there are no fissures on the slope at this stage. The slope, with a depth of soil of only about 300 mm, reaches a saturated state after rainfall infiltration (Wu et al., 2017; Hou et al., 2018), which causes a minor reduction in the strength and stiffness of the slope. The deformation of the slope under the traffic load is small, and the slope remains stable.

Micro-crack penetration stage (Figure 9B): loess contains many micro-cracks inherently, but these micro-cracks do not penetrate in a natural state (Jiang et al., 2014; Wen and Yan, 2014; Xu and Coop, 2016); however, under traffic load on the trailing-edge of the slope, the middle part of the slope surface is subjected to tensile stress, which contributes to the expansion and penetration of the micro-cracks in the soil. Rainfall infiltration reduces the strength of the soil within a certain depth of the slope surface, promoting micro-crack penetration.

Fissure propagation stage (Figure 9C): after rainfall infiltrates the soil along the penetrating micro-cracks, the strength and stiffness of the slope are significantly reduced. The reduction of the stiffness causes an increase in the deformation of the slope, and that of the strength leads to an increase in the shear and tensile failure on the slope surface. With the development of the deformation and failure zone of the slope, the fissure developed from the penetrating micro-cracks gradually propagates downwards, with further increases in the depth, width, and length of the fissure.

Sliding surface formation stage (Figure 9D): a sliding surface is initially formed along the fissure under the continuous action of rainfall and traffic load. The fissure becomes an express channel for water infiltration under precipitation (Jin et al., 2012; Zhou et al., 2014a; Chang et al., 2021). The seepage field near the fissure is greatly affected so that the strength and stiffness of the soil are correspondingly reduced, resulting in a failure zone expanding around the fissure. Meanwhile, with propagation of the fissure, the constraint of the slope near the fissure is gradually reduced, and as a result, the slope becomes more prone to deforming under traffic load. When the area of the failure zone and the magnitude of the deformation reach a certain level, the slope sliding surface gradually forms.

The study of deformation development process of the slope after stabilisation treatment is key to early warning and hazard assessment of landslides (Tu et al., 2009; Springman et al., 2013; Zhang et al., 2019). Previous studies found that the slope deformation-sliding development process is mainly divided into four stages: initial deformation, uniform deformation, accelerated deformation, and sharp deformation (sliding) (Tu et al., 2009; Xu L. et al., 2018; Chen et al., 2018; Sun et al., 2019; Li et al., 2020). However, the deformation-stabilisation development process for slopes with fissures after stabilisation treatment under rainfall and traffic load is much less frequently explored.

Figure 10 illustrates the monitored horizontal displacement of the Yangpoyao slope and the rainfall over the study site. The displacement was monitored using the deep borehole inclinometer labelled in Figure 10A, and the rainfall was monitored by rain gauges at the study site. The displacement was recorded from January 1, 2017 to June 1, 2019 (the data recorded thereafter are not presented since they are generally unchanged compared with those recorded within the first few months of 2019). During the monitoring period, two stabilisation measures were taken to prevent the fissure from further propagating. One of the treatments was that the traffic load on the slope trailing-edge was removed from February 2017, and the other was that a lawn was planted on the slope surface in the spring (especially April) of 2017, which restrained infiltration into the slope.

From the deep horizontal displacements in the middle of the slope (Figure 10A), the displacements are large within the depth range from 0 to 6 m and are small at the depth below 9 m. The displacements change greatly within the depth range from 6 to 9 m, indicating that the depth of the potential sliding surface is mainly within that range. On the other hand, the depth of the potential sliding surface shown in the figure generally matches that predicted in Figures 9C,D, which indicates that the simulation in is reasonable.

From the monitored horizontal displacements of the slope surface (Figure 10B), the slope deformation development process can be divided into three stages: rapid deformation (from January 2017 to April 2017), slow deformation (from May 2017 to November 2017), and stable deformation (from December 2017). There is a significant boundary between the rapid deformation stage and the slow deformation stage, but the boundary between the slow deformation stage and the stable deformation stage is not obvious. The characteristics of each stage are described below.

In the rapid deformation stage, the total horizontal deformation of the slope increased from 0 to 400 mm in the 4 months from January 2017 to April 2017. The monthly deformation increment was 80 mm in January 2017 and increased to 150 mm in February 2017. To suppress the growing deformation and improve the stability of the slope, the traffic load on the trailing edge of the slope was removed in February 2017, as a result of which the monthly deformation increment was reduced to about 90 mm in March 2017. A lawn was planted on the slope in March 2017, preventing infiltration, reducing the monthly deformation increment in April 2017 to 10 mm.

In the slow deformation stage, the total horizontal deformation of the slope increased from 400 to 440 mm with a monthly deformation increment of 6 mm in the 7 months from May 2017 to November 2017. The average rainfall amount in July and August in 2017 reached 110 mm, but because the lawn (planted in March 2017) had grown on the slope by this stage and prevented rainfall from infiltrating into the slope along the fissure, the rainfall scarcely affected the deformation.

In the stable deformation stage, the increment of the total horizontal deformation of the slope within the 2 years (2018 and 2019) was about 30 mm. Since the deformation of the slope is very slow, it can now be recognised as a stable slope without risk of landslide.

Slope stabilisation treatments generally include reducing driving forces, increasing resisting forces, or both (Chen and Liu, 2007; Wang et al., 2014). Many effective techniques have been developed to prevent landslides, such as the pile-anchor combined structure, anchor-lattice beam combined structure and micro-pile structure (Iverson et al., 2000; Chen and Liu, 2007; Tu et al., 2009; Wang et al., 2014; Tang et al., 2015); however, these methods are not completely suitable for the stabilisation treatments of loess slopes with seepage fissures considering the efficiency and cost. Based on the analysis in Fissure-deformation Development Process Under Rainfall and Traffic Load and Deformation Development Process After Stabilisation Treatment, some insights are obtained from the stabilisation measures applied to the Yangpoyao slope.

1) The presence of fissures is one of the important internal factors that affect the trend in deformation development and the stability of loess slopes (Zhou et al., 2014a; Chang et al., 2021). To avoid the development of further fissure-deformation on fissured loess slopes, it is necessary to avoid loading on the trailing edge which results in tensile stress on the slope surfaces. The tensile stress on the slopes will cause the existing micro-cracks in the loess to expand and penetrate, resulting in the formation of fissures on the slopes. The fissures act as channels for rainwater infiltration, greatly weakening the loess and posing serious threats to the stability of such slopes (Jin et al., 2012; Tang et al., 2015). For the fissured slope in the present study, one of the treatments used to avoid tensile stress zones occurring on the slope surface is to avoid traffic load too close to the trailing edge of the slope.