Stability analysis of unsaturated soil slope considering softening and non-softening approach

. In conventional geotechnical engineering practice, peak shear strength parameters are widely used in the design of geo-structures constructed with or within unsaturated soils. However, a reduction in shear strength from the peak to the residual value is typically observed during the shear deformation in typical fine-grained unsaturated soils. Several geo-structures in unsaturated soils undergo a large shear deformation prior to reaching the failure condition. Thus, the factor of safety of such geo-structures will gradually but significantly decrease with the development of the shear deformation. For this reason, the strain-softening behaviour of unsaturated soils should be considered for reliable design of the geo-structures in unsaturated soils. In this study, an unsaturated clay slope under 10 years rainfall infiltration were modelled using softening and non-softening approach using commercial finite element software SIGMA/W. The responses of the studied slope to the long-term rainfall infiltration were analysed. The results of the softening and non-softening analysis were compared. This study provides valuable information with respect to the significance of the strain-softening that are useful in the rational design of slopes in unsaturated soils.


Introduction
The shear strength of unsaturated soils gradually reduces from a peak to a residual value when there is a significantly large shear deformation. Such a phenomenon is typically referred to as strain-softening behaviour. The geotechnical infrastructures in unsaturated soils typically exhibit progressive failure behaviour due to strain-softening. For example, local shear zones develop in an unsaturated soil slope due to various environmental factors (e.g., rainfall infiltration and drying-wetting cycles) [1][2][3]. In the local shear zones, the shear strength of unsaturated soils will be reduced and the extra shear stress which is greater than the post-peak shear strength will be redistributed to the surrounding zones. As a result, large deformation is developed progressively in a typical fine-grained unsaturated soil slope; due to this reason, the factor of safety reduces gradually and reaches a failure condition.
In traditional slope stability analysis, non-softening approach (i.e., only the peak shear strength is considered) is conventionally used. Such an approach underestimates the slope deformation and hence overestimates the factor of safety. For this reason, a softening analysis approach that can consider the shear strength reduction of unsaturated soils during shear deformation is required for rational analyses and design of geotechnical infrastructures associated with finegrained unsaturated soils.
From 1990s, the softening approach has been widely used for the stability analysis of saturated soil slopes [4][5][6][7][8]. In recent years, there were also several studies considering the softening approach in the stability * Corresponding author: vanapall@eng.uottawa.ca analysis of unsaturated soil slopes [9][10][11][12]. These studies highlighted the significance of strain-softening in modelling the behaviours of unsaturated soil slopes.
In this paper, an approximate approach is developed considering the strain-softening of unsaturated soils for modelling the progressive failure of unsaturated soil slope using commercial finite element software SIGMA/W. An unsaturated clay slope under 10 years rainfall infiltration was simulated using softening and non-softening approaches, respectively. The responses of unsaturated soil slope to long-term rainfall infiltration were analysed. The mechanical behaviours of the studied slope were compared between the softening and non-softening analysis. Fig. 1 shows details of a numerical model that was established to study the behaviours of unsaturated soil slope under long-term rainfall infiltration. The studied soil slope has a 20 m height and 30 m length. The rainfall infiltration typically has a significant influence on the soils in a shallow depth [13]. Thus, finer mesh (element size = 0.5m) was used for the region within the depth of 3m below the ground surface, while larger mesh (element size = 1m) was used for the other regions.

Model slope
The initial stresses were established by using the 'Insitu' analysis of SIGMA/W. The initial pore water pressure was determined based on the position of ground water table, which was set at the bottom of the numerical model. This can result in an initial suction of 392 kPa at the crest of the slope.

Boundary conditions
The horizontal displacement was fixed, and the vertical displacement was free at the right and left lateral boundaries. Both the displacements in the horizontal and vertical directions were fixed at the bottom boundary. The water flux was assumed to be zero at the two lateral boundaries and the bottom boundary.
The water flux was assumed to be a step function of time on the ground surface for simulating the rainfall infiltration. The water flux on the ground surface was changed 30 days (i.e., 1 month) intervals. The rainfall for 3600 days (i.e., around 10 years) was considered.
The information of daily rainfall from January 2013 to December 2022 for Ottawa region was collected from 'weatherstats.ca based on Environment and Climate Change Canada data'. The monthly average rainfall was calculated and shown in Fig. 2, which was used to formulate the step function of the water flux on ground surface vs. time.

Material properties
The elastic-plastic model (i.e., elastic, perfectly-plastic model) in SIGMA/W [14] was used in this study. In the elastic-plastic model, SIGMA/W can calculate the shear strength of unsaturated soils using the equation below.
where  is shear strength; (ua -uw) is suction; (n -ua) is net normal stress; c' is effective cohesion; ' is effective internal friction angle;  is volumetric water content that can be estimated using soil water characteristic curve (SWCC); s is saturated volumetric water content; r is residual volumetric water content (assumed as 5% of s by default [14]).
A CH soil reported in the literature [15] was used as the material of the studied slope. A series of single stage suction-controlled ring shear tests were conducted on the CH soil under different net normal stresses and matric suctions [15]. Thus, the peak and residual shear strength of the CH soil corresponding to different combinations of net normal stress and matric suction can be determined from the reported suction-controlled ring shear tests. The peak shear strength parameters (c'p and 'p) and the residual shear strength parameters (c'r and 'r) can be approximately estimated by fitting the test results using Eq. 1, which is shown in Fig. 3a.
For reproducing the strain-softening, the shear strength parameters of the CH soil were assumed to decrease from the peak to the residual value with increasing deviatoric strain, d, as shown in Fig. 3b. The value of d,p and d,r was assumed to be 4% and 20%, respectively in this study. The deviatoric strain was defined as Eq. 2 in SIGMA/W [14].
where x, y and z is normal strain in the x, y and z direction; xy is shear strain.
The SWCC of the CH soil was also reported with corresponding dry unit weight and specific gravity [15]. The SWCC was estimated using Fredlund and Xing's equation [16] in this study. In addition, the permeability function of the studied soil can be estimated automatically by SIGMA/W based on Fredlund and Xing's SWCC equation. The other soil properties were assumed to be common values. All the properties used in this study were summarized in Table 1. 3 Softening and non-softening analysis approach SIGMA/W provides a 'stress redistribution' analysis that is used for solving over-stressing problem (i.e., the current computed state of stress is higher than the soil strength) [14]. The over-stressing typically happens when there is an increase in the pore water pressure while the total stress remains constant in the case of rainfall infiltration or when there is some strength loss due to a change in the soil particle structure. The 'stress redistribution' analysis in SIGMA/W can calculate the unbalanced load (i.e., a portion of the stress is unbalanced by the available shear resistance) in each element. The unbalanced load is applied to the mesh which can create deformation and stress changes. This procedure is iteratively repeated until the unbalanced load is zero or within some tolerable convergence criteria.
A modelling procedure that can approximately consider the strain-softening of unsaturated soils was developed based on the stress redistribution analysis in SIGMA/W, which is described as below: (i) Establishment of subregions: the active zone of the slope was divided into several subregions as shown in Fig. 4. The material properties were defined for each subregion separately. It is assumed that, within each subregion, the shear strength parameters have same values and reduce with the average deviatoric strain of the subregion.
(ii) SIGMA/W coupled analysis: After establishing the initial condition, the boundary condition on the ground surface was changed to the water flux boundary of the first year (i.e., the rainfall infiltration from the 1 st to the 12 th month in Fig. 2). Each month is assumed to have 30 days. A coupled analysis was conducted for 360 days with a time increment of 30 days.
(iii) SIGMA/W stress redistribution analysis due to rainfall infiltration: a stress redistribution analysis was conducted using the coupled analysis in step (ii) as parent analysis. This procedure can adjust the distribution of deformation and stress created in coupled analysis by removing the influence of over-stressing caused by rainfall infiltration.
(iv) SIGMA/W stress redistribution analysis due to strain-softening: the deviatoric strain created in last step at all the nodes within each subregion can be exported. An average deviatoric strain can be calculated for each subregion. Therefore, the values of c' and ' of each subregion can be calculated based on the corresponding average deviatoric strain following Fig. 3b. Then, the c' and ' of each subregion in the numerical model were changed to the new values and a stress redistribution analysis was reconducted to create new deformation and stress. The deformation and stress created in this step was considered as the responses of the studied slope to 1 year rainfall taking account of strain-softening of unsaturated soils.
(v) The boundary condition on the ground surface was changed to the rainfall infiltration of the next year.
Step (ii) -(iv) were conducted using the last stress redistribution analysis as parent analysis. This procedure was repeated until the simulation of the studied slope under 10 years rainfall (i.e., 3600 days in the numerical model) was finished.
The modelling procedure of the non-softening approach is similar with that of the softening approach.
However, in non-softening approach, the c'p and 'p was used as effective shear strength parameters of all the subregions throughout the modelling procedure. Thus, the stress redistribution analysis due to strain-softening can be neglected.  softening approach. As discussed in last section, for modelling the studied slope for each year, a coupled analysis was conducted first to create the deformation and stress due to rainfall infiltration. As a next step, a stress redistribution analysis was conducted to adjust the deformation and stress by removing the over-stressing caused by the rainfall infiltration and strain-softening. For this reason, the displacement-time curve exhibits a staged behaviour (i.e., the dash lines in Fig. 5). The stages of the displacement-time curve are not obvious in the non-softening analysis. However, in the softening analysis, the stages of the displacement-time curve are well-defined with increasing time as soil layer in the slope is strain-softened. A smooth displacement-time curve can be obtained by linking the data points after the stress redistribution analysis of each year (i.e., the solid lines in Fig. 5). The smooth displacement-time curve will be used to analyse the responses of the studied slope to long-term rainfall infiltration. As shown in Fig. 5, within the first three years, the results of the softening and non-softening analysis are almost the same. The slope exhibited increasing horizontal displacement and heave at the top of the slope. This can be attributed to the soil expansion caused by the loss of suction during rainfall infiltration. Fig. 6 shows the suction profiles within the depth of 9 m below the point of interest in Fig. 5. It can be found suction decreased significantly in the first three years.

Analysis results
Figs. 7a and b show the contours of the deviatoric strain in the studied slope after 3 years rainfall. It can be found the deviatoric strains are also the same for softening and non-softening analysis. The deviatoric strain was still small (< 6%) in most regions of the slope. In a shallow depth, the deviatoric strain distributed relatively uniform along the slope. This means the strain-softening did not happen yet in the slope within the first three years, which contribute to the same results for both the softening and non-softening analysis.
From the 3 rd to the 7 th year, the horizontal displacement at the top of the slope increases gradually with a fluctuating increase rate in softening analysis. In the meantime, the vertical displacement at the top of the slope remains almost constant. This means the slope was relatively stable despite the development of some deformation.
Such a behavior can be attributed to the combined effects of the soil expansion caused by rainfall infiltration and the slope sliding caused by strainsoftening. It is important to note that the suction loss continues from the 3 rd to the 7 th year due to the rainfall infiltration. However, the suction loss is less in comparison to the first three years (Fig. 6). This can contribute to an increase in the horizontal displacement and the heave at the slope top. At the same time, strainsoftening initiates and develops progressively during this period (i.e., from the 3 rd to the 7 th year). Due to this reason, an increase in the horizontal displacement and the settlement at the slope top will arise.
Figs. 7c, e and g show the contours of deviatoric strain from the 3 rd to the 7 th year. The results suggest a relatively large deviatoric strain (> 20%) first occurred in a shallow depth (within 0.5m below the ground surface) near the top of the slope (Fig. 7c). After that, the region of large deviatoric strain developed to a larger extent downwards along the slope and to a deeper depth (Fig. 7e). After 7 years rainfall, the region of large deviatoric strain developed to 70% of the slope and reached a depth of 1.5m below the ground surface (Fig.  7g). The maximum deviatoric strain in the slope reached around 40%. The region of the maximum deviatoric strain was almost parallel to the slope. From the 7 th to the 10 th year, the horizontal displacement and the settlement at the top of the slope increases significantly with a growing increase rate in softening analysis. This means the slope became unstable and obvious sliding occurred during this period. Fig. 7i shows the contours of deviatoric strain after 10 years rainfall. It can be found the region of large deviatoric strain has developed throughout the entire slope and has reached a depth of 6m below the ground surface. The maximum deviatoric strain reached more than 70%. of the non-softening analysis deviated from the curve of the softening analysis (Fig. 5). The difference between the results of the softening and non-softening analysis become more pronounced with increasing time.
From the 3 rd to the 10 th year, the horizontal displacement and the heave at the slope top in the nonsoftening analysis increased slightly (Fig. 5). In other words, in the non-softening analysis, the soil expansion caused by suction loss was still a predominant factor that contributed to the deformation of the slope during this period (i.e., from the 3 rd to the 10 th year). The obvious sliding of the slope cannot be captured by the nonsoftening approach.
Figs. 7d, f, h and j show the development of deviatoric strain of the slope from the 3 rd to the 10 th year in non-softening analysis. These results suggest the distribution of deviatoric strain is relatively uniform along the slope throughout the 10 years rainfall. The maximum deviatoric strain is around 17%, which only accumulated in a small region at the toe of the slope. The progressive development of the large deviatoric strain (shown in Figs. 7c, e, g and i) cannot be reproduced by the non-softening approach.  ( ) 13 2 r skeleton a a w sr where 1 is major principal stress; 3 is minor principal stress.
The results suggest that stress path of the softening and non-softening analysis was the same within the first three years. The maximum shear stress increased first and then decreased slightly, while the mean skeleton stress kept decreasing. After 3 years rainfall, the stress path reached the peak shear strength envelope.
From the 3 rd to the 7 th year, the stress path in the softening analysis moved along the peak shear strength envelope. However, no strain-softening occurred during this period. This is consistent with the results of deviatoric strain shown in Figs. 7c, e and g. The large deviatoric strain developed mainly within the depth of 1.5 m from the 3 rd to the 7 th year.
From the 7 th to the 10 th year, the stress path in the softening analysis deviated from the peak shear strength envelope and moved towards the residual shear strength envelope. This indicates the strain-softening of the soil at the depth of 3 m below the ground surface.
In comparison to the softening analysis, the stress path in the non-softening analysis only moved along the peak shear strength envelope without any reduction from the 3 rd to the 10 th year.

Discussion
The modelling procedure developed in this study can reproduce the strain-softening behaviours of unsaturated soils and simulate the progressive failure of unsaturated soil slope. However, there are some limitations.
Firstly, SIGMA/W does not have a spatial function for defining the shear strength parameters of each node of the numerical model. For this reason, the active zone of the slope was divided into several subregions. The shear strength parameters were defined for each subregion separately and were reduced according to the average deviatoric strain of each subregion. This method Secondly, in the coupled analysis of SIGMA/W, the shear strength parameters are assumed to be constant. The reduction in the shear strength parameters should be conducted separately at different time increments based on the stress redistribution analysis. Therefore, the development of large deviatoric strain and strainsoftening during the couple analysis may be underestimated.
Thirdly, the stress redistribution analysis of