Numerical Modeling for the Effect of Soil Type on Stability of Embankment

Dike construction has been widely used because of its potential to protect people and properties from overtopping flows. Water levels may exceed a dike crest and cause overtopping flow during high river discharge. This phenomenon has caused serious damage to the dike body due to the reduction of soil shear strength. The increase of water content within particles and its relationship with the development of breach channel failure in downstream and upstream slopes are affected by a series of geotechnical and hydraulic aspects. Transient seepage and slope stability analyses (FOS) were performed in this study using 2D finite element methods and time-history measurements under the effect of sandy and very silty sand soils. The numerical model of SLIDE 2018 was limited by its inability to incorporate all physical processes governing an overtopping breach failure. Numerical analyses were performed to simulate the development of pore pressures and water content at six positions in the dike’s upstream and downstream slopes in physical experimental tests using the van Genuchten Equation and the limit equilibrium method. The numerical results revealed that fine particles increase the pore water pressure and reduce the FOS. Appropriate dike design and maintenance are dependent on surrounding hydraulic conditions, dimensions, and soil types. Non-cohesive materials with fine particles were preferable.


Introduction
The development of dike construction is highly related to people's lives and property [1]. Dikes are widely used to save various types of lands and avoid the danger of overtopping failure [2,3], drinking and irrigation, energy production and recreation purposes [4] as well as support power generation [5]. Overtopping failure, which occurs due to the transition of overtopping wave into the downstream slope, is considered the main cause of dike's failure [6,7]. The dike's breach failure is dependent on different geotechnical and hydraulic parameters, such as excess rainfall, upstream slope collapse, settlement of the dike crest, inadequate design and construction of the spillway, and global warning [8]. Breach channels inside the dike are initiated due to water infiltration above the crest during overtopping failure. Consequently, the development of the breach channel widens in the vertical and horizontal directions from the downstream slope until it covers the upstream slope. The continuous erosion process transports sediment materials along the foot of the dike. Shapes of breach channels during the erosion process are triangular, rectangular, and trapezoidal [9].
Froehlich [10] observed three types of breach growth models. The shape of breach channel due to overtopping failure in model "A" is initially triangular and continue to evolve until it reaches the valley floor whilst the shape of breach channel is a trapezoid. The initial breach channel in model "B" is trapezoid in shape and continue to evolve until the end of the erosion process. The progression of breach channel width is constant in model C. The breach channel is undercut and the side slope is sufficiently saturated to start failure during the erosion process in downstream and upstream slopes of river dike. Xu and Zhang [11] stated that the use of multiple regression analysis with variables, such as reservoir volume and embankment channel height, is more accurate than applying single regression analysis and demonstrated that erodibility is an important factor that influences the embankment breach parameter. Stability of the dike slope against the erosion of the dike crest is dependent mainly on soil suction and cohesion.
Powledge et al. [12] investigated the influence of compaction efforts on a cohesive (clay) dike. Increasing the compaction effort from 95% to 103% of the standard proctor value will reduce empty pores between soil particles and delay the reduction of matric suction and breach development. Powledge et al. [13,14] performed large laboratory and field tests to evaluate the breaching mechanism. The erosion regime was divided into three zones depending on the velocity of the water flow. The rate of erosion is a function of different parameters, such as embankment configuration, maximum flow velocity and voids in the slope. The researchers observed that the water surface during the breach failure inside the dike is accelerated in the downstream slope with the exception of long channel reaches such as in the case of landslide dams. Coleman et al. [15] observed the reduction of the dike crest and increase of vertical and horizontal erosion inside the breach channel during spatial overtopping tests. Although sediment transport is eroded into the downstream slope, the transmitted area between downstream and upstream slopes and the toe of the upstream slope is nearly removed. The shape of the last stage of the stabilizing erosion process is a curved hourglass, similar to the results of Rozov [16].
Chinnarasri et al. [17] explored the behaviour of homogeneous embankment soil during the progression of the erosion process inside rectangular pilot channels. Inflow to the upstream reservoir, which was measured using an orifice flow meter, increased until outflow occurred through the pilot channel. The researchers indicated that the vertical erosion process first develops in the dike crest during the initiation of breach channel failure whilst the longitudinal profile of the breach channel is parallel to the downstream slope. The lateral erosion process becomes dominant in the upstream slope until the end of the erosion process in the bed of dike construction. The researchers stated that tractive shear forces and side-slope collapses accelerate the widening erosion process. Schmocker and Hager [18] constructed trapezoidal embankment dikes with various grain sizes and embankment geometries to observe the lateral erosion process during overtopping tests under different scale factors. They stated that scale series of large grain sizes tend to accelerate the variability in the breach channel failure during the erosion process, due to the governance of sliding failures observed in the downstream slope. Hassan and Ismail [19] analysed the effect of constant inflow discharges on the response of pore water pressure inside the dike sand embankment, and observed that high inflow increases the water content and thus reduce soil stability. The velocity of water accelerated the erosion in dike slopes.
Hassan and Ismail [20] examined the effects of two dike slopes of 1V:3H and 1V:2.5H on the evolution of breach channel failure in the embankment. A trapezoidal channel was initiated to conduct the physical overtopping tests. Lateral and horizontal erosion processes are fast for steep slopes because of the high water velocity. Numerical and analytical models have been used to simulate breach failures of dike embankment during overtopping failure [21][22][23][24]. These models use a series of differential equations in solving nonlinearity behaviour of seepage mechanism at specified boundary conditions. The accurate prediction of FOS is an indicator used for the stability of dike slope during water flow seepage [25].The validation of wave dissipation of slope breaches is highly required for obtaining run-up and overtopping over the structure [26]. However, a previous study on physical experimental tests and numerical models [27], demonstrated that the influence of some important geotechnical and hydraulic parameters on the development of matric suction, volumetric water content and slope stability remains unclear. Examples of these geotechnical and hydraulic aspects include different soil types for upstream and downstream slopes. These aspects must be considered because they govern the erosion process and side slope stability during overtopping failure. Two types of coarse sand and very silty sand soils were used in this study to determine responses of positive pore pressure, volumetric water content and factor of safety (FOS) using SLIDE 2018 software. The terms "E4" and "E6" indicate numerical analyses of transient seepage and slope stability for coarse sand and very silty sand soils, respectively.

Research Methodology
Physical spatial overtopping tests were described previously by Hassan and Ismail (2018). Measurements of matric suction and volumetric water content in physical tests are conducted using a tensiometer and TDR sensors. Laboratory tests results are subsequently used as an input parameter for the numerical modelling of seepage flow and slope stability analysis for physical experimental tests, as presented in Figure 1. They are used for estimating the hydraulic conductivity and volumetric water content functions during the transient seepage flow from the base of the upstream slope until the end of the dike crest.    The volumetric water content (VWC) functions are calculated directly using the WP4C. The WP4C is used to calculate the SWCC for coarse sand and very silty sand soils using chilled technique system. The SWWC results showed that the saturated and residual water contents are used to evaluate the volumetric water content function. An estimated method of van Genuchten [28] is used to predict the empirical SWCC as follows: where α ψ , n ψ , and m ψ are a function of the AEV, slope of straight line segment in SWCC and the residual water content, respectively. These equations produce high flexibility in representing the SWCC as well as cover a wide range of matric suctions and different soil types due to the m parameter that provides increased stability during parameter optimisation. Sisson & van Genuchten [29] proposed the relationship of m = 1-1/n (n >1, 0 < m < 1) to commit a constant m and to provide a closed form-expression. SLIDE 2018 software is used to establish the volumetric water content functions, and built -in function data of coarse sand and very silty sand soils are utilized to determine curve fitting parameters. Input parameters of soil properties for estimating the volumetric water content function are listed in Table 1. Hence, the hydraulic conductivity function can be estimated once the volumetric water content function and hydraulic conductivity are specified. Required input parameters for estimating the hydraulic conductivity function are presented in Table 2. Measurements were conducted in six positions (groups A-F) in upstream and downstream slopes similar to physical tests as shown in Figure 4.

Figure 4. Locations of tensiometer and TDR sensors in dike embankment
Groups A and B are located in the downstream slope, groups C and D are within the transition area between downstream and upstream slopes, and groups F and E are in the upstream slope. Responses of pore water pressure and volumetric water content for each group were measured during the transient seepage analysis from the toe of the upstream slope until the end of the dike crest. Distinct group distances, measured from the upstream slope towards the downstream slopes for both E4 and E6, are presented in Table 3. Time analyses of FOS based on General Limit Equilibrium (GLE) are 600 and 840 s for E4 and E6, respectively whilst the boundary conditions are similar to the seepage analysis (Table 4).

Soil Water Characteristics Curve Test
Responses of negative pore water pressures (MPa) are plotted against those of volumetric water contents (%) for both soils as shown in Figure 5. The saturated water content, residual water content and air entry value (AEV) were determined on the basis of SWCC. The water content is constant for both soil types in the beginning of tests. The SWCC curve then decreases as the matric suction increases at 0.2 and 1.2 MPa for E4 and E5, respectively. The volumetric water content becomes constant with further increments in matric suctions at 0.9 and 2 MPa because the high permeability of coarse sand prolongs de-saturation. Coarse sand and very silty sand soils presents saturated water contents of 0.19 and 0.27, and AEV values of 0.3 and 1.4 MPa, respectively. The saturated zone of very silty sand is longer than that of coarse sand due to the existence of fine particles. Fine particles typically absorb a large amount of water due to their small and numerous pores. The residual water content of silty sand is higher than that of coarse sand soil because of the faster water evaporation within large particles of the latter soil.
Large pores are at high matric suctions are emptied and small amounts of water are retained. Saturation water contents, AEV and residual water contents coarse sand are lower than those of the other soil type because the higher percentage of large particles (sand and gravel) in the former significantly accelerates water evaporation. Increasing the percentages of coarse particles increases the dewatering potential. Values of saturated water content, residual water content and saturated hydraulic conductivity are used to estimate the hydraulic conductivity and volumetric water content functions using the van Genuchten equation. In addition, a series of laboratory tests was carried out to determine basic soil properties of coarse sand and very silty sand soils as shown in Tables 5 and 6, respectively. These tests are used as input parameters for analyzing the mechanism of seepage flow and slope stability during overtopping moment, to understand the mechanism of seepage flow and slope stability during overtopping moment.

Numerical Processes for the Transient Water Level Inside Dike Embankment
Seepage flow analysis using dike embankments is essential to calculate the change in pore water pressures and volumetric water contents during the overtopping moment in numerical modelling. The seepage flow through saturated-unsaturated soil must be considered because the dike soil state consists of saturated-unsaturated conditions. The scheme for the seepage analysis is illustrated in Figure 6.

Estimation of the volumetric water content
· SWCC test · Typical data values saved in model

Estimation of Hydraulic conductivity function
· Saturated hydraulic conductivity, · Volumetric water content function Discretization of the model geometry

Assignment of boundary conditions
· Transient water level from the toe of the upstream slope until the dike crest where is h Total hydraulic head, k x is Coefficient of permeability in x-direction, k y is Coefficient of permeability in ydirection, Q is Applied boundary flux such as evaporation, rainfall etc, m w is Slope of the SWCC, and θ w is Volumetric water content.
The numerical model of 2D seepage and slope stability analysis (SLIDE 2018) is used to overcome the complexity and non-linearity of this equation [30]. This numerical method simulates the seepage flow analysis inside dike embankments with the assistance of three functions installed inside the program. The three functions include the SWCC, which defines the relationship between the matric suctions and volumetric water contents; permeability function, which indicates the permeability with respect to matric suction and boundary flux, Q. The homogeneous dike embankment is modelled with SLIDE 2018 to simulate physical experimental tests during the overtopping moment. The development of pore water pressures and volumetric water contents for groups of A, B, C, D, E and F is analysed on the basis of the influence of soil types, as presented in Tables 5 and 6. However, the effect of inflow discharge is ignored. External boundaries of trapezoidal dike embankment are drawn on the basis of Table 7. Soil properties of coarse and very silty sands are then assigned as input parameters according to Tables 8 and 9, respectively.   Effective cohesion (kPa) 0 Effective angle of friction ( ○ ) 37 Saturated water content (%) 19 Residual water content (%) 1 Table 9. Required soil properties of very silty sand for numerical modelling

Saturated water content (%) 27
Residual water content (%) 6 Input parameters of the dry density of soil, shear strength parameter, volumetric water content function and hydraulic conductivity function for the two soils are applied to the dike model on the basis of Tables 8 and 9 for coarse sand and very silty sand soils, respectively.. The dry density, effective cohesion and effective angle of friction are used to define properties of coarse sand and very silty sand soils, whilst the saturated and residual water contents are utilized to estimate volumetric water content and hydraulic conductivity functions using van Genughten equation. The dike model is discretized or meshed into finite elements for the result computation, as shown in Figure 9. An initial water level surface is assigned at the base of the dike to represent the dike's water pressure prior to the transient condition. All points inside the embankment are PWP= 0 due to the position of the water level at the dike base. The transient condition is then applied along the upstream slope and dike crest to simulate the overtopping moment, whilst the maximum height of the water level is 30 cm.
where c ′ is the effective cohesion, ϕ ′ is the effective angle of internal friction, u w is the pore water pressure, N s is the slice base normal force, and W x is the slice weight. , , , , and w are geometric parameters, and λ is the inclination of the slice base. Values of effective cohesion and effective angle of friction are used to compute the groundwater analysis on the basis of finite element mesh. Pore water pressures computed from the seepage analysis are then utilized to determine the slope stability analysis as shown in Figure 10. Time analysis showed that, GLE method generated one FOS for all slopes, and F M and F f are equal when Equations 3 and 4 are satisfied. Sequences of slope stability analysis are similar to those of the seepage flow analysis with the exception of adding limits of potential slope failure in the downstream slope and specifying the required time (in seconds) . Figure 11 shows the upper and lower limits of the potential slope surface failure in the downstream slope of the dike embankment. The total time requirement for slope stability analysis is 600 and 840 s for E4 and E6, respectively. The FOS is plotted against the time requirement of the transient water level along the upstream slope and crest after computing the slope failure results.

Assignment of shear strength parameter
· Effective cohesion · Effective angle of friction

Assignment of boundary condition
· Transient water level from the base of upstream slope up to the end of dike crest

Results and discussions
Although conducting the numerical modelling for analysing the mechanism of seepage flow and slope stability is important, the study presents some limitations due to the obscenity of the breach channel analysis. However, GLE is still considered a suitable method for calculating FOS analysis [31] and [32], and interpreting dike slope failure results. Figure 12 presents the developments of pore water pressures and volumetric water content for E4 during the infiltration of water level inside the dike embankment. The behaviour of large particles (sand and gravel) in E4 is controlled by the weight of the particle and associate friction. Other phenomena such as electromagnetic and intermolecular forces are negligible compared with weight. Therefore, the permeability of water flow is easiest in E4 compared with E6.The seepage water flow increases the pore water pressures for all groups due to the increasing percentage of water content inside soil particles. The mechanism of unsaturated soil is gradually converted into saturated zone with increasing water content and decreasing negative pore water pressure. The infiltration of water level at t =0 s is slight in the toe of the upstream slope because high tension in soil particles enhances the stability of dike soil. High shear bed stresses increase the resistance of particles against the boundary condition of the overtopping flow. The increase of the advancement rate of wetting front in the unsaturated zone of particles at latter seconds due to the reduction in dike shear strength, increase the possibility of erosion. The faster transition of negative into positive pore water pressure in group F is similar to the results of Riffai and Nistor [33], (i.e. before t = 25 s) compared with that of other groups and caused by its location near the toe of the upstream slope. Maximum pore water pressures are observed at t = 600 s with 2.21 kPa and 19% volumetric water content. A small dike model leads to small magnitudes of pore water pressures during the transit water level into the downstream slope. The gentle slope helps accelerate the distribution of water level near group F whilst the water saturation process begins to decreases near group E due to the effect of gravity force and unsaturated voids in the middle of the dike. Group E shows lower negative pore pressure in the beginning of the saturation process compared with group D. These differences are due to the low velocity of water seepage flow of group D that fails to replace the negative pore pressure-and -thus increase the positive pore water pressure and volumetric water content inside particles. The pore water pressure inside voids of group C increases during the seepage flow occuring below the dike crest. Pore water pressure of group C is higher than that of group D as the soil saturates. Pore water pressures at t =480s are 0.85 and 0.73 kPa for groups C and D, respectively, mainly caused by the increase of the wetting front of soil voids with the increase of velocity of seepage flow near the downstream slope. Responses of positive pore water pressures of groups A and B occur at the end of the seepage analysis. The volumetric water content of group B slightly increases by 1.6 % at = 300 s. The hydraulic conductivity of soil is higher than the intensity of water saturation that results in slow saturation. Figure 13 shows the responses of pore water pressures and volumetric water content in very silty sand soil. All groups are saturated with water content during the transition of the water level into the dike crest. The soil particles are filled with a high percentage of air flow and remain attached to the water content in the beginning of the saturation process near the toe of the upstream slope. The penetration of horizontal water levels is difficult at the beginning of saturation. Water-filled pores block the path of airflow during the infiltration of water level due to the high porosity of E6. Water molecules are strongly attached to the surface of fine particles. Consequently, channels inside particles for the transportation of water flow become narrow .The velocity of the saturation process for each group is mainly dependent on their positions, whilst the responses in pore water pressure and volumetric water content of groups E and F are faster than those of other groups. Positive pore water pressures occur at t = 480 s with water contents of 25.9% and 27 % for groups C and D, respectively. This result indicates the faster seepage process for soil layers below the dike crest compared with layers near the bed of the dike model. Moisture contents subsequently increase for groups A and B due to higher matric suctions, whilst the hydraulic conductivity of soil is less than phreatic water level. Consistent with findings of physical tests, the results showed that high pore water pressure is generated for coarse sand in group E. The increase of pore water pressure of (E6) decreases gradually after t = 500 seconds, due to the nearly full saturation condition in the upstream slope that increases the difficulty in replacing the matric suction inside soil particles. Figures 14 to 19 show the comparison of the responses of matric suctions and volumetric water content between E4 and E6. Responses of pore water pressures and volumetric water contents are generally higher in E4 than those in E6. The low hydraulic conductivity of E6, due to the existence of fine particles, allows the absorption of a large amount of seepage flow and plugs the movement of water into the soil because of high matric suction. Fine particles in the very silty sand affect the development of the matric suction and volumetric water content. Volumetric water contents are 0.189 and 0.196 at t = 420 s for E4 and E6, respectively. The higher saturated conductivity of E4 than that of E6 increases the penetration of seepage flow and reduces the capability of particles to store water. The large decrease of matric suction for both E4 and E6 of group A occurs at t = 600 and 840 s, respectively, due to the effect of high matric suction trapped in the downstream slope on the concentration of water inside soil voids.
The saturation of E6 is slower (i.e. before t = 500 s) than that of E4 because the high hydraulic conductivity of fine particles significantly reduces the ability of infiltrating intensity to decrease matric suction for group B. However, the water content of E4 is significantly increases by t = 540 s compared with the 16.2% increase of E6, and the unsaturated soil state (matric suction) becomes saturated (positive pore water pressure) during these periods. This phenomenon is due to the increase of seepage flow velocity above the downstream slope that enables the infiltration intensity to overcome the hydraulic conductivity of the soil. The response of pore water pressure of -0.55 kPa of E4 for group C occurs early at t = 180 s, with an abrupt increase (13.4%) in saturation water content. This phenomenon is due to the high saturated hydraulic conductivity of layers in the upstream slope and crest that accelerates the saturation near the bed of the dike model. The increase in water content of E4 slightly changes after t = 200 s compared with that of E6. The increase sizes of soil voids, in E6 due to cohesive materials, results in high water content. Cohesive forces resist the strength of dike materials. As shown in Figure 17, the seepage of water infiltration of E4 in group D is fast in shallow depths near the crest. Therefore, the infiltration intensity is higher than the hydraulic conductivity of soil, and the maximum water content occurs at t = 200 s. Enhancements of pore water pressures of E4 occur faster than that of E6, whilst volumetric water contents are higher in the latter for groups E and F. Inter-particle forces connect and govern fine particles and effect the absorption of water content of the dike soil. These forces decrease the permeability of water infiltration within soil particles; thus, the distribution of vertical and horizontal water levels near the toe of the upstream slope towards group F is impeded. Figure 20 shows the distributions of pore water pressures of E4, and those of E6 during the transition of water level from the toe of upstream slope until the dike crest. The blue line represents the phreatic water level at t = 25, 180, 600 seconds. Positive pore water pressures for the saturation zone are reduce gradually near the water level, due to the effect of the unsaturated zone (two phase zone) whereas the negative pore water pressure begins to increase. In addition, the negative pore water pressure is high near the downstream slope. Water contents highly reduce in the dry zone of unsaturated soil. Similar to the findings of physical experiments, these results indicated that the reduction of soil shear strength parameter (cohesion and friction) or FOS occurs before the initiation of breach channel failure in the dike crest. Therefore, vertical and horizontal erosion processes will develop fast in the middle of the dike and the upstream slope.  [34] results, increasing the water content of soil particles generates high driving forces (weight of failure surface), which become stronger than resisting forces because of the continuous saturation of soil particles, with the decrease of FOS. FOS is less stable in E6 than that in E4 because the lower saturated hydraulic conductivity of fine particles in E6 that prevents easy seepage flow between particles. The dike slope instability in E6 results in steepness of slope during the breach channel failure.

Figure 21. Comparison of FOS between E4 and E6 using LEM
The steepness of the slope similarly occurs in E6, in which the majority of the transition area between downstream and upstream slopes and the area of the upstream slope are highly eroded. The breach channel slope for both soils becomes gentler during the progressive erosion of the upstream slope and gradually remains constant throughout the breach failure. Similar to the results of previous studies [35,36], down cutting of the breach channel and high erosion occurring near the sides of the breach entrance result in a severely undercut breach channel entrance. Undermining of the saturated dike slope in E4 is fast due to the high water turbulence in porous sand materials.
Hence, the rapid change of water level at the dike structure generally leads to the degradation and reduction of soil parameters, and a high threat of structural failure exists under this condition. Constructing a homogeneous embankment with upstream and downstream slopes of 1V:3H is important to achieve an appropriate design, increase the weight of the dike body, and thus resist the water level pressure in front of the upstream slope. Non-cohesive dike materials with fine particles (silt or clay) are preferred given that they absorb the water level and thus delay and reduce the danger of overtopping failure. Further investigation is still needed to simulate the progression of the erosion process on upstream and downstream slopes and its components, such as breach channel outflow, time to peak, and others, using advanced numerical models. An in-depth study can analyze the FOS using 3-d slope stability software programs and compare the results with those of a 2-d slope stability analysis.

Conclusion
The present study analyzed the time responses of pore pressure, water content, and FOS of coarse and very silty sand soils during the overtopping moment across upstream and downstream slopes. This effect increases the water content inside particles and thus reduces the shear strength of both soils. The negative-to-positive conversion of pore water pressure is mainly dependent on the locations of groups. Water infiltration has a significant impact on groups on the upstream slope, resulting in the conversion of unsaturated soil to partially saturated soil. Conversely, the negative pore water pressure is high for groups on the downstream slope due to the shortened saturation process. Pore pressure is high in coarse sand, whilst the water content is high in very silty sand, because the high permeability of coarse particles corresponds to the easy transportation of water molecules. Fine particles in the very silty sand soil absorb large amounts of water and delay water infiltration due to their narrow voids. The gradual decrease of FOS values in both soils with the increase of water content in the dike embankment reduces the soil shear strength. In addition, the FOS of the very silty sand soil is less than that of the coarse soil because the lower saturated hydraulic conductivity of the very silty sand soil reduces its frictional strength.

Data Availability Statement
The data presented in this study are available in the article.