Simulating the Overtopping Failure of Homogeneous Embankment by a Double-Point Two-Phase MPM

Abstract

Embankments are usually constructed along rivers as a defense structure against flooding. Overtopping failure can cause devastating and fatal consequences to life and property of surrounding areas. This motivates researchers to study the formation, propagation, and destructive consequences of such hazards in risk analysis of hydraulic engineering. This paper reports a numerical simulation of failure processes in homogeneous embankments due to flow overtopping. The employed numerical approach is based on a double-point two-phase material point method (MPM) considering water–soil interaction and seepage effects. The simulated results are compared to available laboratory experiments in the literature. It was shown that the proposed method can predict the overtopping failure process of embankments with good accuracy. Furthermore, the effects of the cohesion, internal fiction angle, initial porosity, and maximum porosity of soil on the embankment failure are investigated.

1. Introduction

Embankments are frequently used along rivers against flooding events, the probability of which has increased due to increasing intensity extreme weather events in recent years [1]. Their overtopping failure may cause serious loss of life and property downstream. For example, the Banqiao dam-breaking flood due to Typhoon Nina in 1975 lead to about 62,000 deaths [2]. Therefore, a deep investigation of overtopping failure mechanisms and processes in embankments is important for the risk management and assessment. Among the available study methods, numerical modeling has distinctive advantages such as the relatively low cost and the flexible application to almost arbitrary site configurations. In the past decades, many numerical investigations of the overtopping failure in embankments have been performed. For example, Powledge et al. [3] studied the embankment erosion due to overflow, Tingsanchali and Chinnarasri [4] used a one dimensional model to study the overtopping failure in embankment. Chinnarasri et al. [5] investigated the overtopping-induced progressive damage in embankment. Leopardi et al. [6] reviewed the modeling of embankment overtopping, Pontillo et al. [7] employed a two-phase model to evaluate the dike erosion induced by overtopping, Volz et al. [8] modeled the overtopping failure of a non-cohesive embankment by the dual-mesh approach, Mizutani et al. [9] simulated the overtopping failure of embankment considering infiltration effects, Guan et al. [10] and Kakinuma and Shimizu [11] developed a 2D shallow-water model for embankment breach, Evangelista [12] simulated dike erosion induced by dam-breaking flows using a depth-integrated two-phase model, Larese et al. [2] modeled the overtopping and failure of rockfill dams using the particle finite element method (PFEM). Although many efforts have been made to numerically reproduce the main features of the process, the overtopping failure of embankments is still poorly understood. Moreover, most of the aforementioned numerical methods are mesh-based approaches and have difficulties in accurately modeling the time evolution of the embankment breach, since the real overtopping failure phenomenon is, in fact, an unsteady process involving large deformations, free-surface flows, moving boundary, and water–soil interactions. To deal with these challenges, several meshless methods, such as smoothed particle hydrodynamics (SPH) [13,14], moving particle semi-implicit method (MPS) [15], element-free Galerkin method (EFGM) [16], and material point method (MPM) [17,18], have been developed to model large deformation problems. For examples, Gotoh et al. [19,20] modeled embankment erosion due to overflow using MPS, Li et al. [21] simulated erosion of HPTRM levee based on SPH, Zhang et al. [22] modeled failures of dike due to water level-up and rainfall using SPH, Liu et al. [23] simulated piping erosion process of dike foundation by EFGM. Nikolic et al. [24] presented a discrete beam lattice model capable of simulating localized failure in a heterogeneous fluid-saturated poro-plastic solid. Zhao and Liang [25] applied MPM to model seepage flow through embankments, Martinelli et al. [26] modeled the failure of a sand dike due to seepage flow using two-phase MPM. In this work, the MPM was employed for modeling the overtopping failure of embankments because it has several advantages over other meshless methods such as EFGM, MPS, and SPH: The boundary condition realization is as easy as in the finite element method (FEM) for the use of background mesh in MPM [27]. Time-consuming neighbor particle searching is required in SPH, MPS and EFGM, but the MPM requires only the identification of particles relative to the background mesh. The MPM also avoids the tensile instability that is evident in the SPH [27].

The MPM was first proposed by Sulsky et al. [17] in 1994, which uses the Lagrangian particle and the Euler background mesh to describe large deformation problems. The particles move freely across the background mesh and carry all physical parameters such as stress, density, mass, velocity, and other historical variables. The MPM combines the advantages of the Euler method and the Lagrangian method to avoid the interference of convection terms and mesh distortion. It also can treat contact problems between different objects without any additional interface elements. In addition, the traditional soil constitutive laws such as Mohr–Coulomb, Drucker–Prager, and Cam Clay models can be easily implemented in MPM at the particle level [17,18]. Recent developments in the multiphase MPM formulations considering the soil–water interaction can be categorized into two main groups: Single-point approach and double-point approach. In the first approach, each material point (MP) possesses both the soil and water information. The water is described in the Eulerian framework while the soil skeleton is described in the Lagrangian framework. The single-point approach does not guarantee the mass conservation of water. In the second approach, two sets of MPs are used to represent water and solid maters. The double-point approach can guarantee both the solid and water mass conservation. Mackenzie-Helnwein et al. [28] used the double-point approach to model soil–water mixtures through a drag model for capturing the interphase interaction. Abe et al. [29] employed the double-point approach to model saturated soil without considering the relative acceleration between water and the soil skeleton. Bandara and Soga [30] developed a double-point two-phase MPM formulation based on the mixture theory considering the relative acceleration between water and the soil skeleton. This approach not only can be used for modeling rapid deformation problems but also guarantees both the soil skeleton and water mass conservation. It is therefore very suitable for modeling soil–water interaction problems such as currently considered embankment failure due to erosion.

The present work aims to numerically investigate the failure mechanism and failure process of homogeneous embankment due to overtopping flows based on a double-point two-phase MPM method.

2. Mathematical Formulation

The double-point two-phase MPM formulation [26,31] based on the mixture theory is used to model the water–soil interactions. In this approach, the soil skeleton and the pore liquid in the saturated soil are represented separately by two sets of material points: Solid material points and water material points (Figure 1). This paper focuses on the application of two-phase MPM to simulate embankment overtopping, only some of the basic equations solved are quoted for the convenience of readers in understanding the two-phase MPM formulation, and for detailed description of the algorithms readers are referred to the literature [26,31].

The governing equations of the double-point two-phase MPM include the mass and momentum conservations of both solid and liquid phases with their respective constitutive equations. The mass conservation equations for the solid and liquid phase are respectively given by [31].

$$\frac{D{n}_{\mathrm{L}}}{Dt}=n_{\mathrm{S}}\nabla \cdot v_{\mathrm{S}}$$

(1)

$$\frac{D{\epsilon }_{\mathrm{L}}}{Dt}=\frac{1}{n_{\mathrm{L}}}[n_{\mathrm{S}}\nabla \cdot v_{\mathrm{S}}+n_{\mathrm{L}}\nabla \cdot v_{\mathrm{L}}+(v_{\mathrm{L}}-v_{\mathrm{S}})\nabla n_{\mathrm{L}}]$$

(2)

where $n_{\mathrm{S}}$ is the volumetric concentration ratio of solid, $n_{\mathrm{L}}$ is the volumetric concentration ratio of liquid. For saturated soils, $n_{\mathrm{L}}$ is equivalent to porosity of the soil skeleton and $n_{\mathrm{L}}+n_{\mathrm{S}}=1$. $n_{\mathrm{L}}^{}$ is a discontinuous function at the transition between free surface water and porous medium and it is characterized by two constant values: One in the free surface water and the value of the porosity n in the porous medium. Martinelli [33] proposed the concept of a transition zone in which an interpolated liquid concentration ratio is introduced which gives a smooth transition of $n_{\mathrm{L}}$ between the free surface water and the porous medium. $v_{\mathrm{S}}$ is the velocity vector of the solid phase, $v_{\mathrm{L}}$ is the velocity vector of the liquid phase. ${\epsilon }_{\mathrm{L}}$ is the volumetric strain of the liquid. $\frac{D(•)}{Dt}$ denotes the material time derivative. Equation (1) represents the variation of volumetric concentration ratio of the liquid phase (i.e., porosity), and Equation (2) is also known as the storage equation and represents the volumetric strain rate of the pore liquid.

The momentum conservation of the solid and liquid phase can be respectively expressed as [31]:

$${\overline{\rho }}_{\mathrm{S}}{a}_{\mathrm{S}}=\nabla \cdot {\overline{\mathsf{\sigma }}}_{\mathrm{S}}+f_{\mathrm{L}}^{\mathrm{d}}+{\overline{\rho }}_{\mathrm{S}}g$$

(3)

$${\overline{\rho }}_{\mathrm{L}}{a}_{\mathrm{L}}=\nabla \cdot {\overline{\mathsf{\sigma }}}_{\mathrm{L}}-f_{\mathrm{L}}^{\mathrm{d}}+{\overline{\rho }}_{\mathrm{L}}g$$

(4)

where ${\overline{\rho }}_{\mathrm{S}}$ and ${\overline{\rho }}_{\mathrm{L}}$ are respectively the partial densities of the solid and liquid, computed as the ratio of the mass of each constituent with respect to the reference volume. g is the gravity vector, $a_{\mathrm{L}}$ is the liquid phase acceleration, $a_{\mathrm{S}}$ is the solid phase acceleration. ${\overline{\mathsf{\sigma }}}_{\mathrm{S}}={\mathsf{\sigma }}^{\prime }-n_{\mathrm{S}}{\mathsf{\sigma }}_{\mathrm{L}}$ and ${\overline{\mathsf{\sigma }}}_{\mathrm{L}}=n_{\mathrm{L}}{\mathsf{\sigma }}_{\mathrm{L}}$ are the partial stresses of solid and liquid phases respectively, ${\mathsf{\sigma }}^{\prime }$ is the effective stress tensor, and ${\mathsf{\sigma }}_{\mathrm{L}}$ is the liquid phase stress tensor. $f_{\mathrm{L}}^{\mathrm{d}}$ is the drag force, which represents the water–soil interaction due to the velocity difference between the two phases, it can be calculated by [34]

$$f_{\mathrm{L}}^{\mathrm{d}}=\frac{n_{\mathrm{L}}^2\mu }{\kappa }(v_{\mathrm{L}}-v_{\mathrm{S}})+\beta n_{\mathrm{L}}^3{\rho }_{\mathrm{L}}\left|v_{\mathrm{L}}-v_{\mathrm{S}}\right|(v_{\mathrm{L}}-v_{\mathrm{S}})+{\mathsf{\sigma }}_{\mathrm{L}}\nabla n_{\mathrm{L}}$$

(5)

where $\mu$ is the liquid dynamic viscosity, $\kappa =\frac{D^2}{A}{n}_{\mathrm{L}}^3/{(1-n_{\mathrm{L}})}^2$ is the soil intrinsic permeability with $D$ the average diameter of grains, $\beta =B/\sqrt{\kappa A{n}_{\mathrm{L}}^3}$ is the non-Darcy flow coefficient, the empirical constants $A$ and $B$ are respectively set to 150 and 1.75 [34].

In order to fully describe the behavior of saturated soils, the constitutive equations for both phases are required. Assuming the validity of Terzaghi’s effective stress concept, the general form of stress–strain relationship for soil skeleton is given by

$$\frac{D{\mathsf{\sigma }}^{\prime }}{Dt}=D\frac{D\epsilon }{Dt}$$

(6)

where D is the tangent stiffness matrix, ${\mathsf{\sigma }}^{\prime }$ is the effective stress vector, and $\mathsf{\epsilon }$ is the total strain. For the liquid phase, the volumetric stresses is updated by

$$\frac{D{\mathsf{\sigma }}_{\mathrm{Vol},\mathrm{L}}}{Dt}=K_{\mathrm{L}}\frac{D{\epsilon }_{\mathrm{Vol},\mathrm{L}}}{Dt}$$

(7)

where $K_{\mathrm{L}}$ is the liquid bulk modulus, ${\epsilon }_{\mathrm{Vol},\mathrm{L}}$ is the volumetric strain of the liquid. In addition, for pure liquid or fluidized mixture, the deviatoric part of the stress tensor is calculated by

$${\mathsf{\sigma }}_{\mathrm{dev},\mathrm{L}}=2\mu \frac{D{\epsilon }_{\mathrm{Vol},\mathrm{L}}}{Dt}$$

(8)

In the proposed double-point two-phase MPM, the saturated soils can be considered as a solid-like or liquid-like state according to the porosity. As shown in Figure 2a, in the case of low porosity, the solid skeleton grains are in contact and the soil behaves like a solid, its response thus can be modeled by a soil constitutive model. For a high-porosity soil, as shown in Figure 2b, the grains are not in contact and float together with the liquid phase and the soil response is modeled by the Navier–Stokes equation. In the present work, a pre-given maximum porosity n_max is used to distinguish the two aforementioned states. When the soil porosity is less than the maximum porosity ($n_{\mathrm{L}}=1-n_{\mathrm{S}}<n_{\max }$), the mean effective stress decreases as the porosity increases and the mean effective stress vanishes once the grains are not in contact. When the porosity is larger than n_max, fluidization occurs.

3. Numerical Examples

In this section, the proposed double-point two-phase MPM is first validated by an example of flow through a porous block. Then a numerical investigation of failure process in a homogeneous embankment due to flow overtopping was performed and the effects of the cohesion, internal fiction angle, initial porosity, and maximum porosity of soils on the embankment failure are investigated. All calculations are conducted with the software Anura3D_v2016 (www.anura3d.com), which implements the double-point two-phase MPM formulation.

3.1. Flow through a Porous Block

In order to validate the capability of the current MPM algorithm to model the water–soil interaction, a numerical simulation of the flow through a porous block was performed and the calculated results were compared with the experimental data obtained by Liu et al. [35]. The experiment as sketched in Figure 3 was carried out in a transparent water tank with a size of 0.892 m × 0.44 m × 0.58 m, the porous block, which was made of gravel, was located at the center of the water tank (x = 0.30–0.59 m), and it was 0.29 m long, 0.44 m wide, and 0.37 m high. The water on the left side of the water tank was separated from the porous block by a movable gate with a thickness of 0.02 m, and the height was set to 0.35 m, and the porous block was confined in the initial region to ensure that the porous media was not allowed to move. The material parameters used in this simulation are listed in

Table 1. The material parameters for flow through the porous block.

Material	Parameter	Value
Porous block	Density (kg/m3)	2700
	Young modulus (kPa)	1000
	Poisson’s ratio	0.3
	Initial porosity	0.49
	Mean diameter (mm)	15.9
Water	Density (kg/m3)	1000
	Bulk modulus (kPa)	2.15 × 104
	Dynamic viscosity (kPa·s)	1 × 10−6

. The computational domain was divided into 4104 tetrahedral elements. Each mesh element in the porous block part contained four material points, while each element in the water contained 10 material points.

The simulation process was divided into two stages, the first stage was the initialization of stress, also called the gravity loading process. At this stage, the lateral displacement of the water was constrained and the stress initialization was done by increasing the gravity until the solution converged to a quasi-static equilibrium, which was checked by evaluating the normalized kinetic energy and the normalized out-of-balance force falls below a predefined tolerance as given in [31]. In addition, a local damping coefficient of 0.75 was applied to accelerate the convergence to quasi-static equilibrium allowing a considerable reduction in the computational time. In the second stage, the lateral constraint of water was removed to simulate the instantaneous moving process of the gate, and then water flow impacted the porous block. Figure 4 shows some snapshots of the flow through the porous block at time of 0.35, 0.75, 1.15, and 1.95 s and the calculated free-surface profiles agree well with that of the experimental results [35], and comparisons of free surface profiles are shown in Figure 5. It is proven that the double-point MPM can simulate the process of flow through the porous medium considering the water–soil interaction.

3.2. Overtopping Failure of Homogeneous Embankments

For a better understanding of the embankments failure mechanism and process due to the overtopping flow, in this section, MPM modeling of the overtopping erosion of a homogeneous embankment considering free water surface flow, water–soil interaction, and seepage effects was performed and the simulated results were compared with the experimental data obtained by Evangelista et al. [12]. The effects of the cohesion, internal fiction angle, initial porosity, and maximum porosity of soil on the embankment failure were investigated. As sketched in Figure 6a, this experiment was carried out in a Perspex horizontal channel with a 0.40 m wide rectangular cross-section and transparent walls. A 3.00 m long tank with a moving gate was used to store water, in which H represents the water level. The gate was suddenly lifted to produce the dam-break flow, which propagates over the downstream channel bottom. The first segment of the rigid bed was 1.00 m long and had a mild slope (5%), the subsequent one was a 2.00 m long, 0.05 m thick horizontal segment, on which the embankment was built.

According to Evangelista et al. [12], the friction at the lateral interface of the channel has no significant wall effects and the embankment undergoes an almost plane process because the flow overtops the entire embankment width. For these reasons, in our simulation, the overtopping failure process was simplified approximately as a two-dimensional problem. In addition, the gate opening was considered instantaneous following the work of Lauber and Hager [36]. The computational domain was divided into a set of tetrahedral elements with the size of 0.05 m. Each mesh element in the embankment part contained four material points, while each element in the water contained 10 material points (see Figure 6b). The water was modeled as a weak compressible fluid, and the Mohr–Coulomb model was used as the constitutive model of the soil. The material parameters used in this simulation are listed in

Table 2. The material parameters of the two phase materials.

Material	Parameter	Value
Sand	Density (kg/m3)	2680
	Young modulus (kPa)	1000
	Poisson’s ratio	0.3
	Internal friction angle (°)	36
	Cohesion (kN)	0
	Initial porosity Maximum porosity	0.40 0.45
	Mean diameter (mm)	0.2
Water	Density (kg/m3)	1000
	Bulk modulus (kPa)	2.15 × 104
	Dynamic viscosity (kPa·s)	1 × 10−6

. Three different cases with initial water level of H = 0.2, 0.25 and 0.395 m were conducted.

For the first case (H = 0.2 m), Figure 7 shows several snapshots of the simulated water surface and the embankment profiles along the channel at different times of 0.35, 0.8, 1.4, 2.0, and 3.0 s. Comparisons between the numerically simulated results and the experimental ones [12] are illustrated in Figure 8. Specifically, the embankment (black color) profiles and the water surface (blue color) at times 1.4 and 3.0 s, respectively, after the gate removal are plotted: Solid lines represent the numerical results and dotted lines indicate the experimental ones. In this case, the initial water level in the water tank was almost equal to the top elevation of embankment, thus the gate opening–induced dam break flow did not have enough energy to overtop the dike embankment, and therefore, it induced erosion only on the upstream slope of the embankment.

For the second case (H = 0.25 m), Figure 9 shows several snapshots of the simulated water surface and the embankment profiles along the channel at different times of 0.35, 0.8, 1.4, 2.0, and 3.0 s. Comparisons between the numerically simulated results and the experimental ones [12] are illustrated in Figure 10. Specifically, the embankment (black color) profiles and the water surface (blue color) at times 1.4 and 3.0 s, respectively, after the gate removal are plotted: Solid lines represent the numerical results and dotted lines indicate the experimental ones. In this case, the initial water level in tank was higher than the top elevation of embankment, thus the gate opening–induced dam breaking flow overtopped the embankment and eroded it along its entire profile. The upstream slope erosion was more prominent due to the higher potential energy of the dam breaking flow. This dam breaking flow also reached the crest of the embankment this time, progressively eroding it as it passed.

For the third case (H = 0.395 m), Figure 11 shows several snapshots of the simulated water surface and the embankment profiles along the channel at different times of 0.35, 0.8, 1.4, 2.0, and 3.0 s. Comparisons between the numerically simulated results and the experimental ones [12] are illustrated in Figure 12. Specifically, the embankment (black color) profiles and the water surface (blue color) at times 1.4 and 3.0 s, respectively, after the gate removal are plotted: Solid lines represent the numerical results and dotted lines indicate the experimental ones. In this case, the potential energy of the gate opening induced dam break flow was about two times that in the previous case, thus leading to more prominent erosion in the embankment. The erosion was quite homogeneous along the entire profile, and the embankment profile shape was rounded during the process. During the process, steep portions of the embankment slopes underwent a combination of surface erosion and sliding failures due to the loss of supporting material.

In general, the erosion process of the embankment was greatly affected by the upstream water level elevation. The higher the initial water level, the more serious was the rate of damage of the embankment. By analyzing the embankment erosion process at low, medium, and high water levels, it can be seen that when the water level is low or medium (H = 0.2, 0.25 m), the water flow coming to the embankment takes a long time, and the erosion is mainly concentrated on the upstream slope surface. When the water level is high (H = 0.395 m), the water flow has a larger energy and rises rapidly along the embankment and overtops the crest to the downstream. Under the erosion and seepage of the water flow, the entire embankment is damaged seriously. Overall, the numerical simulation can relatively reproduce the overtopping erosion phenomenon, which validates the capacity of the proposed double-point MPM to model the overtopping failure of the homogeneous embankments with considering free water surface flow, water–soil interaction, and seepage effects.

In the previous simulation, the homogeneous embankments were considered as non-cohesive. However, we must remember that many embankments were built from cohesive materials, which behave quite differently than non-cohesive materials. In order to study the effect of cohesion on the embankment erosion process, the previous third case (H = 0.395 m) was selected, the simulations were conducted with three cohesions of c = 0, 5, and 10 kPa, the other conditions remained unchanged. Simulated results are shown in Figure 13, the erosion rate of the upstream and downstream slopes was more serious when c = 0 kPa. However, when c = 5 and 10 kPa, the rate of embankment erosion was very small and almost negligible. It can be seen that the increase of soil cohesion can effectively improve the resistance. As the cohesion increases, the erosion caused by the water flow is weakened, but when increased to a certain extent, the effect is weakened.

The effect of the internal fiction angle on the embankment erosion process was also investigated. Three different internal fiction angles ($\phi$ = 30°, 36°, 41°) were used in the case of the initial water level H = 0.395 m, and the other parameters were the same as in the previous simulations. The calculated embankment profiles for different internal fiction angles at t = 1.4 s and t = 3.0 s are shown in Figure 14. Results show that the larger the internal friction angle is, the weaker the erosion will be. It is the internal friction angle that influences the contact between soil particles. For the small internal friction angle, it has weak resistance to overtopping flow, whereas a large friction angle can retard the flow.

In order to investigate the effect of the initial porosity on the embankment erosion process, three different initial porosities (n = 0.2, 0.3, 0.4) were selected to model in the case of the initial water level H = 0.395 m, and the other parameters were also the same as in the previous simulations. The calculated embankment profiles for different porosities at t = 1.4 s and t = 3.0 s are shown in Figure 15. It can be seen that in the early stage (t = 1.4 s) of the embankment failure, the erosion was not very obvious, only surface erosion caused by the water flow occurred and there was no infiltration. However, the profiles of the embankments at t = 3.0 s were different. The larger the porosity of the soil, the more obvious the embankment erosion. It indicates that the water flows began to infiltrate into the embankment in the middle stage, and the embankment erosion was induced by the water flow and seepage, and the seepage failure mainly occurred in the middle and later stages. With the gradual infiltration of the water flow, the water flow infiltrated completely in the later stage, and the pores of the soil were filled with water, the soil became saturated. When the effective stress of the soil skeletons is 0, the soil may be liquefied and moved with the water flow. Hence, decreasing the porosity of the soil can slow down the infiltration rate and reduce the damage of the embankments.

In the proposed double-point two-phase MPM, a pre-given maximum porosity n_max was used to distinguish the saturated soils as a solid-like or liquid-like state. When the porosity was larger than n_max, fluidization occurred. The grains are not in contact and float together with the liquid phase. Here, three cases of different porosities (n_max = 0.45, 0.50, 0.55) were conducted to study the effect of maximum porosity on the embankment erosion process, and the other conditions were the same as in the previous simulations. The calculated embankment profiles for different cases at t = 1.4 s and t = 3.0 s are shown in Figure 16. Results show that the maximum porosity plays an important role on the dam failure due to overtopping flow. The larger the maximum porosity is, the more difficult it is for soil particles to reach the fluidization state. Hence, it is very important to get the right maximum porosity in the simulations, and more accurate simulation results can only be obtained with the appropriate value.

4. Conclusions

In this paper, a double-point two-phase material point method was employed to model problems involving free water surface flow, large deformations, and water–soil interaction. The proposed method was first validated by the example of the flow through a porous block problem. Then a numerical investigation of failure process in homogeneous embankments due to overtopping flow was performed. By comparing the water flow and embankments profiles with the experimental data, it was shown that the current double-point two-phase MPM can predict the overtopping flow pattern and the development of embankment erosion with good accuracy, and the higher the initial water level, the more obviously and faster the embankment is eroded. Furthermore, the effects of the cohesion, internal fiction angle, initial porosity, and maximum porosity of soil on the embankment failure were investigated. The results show that these parameters have important influence on the erosion process. With the increase of the cohesion, the resistance of the embankment can be improved to a certain extent, but the effect of the cohesion is gradually weakened. The larger the internal friction angle is, the weaker the erosion will be. In addition, the larger the initial porosity of the soil, the more obvious the embankment erosion. The larger the maximum porosity is, the more difficult it is for soil particles to reach the fluidization state. These investigated results indicate that the double-point two-phase MPM is capable of predicting and reproducing the failure process of embankments due to the overtopping flow. The proposed method is an alternative promising tool for investigating complex failure mechanism in water–soil interaction.