Effect of the Water-Air Coupling on the Stability of Rainfall-Induced Landslides Using a Coupled Infiltration and Hydromechanical Model

The water and air ﬂ ow in soil pores in response to rainfall in ﬁ ltration is an important factor a ﬀ ecting the stability of rainfall-induced landslides. The objective of this paper is to investigate the e ﬀ ect of water-air ﬂ ow on the deformation and stability of rainfall-induced landslides based on a coupled in ﬁ ltration and hydromechanical model. The model is formulated by the water and air ﬂ ow equation, the mechanical equilibrium, and porosity equations. The rainfall in ﬁ ltration is introduced as a ﬂ ux boundary, which is determined according to the comparison of the in ﬁ ltration capacity and rainfall intensity. The numerical model and solution approach are validated by simulating the Liakopoulos drainage test and a rainfall in ﬁ ltration experiment with satisfactory results. Taking the Tanjiawan landslide as a case study, the water and air ﬂ ow in response to rainfall in ﬁ ltration and its e ﬀ ect on deformation and stability are examined. The results show that the pore air is gradually trapped and compressed due to rainwater in ﬁ ltration. The entrapped air has a slowing e ﬀ ect on the rainfall in ﬁ ltration and a pushing-out e ﬀ ect on the front sliding body, which is an important driving force for the evolution of the slope deformation and stability due to rainfall in ﬁ ltration.


Introduction
Rainfall is one of the most important factors to trigger the landslide hazards. Based on previous experimental and numerical studies [1][2][3][4], rainwater infiltration into a soil slope may increase underground water level and soil weight and deteriorate soil mechanical properties, leading to slope instability. However, some landslides are observed even on the condition that only a thin saturated layer formed near the surface due to rainwater infiltration [5]. The effects of the changes in pore water and the degradation of mechanical properties are limited. The rainwater infiltrates into the pores in unsaturated soil that occupied by water and air, and the flow of air has a significant effect on infiltration of water. The subsequent fluid flow and soil deformation in response to rainwater infiltration and its effect on the slope deformation and stability are complicated processes involving coupled infiltration and water-air flow and soil deformation with complex initial and boundary conditions [6].
The rainfall infiltration is one of the most important boundaries affecting the transient state of subsurface water and air flow. Most of the rainwater infiltration and runoff models focus on the shallow water flow on the slope surface and ignore the influence of the subsurface fluid flow [7][8][9]. However, runoffs can still be observed with very low rainfall intensity [10]. Rainfall may not infiltrate into slope through some area of the surface. The movement of rainwater and the redistribution of infiltration may change the hydrodynamic loading path on slopes and affect the acceleration of the slope deformation and instability.
Many rainfall infiltration analysis methods for unsaturated soil slopes consider pore-air pressure to be equal to zero by using Richards' equation [11] to model the subsurface fluid flow [12][13][14][15][16]. These models ignored the coupling with airflow in pores assumed that the airflow is much freer compared with water flow. Many researchers [17][18][19][20][21][22] have shown that the airflow is restricted in the pore space, and it can affect the rainwater infiltration and slope stability. Some results showed the air transport in soil pores has a slowing effect on the movement of the rainwater and is unfavorable to slope stability [23]. In recent years, numerical models for coupled water and air flow addressing the significance of the airflow on rainfall-induced landslides have gained more popularity. Some researches [24] considered a coupled water-air two-phase flow model in combination with a limit equilibrium method to assess the slope stability. Other studies [25][26][27][28][29] applied a coupled water-air-soil three-phase flow-deformation model to simulate the water and air flow response to rainfall and to assess the deformation and stability evolution considering the contributions from both pore-air and pore-water pressures. Theoretically, the coupled approach should produce results that are more realistic. However, the assessment of how the slowing of rainwater infiltration affects the slope stability still remains unresolved.
The objective of this paper is to improve an understanding on the coupled water-air flow in response to rainfall infiltration and its effect on the evolution of deformation and stability of rainfall-induced landslides. A coupled infiltration and water-air flow and soil deformation model was adopted to examine the effect of airflow and rainfall infiltration on the dynamic deformation and stability of the Tanjiawan landslide. The results show the entrapped air due to rainfall infiltration is an important driving force for the slope deformation and stability.

Model Description
2.1. Governing Equations. The governing equations of the soil-water-air coupling system can be classified into four groups: the mechanical equilibrium equation for the soil-water-air mixture (Equation (1)), the mass conservation equation for the solid phase (Equation (2)), the water flow equation (Equation (3)), and the airflow equation (Equation (4)).
where σ is the total stress, S r is the water saturation, n is the porosity, b is the body force, k is the intrinsic permeability, g is the gravitational acceleration vector, and t is the time. ρ s , ρ w , and ρ a are the densities of solid, water, and air, respectively. p w and p a are the pressure of pore air and pore water, respectively. μ w and μ a are the viscosity of water and air, respectively. k r w and k r a are the relative permeability of water and air, respectively. Q w and Q a are the source of water and air, respectively.

Constitutive
Relations. The solution of the coupled flowdeformation model governed by Equations (1)-(4) can be obtained numerically with appropriate constitutive relations listed in Table 1. 2.3. Factor of Safety. Factor of safety F s is a wide-used indicator to access the slope stability. In this paper, the Bishop integral method is adopted to obtain the factor of safety with a given sliding surface. The formula is given as the ratio of

Constitutive relations Expression Description
Effective stress σ′ = σ − α B 1 − S r ð Þp a I − α B S r p w I σ′ is a bishop-type skeleton stress tensor, α B is the Biot coefficient, and I is the identity tensor Soil-water characteristic curve (SWCC) [30] S e = 1 + p c /p 0 ð Þ 1/1−λ −λ p 0 and λ are the model parameters S e = S r − S rr /1 − S rr Relative permeability of water [31] k Relative permeability of air [32] Intrinsic permeability [25] k = k 0 n/n 0 ð Þ 3 1 − n 0 /1 − n ð Þ 2 n 0 is the initial porosity, k 0 is the initial intrinsic permeability Density of water [29] ρ w = K w ρ T /K w − p w + 101325:0 ρ T = 998:2065 kg/m 3 , K w = 2:15 × 10 9 Pa Density of air ρ a = 7:3756 × 10 −6 p a 2 Geofluids resistance forces to sliding forces along the sliding surface, written as where c i is the cohesion, φ i is the fiction angle, τ n and σ n ′ are the shear and normal stress on the sliding surface, respectively, and l 0 and l 1 are the starting and ending points of the sliding surface, respectively.

Numerical Implementation
2.4.1. Rainfall Infiltration Boundary. The rainwater runoff is calculated based on the assumption that surface runoff equals to the difference between the rainfall intensity and the infiltration capacity at a specific space and time, yielding where v r is the rainwater runoff velocity, I n is the rainfall intensity, and v in is the maximum inflow rate (outflow is negative). On the rainfall infiltration boundary, the infiltration water flow can be obtained as The maximum infiltration rate v in is obtained from the subsurface flow velocity of the water-air mixture v m , given as The water saturation and pore-air pressure are set as the basic variables, and the solution of Equations (3) and (4) is replaced by that of Equations (9) and (4). The iterative method for solving the two-phase flow equations was proven to have good performance of convergence and numerical stability [33]. 1. For the time marching (t n , t n+1 ], let {n n , S r,n , p a,n , u n } be available. 2. Set k =0 and initialize the guess n 0 n+1 = n n and S 0 r,n+1 = S r,n . 3. Find {S k+1 r,n+1 , p k+1 a,n+1 } by solving water and air flow equations (10) and calculate n+1 by solving Eq. (11). 5. Find by solving Eq. (12). 6. Check the convergence criteria: ELSE k = k +1, GO TO 3. 7. Check the time-marching procedure: IF t n + Δt < t max , THEN GO TO 8; ELSE GO TO 9. 8. n = n +1, GO TO 1. 9. Output the results.
Algorithm 1: The iteration procedure for solving the displacement-saturation-pressure-porosity system.

Geofluids
The coupled flow-deformation equations are spatially discretized by the Galerkin-type FEM method. The discre-tized forms of the flow, equilibrium, and porosity equations can be expressed as

Geofluids
The expressions for the coefficient matrices C and K and the vector f are presented in Appendix A.
The time discretization is carried out based on finite difference method. At each time step, the variations of the basic variables are assumed to be linear. where The time derivatives are expressed as follows: 2.4.3. Iterative Procedure at the Discretized Level. The iterative procedure of the displacement-saturation-pressure-porosity system is carried out in sequential coupling scheme, shown in Algorithm 1. The main reason for choosing this sequential iteration is the difference of modulus of the elasticity and permeability, which may in range by several orders of magnitude. The iterative processes will continue until the values of the displacement vector, pore-air pressure, water saturation, and porosity change little in a given time step. The updated variables are then adopted as input data in the iteration at next time step.

Validation
3.1. The Liakopoulos Drainage Test. The numerical model was validated against the experimental results from Liakopoulos drainage test [34] on Del Monte sand. The test has been adopted as a standard validation for the flowdeformation problem [19,35,36]. The test can be simplified as an axisymmetric problem, shown in Figure 1(a). The soil column was discretized into 40 meshes with 63 nodes. The initial and boundary conditions and material parameters are obtained from Hu et al. [27]. The initial time step size was 10 s, and the size was updated automatically according to the convergence of the calculation. The tolerances were specified as Tol 1 = 0:01, Tol 2 = 0:001, Tol 3 = 0:001, and Tol 4 = 0:001. The simulated water outflow rate is compared with the experimental data, shown in Figure 1(b). It can be observed that the infiltration rate agrees well with the experimental results. The pore-water pressure is monitored in the simulation and plotted in Figure 1(c). The figure shows the pore-water pressure calculated from the numerical simulation decreases a little faster than the experimental data in the early 10 min, and the predicted and experimental results are in good agreement afterward. The settlement deformation at 120 min is 1.58 mm shown in Figure 1(d), which coincides with the results of previous studies [28]. The good agreements show the feasibility and applicability of the numerical

A Laboratory Test
Reproducing Rain Infiltration into a Soil Column. A laboratory test was performed to reproduce the coupled air-water flow processes induced by rainfall infiltration in this section (Figure 2(a)). The experimental results were used to investigate the performance of the numerical formulation. The soil column was 0.9 m in height and 0:5 m × 0:5 m in horizontal cross-section, shown in   7 Geofluids Figure 2(b). The rainwater infiltrated into the soil column from the top, and the runoff was collected and weighed to calculate the water infiltration. Six pressure sensors and four soil-moisture sensors were laid out along the depth. The gravel was placed at the bottom to observe the water outflow, and one pressure sensor was used to monitor the outlet pressure.
The rain infiltration into the soil column can be simplified as a plane problem, and the FEM mesh is shown in Figure 2(c). An initial saturation of 40% and porosity of 0.36 have been measured before the test. The sidewalls are sealed for both water and air flow and fixed with zero normal displacement. The rainfall infiltration and the atmosphere pressure are applied on the top surface. The pressure    Geofluids on the bottom is equal to the measured outlet pressure. The soil sample is loamy sand from the Three Gorges region. The model parameters could be calibrated by experimental tests on the soil samples, given as p 0 = 11:3 kPa, λ = 0:25, and S rr = 0:1. Figure 2(d) shows the rainwater infiltration into the soil column versus time. It is observed that the simulated cumulative water infiltration agrees well with the experimental results, and the infiltration rate decreases to a low value in the first few minutes. Figure 2(e) plots the comparison of the simulated water saturation with the experimental results. The good agreements illustrate that the coupled model can describe the overall behavior of the subsurface flow in the rainfall infiltration processes. The saturated zone is formed at the top and moves slowly toward the bottom unsaturated zone. The rainfall infiltration process is greatly retarded because the lateral and bottom surfaces are sealed to prevent air and water flow. The wetting front reaches a depth of 20 cm after a long time of 80 h. The rainwater infiltration exhibits a "trapping" effect on airflow to increase the poreair pressure. The trapped air delays the propagation of the wetting front.

Geological and Engineering
Background. The Tanjiawan landslide is located on the right bank of the Zhaxi River in the Three Gorges Reservoir. It is approximately 43 km away from the Three Gorges Dam. Figure 3 shows the layout of the Tanjiawan landslide. The landslide is nearly in west-east direction. Both sides of the landslide are gullies, the back is located at the foot of the bedrock, and the front is a high scarp with a deep valley, reaching the Zhaxi River. Figure 4 Figure 6 shows the boundary conditions for the numerical simulation. The upper and right boundaries are subject to reservoir water and specified as the Dirichlet boundary with S r = 1:0 and p a = p w below the water level and p a = p atm above the water level. The rainfall infiltration with monitored rainfall intensity is applied on the upper surface. The bottom boundary is considered rigid and sealed for water and air flow. The right and left boundaries are fixed with zero normal displacement. The initial distribution of the water saturation and pore-air pressure are obtained by simulating a two-year water and air flow under the preceding boundary from full saturation. The initial geo-stress is determined by applying the gravity loads and fluid pressure based on mechanical equilibrium method. Figure 7(a) illustrates the typical flow vector of a water-air mixture at Feb 2008. The slope surface can be macroscopically divided into two parts, the inflow zone on the back and the outflow zone on the front, according to the flow direction of water-air mixture. Figure 7(b) displays the variation in the inflow and outflow flux of the rainwater in this period. As the water saturation near the upper surface increases, the inflow and outflow fluxes generally decrease. In the rainy season (Apr~Jul 2007, Apr~Oct 2008, and Mar~May 2008), the decreasing rate becomes much slower, and even slight increases are observed. The inflow and outflow fluxes depend on the rainfall intensity and the initial water distribution near the slope surface. that the pore air is gradually trapped and compressed, and its escape is restricted. A large zone of low saturation is generated and maintained in such a long period. The wetting front moves more slowly as the zone of low saturation becomes smaller, indicating that the trapped pore air has a slowing effect on the movement of the wetting front. After April 2008, a high saturation zone is formed near the sliding surface. At the hole V2 in the outflow zone where water-air mixture flows outwards, the water saturation is larger at the deeper location. Figure 9 compares the pore-air pressure evolutions in holes V1 and V2 at different depths. The pore-air pressure near the slope surface is approximately zero. The pressure increases as the depth increases. At hole V1 in the inflow  11 Geofluids zone, the pore-air pressure at depths of 10 m and 22 m increases in the first and second rainy seasons, and the peak of the pore-air pressure occurs at about 40 d later than the end of the rainy season. This indicates that the increase of pore-air pressure is related to the rainfall, and the wetting front moves very slowly and continues extending even after the end of the rainy season. The air in soil pores has difficulty in escaping from the surface in the inflow zone. At hole V2 in the outflow zone, the value of pore-air pressure is almost synchronized with that of rainfall intensity. Almost no delay effect is observed. Figure 10 shows the variation of the ground displacement and a comparison with monitored data at monitoring points ZG331, ZG332, and ZG333. The simulated results agree well with the monitored data with acceptable differences. The entrapped air in the slope flows out below the center of the slope surface, where the saturation is high. The entrapped air has a direct pushing-out effect on the center of the sliding body, where the maximum displacement is observed (ZG331), with a magnitude of approximately 60 mm. The variation in deformation of the Tanjiawan landslide is closely associated with the rainfall variation, and the acceleration stage of deformation occurred in the rainy seasons of 2007 and 2008. An important reason may be that the sliding body is saturated and the pore air is compressed mainly in the first rainy season, and the pore-air pressure increases little afterward. Figure 11 plots the evolution of the factor of safety and the corresponding rainfall intensity in the three-year period. It is observed that the factor of safety rises slightly at the beginning of the rainy season due to an increase of the slope weight as a result of rainwater inflow. With continuous rainfall infiltration into the slope, the pore-air pressure increases in the first rainy season, where the factor of safety decreases significantly. The pore-air pressure occurs a second growth, and the factor of safety decreases fast after the end of a heavy rain month in the second rainy season. The change of factor of safety is relatively stable, and the pore-air pressure presents no further increase in the third rainy season. It illustrates that the flow-deformation process considering pore airflow lasts for a longer time, and the increase of pore-air pressure is unfavorable to the stability of rainfall-induced landslides.

Water and Air Flow Response to Rainfall Infiltration.
In the process of rainfall infiltration, the slope surface is the main region for water and air inflow and outflow. According to the flow direction of water-air mixture, the slope surface can be divided into an inflow zone located in the upper part of the slope and an outflow zone located in the lower part. 12 Geofluids air mixture is less affected by runoff water pressure in the rainfall infiltration process, and it can be used for identifying whether rainfall infiltration produces runoff. The boundary conditions affect the airflow in pores and infiltration of rainfall. A layer of high saturation zone is formed near the slope surface during the infiltration of rainwater. The pore air is hard to escape from the slope surface due to the low hydraulic conductivity of pore air in the high saturation zone. The pore air is compressed and entrapped by the advance of wetting zone and delays the movement of the wetting front, leading to an increase of pore-air pressure. The entrapment of pore air is a path-dependent process which depends on the formation of the air sealing of slope surface. The variation of real-time rainfall conditions and the hydraulic properties of unsaturated soil can provide different hydraulic path on the slope and produce diversity of rainwater infiltration and the increase of poreair pressure.

Airflow
Effect on the Slope Stability. The compressed pore air can rapidly make the top pore-water pressure apply on the middle and front of slope via increasing the pore-air pressure. The pore-air pressure on the front of slope is close to the top pore-water pressure with less pressure loss because the viscosity of air is much less than that of liquid water. The increased pore-air pressure provides and maintains a pushing-out force on the front sliding body for a long-term timescale, which is unfavorable to the stability of the slope, shown in Figure 12. Although the air pressure is as small as only hundreds of pascals, the direction of the pore-air pressure force is directly pushing the sliding body out, and the cumulative force of the pore-air pressure is 13 Geofluids significant and nonnegligible as it acts on a large surface area with more than tens of thousands of square meters. The loading effect of airflow relies on the compression of pore air that may be dominant in some shallow landslides with long duration rainfall.
The transient increase of the slope stability at the beginning of rainfall may be due to a slight increase of the slope weight and the low pore-air pressure that has not yet risen. Rainwater infiltration forms a thin transient high saturation area near the surface of the slope. In the saturated area of the slope, rainwater will form static pore-water pressure and act on the solid particle framework in the form of seepage volume force. The effect of static seepage force in the high saturation area on slope stability depends on the magnitude of pore-water pressure, which is small in general circumstance.
Pore-water flow forms a high saturation zone near the sliding surface with accumulation of time, especially for deposit slopes where the permeability of sliding bed is far less than that of sliding body. The higher the soil saturation, the smaller the matrix suction and the lower the shear strength. The degradation of the material properties of the sliding body may cause the redistribution of slope stress and trigger slope deformation and instability. The rainfall infiltration may generally act as a generalized loading and a degradation factor of material properties on the slope deformation and stability.
This study provides an understanding of the macroscopic behavior of the pore fluid flow during rainfall infiltration and an insight into the airflow effect on the dynamic behavior of slope stability. It should be noted that failure of a slope induced by rainfall may be a result of many factors, such as heterogeneity and uncertainty in soil properties [37,38], surface erosion [39], and preferential flow [40,41]. Further researches on the modeling and effect of these factors can be performed to improve the understanding of rainfall-induced landslides.

Conclusions
In this paper, the evolution of deformation and stability of landslides in response to rainfall infiltration is investigated using a numerical model of coupled infiltration and hydromechanical processes. The results show that the slope surface can be macroscopically divided into the inflow zone on the back and the outflow zone on the front according to the flow direction of water-air mixture. The rainwater macroscopically infiltrates into the slope through the inflow zone and compresses the air within the pores. The pore-air pressure then increases as the advance of wetting front. The entrapped air reduces the infiltration rate and provides a pushing-out force on the sliding body. The middle and front of the sliding body are the direct parts bearing the increased air pressure, and the deformation increases with an increase of pore-air pressure. In the processes of rainfall infiltration, the factor of safety decreases rapidly even after the end of the rainy season, illustrating that the coupled water and air flow and soil deformation in response to rainfall infiltration is a very slow process and could influence the timescale of the stability evolution of rainfall-induced landslides.
where B is the geometric matrix.
(4) For the porosity equation

Data Availability
The data used to support the findings of this study are included within the article.

Conflicts of Interest
The authors declare no competing interests.