Study on Initiation and Mechanism of Seepage-induced Slope Failure in Jianghan Plain, Central China

The problem of seepage erosion is widespread in occurrence, but is rarely recognized as one of landslide causes. This paper describes and explains one event of slope erosion in Jianghan plain and is focused on the role of seepage erosion. The laboratory test is performed to reconstruct the hydraulic processes leading to slope failure with the application of different hydraulic loads. The observation and measurement reveal that slope alteration experienced three stages under seepage, and hydraulic gradient played a key role during seepage erosion phase. Using these results, a CFD-DEM model is developed to investigate how the seepage-induced soil particle migration would affect the slope evolution. The geotechnical properties control the failure processes and cause a significant change to the toe erosion and loosened zone. Although the instability of slope is mainly attributed to seepage gradient force, the occurrence of seepage erosion affects the failure mode. This work provides new guidance for hydrodynamic landslide prevention.


Introduction
Seepage erosion has been recognized as a problem of global significance. Many cases of instabilities induced by seepage erosion have occurred in different geological settings and materials (Hagerty et al. 1981;Abam 1993).The dislodgement of soil particles due to seepage-induced internal erosion is suggested to potentially play a prominent role in the failures of embankments and dams (Fell et al. 2000).The most commonly adopted terms are piping (Glynn et al. 2012;Wilson et al. 2013) or backward erosion (Robbins et al. 2019;Semmens et al. 2021).In contrast to the well understood failure modes of levee/dam breaching, the processes of seepage-induced landslide are unknown.
Old and recent fluvial terraces or stream banks are locally susceptible to the occurrence of such processes due to sedimentological, morphological and hydrological reasons. (Jones 1987;Wilson et al. 2007).During the early summer of 2016, one incident about seepage erosion in slope toe was observed along a fluvial terrace in Jianghan Plain, China (Fig. 1). The shallow sediments in JHP are mainly composed of an upper silty-clay layer functioning as an aquitard and an underlying sand functioning as an aquifer (Du et al. 2018). Large natural or artificial ponds can form behind slope or in topographic depression. In comparison with the previous rapid slope collapses in other plain areas (Swanson et al.1989;Hagerty 1991;Crosta and di Prisco 1999;Fox and Wilson 2010), the 2016 event initiated at a lower gradient, was relatively slow, and generated uneven settlements, causing the tensile fractures development on the upper slope.
Subsurface flow can cause hillslope instability through three different but interrelated mechanisms: increased soilwater pressure (Rinaldi et al. 2004), seepage gradient forces (Iverson and Major 1986), and seepage erosion ).More precise expression is a combination of different mechanisms that may or may not coexist. The final consequence of movements of eroded material by the concentrated flow to its downstream end where it intersected the ground surface referred to as a sand boil (Pratama et al. 2020). It is unclear if sand boils were related to failure, or simply very localized surface features in JHP, as the visual impact of features produced by seepage erosion may be much less than the impression created by massive slumps. Due to the lack of observation and record of erosion after landslide outbreak and the limited understanding of the specific process of advancing geotechnical instability in previous research, it is difficult in differentiating the role of this mechanism.
Previous studies determined that erosion processes in slope involve more complex interactions between intrinsic and external factors. With the increase of pore pressure inside the slope, the effective stress acting on a single soil particle decreases, and the resistance to seepage pressure and hydrodynamic force increases the erodibility of the slope. When the slope could resist the seepage gradient forces and the seepage is concentrated, particle mobilization initiates (Chu-Agor 2008a).Different seepage responses may occur: one is heaving which involves the upward movement of a relatively impervious surface layer resulting from subsurface water pressures in excess of its weight, and the other is sand boiling which involves the movement of subsurface sand to the surface by flowing water (Wolff 2002).However, the exact point at which a macro-pore becomes a pipe is elusive and variable, depending particularly on soil properties (e.g., soil composition, permeability) and local geologic conditions (Jones 2010;Nieber and Sidle 2010;Midgley et al. 2013). Bonelli (2018) found that the lower cohesion materials beneath the blanket are vulnerable to being eroded by the concentrated seepage forces.
In addition to the intrinsic factors, hydraulic condition can be considered as short-term trigger in facilitating landslide. Hydraulic gradient (i), a dimensionless ratio between the difference in hydraulic head and distance between two points, shows a strong correlation with the occurrence of erosion (Robbins et al. 2020). When there are ponds at the uphill boundary and the slope is covered by a thin layer of cohesive soil, vertical infiltration into the soil matrix is ignorable (Shen et al. 2020).As the water level rises, a large hydrostatic imbalance is formed under the clay confining, providing an opportunity for lateral subsurface flow (McDaniel et al. 2008;Wilson et al. 2019).To correct this imbalance, water seeps from the high-pressure water-side to the lower-pressure land-side through the soil matrix. In this case, if the fluid strongly affects the interparticle frictional strength, particle entrainment may occur (Nouwakpo et al. 2010).At boiling point on the ground, the hydraulic gradient quickly dissipates, depositing eroded material around the mouth of the pipe. This article investigated the evolution of seepage erosion to explain the observed field phenomena through some small-scale physical models. The entire morphing slope and loading conditions has been quantified using sensors, to provide an understanding of specific failure criterion. A simulation of grain-scale fluid-grain interactions through a coupled CFD-DEM numerical approach is proposed to model the seepage-induced fine particle migration within continuous grading soils, combined with experimental measurements to test the influence of the geotechnical properties and hydraulic gradients on initiation, continuation and progression of slope instabilities. The results could reveal unexpected links with seepage-induced internal erosion and hydrodynamic landslide, and assist in more accurate evaluation. This research is placed in the context of Jianghan Plain but is generally applicable to other sites as well because hydrogeological conditions similar to Jianghan Plain are quite common globally.

Field Site Description and Sampling
The research site is in the hinterland of Jianghan Plain in Hubei Province, China (Fig. 1). Song (2020) provides an overview of geographic map in location of Jianghan Plain. Total station equipment and exploration drill holes data from previous studies was used to conduct a site survey in June, 2016. Figure 2 shows a simplified stratigraphic section through the topography. The topography is relatively gentle, for the slope angle of 15°-20° and the slope height of 17m. The overlain clay with 1-3 m thick makes the aquifuge. The thickness of confined aquifer is about 10-20m, which mainly consist of sand. The highest water level of the pond was almost level with the top of the slope due to recent precipitation events, leading the high local hydraulic head, which was prone to developing subsurface erosion. As observed in the study area , the concentrated seepage carried eroded material to the ground surface at slope toe where the maximum flow can reach 5L/s and the mature trees nearby tilt to a certain extent. A large crevice was revealed at the back edge of the slope, with a maximum width of 0.50m and a maximum depth of 0.45m.
A series of laboratory tests were conducted to analyze the physical properties of sampling soil which was collected from depth of 3.0-5.0m in slope. As shown in Table 1, the specific gravity obtained by pyknometer method is 2.65, the permeability coefficient obtained by the conventional permeability test instrument is 1.28×10 -5 , and the internal friction angle obtained by triaxial compression test is 27. Figure 3a provides a description of the grain size distribution of soil guided by sieving test. The nonuniformity coefficient ( ) is 6.667> 5, and the curvature coefficient ( ) is 1.667∈ [1,3], which reveal that the particle size distribution is nonuniform and gradation is continuous.
Internal erosion includes four distinct mechanisms: regressive erosion, erosion along a concentrated leak, interfacial erosion, and suffusion. When fine particle size is smaller than the pore size, the topsoil particles affected by the seepage force may break away at any time, and the direction of destruction determines the development direction of the concentrated seepage channel. For the continuous graded soils, empirical formulas were used to measure the susceptibility of suffusion (i.e., the internal stability of soils).The average pore diameter 0 is given by where ( = 20) is the particle diameter corresponding to % of the soil mass. The calculated average pore diameter ( 0 ) of 0.035 mm is a bit greater than the maximum flowable particle diameter ( 5 ) of 0.030 mm. Sandy substrates contain little cohesive silt particles and coarse sandy gravel, however the cohesionless processes still dominate (Baas et al. 2013). Very fine particles have chance to move in the connected pore channels under a certain hydraulic action, and dynamically similar to silt-free sand condition. The soil particles account for a higher proportion and will be drove to roll up by water flow, and the apparent bulk density of water will greatly increase as well in sand boil zone. Considering the constraint of the surrounding media on the soil particles and the slope conditions throughout the sand boil, the critical hydraulic gradient of non-cohesive soil can be estimated to be where is the material internal friction angle and is the slope angle. The value of critical hydraulic gradient ( ) is calculated to be 0.29 from Equation 2, which is lower than the measured hydraulic gradient at slope toe of 0.31. The progression of erosion is controlled by the local hydraulic gradient. These results indicate that erosion along a concentrated leak may occur, but how internal erosion further induces the instability and failure of the slope remains to be studied by model test.

Slope alteration experienced three stages under seepage
The development of distress features evolve in stages. To achieve more comprehensive insights into the characteristics and controls of failure processes, a small-scale model was built in this experiment with the similarity ratio (n) of 50. By controlling the inlet flow and keeping the water level constant, the lateral seepage in slope to propagate until the soil is fully saturated was simulated. Figure 4 shows the experimental configurations. The seepage simulation system consists of a water feed device with controller, water reservoir and porous disc. The model chamber was open at the toe of the slope, so that seepage into the slope could be drained out. A miniaturized slope of a 18°gradient was formed inside the model chamber. The dual soil slope consisted of a clay top layer of 0.06 m and sand layer at the bottom. Soil materials were collected from the site. Using the equivalent substitution method, the test material parameters were calibrated to ensure same similarity constant of density, internal friction angle and moisture content. The correlation coefficient of key physical quantities in the model test is shown in Table  2.The particle size curve after replacement is presented in Fig. 3b. The purpose of model test was to observe slope failure and to investigate the hydrologic processes in slope. Six earth pressure sensors (diameter 10mm, range 0.2MPa, accuracy 0.5%FS) were placed in two rows at specified depths to record the pressure between soil layers. Displacement monitors were mainly composed of steel wire (diameter 1mm) and bar with graduated mark (range 150mm, accuracy 0.3%FS), distributed at 20cmm, 50cm, 80cm from the slope toe. Three groundwater monitoring tubes were located at the front, middle and rear of the slope respectively to measure water level (Fig. 5).
In this test case, a series of experiments with gradually increasing hydraulic load was performed, referred to as gradual. The initial hydraulic head difference(△H) of 0 cm was increased in steps of 10 cm, and kept constant (equilibrium) for 3h.If no erosion was observed, this process was repeated again until the hydraulic head was increased up to a maximum, and the test was stopped. In practice however, the water level of ponds can vary quickly due to rainfall discharge effect. As a complement, an experiment was performed with a quick increase in hydraulic load to the value of 30 cm. Figure 6 shows the front view of the failure features at the gradual test. There are three stages during the failure processes of slope.
(a) Wetting. For the hydraulic loads HL=10cm, the groundwater gradually migrated laterally up to the slope during the experiment. After 30min, seepage flow initiated and the slope toe saturated first. At the surface of the slope, there appeared some small linear cracks in clay layer (Fig. 6b-c). The whole slope was basically stable without obvious signs of sliding deformation. (b) Seepage erosion. For the hydraulic loads HL=20cm, the seepage zone extended to the surface of the slope. Due to the poor permeability of the overlying clay layer, the water in the sand layer was difficult to discharge, and the accumulated seepage water pressure could not be dissipated, resulting in the heave at slope toe and some small tensile cracks in the middle of the slope body. The initial sign of this stage is heave induced radial cracking (Fig. 6e). In the continuous seepage process, radial fractures broadened. For the hydraulic loads HL=30cm, a large number of soil particles were carried by water and lost. A sand boil formed at the foot of the slope eventually (Fig. 6f). (c) Major failure. It was interesting to observe the formation of tension cracks similar to those in previous documented laboratory tests of the lysimeter (Fox et al.2006Chu-Agor et al. 2008b),as shown in Fig.  6(h).The beginning of this stage was marked by erosion at slope toe (Fig. 6g). In this stage, the main failure forms were as follows: erosion at slope toe → tension crack→retrogressive toe sliding followed by movements with relative motion between soil. The tension crack connected with a potential shear zone, and the soil mass at the front edge continued slow sliding (Fig. 6h), finally forming a landslide.
Our data reveal that the initiation of seepage erosion depends greatly on the hydraulic loads. The mechanical behaviors of slope were similar in both quick and gradual loading experiments with same hydraulic head difference △H=30cm, as shown in Fig. 7. At equilibrium under steady flow (C2-D1 phase), water is under pressure due to the weight of the upgradient water and the confinement of the water under impermeable layer (Fig. 7b). A body of soil near the slope toe would begin to move, due to the seepage forces overcoming the weight of gravity (Fig. 7d). The movement of the clay layer will lead to heave and cracking through which concentrated flow. This result can be seen clearly in Fig. 6e. For sand boiling to occur, the hydraulic gradient of the seepage must be substantial enough. The D2-D3 phase indicated the head loss required to transport sand vertically through a sand boil (Fig. 7b).The sampled sand ejected from the sand boil consisted of fine-grained sand with of 0.114 mm. The sand obtained from the throat was considerably coarser and more poorly graded than the material deposited, with of 0.319 mm. In addition to the hydrodynamic head loss due to fluid flow, the presence of coarser material in the sand boil throat may result in higher head loss due to the drag force exerted. A second observation is that the location of sand boil is non-axisymmetric. This can be attributed to the spatial variability of soil properties in sand.

Basic modeling principles
Seepage flow phenomena of the miniaturized landslide at the particle-scale were simulated using a CFD-DEM model. The numerical modeling aimed to assist in understanding how the seepage-induced soil particle migration in continuous graded soils would affect the slope evolution.
In the model, an individual grain is approximated by an idealized shape, e.g., a sphere. The movement of the solid granular medium (particles) is traced by the discrete element method (DEM) using Newton's second law (Cundall and Strack 1979).The governing equations for translational and rotational motion of the particle can be expressed as where and are the translational and angular velocities of a particle respectively, is the particle mass, is the gravitational force, is the contact force on particle , is the particle-fluid interaction force on particle from the fluid phase, is the moment acting on particle , is the moment of inertia. With the velocity and position of a specific particle, a contact is created as two particles overlap, and the contact points and overlapping distance between the two particles are defined. Figure 8 shows the mechanical contact model between contacting bodies in sand and clay soil model. A contact force between two particles consists of a normal force and a tangential force. The relationship between displacements and contact forces in non-cohesive granular material was implemented by the linear spring-dashpot model (Zou et al. 2013).The parallel bond model which can transfer forces and also moments was applied to analyze each particle in cohesive material (Shen et al. 2016).
The fluid calculation was solved by the finite volume method (Meakin et al. 2008) and realized by the seepage module Computational Fluid Dynamics (CFD) in Particle Flow Code (PFC) software (Desai and Holloway 1971). In the un-resolved model, the flow around each particle is volume-averaged within a local domain through solving the locally averaged Navier-Stokes equations (Anderson and Jackson 1967).For incompressible fluid with constant density, the conservation of mass and momentum equations can be written as follows: where ∇ •is the divergence operator, ̅ is the averaged fluid velocity vector, ̅ is the averaged fluid pressure, is porosity, is the fluid density, is the gravitational acceleration vector, ̅ is the averaged interaction force vector, and is the fluid stress tensor = [∇̅ + (∇ ̅) − 2 3 (∇̅) ]for Newtonian fluids, is the identity tensor. The fixed coarse grid fluid scheme (Pirnia et al.2020) was used to solve the seepage flow problem in equations (5) and (6).The averaged interaction force imposed by particles to the fluid is given by where ∆ is the mesh volume of a CFD cell, is the particle number inside the cell, and ̅ is the interaction force imposed by the fluid to the th particle in the cell.
The interaction force between a particle and a fluid consists of two parts. One is hydrostatic force which contains the fluid pressure gradient force, or the buoyant force. And the other is hydrodynamic force which includes the drag, viscous, lift, Basset and virtual mass forces and so on (Crowe et al. 1998). Since the drag force is dominant in a low Reynolds number flow, the other terms in hydrodynamic force are neglected. Therefore, the dominant interaction forces between particles and fluid on an individual particle include drag force, pressure gradient force, and viscous shear force: where is the volume of particles. The drag force imposed on a particle was treated using well-established empirical models (Zhu et al. 2007) based on the correlation of Di Felice (1994).
At each time step, DEM cycle will give information, such as the current positions and velocities of individual particles, to CFD cycle for the evaluation of porosity and fluid drag force in a computational cell (Tsuji et al. 1993). The particle-fluid coupling simulation is performed to apply a seepage (drag) force on each particle at each iteration, and explicitly updates the particle dynamics through Newton equations (Tao et al. 2017). The porosity for different fluid cells changes accordingly, and then affects the fluid parameters of the grid.

Numerical modeling
In the basic case simulation, we prescribed the domain geometry according to the geologic section (Fig. 9a). The length, width and height of the numerical model are set to 100 cm, 1 cm and 35 cm respectively. Circles with numbers represent typical location for measuring displacement and stress of the specific particles in the central YZplane at Y=0.5cm.The slope model is composed of 25377 particles with different sizes where generated particles are randomly distributed. Table 3 lists the parameters of particle properties based on the typical values suggested by Hosn et al. (2018), which were adjusted to be basically consistent with macroscopic properties. To balance computational efficiency and accuracy, the fluid is discretized into a two-dimensional structured mesh of triangular elements (Fig. 9b). For more details on the technique of coupling two-dimensional CFD with three-dimensional DEM, the readers can refer to Xu and Yu (1998); Zhou et al. (2010). In this study, the total number of CFD cells is 360 and the time step is 5×10-4 s.
To validate that the numerical model gives physically realistic results, hydraulic prerequisites with different constant boundary water levels were simulated for triggering landslides. Similar to previously mentioned experiments, the hydraulic pressure difference (△P) is created in the domain along the x-direction to induce a rightward seepage flow. The fluid calculation time was used to simulate the erosion time. Figure 10 compares the experimental measured and simulated displacement in x-direction for the specific location near the slope toe when HL=30cm. In the numerical simulations, the final displacement is 3.89cm which is close to that in the experiments ( =3.65cm, 3.50cm). Initially, the horizontal movement increases over time, and then the curve show a sudden deceleration in the increasing displacement because the sand boiling happens. The influence distance between the sand boil and the monitoring point leads to a slight discrepancy in the final displacement between the numerical model and the physical experiment.
Figure shows the vertical displacement field of particles in numerical simulations with different water boundary. It can be easily seen a few of eroded particles (Fig. 11b) and small cavities on the slope surface (Fig. 11c) in the regions of dense red particles. As the hydraulic gradient increased, the heave upwards zone at slope toe and loosened zone moving downwards at rear of slope extend, because the soil mass rearrange in response to the increased seepage force. These results indicate that the model gives similar behavior to observed sliding deformation in experiments.

Suffusion enhance macro-void thereby decreasing slope stability
To test the effects of the internal erosion on initiation of failure, the parameters of the particle size distribution were varied within the simulation domain. Three analysis cases were designed using the random distribution particle generation method, which are listed in Table 4. The particles for Case 1 were small and well sorted, and the grains for Case 3 were relatively coarse but the sorting was poor, which are consist with the characteristic of lacustrine and alluvial sediments respectively in the Jianghan plain region (Gan et al. 2018).
The potential of fine-grained particles migration through the pores of the coarse skeleton differed from different particle size distributions. From case A to case C, the number of eroded fine soil particles increased as the coarser skeletal grains increased, because the coarser grains contain more void space that allows finer particles to continue moving. When enough material are removed, small cavities which begin along weakness in the ground grow into pipes to promote migration of particles ( Fig.12b-c), which played a key role in whether the slope could make the transition from seepage erosion phase to the major failure phase. In contrast to other cases, the heave upwards zone at slope toe and tension crack formation never observed in case A (Fig.12a). Figure 13a-b show the development of stress in the slopes at monitoring points. For relatively uniform particles (Case A), the channel of particle movement was limited, and the intruding movement was poor. In this case, the high stress in x-direction kept the relatively stable state, which indicates the eroded particles became clogged in a small region. There was no significant particle movement or erosion happened. When water enters the slope, the infiltration process adds weight to the soil block, increasing the driving forces. This tendency can be seen clearly where the stress in x-direction rapidly increased first. The permeability coefficient of the overburden clay is more than two orders of magnitude smaller than that of the sand, thus the water head difference basically acts on the bottom of the overburden clay during high water event (HL=30cm).Seepage gradient forces played a bigger role especially at the toe part of the slope. Water seepage could promote the upward movement of particles, leading the particles more easily in suspension state. Thus the stress in z-direction maintained at a certain small value close to 0 and the stress in x-direction decreased. The later fluctuation of the stress for Case B and Case C means that the complex interaction between the seepage water and the soil forces determines the failure development, which is a nonlinear dynamic process. As a result of an imbalance in seepage-induced uplift forces, concentrated seepage passes through defects in the clay layer.
As can be seen from Fig.13 c-d, changes in the internal erosion ratios of soil differ for different cases. During seepage erosion phase for Case B and Case C, particles were carried out by water flow (Fig.13d) from the uplifted clay overburden (Fig.13c). Thereafter, the mass resumed movement forward at an almost constant velocity. The vertical movement of particles for Case C firstly reached the maximum (i.e., the seepage channel had penetrated) when t=0.65s, whereas for Case B, t=1.05s. It can be further concluded that a wider size range provides greater potential for fine particles to move. The mobilization and loss of particles by seepage had a marked effect on continuous flow sliding and undercutting. Due to the continuous sliding of the front slope, the uneven settlement crack expanded into tension crack and became the potential slip surface. Furthermore, the influence of grain size is an important variable. A wider size range means more coarse grains. The migration and loss of grains with larger diameter will lead to deeper settlement of the slope. In this way, the slope slide is easier to form.

Conclusions
The discovered mechanism of both slope erosion and geotechnical instability, particularly the role and importance of seepage erosion, provided critical insights into control and evaluation of hydrodynamic landslides. Understanding the factors for the initiation of seepage erosion and triggering of slope failure is an important step in recognizing and preventing future disasters in similar settings. The laboratory experiment and numerical modeling presented herein demonstrate that (a) The failure processes of slope can be divided into three stages, that is, wetting, seepage erosion and major failure. In practice, indicators of the threatened slope, such as sand boils and tension cracks, are easy to be observed and recorded. The seepage erosion phase often lasted several days before reach a limiting condition, allowing time to draw the reservoir water down to prevent the evolution of next stage. Besides, research of major failure phase shows that local instabilities develop rapidly, thus remedial works need to carry out in time. (b) The instability of hillslopes is generally attributed to seepage gradient force, but seepage erosion also affects the development processes of slope deformation that hydraulic gradient and erosion time originally dominate. The properties of the soil material itself causes a significant change to the toe erosion and loosened zone. The transport capacity of fine grains through coarser soil matrix increases as grain size range increases and then sand boil occurs earlier. Meanwhile, as the grain size range increases, the larger size of eroded material contributes the internal subsidence of slope. (c) Slope failure mode is analyzed with respect to micromechanical features of stress and displacement state. The concentrated rainfall period and enriched alluvial deposition in jianghan Plain are the triggering factors of this investigated case. The tension cracking at the rear of slopes, include erosion at slope toe, are the precursors of failure. The tension crack connects with a potential shear zone, and the soil mass at the front edge continue slow sliding by seepage, finally forming the overall slope failure.

Figure 1
Location map of the studied area.

Figure 2
Cross section of the uvial terrace slope as obtained from eld surveys and photographs of the 2016 Jianghan Plain event. Soil slope with a type of internal erosion during high water events along with example pictures illustrating discontinuous gullies, and a sand boil at the toe of the slope.  The Schematic diagram of the model test apparatus.

Figure 5
Geometry of the model slope and con guration of monitoring sensors.

Figure 6
Experimental photographs of the slope at each moment of failure:(b)settlement crack (after 78min); (b)linear cracks (after 128min) ;(e) heave induced radial cracking at foot of the slope (after 282min); (f)sand boil at foot of the slope (after 430min); (i)erosion at slope toe(after 457min); (g)tensile crack formation at the rear of slope(after 489min); (h) slow retrogressive toe sliding with visible deformation as a whole. Pictures shown are for the gradualloadingexperiments.    Time-varying experimental measured and simulated slope displacement in x-direction (HL=30cm, P=1000pa).

Figure 11
The displacement eld of particles in z-direction with different water boundary (t= 1.2s).

Figure 12
The displacement eld of particles in z-direction with different particle size distribution (HL=30cm, t= 1.2s).

Figure 13
Effect of particle size distribution on development of deformation and stress of two monitoring points Circle 5 and Circle 6.