Dam Breach Simulation with the Material Point Method

Dam embankment breaches caused by overtopping or internal erosion can impact both life and property downstream. It is important to accurately predict the amount of erosion, peak discharge, and the resulting downstream flow. This paper presents a new model based on the material point method to simulate soil and water interaction and predict failure rate parameters. The model assumes that the dam consists of a homogeneous embankment constructed with cohesive soil, and water inflow is defined by a hydrograph using other readily available reach routing software. The model uses continuum mixture theory to describe each phase where each species individually obeys the conservation of mass and momentum. A two-grid material point method is used to discretize the governing equations. The Drucker–Prager plastic flow model, combined with a Hencky strain-based hyperelasticity model, is used to compute soil stress. Water is modeled as a weakly compressible fluid. Analysis of the model demonstrates the efficacy of our approach for existing examples of overtopping dam breach, dam failures, and collisions. Simulation results from our model are compared with a physical-based breach model, WinDAM C. The new model can capture water and soil interaction at a finer granularity than WinDAM C. The new model gradually removes the granular material during the breach process. The impact of material properties on the dam breach process is also analyzed.


Introduction
Most flood routing of water through water control structures (dams) before the middle 1960s was computed manually. Then, routing software on computers began to replace manual methods. In the 1990s, approximately one hundred sites were selected for more in-depth evaluation, data collection, and data analysis. Additional tests were conducted in the USDA-ARS outdoor Hydraulic Engineering Research Unit (HERU) Laboratory near Stillwater, Oklahoma. Kansas State University cooperated with the USDA to develop software for water control structure design and simulation analysis. Advances in scientific computing capability have enabled more advanced tools to be developed to benefit all types of hydraulic science research.
Interest in the potential flooding due to overtopping has existed for years. There are many simulation processes that can be accurately modeled using computer simulation. Fluid simulation methods are divided into Eulerian and Lagrangian methods. The most popular fluid simulation approach using the Eulerian method is Computational Fluid Dynamics (CFD). CFD for dam break simulation has been extensively investigated [1][2][3]. CFD supports three discretization methods: finite difference, finite volume, and finite element methods. These three discretization methods use distinct mathematical models and governing equations for computation. In our previous study [4,5], we explored the Finite Volume Method(FVM) for dam breach simulation using the open-source toolbox OpenFOAM. Cao and Neilsen [5] simulated the flooding due to overtopping without considering the deformation of the solid dam structure.
Unlike Eulerian methods, there is no numerical dispersion term in the Lagrangian method. Therefore, Lagrangian techniques are useful for simulations when large topological changes occur in the fluid interface. Smoothed Particle Hydrodynamics (SPH) is a well-known method in the computer simulation field. Instead of considering the dam embankment structure as a solid body material as in [4,5], the granular material, like soil, can be represented as either a continuum or a set of individual particles using SPH. SPH is used to simulate large surface flows [6,7] and dam break flows [8][9][10].
In addition, the Discrete Element Method (DEM) is a very popular particle-based system to handle fluid simulation. Often, CFD is coupled with the DEM in many studies and engineering fields including rock slides [11,12], granular flow in water reservoirs [13,14], fluid-particle interaction in dam breaks [15,16]. Other than the CFD-DEM, the coupled SPH-DEM is encountered in many multi-species simulation studies as well. Rungjiratananon et al. simulated sand-water interaction in real time using a hybrid SPH-DEM [17]. Wu et al. [18] and Canelas [19] used the SPH-DEM to simulate multi-phase free surface flow. Lenaerts and Dutré [20] also proposed a dynamics framework for the simulation of both fluid and porous granular material using SPH-DEM.
Another approach, the Material Point Method (MPM), was initially introduced by Sulsky [21] and used to tackle more complex problems. The MPM computes forces using a fixed Eulerian grid, while it also stores information on Lagrangian particles. The fixed grid handles the topology changes like collisions between the material. The MPM stores the information using two distinct representations, and it must be transferred between them. The MPM also avoids re-meshing and storing connectivity between particles. The large deformation landslide problem using the MPM has been widely investigated. Soga et al. [22] introduced MPM to analyze large deformation of landslide mass movements and postfailure behavior. Troncone et al. [23,24] used the MPM to simulate landslide triggered by an increase in the groundwater level or water pressure. Yerro et al. [25,26] applied the MPM to analyze landslides for brittle soils or unsaturated soils. Furthermore, there is much potentiality of the MPM for analyzing geotechnical problems including the landslide that occurred at Oso in USA [27], the landslides of Senise and Maierato in Italy [28,29], and the Sainte-Monique landslide [30]. Dong et al. applied the MPM to solve submarine landslide problems [31,32]. In this paper, we use MPM for overtopping dam breach simulation as shown in Figure 1. To handle multi-species interactions such as water and porous soil [33], Bandara et al. [34,35] introduced models for soil deformation and porous fluid flow using the MPM. Tampubolon et al. [36] simulated the interaction of sand and water mixtures using the MPM and obtained encouraging results. For the porous material property, Klár et al. [37] used the improved Drucker-Prager plastic flow model with volume correction. For the MPM implementation, Arduino et al. [38] and Jassim et al. [39] examined that momentum exchange using the two-grid MPM for the multi-species interaction.
For dam safety simulation, the two most common approaches for the discretization of solids or liquids are Eulerian or Lagrangian. For Eulerian-based methods, the quantities of interest are in fixed locations or fixed grids such as in CFD or CFD-DEM mesh modeling. For the Lagrangian-based methods, the quantities of interest are attached to the materials including the SPH, SPH-DEM, and MPM particle-based methods. Among these modeling methods, in order to simulate dam erosion with multi-species interactions between soil and water, the proposed new dam simulation models use the material point method, which combines the aspects of both types of discretization.
In particular, we construct a new model for the multi-species material point method to allow interaction between soil and water. We allow wet soil transitions from cohesive grains to flowing sediment as water saturation increases. We create an overtopping dam erosion simulation model using multi-species particles and a two-grid material point method. For validation, the MPM simulation model of overtopping dam erosion is compared with available experimental data and the results of other physical-based models like WinDAM [40].

Mathematical Background
Tampubolon and Angeles [36], Atktn and Craine [41], and Borja [42] all considered multi-species using mixture theory. Therefore, the soil and water were modeled as a multi-species continuum using mixture theory.

Conservation Laws
The conservation of mass and the conservation of momentum are calculated individually, and each species obeys the following conservation laws with respect to its own motion. Conservation of mass, represented by the following equation, means that the quantity of mass in the system remains constant over time.
Likewise, conservation of linear momentum is represented by the following equation.
Here, the superscript s is used to represent soil, w is used to represent water, total mass density is ρ = ρ s + ρ w , and total momentum is the sum ρv = ρ s v s + ρ w v w . The velocity v = 1 ρ (ρ s v s + ρ w v w ). After summing the two species, the standard conservation of mass is obtained in (1), and we also note that the conservation of linear momentum for the individual species implies the conservation of linear momentum for the mixture (2).

Deformation Gradient
The deformation gradient represents how deformed the material is locally. It is used to measure how the material has locally rotated and deformed due to its motion [36]. Plasticity is represented by factoring the deformation gradient into elastic and plastic parts as F = F E F P [37].
By factoring the deformation gradient in this way, the deformation history is divided into two pieces, plastic deformation F P and elastic deformation F E . Figure 2a shows the single-phase deformation gradient F and the material properties. Figure 2b displays the material property for multi-phase deformation between water and soil, and the overlapping region is captured through momentum exchange.

Soil Model
Tampubolon and Angeles [36] modified the model of Klár et al. [37] to include cohesive stresses for a porous material and water mixture. Therefore, the amount of cohesion varies with the saturation of water in the mixture. The soil partial stress σ s is defined in terms of the hyper-elastic potential energy density ψ s as: F s is the deformation gradient of soil motion. The F s,E T and 1 det(F s ) terms arise because of the potential energy density in terms of the deformation gradient. As done in finite strain elastoplasticity [43], F s = F s,E F s,P defines plastic flow for porous soil. F s,E represents the compression and shearing, while F s,P represents the sliding and separation. The Drucker-Prager [44] plastic flow and yield condition is used to determine the elastic and plastic deformation gradient. The Drucker-Prager yield criterion is a pressure-dependent model, and the criterion was introduced to deal with the plastic deformation of soils. Bonet and Wood [43] provided the background for elastoplastic constitutive modeling. Dry porous material can be modeled effectively with the assumption made by Drucker-Prager since the yield condition is defined by the constraint that shear stress should be no larger than the compressive normal stress in all directions [36]. For the elastic part of the constitutive behavior for the soil phase, the elastic potential energy density is defined in terms of the logarithmic strain as: Cohesive effects can be modeled by modifying the elastic stress yield condition to be: where F S = UΣV is the singular decomposition of F S and = log(Σ), d = 2, 3 is the spatial dimension, and c C ≥ 0 increases with the amount of cohesion in the material.

Water Model
Water is modeled as incompressible [45] with partial stress: p w is the water pressure. J w is the ratio of the current to initial local volume of material in the water phase.
where k w is the bulk modulus of the water and γ is a term that more stiffly penalizes large deviations from incompressibility.

Momentum Exchange
The momentum exchange term P s , P w for water and porous soil interactions can be considered as a combination of dissipative and reversible interactions [42]. Bandara and Soga [35] introduced the formulation, and they assume: where c E = n 2 ρ w ĝ k is the drag coefficient, n is the soil porosity,k is the soil permeability, and g is the gravitational acceleration. φ w = ρ w ρ is the water volume fraction, and p w is the water pressure.

Cohesion and Saturation
The saturation is estimated by using the percentage of water in the mixture to the total density, which is φ w = ρ w ρ w +ρ s . The soil cohesion is modeled as unitless and different from the typical classical Mohr-Coulomb failure criterion. The soil cohesion varies as a function of water saturation. The cohesion of soil is zero when it is completely dry (φ w = 0). Robert and Soga [46] found an increase to a maximal value c C when the φ w is set to 0.4. In the experiment, a model was created in which the water interacts with wet soil, and also, the wet soil can keep its shape. Soil cohesion is set to the maximal at the beginning, and the cohesion decreases linearly with increasing saturation beyond this point, which means if full saturation for the mixture is reached, the cohesion equals zero and the ratio φ w = 1.

Discretization
The Material Point Method (MPM) computes the forces using a fixed Eulerian grid, but also keeps the information about the Lagrangian particles. Information can be transferred through the grid and particles since the MPM use two distinct representations. There are other approaches in the engineering literature [35,38,39] using the MPM to solve the multi-species problem. Discretization of the continuum equation using two sets of the grid according to Tampubolon and Angeles [36] is followed in the model. One is associated with soil, and the other is associated with water. The superscripts α = s, w are introduced to indicate the corresponding species.

Transfer to Grids
As the process demonstrated in Figure 3, the primary representation of the state is stored on particles in the MPM in Figure 3a. The MPM transfers the mass and momentum of each species to the grid. The mass and momentum of each species are used to compute velocity on its corresponding grid in Figure 3b. The related velocity spatial derivatives are applied, including mass, position, velocity, and affine momentum. The information needs to be transferred between grid and particle representations. Each particle and grid node is assigned a weight that determines how firmly the particle and node interact. If the particle and grid node are close together, the weight should be large. Otherwise, the weight should be small. As demonstrated in Figure 3, information like velocity and momentum exchange can be transferred between grid nodes and particles according to the weights. The mass for grid nodes and particles is initialized, and weight, weight gradient, and kernel are calculated using cubic b-splines kernel for the particle. The first step is the transfer of state particles to the fixed Cartesian grid, then distributing the mass of each particle to its surrounding grid nodes, summing up the mass for grid nodes, and multiplying the surrounding weight and particle mass, as addressed in the equation: Weights are computed based on a kernel as w a,n i,p = N(x a,n p − x a,n i ), where x a,n p and x a,n i are the locations of the particle and grid node locations based on each species a = s, w.
v a,n i = 1 m a,n i ∑ p w a,n ip m a p (v a,n p + C a,n p (x a,n i − x a,n p )) (11) For this step, the velocities for the particle and nodes are initialized, and the grid nodes' velocity is calculated using velocity transfer using C a,n p = B a,n p (D a,n p ) −1 where D a,n p is the cubic b-splines kernel and B a,n p is the affine momentum matrix followed by the APIC transfer [47].

Update Grid Momentum
As shown in Figure 3, the MPM computes the forces using the elastic and plastic deformation gradient and solves for coupled water and soil grid velocities. The explicit grid node force is determined using the particle volume times energy differential, the particle deformation gradient transpose, and the weight gradient.
The forces in the soil and water phases are computed as: Stomakhin [48] mentioned thatx a i can be considered as the position of the grid node i corresponding to species a that is deformed from its original position x a i by the amount of tv a,n+1 i . Then, the grid nodes' velocity is calculated for both soil and water materials using drag coefficient c E and the updated grid nodes' explicit forces, as shown in Figure 3b,c.
(M + tD)v n+1 = Mv n + t(Mg + f (x(v n+1 ))) (14) where M is the mass matrix, v is the velocity matrix for both soil and water particles, and g is gravity.

Update Particles
As illustrated in Figure 3c,d, the MPM updates all particle states, including the cohesion based on saturation, as well as plasticity return mappings. The water determinant J w is monitored. Since the effects of plasticity were not taken into consideration during the implicit solving for momentum,F s,E,n+1 evolves with the grid during the grid momentum update.
For each grid, an indicator function is set for the overlapping region between the soil and water constituents. Then, the soil particle cohesion is computed using the sum of surrounding grid nodes' cohesion.
The cohesion linear function of water saturation is used.
Next, the Drucker-Prager projection is applied, and volume correction treatment is introduced. The plasticity is defined in terms of the singular value decomposition of the deformation gradient. F s,E,n+1 p = U p Σ p V T p , and p = lnΣ p . For the artifact, an extra scalar attribute v s,n cp is added at each time stamp [36]. For the next step, the grid node velocity is changed for collision and friction based on the explicit case. Then, the velocities for particles using the surrounding grid node velocities are updated. The velocity of particles is updated according to: The position of particles is updated according to: Finally, the algorithm updates the deformation gradient and position for both soil and water particles. A summary of all symbols used is given in the glossary in Table 1.

Implementation and Simulation Results
The project was implemented in C++ and Visual Studio 2017, and the implementation idea was inspired by Xia's project [49]. The simulation times shown in Table 2 are measured on a PC with Intel(R) Core i5-3570K CPU, 16 GB RAM, and an NVIDIA GeForce GTX 470 GPU.
Parametric studies were conducted and are shown in Figures 4-6 and Table 2. Figure 4 shows the different selection of the cohesion parameter c Cp from the range 0.002 to 0.008, which will have an influence on the overtopping dam breach erosion simulation. The result shows that with the higher cohesion demonstrated in Figure 4(4), the dam structure is more stable during the overtopping breach for the same timestamp. Especially on the right edge of the dam in Figure 4(4), the dam structure contains less shape deformation than the low cohesive dam structure at the same timestamp compared with Figure 4(3). The soil with a higher cohesion parameter for the earthen embankment homogeneous cohesive dam will maintain the dam structure from the overtopping dam breach.   Figure 5 demonstrates two initial reservoir setups that will have an impact on the overtopping dam breach erosion along the downslope embankment. Figure 5(a1-a4) shows the empty initial reservoir overtopping the breach process at each timestamp. For Figure 5a, the quantified quantified parameters are shown in Table 2, Example 1. The initial reservoir level is zero meters, and the dam dimensions are a crest width of 6 m and a dam height of 3.75 m. Figure 5(b1-b4) shows the full reservoir overtopping breach process at each timestamp. For Figure 5b, the quantified parameters are shown in Table 2, Example 2. For Example 2, the initial reservoir level is 3.25 m, and the dam dimensions are a crest width of 6 m and a dam height of 3.75 m. The result clearly shows that the dam in Figure 5b experiences more erosion at the downward levee during the breach process. Furthermore, the WinDAM C [40] 05-HR-OBA example was applied to the proposed MPM simulation, and the comparison result is shown in Figure 6. WinDAM C [40] is a module used to analyze overtopping earth embankments. Figure 6a is the output from WinDAM C. Figure 6b is the simulation results using the proposed MPM simulation model. The parameter setting and condition for the cohesive dam overtopping failure cases are shown in Table 2. The dimension and dam structure specs are set up in the column 05-HR-OBA in Table 2. WinDAM C simulates the result using the same 05-HR-OBA example. As shown in Figure 6, from Figure 6a's WinDAM output, the dam breach erosion process gradually deforms the right top edge of the dam crust. The physical-based erosion model used by WinDAM is the Hanson/Robinson stress-driven model [50]. Hanson/Robinson's model assumes a vertical cut in a cohesive soil is located at a distance of half the bank height back from the edge. Furthermore, the head cut is made to advance through a series of mass failures driven by erosion at the base. Therefore, the model will erode the cohesive soil on the base surface and form a wedge-shaped notch. In the proposed MPM simulation result in Figure 6b, the head cut during dam overtopping breach occurred progressively. As the saturation of the water and soil mixture reaches the critical point, the right downward levee starts to erode. The right bottom is enlarged due to the erosion and gradual soil movement, and the deformation shape of the dam head cut is not necessarily a wedge shape. Compared with the output from WinDAM C, it shows that the erosion deformation on the levee happens dramatically since WinDAM C uses the Hanson/Robinson stress-driven model, and in the proposed MPM simulation result, the erosion process occurred smoothly due to the MPM calculation based on the interaction between each particle and the surrounding grids.  [52] model formulates the breach evolution along the dam's axis, as well as the head cut formation and migration in the longitudinal section. Hanson et al. [51] provided experimental data for future numerical model development. All these physical models contain mathematical equations and only work on the physical erosion model, and the simulation and visualization effects of erosion processes are not well considered. Nia et al. [8], Pu et al. [9], and Zhao et al. [54] used SPH or the MPM to simulate the dam break process, but they only talked about single-phase water. Tampubolon and Angeles [36] and Yan et al's [55] study about the fluids' and solids' interaction mainly focused on the visualization results. Tampubolon and Angeles' [36] MPM model and simulation achieved excellent results, but they lacked a comparison with real dam breach examples.
Based on the simulation results, the proposed MPM is able to provide extensive details during the erosion process on the levee and gives hydraulic scientists robust animated visualization results for the dam safety study. The proposed MPM also compares the existing physical model results using WinDAM. The MPM simulation model was developed to predict the overtopping breaching process of a cohesive dam. Unlike other large-scale physically-based dam breach models [52,53], the proposed MPM clearly illustrates the flow of soil sediment and water mixture. The performance of the 05-HR-OBA case shows that the proposed model gives rich details when the deformation is happening. Furthermore, the proposed MPM model can perform a parametric study, a synthetic study on cohesion, and compare a physical model.

Conclusions and Future Work
We propose a dam breach simulation framework using the material point method. Based on the simulation results, the proposed MPM can provide extensive details during the erosion process on the downward levee. We also compare the proposed simulation results with the existing physical-based model, WinDAM. For the current implementation, the model we use for the momentum exchange and energy density function is a simplified version, and it can be improved. Furthermore, the simulation time for the proposed MPM method is 60 s, but the simulation time for WinDAM C is around 16 h. For the current stage, the proposed MPM simulates the process in a shorter time. In a future study, we could implement the proposed simulation on a large scale and also try to implement the complete version of the momentum exchange and energy density function. By optimizing the experiment's implementation, we hope we can achieve more experimental accuracy and reduce the run time at the same time.
Overall, our proposed simulation framework provides rich details during the dam breach process. Unlike simply considering the dam break removal part as a triangle in a physics-based model like WinDAM, the proposed framework can capture the water and soil interaction and gradually remove the soil during the breach process when simulating using the existing dam breach example.

Acknowledgments:
The authors gratefully acknowledge the support of the USDA Hydraulic Engineering Laboratory in Stillwater, Oklahoma, USA, for providing the data and WinDAM software used in this study.