Stability Analysis of Road Embankment Slope Subjected to Rainfall Considering Runoff-Unsaturated Seepage and Unsaturated Fluid–Solid Coupling

Rainfall is an important triggering factor influencing the stability of soil slope. Study on some influences of the rainfall on the instability characteristics of unsaturated soil embankment slope has been conducted in this paper. First, based on the effective stress theory of unsaturated soil for single variable, fluid–solid coupling constitutive equations were established. Then, a segment of red clay embankment slope, along a railway from Dazhou to Chengdu, damaged by rainfall, was theoretically, numerically, and experimentally researched by considering both the runoff-underground seepage and the fluid–solid coupling. The failure characteristics of the embankment slope and the numerical simulation results were in excellent agreement. In the end, a sensitivity analysis of the key factors influencing the slope stability subjected to rainfall was performed with emphasis on damage depth as well as infiltration rainfall depth. From the analysis in this paper, it was concluded that the intensity of rainfall, rainfall duration, and long-term strength of soil have most effects on slope stability when subjected to rainfall. These results suggest that the numerical simulation can be used for practical applications.


Introduction
Rainfall in mountains and hills regions, such as in Southwest China, is often heavy and prolonged; thus, disasters associated with slope instability induced by rainfall usually happen, resulting in severe impacts to the transportation systems, human beings' safety, and wellbeing [1][2][3]. Therefore, it has become a significant and urgent subject to take rainfall into account as a considerable factor in slope stabilizing [4,5].
Slope instability induced by rainfall is affected by both water flows on slope surface and water infiltration. Generally, water flow on slope surface produced by rainfall mainly includes two basic types [6,7], of which one so-called "super-infiltration surface runoff" occurs when the rainfall intensity exceeds the infiltration capacity of the slope, and the other occurs after the rainfall along the instauration shallow regions and finally transfers into surface flow when the regions are saturated. In geotechnical engineering, problem caused by 'super-infiltration surface runoff' is not uncommon, so more attention has been paid on it. The process of water infiltration into the slope induced by the rainfall consists of two stages [8]: the first is the "complete-infiltration stage," which happens when the rainfall intensity and the duration are low. The rain water is able to completely infiltrate into the soil matrix and there is little residual water to form surface runoff. The second is pressured-infiltration stage which forms after the soil has become saturated. These two stages are not strictly separated or sequent because the infiltration and runoff are determined by many factors, such as the soil percolation characteristics, rainfall intensity, rainy duration, and groundwater table. Most of the current analysis methods for rainfall infiltration usually simplify the surface runoff, in which the assumption is taken that the surface runoff completely drains away instantly, and the depth of surface standing water is zero. By this way, the boundary conditions and rainfall influences on slope infiltration are taken into consideration [6,[9][10][11][12]. As a consequence, these analyzing methods cannot completely reflect the actual processes and mechanisms of rainfall infiltration and soil yielding when subjected to runoff [6,13]. Accurately estimating the runoff in surface areas and also considering the effects of depth of standing water on the infiltration could be helpful to accurately estimate the quantity of the infiltration, as well as the variation of the pore water pressure. Hence a more precise yield analysis of slope stability and slope reinforcement can be presented. At present, the numerical analysis methods of slope stability under rainfall condition have the limit equilibrium method [6,[9][10][11]14], finite element method [15][16][17], numerical simulation of fluid solid coupling of unsaturated soil [18,19]. The stabilization enhancement of soil slope can be achieved by various means such as consolidation and compaction, soil replacement, physical and chemical improvement, earth reinforcement, relief well and other anti-water infiltration and retaining engineering facilities [20][21][22][23][24].
With the rapid development of high-speed rail and highways in mountain regions of China, attentions focused on slope instability induced by rainfall, especially for natural soil slope and road cuts, have become higher and higher. This type of slope failure is usually characterized by shallow failure whose destroying depth is usually less than 3 m above the groundwater table, and generally with small volume on steep soil slopes of 30°-50° [25,26]. In addition to China, in many other regions, such disasters induced by rainfall are also common and confusion, such as Hongkong [25], Taiwan [27], and Korean [28]. Similar researches for the triggering mechanism and treatments of the rainfall inducing slope failure have been performed. Some of these achievements have been combined into the local engineering practices [29][30][31]. In regions of rainfall-induced slope failure, pore water pressure and infiltration force may form in the region of slope surface, and may cause shallow surface failure, then gradually induce complete failure of the entire slope. This problem has gradually gained wide recognition.
Relative engineering projects and researches [6,9,11,32] have been conducted by simplifying boundaries of rainfall infiltration. The temporary distribution of matrix suction has been calculated using the unsaturated seepage of underground water, and the stability assessment method of limitedequilibrium theory based on the shear-strength theory of unsaturated double stress invariants ( − u a ) and (u a − u w ). However, the solid-fluid coupling effect during infiltration is sparsely been considered. In addition, the failure of slope is gradual due to the fact that the stress deviates from the average of the soil. Thus, failure first occurs in regions of high shear stress, and then causes the nearby regions to fail as a result of shear stress transfer to the un-softened soil. This process continues until a new equilibrium forms, ending in local failure or complete failure. For road embankments subjected to rainfall, failure of the slope surface often occurs and then causes cracks or subsidence to the entire body, but in short time, total and complete, instability rarely happens. Even so, 1 3 cracks or subsidence are also not allowed, especially for highspeed railways.
In this paper, the fluid-solid coupling equations of twophase unsaturated air-water have been established based on the effective stress theory of single variable unsaturated soil. By using the well-known FLAC software [33] and its own FISH program and C++program connector, the solid-fluid coupling program of ground surface flow and underground infiltration has been implemented. At the same time, the road foundation stability under rainfall has been analyzed using the coupling method.

Constitutive Model of Runoff-Underground
Seepage Coupling

Flow Model of Surficial Water
For simplification, the rainfall intensity is assumed to be a constant value. In order to consider the influence of slope and rainfall direction to the slope surface, in practice, the Chezy formula or Manning Equation can be taken to substitute the momentum equation, and with unit-width flow, that is g = uh on the slope surface. Subsequently, the governing formulas of Saint Venant Slope Current can be revised as (1) [34] where n is the coefficient of roughness as defined by Engman [35]; x is coordinate value of x-axis; g is the acceleration of gravity; and u, h, and q stand for the velocity of flow and location of x on the slope, the water depth as well as the rate of flow, respectively; , , and I, respectively, express the slope angle, the rainfall direction (which is determined by the intersection angle to the vertical direction, with the counterclockwise direction as the positive direction), and rainfall strength, respectively, as shown in Fig. 1.

Flow Model of Underground Water
By the continuity equations of the underground water movement, the basic differential equations of two-phase water-air flow in embankment can be deducted from the principle of the mass conservation law, which is detailed as (2) where S a and S w , respectively, are saturations of air and water, and they satisfy the relationship that S a + S w = 1; H a and H w are the air head and the water head, respectively; q a and q w refer to the sources of air and water, respectively; and a and w refer to density of air and water, respectively; x i is the Cartesian coordinates, i = x, y, and z.

Implementation of Surface Runoff and Underground Water Infiltration Coupling
In conjunction with the governing equations of surface water and underground water (1)-(3), the coupling is realized by way of infiltration between surface and underground water. The alternating iterative method was taken to analyze the coupling of surface and underground water. The flow chart of the coupling process is provided in Fig. 2.

The Fluid-Solid Coupling Model of Unsaturated Soil
The particles of the soil are assumed to be incompressible, and the characteristics of fluid-solid interaction can be stated as follows: (1) the effective stress increment is Terzaghi effective stress, and the pore pressure is replaced by fluid pressure increment of average saturated weight; (2) the volumetric strain is produced by effective stress according to Terzaghi Theory; (3) the volumetric strain affects the changing of fluid pressure; (4) Bishop effective stress [36,37] is applied to check the plastic yielding of the constitutive model. Based on Biot consolidation theory, unsaturated soil fluid-solid coupling model must satisfy the fluid equilibrium equation, momentum equilibrium equation, compatibility equation, as well as other requirements such as conduction, capillarity, fluid and mechanical constitutive equation laws, respectively.

Conduction Law
The governing equations of water and air migrations can be expressed as (3) by Darcy's Law:  where k ij is saturated flow coefficient, which is a tensor defined as the ratio of intrinsic permeability to dynamic viscosity; r is relative permeability of the fluid, an empirical function of saturation S w ; is the dynamic viscosity; P is the pore pressure; is the fluid density; g is the acceleration of gravity.

Capillary Law
Capillary pressure is the difference of pore pressure in the fluid, expressed as P c = (P g − P w ), which is an empirical function of effective saturation.

Flow Equilibrium Law
For the slightly compressible fluid, the equilibrium relation is as follows: where is the change of fluid volume (that is the change of the volume of unit-volume fluid in the porous material); q v is the volume strength of fluid source.

Fluid Constitutive Law
where K w and K g are bulk modulus of fluid and air, respectively; is the volumetric strain.

Momentum Equilibrium Equations
where is the volume density, and = d + n(s w w + s a a ), in which a and w are air and fluid densities, respectively; d is the dry density of the soil; and ̇u is the velocity.

Mechanical Constitutive Law
The incremental constitutive model of the pore matrix is expressed as where Δ � ij is the change of the effective stresses, whose expression is [36,37]; H is the functional form of the constitutive law. The Coulomb-Mohr Strength Model was adopted in this paper; is the parameter of stress-paths.

Geometric Compatibility Equation
According to the assumption of small-strain and the assumption that the compressive strain is set as positive, the geometric equations of soil are expressed as

Basic Conditions of the Slope
To ensure whether the influence of fluid-solid coupling on the embankment should be considered, the finite element meshed model shown in Fig. 3 is established. Due to symmetry of the embankment, only the left part of the embankment is illustrated, in which, the height of the slope is 10 m, the slope ratio is 1:1.5, and the width of the slope crest is 10 m too. Only the infiltration produced by the vertical rainfall on the slope surface has been considered, and the flow of air is also allowed. In this model, the fillers of the embankment consist of red clay soil (whose particles were crushed to below 2 mm), the compaction coefficient is 0.9 (defined as the ratio of dry density to the maximum dry density). The soil parameters are listed in Table 1, where the parameters of unsaturated soil were inversely deduced based on RBF neural network model [38], the coefficient of roughness is 0.035, and the initial saturation S w is 0.5. The rainfall intensity is 100 mm h −1 ; the duration of rainfall is 3 h, so the rainfall total sums up to 300 mm. The analyses are carried out respectively with and without consideration of the fluid-solid coupling.

"Factor of Safety (FOS)" Distribution
Figures 4 and 5, respectively, show the zones where the FOS is less than 1.0, which are actually the potential failure ones in the slope. In Fig. 4, the fluid-solid coupling (in which the unsaturated infiltration calculation is ahead of the mechanical calculation) is not considered, whereas, in Fig. 5, the fluid-coupling effect was considered.
The FOS is defined as follows: where c s and s are cohesion strength and internal friction angle, respectively, under the corresponding water saturation of "1"; ′ 1 and ′ 3 are the maximum and the minimum  principal effective stresses based on the Bishop effective stress principle [37]. Comparing Fig. 4 with Fig. 5, it can be indicated that the slope did not fail when fluid-solid coupling was without considered. Under the condition that the fluid-solid coupling was considered, shear failure zones have extended 0.15 m in depth from the slope toe to 0.33 m from the slope top, forming a concave band whose ends bulge outwards. The bulging zones have formed a range between 0.15 and 1.43 m from the toe, and another range between 0.33 and 1.04 m from the top. The shear failure zones substantially parallel the slope surface. This indicates that the seepage pressure produced by infiltration fluid and the pore pressure changes caused by the volume change is obvious influencing the failure depth and range, so the fluid-solid coupling cannot be neglected. Figure 6 shows the variation of the FOS at section "1-1." For plane Section "1-1," (shown in Fig. 3), due to the shallow failure performance, only a horizontal depth of 1.5 m away the slope surface is considered in this study. From Fig. 6 at the section "1-1," one can find that the FOS of the slope is all larger than performing in no failure occur when fluid-solid coupling is not considered, whereas, a shear failure zone with a width of 0.21 m developed when the fluid-solid coupling influence has been included.

A Project Illustration Using the Coupling Model
Along the Dazhou-Chengdu Railway Line (connection two main cities in east and central Sichuan Province, southwestern China) between K240 + 10 and K240 + 70, due to a prolonged rainfall lasting from Aug. 7th-8th 2005, plastic failure occurred on the slope, resulting in the road shoulder's sinking and cracking. The field investigation revealed that at the bottom of the roadbed there was a lateral slope of about 5°. The length of the right side was 21.5 m and that of the left side was 15.3 m. The slope gradient was 1:1.75; strip rocks were used as the slope protection. Seven arches were formed with it. Each arch was about 6.4 m wide. Five of the arches were squeezed and bulged. The relatively worse soil plastic deformation was notable in two of those arches (Fig. 7). Through field sampling analysis, the soil in the slope body was determined to be red clay fill. The compaction coefficient had a value of 93%, and the moisture content was 11.5%. Based on the local meteorological data, the rainfall intensity was 5 mm h −1 and lasted about 50 h. The impact of the stone arch slope protection is not considered. Assuming the slope is bare, and adopting the simplified calculating model as shown in Fig. 7, it is found that the zones of plastic failure calculated by numerical calculation conform to the bulging zones observed in the slope. This verifies the feasibility and validity of the coupling analytic and calculating theory determined previously above for the analysis of slope stability during rainfall.

Investigation of the Effect of Each Factor
In order to demonstrate how each factor impacts the stability of the slope, calculations have been conducted by varying individual factors, such as slope gradient, slope height, percolation properties, rain intensity, compaction coefficient, and long-term strength. Two key indexes " L/L 0 " (ratio of failure length to the total length of the slope surface) and failure depth (depth perpendicular to slope surface) are selected to reveal the influences. The soil parameters of the red clay fill are listed in Table 1. The value range of each factor is determined according to relative engineering practices and the authors' experiences. The values listed have basically concluded the situations of the realistic engineering in similar projects. In the end, the factors which mainly influence the results have been selected on a comparative basis.

The Effect of Slope Gradient
The slope gradient was classified into six cases as 1:1. In Fig. 8, it is also shown how the failure depth of the slope at section "1-1" varies with slope gradient variation after the 3-h-long 100 mm h −1 rainfalls. With the slope steepening, the failure depth increases and slowly changes; when the slope gradient remains between 32.83° (1:1.55) and 34.59° (1:1.45), the damaged depth varies greatly and is within a range of 0.12-0.31 m; hereafter the failure depth almost remains constant with further increases in the slope gradient. So in order to decrease the influence of rainfall on the slope and maintain the stability of the slope affected by the rainfall, it is advisable to control the slope gradient below 32.83° (corresponding slope ratio 1:1.55).

The Effect of Slope Height
To analyze the influence of the slope height on the slope stability during rain, the slope height was classified into five cases as 6, 8, 10, 12, and 14 m, and the other conditions are set as described in Sect. 4.1. In Fig. 9, it is shown how L∕L 0 varies with the slope height after the 3-h-long 100 mm h −1 rainfall. L∕L 0 increases continuously with the slope height's increase. When the slope height is near 6 m, L∕L 0 reaches its minimum 0.91 and when the height is 14 m, L∕L 0 reaches its maximum 0.97. The failure depth of the slope at section "1-1" is also shown in Fig. 9. It varies with the slope height changes after the 3-h-long 100 mm h −1 rainfall. The failure depth remains nearly constant, at 0.29 m, whatever the slope height.

The Effect of Horizontal Percolation Properties
The vertical permeability was set as a constant K 11 = 2.81 × 10 −7 m/s, and the ratio of the horizontal permeability parameter to the vertical one was classified into five cases: K 11 ∕K 22 =1, 2, 4, 6, 8. The other conditions are the same as described in Sect. 4.1. The results are summarized in Fig. 10. Figure 10 shows how L∕L 0 varies with K 11 ∕K 22 after the 3-h-long 100 mm h −1 rainfall. From Fig. 10, it is found that L∕L 0 changes little with K 11 ∕K 22 's increase. When the K 11 ∕K 22 is 4, the ratio of L∕L 0 reaches its maximum, 0.992, and when the K 11 ∕K 22 is 1, L∕L 0 reaches its minimum 0.952. Figure 10 shows how the failure depth of the slope at section "1-1" varies with K 11 ∕K 22 after the 3-h-long 100 mm h −1 rainfall. This figure indicates that the failure depth increases and slowly changes with the horizontal permeability parameters' increase; when the ratio of K 11 ∕K 22 stays within 4 and 6, the failure depth varies greatly with a range of 0.29-0.44 m. Thus, to decrease the effect of rainfall on the slope's stability, it is advisable to control the ratio of K 11 ∕K 22 less than 4.

The Effect of Rain Intensity
The rainfall strength was classified into five cases:R = 100 mm h − 1 , t = 3 h; R = 80 mm h − 1 , t = 3.75 h; R = 60 mm h − 1 , t = 5 h; R = 40 mm h − 1 , t = 7.5 h; and R = 20 mmh − 1 , t = 15 h, and the other conditions are the same as described in Sect. 4.1. Figure 11 shows the results that how L∕L 0 varies with changing rainfall intensity under a constant rainfall total. From this figure, it is revealed when R = 100mmh − 1 and t = 3 h, L∕L 0 reaches its minimum 0.952 and when R = 40 mm h −1 , t = 7.5 h, the ratio of L∕L 0 reaches its maximum 1.0. Figure 11 shows how the failure depth of the slope at section "1-1" varies with the rainfall intensity changing under a constant rainfall in total. According to Fig. 11, one finds that with the rainfall strength decreasing and time extending, the failure depth at section of '1-1' increases; however, when the rainfall intensity R stays between 60 and 100 mm h −1 , the failure depth only varies a little and reaches about 0.29 m; when R falls less than 60 mm h −1 , the failure depth varies much more, and when R reaches 20 mm h −1 , the failure depth is up to 0.55 m.

The Effect of the Initial Saturation of Soil
The initial saturation was classified into five cases: S = 0.35, 0.45, 0.55, 0.65, and 0.75, and the other conditions are the same as described in Sect. 4.1. Figure 12 shows how L∕L 0 varies with the initial saturation changing after the 3-h-long rainfall intensity of 100 mm h −1 . From the figure, one can find that L∕L 0 does not keep decreasing with the initial saturation increasing, and when the initial saturation reaches 0.55, L∕L 0 reaches its minimum 0.95; hereafter, L∕L 0 gradually increases, with small changes. Figure 12 also shows how the failure depth of the slope at the section "1-1" varies with the initial saturation after

The Effect of Compaction
In the process of on-site compaction, the slope edge normally could not be tightly compacted due to the slope edge lacking lateral restraints. To analyze the influence of the compaction of the slope edge on the slope stability during raining, the zone 3 m wide from the slope edge to the slope body was assumed weakly compacted and its compaction coefficient was classified into four cases as 87, 90, 93, and 95%, and other zones of the slope were well compacted and their compaction coefficient was 98%. The other conditions are the same as described in Sect. 4.1. Figure 13 shows how the ratio of L∕L 0 varies with the compaction coefficient changing after the 3-h-long 100 mm h −1 rainfall. From this figure, it is revealed that L∕L 0 does not keep decreasing with the compaction increasing. When the compaction coefficient reaches 0.95, L∕L 0 reaches its minimum, 0.96; hereafter, L∕L 0 gradually increases, with small changes. Figure 13 also shows how the failure depth of the slope at section "1-1" varies with the compaction coefficient changing after the 3-h-long 100 mm h −1 rainfall. The failure depth gradually decreases with increasing compaction coefficient; when the compaction coefficient is 0.87, the failure depth reaches its maximum of 0.39 m and when the compaction coefficient is 0.95, the damaging depth reaches its minimum of 0.17 m.

The Effect of the Long-Term Strength
To analyze the influence of the long-term strength on the slope stability during raining time, the short-term strength parameters were C s =28.76 kPa and s = 23.36 • , and the long-term strength parameters were C s =11.22 kPa and s = 34.01 • are provided for stability analysis. Other parameters are shown in Table 1, and the rest of the information is the same as in Sect. 4.1. Figure 14 shows the potential failure zones with and without considering the effect of the long-term strength. The curves reveal that if the influence of the long-term strength is neglected, the potential failure zones only develop in the zones from A to D and B to D, but if the effect of the long-term strength is considered, the potential failure zones will run through from the slope toe to the slope top, almost parallel to the slope surface, and they will develop in the zones of 0.3 m deep below the slope surface.

A Brief Discussion on Sensitivity
The slope stability is determined by many factors. The higher the value of L∕L 0 , the more serious the failure; the deeper the failure, the more serious the failure. The results of the analysis indicate that three factors, rainfall intensity, its duration, as well as the long-term strength are the ones most influencing he slope stability. Thus, for a real engineering practice, attention is suggested to be more paid on the influences of these three factors.

Conclusions
In this paper, the constitutive model of two-phases coupling of unsaturated water and air has been established based on the single variable effective stress theory of unsaturated soil. With the consideration of the fluid-solid coupling of surface runoff and underground infiltration, the theoretical calculating model has been numerically verified by the engineering case analysis of rainfall-caused slope failure within the part of K240 + 10 and K240 + 70 on the Dazhou-Chengdu Railway Line. The numerical calculation and analysis have been conducted for each main factor that may significantly influence the slope stability. The following conclusions have been obtained: 1. The consideration of fluid-solid coupling of surface runoff and underground infiltration is necessary for the analysis of slope stability. Only in this way can the process of rainfall affecting the slope, as well as the gradual failure process of slope can be revealed accurately. The proposed coupling model is reliable to analyze such coupling effect for slope engineering under rainfall. 2. In order to decrease the effect of rainfall on the slope stability, the whole slope surface should be protected, the slope ratio of 1:1.55 is suggested, the ratio of K 11 ∕K 22 should be less than 4, the reinforcing depth in the slope should reach about 0.5 m and the compaction coefficient of the embankment slope should be increased as possible; 3. Whatever the slope height and the soil's initial saturation are, the slope stability under rainfall is little affected by them. Under the conditions of the same rainfall totals, when rainfall strength is below 60 mm h −1 the change of the rainfall intensities has a large influence on the slope stability, whereas when the rainfall intensities are larger than this value, the influence of the variation of the rainfall intensities on the slope stability is relatively much less; 4. Putting together the influences of all factors, the rainfall intensity, the duration, as well as the long-term strength of the soil has the most notable influence on the slope stability under raining conditions.