Characteristics of failure mass and safety factor during rainfall of an unsaturated slope

Water is an important factor triggering destabilization of soil masses in slopes, particular in unsaturated condition. Rainfall and snow melting are some of the main sources of natural water which frequently induce unsaturated slope failure in arid/semi-arid region or in active layers (i.e. frequently frozen/thawed layers) in permafrost. If the slopes are in coastal environment or adjacent to a stream, the soil materials from the sliding masses might be carried away by waves and currents, which can further lead to soil erosion problem. Therefore, characterization of the failure mass and safety factor is critical both for evaluation of slope failure hazard but also for erosion problem. This paper presents a method to study unsaturated slope stability including characterizing the size of failure and the safety factor by analyzing a theoretical unsaturated homogeneous slope subjected to rainfall infiltration. The slope is modelled numerically by the finite element method. A shear strength reduction technique is modified for unsaturated soils and employed to evaluate the factor of safety and the extent of the sliding mass. The results show that the size of failure and the factor of safety vary over the rainfall event. Notably, the extent of the sliding mass for unsaturated slopes depends strongly on the infiltration depth. Large global sliding surface tends to be the dominant mechanism if the infiltration depth is very shallow or very deep. Conversely, local restricted sliding surface becomes more dominant if the infiltration is approximately at a certain critical value.


Introduction
Water is known to be one of the main factors contributing to triggering the destabilization of soil masses in unsaturated slopes.Some of the main sources of natural water are rainfall and snow melting.Rainfall has been known to induce slope failure in unsaturated soils in arid/semi-arid regions such as Hong Kong, Australia and the Middle East while snow melting can destabilize the unsaturated active layers in permafrost in Arctic regions such as Svalbard, Alaska and northern Russia.If the slopes are in coastal environments (e.g. in Svalbard or adjacent to a stream, the soil material from the sliding masses might be carried away by waves and currents, which can further lead to soil erosion problem.Therefore, characterization of the failure mass and the safety factor is an important task in order to evaluate the time and the volume of failure.This can subsequently be used for evaluating the hazard of a slope failure and for estimating the erosion rate, particular in the Arctic coastal areas.
This study limits to unsaturated slope failures caused by rainfall infiltration as the source of natural water.Even though catastrophic slope failures induced by rainfall infiltration is one of the major natural hazards facing many regions in the world, understanding of the triggering mechanisms still requires further research.This is particularly important and challenging for slopes in partially saturated soils because the wetting process can lead to a reduction in suction and, possibly, a rise in positive pore water pressure.Both these effects directly contribute to a reduction of the soil shear strength and, hence, destabilise the slope.
This study investigates an unsaturated homogeneous slope subjected to a rainfall event of constant intensity.The slope is modelled numerically by the finite element method (FEM) with application of a shear strength reduction technique (SRT) modified for unsaturated soils.The extent of the sliding mass and the factor of safety are evaluated at different times during and after the rainfall period.The development of the infiltrated depth is observed in parallel with the development of the slip surface in order to investigate the correlation between these factors.

Geometry and boundary condition
The slope model used in this study is detailed in Figure 1.The slope is 10 m height, having an angle of 2:1 (or 26.6°) and resting on a base of 20 m thickness.The initial stress distribution is in equilibrium with gravity (9.81 m 2 /s) while the initial pore water pressure (p w ) distribution is in hydrostatic equilibrium with the initial water table position (5 m below the ground level at the toe of the slope).The hydrostatic condition results in a maximum initial suction (s) of approximately 150 kPa at the crest (AB in Fig. 1) which falls in the low end of the prevalent surface suction range of semiarid or arid environments such as in Australia [1].
A mesh of 1515 quadrilateral elements is employed.This mesh has been shown to produce a reliable estimate of the factor of safety (FoS) of the slope in saturated condition which remained essentially unchanged for finer discretisations [2].The vast majority of elements are approximate squares or parallelograms (Fig. 1).Each square element, and the most of parallelogram elements, has an approximate area of 1 m 2 .The uniformity of the mesh is important for the evaluation of the size of failure as discussed later.Approximately 1% of elements, in the centre of the mesh, are degenerated into smaller polygons (Fig. 1).The degenerated elements are located at considerable distances from the possible locations of slip surfaces and are, hence, unlikely to have any substantial impact on the evaluation of the size of the failure mass.The analysis was conducted using the finite element code Code_Bright [3].

Hydraulic and mechanical models
The van Genuchten (1980) [4] water retention curve and van Genuchten and Nielsen (1985) [5] permeability function are selected because they are capable of representing the soil hydraulic behaviour in a rather simple and numerically stable way.They are also wellaccepted in for modelling flows in unsaturated soils.These constitutive models are presented briefly in equation 1-5.More details about these relationships can be found in Le et )) ( exp( The soil water retention curve proposed by van Genuchten (1980) [4] is given by equation 1, which relates the effective degree of saturation (S e ) (calculated as a function of the current (S), the maximum (S s ), and the residual (S r ) degree of saturation) to suction s=-p w through the air entry suction parameter s e .In equation 2, the parameter s e is in turn related to the porosity (φ) through parameter η which controls the rate at which s e deviates from its reference value s eo when φ deviates from its reference value φ o [10,11].Similarly, equation 3 describes the dependency of the saturated permeability (k s ) on its reference value k so and the difference between φ and its reference value φ o , as proposed by Kozeny (1927) [12].Equation 4 is the van Genuchten and Nielsen (1985) [5] permeability curve linking the relative permeability k r to S e through the parameter m, which can be geometrically interpreted as the curve gradient.The unsaturated permeability k u is the product of the saturated and relative permeability (i.e.k u = k s k r ).Finally, the unsaturated flow q is calculated using the generalized Darcy's law (equation 5).The symbols u w , ρ w , g and z indicate the pore water pressure, the water density, the gravitational acceleration and the elevation coordinate, respectively.
The selected values of the required soil parameters are shown in Table 1.Rainfall intensity I r is assumed to be constant at 43.2 mm/day.Constitutive parameters (e.g.m, k so , s eo ) are selected to be at around the middle of the likely range for each parameter in order to avoid overly large or small input parameters and, hence, unrepresentative or unrealistic results.For example, an very low permeability (i.e.very small k so ) would result in the rainfall having no effect on slope stability because water would not have sufficient time to infiltrate into the unsaturated soil domain.A linear elastic model with an extended Mohr-Coulomb (MC) failure criterion (equation 6) describes the mechanical behavior of the unsaturated soil [13]: (6) In equation 6, the shear stress at failure (τ) depends on the net normal stress (σ) and suction (s) through the friction angle (φ'), cohesion (c') and a parameter controlling the increase in shear strength with suction (φ b ).The component of strength contributed by suction (i.e. the 3 rd term in equation 6) decreases with decreasing s and becomes zero for a fully saturated soil.The values of the strength parameters (c', φ' and φ b ) assumed in this study are typical of a clay and are based on the values reported by Bishop et al., (1960) [14] for boulder clay and by Gan et al., (1988) [15] for a compacted glacial till.The elastic parameter values, i.e.Young modulus (E) and Poisson ratio (ν), are also typical of a clay and chosen within their respective typical ranges of variation [16].The values of all mechanical parameters are summarized in Table 2. 3. Evaluating the FoS and the size of failure for unsaturated slopes

Shear strength reduction technique for unsaturated soils
The shear strength reduction technique (SRT), first suggested by Zienkiewicz et al. (1975) [17], is widely adopted to estimate the slope factor of safety (FoS) of saturated slopes analyzed by the FEM.The FoS is defined as the factor by which the strength parameters must be divided to bring a slope into the state of being "barely stable" [18].Since the FEM allows tracking of the slope deformation up to and including early shear failure, this definition conveniently outlines a technique to calculate the FoS by gradually and systematically reducing cohesion, c' and the tangent of the friction angle, tanφ' until "early failure" is achieved: In the unsaturated soil context, in addition to the conventional shear strength parameters c' and tanφ', the Mohr-Coulomb criterion is also characterised by tanφ b which controls the expansion of the failure envelope with suction increase.The product of s and tanφ b has a similar role to the cohesion and is often combined with c' to create an "enhanced cohesion".Hence, it is reasonable to apply the strength reduction factor to tanφ b in the same manner as to the other two shear strength parameters.In other words, the FoS for unsaturated soil is defined as: During the SRT analysis, an algorithm is also executed to monitor and check the displacement history of all the nodes in the FE mesh.All the nodes that exhibit sudden increase in displacement during the SRT process, and subsequently abnormally large displacement at the end are recorded.The end of the analysis is the point at which the reduction factor could not increase further and is considered to be equal to the FoS.Details of different criteria which are used for this monitoring algorithm have been discussed in Le et al., (2015) [6].

Simulation process
The simulation process involves three stages corresponding to (1) establishing the initial stresses, strains and p w distributions under the action of gravity and the initial hydrostatic condition; (2) establishing stresses, strains and p w distributions at times of interest over the rainfall period; and (3) estimating of the FoS and the size of failure at these times of interest via the SRT.The three stages are presented in detail below: Stage 1: The p w in equilibrium with the hydrostatic condition (i.e.stable horizontal water table at 5 m below the slope toe, EF in Fig. 1) is fixed at every node while the slope is assumed weightless.Subsequently, the unsaturated slope is gradually equilibrated with the standard gravity (9.81 m/s 2 ).Since this gravity loading process is entirely artificial (i.e.not based on an actual physical event) and has the only purpose of setting up a physically reasonable initial stresses/strains distributions, any deformations, and hence porosity changes, caused by this process are reset to zero after reaching equilibrium.
Stage 2: Rainfall is imposed on the crest, the inclined face and the base surface (ABCD) over a period of 10 days (Fig. 1) using the seepage boundary condition.After 10 days, the rainfall is terminated (i.e.removing the seepage boundary condition) but the simulation is continued for another 355 days to investigate recovery of slope stability.
The seepage boundary condition allows water to infiltrate the soil domain at a constant flux equal to I r as long as the p w at this boundary is less than zero.If the p w becomes equal to or larger than zero, the seepage boundary condition automatically switches to a zero p w boundary to avoid ponding of water on the infiltration boundary.
Boundaries OA, OG and GD (Fig. 1) are assumed to be impermeable which allows the infiltrated water to accumulate inside the soil domain and raises the water table during the rainfall period.During this stage, the p w is no longer fixed at every node (as in Stage 1) to enable the natural response of the hydraulic regime to rainfall infiltration.The mechanical deformation induced by rainfall infiltration is fully coupled with the variation of p w .Stage 3: The SRT is employed to estimate the FoS and the size of failure at four times during (at 0, 0.5, 5 and 10 days) and four times after (at 15, 20, 100 and 365 days) the rainfall event.At each time of interest, the FoS is calculated in an independent SRT analysis in which the stresses and strains, obtained at the corresponding time from the simulation in Stage 2, are imposed as the initial stresses and strains conditions at each element.On the other hand, the p w , also obtained from Stage 2, is fixed at every corresponding node throughout the SRT analysis.The restrain of p w at a constant value is to ensure a fully drained process and hence prevent any p w built up as the result of stress changes caused by the reduction of strength.This is desirable because the SRT process is purely a numerical technique to estimate the FoS and does not represent a physically based phenomenon.

Characterization of the sliding mass
In a conventional slope stability study, the FoS is usually the only variable used to evaluate risk.However, the hazard posed by a slope failure is normally judged, in reality, not only by this factor but also by the size of failure or the extent of the sliding mass.A large failure with an extensive sliding mass can cause serious damage while a small one might not cause any damage upon failure.Hence, a large slope failure is regarded from a practical standpoint as more dangerous than a small one even if the former might have a higher FoS than the latter.In addition, the estimation of the size of failure is very important to study soil erosion problem as it directly determines the amount of materials that might be transported way.In reality, the size of failure is often characterized by the volume of soil material sliding down.This study uses a two-dimensional model, hence the size of failure is characterized equivalently by the sliding area (A s ) and the maximum depth from the slope surface to the failure surface to (D s ).
In this study, the algorithm to monitor the displacement history at every node also allows the extent of the slope failure to be analysed.The value of A s is equated to the total number of nodes which has displaced suddenly at some point during the SRT analysis and has displaced in total by a significant amount by the end of the SRT analysis.The reason is that each node lying on the sliding mass represents an area portion surrounding the node which is, in most cases, comprised of a quarter of each of the four elements sharing the node.Since the mesh consists mostly of square or parallelogram elements of similar areas (1 m 2 ), the portion representing by each node is also equal to 1 m 2 .Some nodes locating close to or on the sliding mass boundaries represent portions of the sliding mass that are smaller than 1m 2 .Therefore, the total number of nodes on the sliding mass slightly overestimates the actual value of A s .However, if A s increases, the total number of nodes locating on the mass also increases and vice versa.The variation of the total number of nodes on the sliding mass therefore reflects the variation of A s , and so it is considered satisfactory in this study to equate these two quantities.If the case studied was a practical slope, it is recommended that the estimation of A s is improved either by using much finer mesh to reduce the discrepancy or by identifying nodes on or close to the boundaries and assigning smaller area portions to these nodes.The estimation of D s involves calculating the distances from all the nodes on the sliding mass to the slope face in order to find the maximum value among those distances.Illustration of a sliding mass (after 100 days) imposed on the mesh is shown in Figure 2. The number of nodes within the shaded area is used to approximate A s and D s as discussed.The failure mode is also been categorized into 'global' when the failure surface extends the whole length of the slope face or 'local' when the failure surface limits to a portion of the slope face.The p w contours and corresponding failure masses (indicated by displacement contours) during and after the period of rainfall are presented in Figures 3 and 4, respectively, for the homogenous unsaturated slope considered in this study.The values of FoS, A s , D s and the corresponding failure modes are tabulated in Table 3.
The use of the FEM allows natural and intuitive failure mechanisms to develop where the slip surfaces tend to pass through the weakest soil regions in the area of potential sliding (Figs. 3 and 4).It is important to note for unsaturated slopes that, unlike in dry or fully saturated slopes, the occurrence of the weakest region is significantly affected by the infiltration pattern which controls the reduction of suction and possibly the development of positive p w close to the slope surface.Therefore, the soil area close to the slope face which has been wetted by the rainfall becomes significantly weaker than unsaturated soils in the area underneath the wetting front.The results show that the FoS and the value of A s vary considerably during the rainfall and during the postinfiltration period.At the beginning of the rainfall when the infiltration depth is insignificant or very shallow (i.e.0 -0.5 day), the size of failure tends to be large, however the FoS are quite high (see Fig. 3a-d and Table 3).As the rain infiltrates deeper into the slope face and the wetting front has migrated a certain distance from the surface (i.e.0.5 -5 days), the size of failure becomes smaller and the FoS drops at the same time (Fig. 3c-f and Table 3).Interestingly, as the rainfall continues and the wetting front migrates even deeper (i.e. 5 -10 days), the FoS also continues to drop but the size of failure actually increases again.After the rainfall has stopped (i.e. after day 10), the size of failure keep increasing over time as the p w readjusts to a new steady state during the post-infiltration period (Fig. 4 and Table 3).The FoS also increases after rainfall due to the recovery of suction in the soil region close to the slope face, then the FoS becomes almost stable from 15 days onward.Both the value of A s and the FoS do not increase back to the initial levels (i.e. at time 0) after the rainfall stops due to the accumulation of water during the rainfall.This is the result of the impermeable condition assumed for the vertical and bottom boundaries (Fig. 1).Therefore, the new steady state at the end (day 365) has higher water table than the initial state.
Notably from Figure 3e, water infiltration into the slope face after 5 days seems to have reached a certain critical depth.At this depth, the wetted soil region is just sufficiently large to contain a shallow local failure mechanism (Fig. 3f).The values of A s and D s at 5 days are, in fact, the lowest and are significantly smaller than their corresponding values at earlier times or later times.The FoS reaches the lowest values at the end of the rainfall (10 days) which is caused by the combination of suction loss in the unsaturated zone and water table rise in the saturated zone.However, the size of failure at 10 days is larger than at 5 days (Fig 3 . f and h) because the wetting front has migrated further which increases the size of the wetted region which contains the sliding mass.
Figures 5 and 6 compare the variation of the FoS, A s and D s over the rainfall period.The variation of failure depth is qualitatively similar to in the variation of failure area.The fact that the size of failure drops to the smallest value during the rainfall (5 days), but the FoS becomes E 2016 the smallest at the end of the rainfall (10 days) suggests that the hazard evaluation of unsaturated slope failure should combine these two types of information (e.g.[9]).
Note that the results for the FoS in this study have been verified using the limit equilibrium method.This is performed using SLOPE/W (GEO-SLOPE International, Ltd, Alberta, Canada) to analyse the unsaturated slope with p w distribution obtained from SEEP/W (GEO-SLOPE International, Ltd, Alberta, Canada).The results show good agreement between the FoS obtained by the FEM in this study and that obtained by GeoSLOPE commercial software.The localized failure mechanisms suggests that the method applied in this study has the potential to be use for investigating the retrogressive failure mechanisms triggered by water infiltration or similar factors (e.g.temperature rise in the cold region).If for example, the FoS corresponding to the localized failure at 5 days is smaller than 1 (e.g. the slope actually fails), then continuous rainfall can subsequently infiltrated into the exposed slip surface and cause subsequent failure once it has infiltrated deeper than the critical depth again.The process can continue if rainfall continues which lead to successive or retrogressive failure.

Conclusion
This paper has presented a method to study the failure mechanism of an unsaturated slope in which the failure has been characterized by both the factor of safety and its size in terms of area and depth.The conventional shear strength reduction technique was modified and adapted for unsaturated soils.The size of failure and factor of safety were evaluated by tracking the displacement history and all the nodes and recorded the nodes that are located on the sliding area.
The results show that the size of the failure and the factor of safety vary with time during the rainfall and also for a period after the rainfall has stopped.The size of failure reaches the smallest value at the transient state during the rainfall while the factor of safety drops to the lowest at the end of the rainfall.Notably, the extent of the sliding mass depends on the infiltration depth.Large global sliding surface tends to be the dominant mechanism if the infiltration depth is very shallow or very deep.Conversely, local restricted sliding surface becomes dominant if the infiltration is approximately a certain critical value which leads to the failure mass to confine to the wetted region.

Figure 1 :
Figure 1: Slope dimensions, geometry and boundary conditions (scale in metre)

Figure 2 :
Figure 2: Estimation of slip surface area (A s ) and depth (D s ) after 100 days

Figure 5 :
Figure 5: Variation of FoS and A s with time

Table 1 :
Values of hydraulic parameters

Table 2 :
Values of mechanical parameters

Table 3 :
Values of studied variables at different times