Simulating the Hydraulic Heave Phenomenon with Multiphase Fluid Flows Using CFD-DEM

: In geotechnical engineering, the seepage phenomena, especially regarding the hydraulic heave, is one of the most dangerous failure mechanisms related to infrastructural stability. Hence, a fundamental understanding of this occurrence is important for the design and construction of water-retaining structures. In this study, a computational ﬂuid dynamics (CFD) solver was developed and coupled with discrete element method (DEM) software to simulate the seepage failure process for the three phases of soil, water, and air. Specimens were constructed with two layers of gap-graded particles to give di ﬀ erent permeability properties in the vertical direction. More signiﬁcant heave failure was observed for the sample with higher permeability in the upper layer. Special attention was drawn to the particle-scale observations of the internal structure and drag force to study the erosion mechanism. The soil ﬁlled with air bubbles produced a higher drag force in the region below the retaining wall and showed a larger loss of ﬁne particles than the saturated soil, particularly in the initial stages. The results indicate that the impact of air bubbles would accelerate the development of the heave or boiling phenomenon and inﬂuence the stability of the system at an early stage. failure at the initial stages of the seepage process. Compared with the saturated condition, samples with a great number of air bubbles caused a higher ratio of drag force and soil weight adjacent to the wall, increasing the velocity of particles and rapidly inducing severe seepage failure. The results suggest that the existence of air bubbles may exacerbate the process of heave or boiling failure and a ﬀ ect the infrastructure stability at an early stage.


Introduction
Seepage failure is a phenomenon that occurs when the seepage pressure exceeds the dead load of a considered soil bulk. McNamee [1] classified seepage failure into two main categories-piping and boiling or heaving. Based on the incidents reported in the International Commission on Large Dams (ICOLD), Foster et al. [2] found that approximately of 46% failure occurred due to piping through the embankment. Regarding the hydraulic heave, this was first addressed by Terzaghi [3], who showed that it can be captured at retaining wall excavation sites and can impact the stability of the system. Thus, a fundamental understanding of seepage impacts on the erosion mechanism is extremely important in the design and construction of geotechnical practice problems.
The flume test is one of the commonly used experimental approaches for studying the piping phenomenon. Adopting this technique, the piping resistance was investigated relating to the specific gravity, stress state, anisotropic condition of the soil, and the axisymmetric seepage flow [4][5][6][7][8][9]. Tanaka and Verruijt [10] analyzed the seepage failure behind the sheet piles and found that the critical hydraulic gradient depended on the anisotropy of the permeability of the soil. The seepage impact produced a significantly larger local hydraulic gradient around the cutoff walls, as seen from the observations of temporal and spatial progressions of seepage-induced erosion [11,12]. Moreover, the literature also suggested the humidity condition as an important factor in the piping resistance. Wan and Fell [13] found that the erosion resistance was increased for samples in a wet condition with the optimum water content.
In recent decades, numerical analysis has become one of the commonly used approaches for studying the mechanism of erosion. The literature contains serval methods used for investigating multiphase media, including the Eulerian and Lagrangian approaches. The Lagrangian approach

Methodology
This study adopted the open-sourced LAMMPS Improved for General Granular and Granular Heat Transfer Simulations (LIGGGHTS) software to simulate the discrete particles [32], and the Open Field Operations and Manipulation (OpenFOAM) software to mimic the fluid flow [33]. OpenFOAM adopts the finite volume (FV) method to discretize the complex geometries into cell control volumes, and the following continuity and Navier-Stokes equations are assumed to satisfy the conservation of mass and momentum in each control volume: Water 2020, 12, 1077 where u f and α f denote the velocity and porosity of the fluid flow, respectively; ∇ is the gradient operator; and u p refers to the averaged particle velocity. The stress tensor for the fluid phase is calculated as τ = ν f ∇u f , in which ν f is the fluid viscosity and K sl denotes the momentum exchange with the particulate phase.
Regarding the solid particles phase, the linear and angular velocities for each particle are estimated on the basis of Newton's second law of motion. The governing equations of translational and rotational motions of a particle with mass and momentum are as follows [34]: where F c ij and M ij are the contact force and torque for particle i with particle j or with the boundary; F f i and F g i refer to the particle-fluid interaction force and the gravitational force acting on particle i, respectively.
Details of the CFD approach, momentum estimation, and coupling procedures can be found in [35]. The interFoam solver was developed in the CFD-DEM engine to couple the two software programs and perform simulations with the two fluid-phase conditions. For the two fluid phases, the density of the fluid ρ f is calculated as ρ f = α 1 ρ 1 + α 2 ρ 2 , in which α 1 and α 2 are the phase fraction parameters, and ρ 1 and ρ 2 are the corresponding density values of the fluid phases, respectively. The interFoam solver defines the gas phase with a phase fraction parameter that equals 0. If the coefficient of the phase fraction is 1, then this means that the cell is filled with liquid. For the free surface (interphase) boundary, the coefficient should lie between 0 and 1. Regarding the interface, the pressure difference (∆p) is estimated based on the Young-Laplace equation, ∆p = σkn (5) in which σ is the surface tension coefficient andn is the surface norm, while k is the curvature of the surface. In this study, the surface tension coefficient is set as 0.072 N/m, which is used as the value for distilled water. A simulation of a water droplet is taken as an example to validate the software, as shown in Figure 1. In the graph, the red color represents the fluid phase of the water, while the blue color denotes the fluid phase of air. The colors become blurred around the interphase, due to the magnitudes of the volume fraction, which ranges between 0 and 1.
where c ij F and ij M are the contact force and torque for particle i with particle j or with the boundary; f i F and g i F refer to the particle-fluid interaction force and the gravitational force acting on particle i , respectively.
Details of the CFD approach, momentum estimation, and coupling procedures can be found in [35]. The interFoam solver was developed in the CFD-DEM engine to couple the two software programs and perform simulations with the two fluid-phase conditions. For the two fluid phases, the density of the fluid   f  is calculated as , in which 1  and 2  are the phase fraction parameters, and 1  and 2  are the corresponding density values of the fluid phases, respectively. The interFoam solver defines the gas phase with a phase fraction parameter that equals 0. If the coefficient of the phase fraction is 1, then this means that the cell is filled with liquid. For the free surface (interphase) boundary, the coefficient should lie between 0 and 1. Regarding the interface, the pressure difference   p  is estimated based on the Young-Laplace equation, Δ =p σkn (5) in which  is the surface tension coefficient and n is the surface norm, while k is the curvature of the surface. In this study, the surface tension coefficient is set as 0.072 N/m, which is used as the value for distilled water. A simulation of a water droplet is taken as an example to validate the software, as shown in Figure 1. In the graph, the red color represents the fluid phase of the water, while the blue color denotes the fluid phase of air. The colors become blurred around the interphase, due to the magnitudes of the volume fraction, which ranges between 0 and 1.  The mesh size was uniformly distributed in each direction within the boundaries of (0.2 × 0.02 × 0.1 m). Figure 2 displays simulations of water droplets with different numbers of Water 2020, 12, 1077 4 of 16 mesh elements. The water droplet process was hard to capture when the coarse mesh cells were adopted. With reduced mesh cell size, the shape of the water droplets was more accurate. The results demonstrate that the developed software can successfully mimic the two fluid phases. Tests were conducted on an Intel Core i7-8700 K processor, with the simulation time presented in Table 1.
Water 2020, 12, x FOR PEER REVIEW 4 of 16 The mesh size was uniformly distributed in each direction within the boundaries of . Figure 2 displays simulations of water droplets with different numbers of mesh elements. The water droplet process was hard to capture when the coarse mesh cells were adopted. With reduced mesh cell size, the shape of the water droplets was more accurate. The results demonstrate that the developed software can successfully mimic the two fluid phases. Tests were conducted on an Intel Core i7-8700 K processor, with the simulation time presented in Table 1.

Observations of the Hydraulic Heave Phenomenon
The numerical simulations started with sample preparations, using the spherical particles in LIGGGHTS software. Figure 2 shows that tests in which the number of mesh elements was larger than 80, captured the air bubbles inside the water droplets during the simulation process. According to Table 1, the simulation time with 200 mesh elements was approximately 44 times the result with 80 mesh elements, without considering the solid phase. Hence, in this section, the size of the mesh cells in the CFD boundary equaled that adopted in Section 2 with 80 mesh elements. This size was used to model the multiphase fluids accurately and to reduce the computation effort.

Impact of Stratified Soils
Stratified soils are widely observed in natural lacustrine and marine deposits, and in engineering structures such as dam embankments. Gupta et al. [36] found that the permeability of soils constructed from clay in the top layer and sand in the ground would decrease with the increase of the top layer thickness. However, when sand is present in the upper layer and clay is present in the bottom layer, the permeability of the soil skeleton increases with the increasing thickness in the top.

Observations of the Hydraulic Heave Phenomenon
The numerical simulations started with sample preparations, using the spherical particles in LIGGGHTS software. Figure 2 shows that tests in which the number of mesh elements was larger than 80, captured the air bubbles inside the water droplets during the simulation process. According to Table 1, the simulation time with 200 mesh elements was approximately 44 times the result with 80 mesh elements, without considering the solid phase. Hence, in this section, the size of the mesh cells in the CFD boundary equaled that adopted in Section 2 with 80 mesh elements. This size was used to model the multiphase fluids accurately and to reduce the computation effort.

Impact of Stratified Soils
Stratified soils are widely observed in natural lacustrine and marine deposits, and in engineering structures such as dam embankments. Gupta et al. [36] found that the permeability of soils constructed from clay in the top layer and sand in the ground would decrease with the increase of the top layer thickness. However, when sand is present in the upper layer and clay is present in the bottom layer, the permeability of the soil skeleton increases with the increasing thickness in the top. This would have a decisive effect on the rate of fluid flow through porous media and an impact on the seepage failure mechanism. Hence, the simulations adopted two different groups of gap-graded particles to study the impact of permeability in the vertical direction on the development of the hydraulic heave phenomenon. One group had radii of 0.25 and 1 mm for the fine particles and coarse particles, respectively, while the other one had radii of 0.5 and 1 mm for the fine particles and coarse particles, respectively.
This section presents the tests that were performed in the saturated condition with three different specimens. Sample A refers to a uniform condition constructed with only one binary-sized particle, with radii of 0.25 and 1 mm. The other specimens were prepared with two layers of gap-graded particles, denoted as samples B and C, respectively. The condition of sample B was k 1 > k 2 , whereas the condition of sample C was k 1 < k 2 , where k 1 and k 2 denote the permeability coefficient of the top and bottom layers, respectively. This work determined the interface of the layers as the bottom of the retaining wall. Figure 3 depicts the initial configuration of the specimen, which measured 0.02 m in the y direction. Regarding a larger excavation, the seepage failure of the soil in front of sheet pile walls could be treated as a two-dimensional problem. Table 2 presents the information of samples.  The critical characteristic of heave failure is the migration of particles inside the soil skeleton. Hence, Figure 4 portrays the particle trajectory to observe the seepage impact. Particles in the graph represent those located on the upstream side, which experience the largest displacement during the erosion process. The results indicate that a higher hydraulic gradient would result in severe seepage failure. Compared with the results of other samples, sample A produced the least transportation of particles, demonstrating that the erosion resistance increased with the increase of fine particles. Regarding the stratified specimens, sample B induced a higher detachment of particles than sample C. The system showed faster particle speed for the region adjacent to the retaining wall. Moreover, the graphs show that sample A induces a relatively lower particle velocity value. Similar to the  After the samples were prepared, the system was submerged in fluid to observe the heave phenomenon. El Shamy and Zeghal [37] examined the impact of shear moduli on the steady-state response of pore fluid and particles subjected to upward fluid flow. For the moduli of 2.9 × 10 10 and 2.9 × 10 6 Pa, they found that the results of the tests were consistent and only observed minor local variations. Hence, Young's modulus was selected as 5 MPa to reduce the computational effort. Table 3 gives the coefficients of parameters used in the simulations. Tests were conducted by setting different hydraulic conditions at the two sides of the retaining wall. The simulations defined the left side as the downstream side, with the water table set to 0.025 m. Regarding the upstream side, this was selected as two different hydraulic heads (h) to study the impact of the hydraulic gradient, as shown in Figure 3.
The critical characteristic of heave failure is the migration of particles inside the soil skeleton. Hence, Figure 4 portrays the particle trajectory to observe the seepage impact. Particles in the graph represent those located on the upstream side, which experience the largest displacement during the erosion process. The results indicate that a higher hydraulic gradient would result in severe seepage failure. Compared with the results of other samples, sample A produced the least transportation of particles, demonstrating that the erosion resistance increased with the increase of fine particles. Regarding the stratified specimens, sample B induced a higher detachment of particles than sample C. The critical characteristic of heave failure is the migration of particles inside the soil skeleton. Hence, Figure 4 portrays the particle trajectory to observe the seepage impact. Particles in the graph represent those located on the upstream side, which experience the largest displacement during the erosion process. The results indicate that a higher hydraulic gradient would result in severe seepage failure. Compared with the results of other samples, sample A produced the least transportation of particles, demonstrating that the erosion resistance increased with the increase of fine particles. Regarding the stratified specimens, sample B induced a higher detachment of particles than sample C. The system showed faster particle speed for the region adjacent to the retaining wall. Moreover, the graphs show that sample A induces a relatively lower particle velocity value. Similar to the observations of Tanaka [38], sample B produced a significant hydraulic heave phenomenon, depending on the higher magnitude of the particle velocity. The results indicate that the top layer, which has a larger permeability coefficient, would result in a severe heave phenomenon. The system showed faster particle speed for the region adjacent to the retaining wall. Moreover, the graphs show that sample A induces a relatively lower particle velocity value. Similar to the observations of Tanaka [38], sample B produced a significant hydraulic heave phenomenon, depending on the higher magnitude of the particle velocity. The results indicate that the top layer, which has a larger permeability coefficient, would result in a severe heave phenomenon.

Impact of Unsaturated Condition
The above investigations demonstrate that the stratified soil skeleton would influence the occurrence of the heave phenomenon. This section further studies the erosion mechanism by considering air bubbles inside the system. The air bubbles, with the size of 0.005 m, were uniformly generated inside the specimen within the height of 0.0325 m in order to change the water content. The antecedent moisture content was estimated as the total volume of water inside the skeleton divided by the volume of the system. Tests were conducted under h = 0.18 m with moisture contents of 90% and 95%, respectively. Figure 7 gives the multiphase fluid flow behavior for different moisture contents. The results indicate that the fluid flow alters the shape of air bubbles and discharges the air bubbles from the downstream side. Regarding sample C, due to the coarse soil skeleton, the discharge rate of the air bubbles was larger in comparison with the results of sample B, as seen by the behavior at 1.0 s.

Impact of Unsaturated Condition
The above investigations demonstrate that the stratified soil skeleton would influence the occurrence of the heave phenomenon. This section further studies the erosion mechanism by considering air bubbles inside the system. The air bubbles, with the size of 0.005 m, were uniformly generated inside the specimen within the height of 0.0325 m in order to change the water content. The antecedent moisture content was estimated as the total volume of water inside the skeleton divided by the volume of the system. Tests were conducted under h = 0.18 m with moisture contents of 90% and 95%, respectively. Figure 7 gives the multiphase fluid flow behavior for different moisture contents. The results indicate that the fluid flow alters the shape of air bubbles and discharges the air bubbles from the downstream side. Regarding sample C, due to the coarse soil skeleton, the discharge rate of the air bubbles was larger in comparison with the results of sample B, as seen by the behavior at 1.0 s.

Impact of Unsaturated Condition
The above investigations demonstrate that the stratified soil skeleton would influence the occurrence of the heave phenomenon. This section further studies the erosion mechanism by considering air bubbles inside the system. The air bubbles, with the size of 0.005 m, were uniformly generated inside the specimen within the height of 0.0325 m in order to change the water content. The antecedent moisture content was estimated as the total volume of water inside the skeleton divided by the volume of the system. Tests were conducted under h = 0.18 m with moisture contents of 90% and 95%, respectively. Figure 7 gives the multiphase fluid flow behavior for different moisture contents. The results indicate that the fluid flow alters the shape of air bubbles and discharges the air bubbles from the downstream side. Regarding sample C, due to the coarse soil skeleton, the discharge rate of the air bubbles was larger in comparison with the results of sample B, as seen by the behavior at 1.0 s.
Taking t = 0.1 s as an example, Figures 8 and 9 show the particle velocity fields of the two specimens to investigate the impact of air bubbles on the transportation of particles. For each simulation, the peak magnitude of the particle velocity occurred around the retaining wall. Moreover, the particle velocity was increased with an increase in the content of air bubbles. The results were similar to the work of Sakai and Maeda [31], which indicated that soil with lower antecedent moisture content would induce a severe seepage phenomenon at early stages. Taking t = 0.1 s as an example, Figures 8 and 9 show the particle velocity fields of the two specimens to investigate the impact of air bubbles on the transportation of particles. For each simulation, the peak magnitude of the particle velocity occurred around the retaining wall. Moreover, the particle velocity was increased with an increase in the content of air bubbles. The results were similar to the work of Sakai and Maeda [31], which indicated that soil with lower antecedent moisture content would induce a severe seepage phenomenon at early stages.   Taking t = 0.1 s as an example, Figures 8 and 9 show the particle velocity fields of the two specimens to investigate the impact of air bubbles on the transportation of particles. For each simulation, the peak magnitude of the particle velocity occurred around the retaining wall. Moreover, the particle velocity was increased with an increase in the content of air bubbles. The results were similar to the work of Sakai and Maeda [31], which indicated that soil with lower antecedent moisture content would induce a severe seepage phenomenon at early stages.   Taking t = 0.1 s as an example, Figures 8 and 9 show the particle velocity fields of the two specimens to investigate the impact of air bubbles on the transportation of particles. For each simulation, the peak magnitude of the particle velocity occurred around the retaining wall. Moreover, the particle velocity was increased with an increase in the content of air bubbles. The results were similar to the work of Sakai and Maeda [31], which indicated that soil with lower antecedent moisture content would induce a severe seepage phenomenon at early stages.   Figure 10 depicts the spatial distribution of particle displacement to show the heave phenomenon at t = 1.0 s. The results are consistent with the particle velocity behavior. The graphs show that sample B containing more bubbles produces a great number of particles with higher displacement. According to Figure 10b, the influence of moisture content on the hydraulic heave phenomenon was negligible for sample C. The results indicate that the air bubbles in the fine soil skeleton would have a more significant impact on the seepage failure than those in a coarse soil skeleton. Figure 10 depicts the spatial distribution of particle displacement to show the heave phenomenon at t = 1.0 s. The results are consistent with the particle velocity behavior. The graphs show that sample B containing more bubbles produces a great number of particles with higher displacement. According to Figure 10b, the influence of moisture content on the hydraulic heave phenomenon was negligible for sample C. The results indicate that the air bubbles in the fine soil skeleton would have a more significant impact on the seepage failure than those in a coarse soil skeleton.

Discussions on the Development of Heave Failure
Section 3 shows that the fluid flow has a significant impact on the particle velocity and the displacement field. This section aims to interpret the impact of antecedent moisture content on the occurrence of heave failure from particle-scale observations of internal structure and force. Based on the spatial distribution of the particle velocity and particle displacement, sample B produced a relatively significant seepage phenomenon under the impact of multiphase fluid flows, which was selected to investigate the erosion mechanism with air bubbles. Figure 11 portrays the temporal and spatial distribution of the particle displacement for the specimen with 90% moisture content in order to study the development of the hydraulic heave phenomenon. At t = 0.1 s, the seepage impact formed flow paths around the boundary, inducing the occurrence of the boiling phenomenon. An obvious heave phenomenon was observed from t = 0.3 s onwards.

Discussions on the Development of Heave Failure
Section 3 shows that the fluid flow has a significant impact on the particle velocity and the displacement field. This section aims to interpret the impact of antecedent moisture content on the occurrence of heave failure from particle-scale observations of internal structure and force. Based on the spatial distribution of the particle velocity and particle displacement, sample B produced a relatively significant seepage phenomenon under the impact of multiphase fluid flows, which was selected to investigate the erosion mechanism with air bubbles. Figure 11 portrays the temporal and spatial distribution of the particle displacement for the specimen with 90% moisture content in order to study the development of the hydraulic heave phenomenon. At t = 0.1 s, the seepage impact formed flow paths around the boundary, inducing the occurrence of the boiling phenomenon. An obvious heave phenomenon was observed from t = 0.3 s onwards. . .

was selected and taken
as an example to show the detailed redistribution of particles under seepage impact, as shown in Figure 12. The graphs display the obvious detachment of particles and suggest that the erosion The above results indicate that the fluid flow would influence the transportation of particles. A central plane inside the specimen with dimensions of (0.012 × 0.005 × 0.01 m) was selected and taken as an example to show the detailed redistribution of particles under seepage impact, as shown in Figure 12. The graphs display the obvious detachment of particles and suggest that the erosion phenomenon starts with particles adjacent to the wall on the downstream side and then regresses to the upstream side, resulting in the formation of flow channels. The results are consistent with the particle velocity behavior. The fluid flow results in a more serious loss of particles with increasing content of air bubbles inside the soil skeleton.
(c) (d) Figure 11. The spatial distribution of particle displacement for the sample with 90% moisture content at t = (a) 0 . .

was selected and taken
as an example to show the detailed redistribution of particles under seepage impact, as shown in Figure 12. The graphs display the obvious detachment of particles and suggest that the erosion phenomenon starts with particles adjacent to the wall on the downstream side and then regresses to the upstream side, resulting in the formation of flow channels. The results are consistent with the particle velocity behavior. The fluid flow results in a more serious loss of particles with increasing content of air bubbles inside the soil skeleton. The spatial distribution of the soil skeleton demonstrates that the seepage impact induces the redistribution of particles and results in different particle interactions. According to the definition of the internal structure (fabric) proposed by Brewer [39], the results imply that the fluid flow alters the internal structure of the system due to the change of porosity and rearrangement of particles during the erosion process, and influences the erosion resistance and shear resistance of the system. Since The spatial distribution of the soil skeleton demonstrates that the seepage impact induces the redistribution of particles and results in different particle interactions. According to the definition of the internal structure (fabric) proposed by Brewer [39], the results imply that the fluid flow alters the internal structure of the system due to the change of porosity and rearrangement of particles during the erosion process, and influences the erosion resistance and shear resistance of the system. Since the redistribution of particles mainly appears around the boundary, four different representative volume elements (RVE), with the dimensions of (0.01 m × 0.01 m × 0.01 m), were selected to study the characteristics of the internal structure during erosion. These elements were located underneath the boundary, and on the left and right sides of the retaining wall, as shown in Figure 13. . .
, were selected to study the characteristics of the internal structure during erosion. These elements were located underneath the boundary, and on the left and right sides of the retaining wall, as shown in Figure 13.  Figure 14 presents the evolution of the void ratio for the different elements. Similar behavior was observed for the specimen with different antecedent moisture contents. Regarding the element near the free surface on the downstream side, the void ratio was higher at the initial stage due to the coarse soil skeleton. During erosion, the fluid flow triggered the movement of fine particles to fill this element and induced a reduction of the void ratio due to the occurrence of heave or boiling failure. The detachment of particles would affect the loss of particles on the upstream side and form flow Figure 13. The representative volume elements inside the system (unit: mm). Figure 14 presents the evolution of the void ratio for the different elements. Similar behavior was observed for the specimen with different antecedent moisture contents. Regarding the element near the free surface on the downstream side, the void ratio was higher at the initial stage due to the coarse soil skeleton. During erosion, the fluid flow triggered the movement of fine particles to fill this element and induced a reduction of the void ratio due to the occurrence of heave or boiling failure. The detachment of particles would affect the loss of particles on the upstream side and form flow channels, which can increase the porosity on the upper parts of flow paths. Figure 13. The representative volume elements inside the system (unit: mm). Figure 14 presents the evolution of the void ratio for the different elements. Similar behavior was observed for the specimen with different antecedent moisture contents. Regarding the element near the free surface on the downstream side, the void ratio was higher at the initial stage due to the coarse soil skeleton. During erosion, the fluid flow triggered the movement of fine particles to fill this element and induced a reduction of the void ratio due to the occurrence of heave or boiling failure. The detachment of particles would affect the loss of particles on the upstream side and form flow channels, which can increase the porosity on the upper parts of flow paths. To further analyze the erosion mechanism, the temporal changes of particle size distribution inside the elements were tracked and recorded, with the results displayed in Figures 15 and 16. The particles of 0.00025 and 0.0005 m were not contained in elements RVE2 and RVE4, respectively. The impact of air bubbles was slight for the element near the base of the system. The loss of fine particles was more significant around the wall. The graphs imply that the seepage phenomenon starts at the free surface on the downstream side (RVE3) and induces a drop in the content of the 0.0005 m particles due to fluid flow, which results in a rapid increase in the content of fine particles. Compared with the results of other samples, the specimen with 90% antecedent moisture content produced a serious reduction of fine particles in the region below the retaining wall, which caused further movement of particles on the upstream side and increased the content of coarse particles in RVE2, particularly at the initial stages. The results indicate that the impact of multiphase fluid flows would induce rapid occurrence of the heave phenomenon and influence the stability of the system at an early stage. To further analyze the erosion mechanism, the temporal changes of particle size distribution inside the elements were tracked and recorded, with the results displayed in Figures 15 and 16. The particles of 0.00025 and 0.0005 m were not contained in elements RVE2 and RVE4, respectively. The impact of air bubbles was slight for the element near the base of the system. The loss of fine particles was more significant around the wall. The graphs imply that the seepage phenomenon starts at the free surface on the downstream side (RVE3) and induces a drop in the content of the 0.0005 m particles due to fluid flow, which results in a rapid increase in the content of fine particles. Compared with the results of other samples, the specimen with 90% antecedent moisture content produced a serious reduction of fine particles in the region below the retaining wall, which caused further movement of particles on the upstream side and increased the content of coarse particles in RVE2, particularly at the initial stages. The results indicate that the impact of multiphase fluid flows would induce rapid occurrence of the heave phenomenon and influence the stability of the system at an early stage.  The redistribution of particles would also impact the rearrangement of particle interactions. The probability distribution of the contact orientation was examined with different saturated conditions. Figure 17 defines the contact orientation in this study. The RVE1 element showed significant differences in the contents of the particles. This element was selected to study the behavior of contact orientation during erosion, with the results displayed in Figure 18. Compared with the initial results, the seepage impact increased the probability of contact orientation within the span of   The results indicate that the fluid flow in the soil skeleton with lower moisture content induces the rotation of inter-particle contacts in a relatively vertical direction, particularly in the initial stage, which may exacerbate the occurrence of the boiling phenomenon. The redistribution of particles would also impact the rearrangement of particle interactions. The probability distribution of the contact orientation was examined with different saturated conditions. Figure 17 defines the contact orientation in this study. The RVE1 element showed significant differences in the contents of the particles. This element was selected to study the behavior of contact orientation during erosion, with the results displayed in Figure 18. Compared with the initial results, the seepage impact increased the probability of contact orientation within the span of (0 • , 40 • ). The results indicate that the fluid flow in the soil skeleton with lower moisture content induces the rotation of inter-particle contacts in a relatively vertical direction, particularly in the initial stage, which may exacerbate the occurrence of the boiling phenomenon. The redistribution of particles would also impact the rearrangement of particle interactions. The probability distribution of the contact orientation was examined with different saturated conditions. Figure 17 defines the contact orientation in this study. The RVE1 element showed significant differences in the contents of the particles. This element was selected to study the behavior of contact orientation during erosion, with the results displayed in Figure 18. Compared with the initial results, the seepage impact increased the probability of contact orientation within the span of   The results indicate that the fluid flow in the soil skeleton with lower moisture content induces the rotation of inter-particle contacts in a relatively vertical direction, particularly in the initial stage, which may exacerbate the occurrence of the boiling phenomenon.

Spatial Distribution of Drag Force
Section 4.1 shows that the unsaturated condition would influence the internal structure of the system. This section further studied the impact of air bubbles on the seepage failure mechanism based on the force field observations. The simulations investigated the drag force as the interaction force between fluid and particles using the Di Felice model with the following equations [40]:

Spatial Distribution of Drag Force
Section 4.1 shows that the unsaturated condition would influence the internal structure of the system. This section further studied the impact of air bubbles on the seepage failure mechanism based Water 2020, 12, 1077 13 of 16 on the force field observations. The simulations investigated the drag force as the interaction force between fluid and particles using the Di Felice model with the following equations [40]: Re where C d is the drag coefficient and Re p is the Reynolds number; d p is the diameter of the considered particle; and F dp and V p are the drag force and volume of each particle, respectively. The transportation of particles occurred when the seepage force exceeded the weight of the soil bulk. Hence, the ratio between the drag force and corresponding soil weight, shown as F d /W, was investigated to analyze the occurrence of seepage failure. According to Figures 15 and 16, the impact of air bubbles on the movement of particles was more significant in the initial stages. Thus, the temporal-spatial distribution of drag force was studied within the simulation time of 0.4 s. Figures 19-21 portray the results of F d /W in the central region of the system for sample B. For each simulation, the peak magnitude of F d /W was observed underneath the retaining wall at the initial stage. Moreover, the values of F d /W were increased around the boundary, with a reduction of the antecedent moisture content of the system, which caused higher particle velocities and aggravated the movement of particles, resulting in a significant seepage failure phenomenon.
Water 2020, 12, x FOR PEER REVIEW 14 of 16 antecedent moisture content of the system, which caused higher particle velocities and aggravated the movement of particles, resulting in a significant seepage failure phenomenon.   antecedent moisture content of the system, which caused higher particle velocities and aggravated the movement of particles, resulting in a significant seepage failure phenomenon.

Conclusions
This paper studied the occurrence of heave failure under multiphase fluid flow using the CFD-DEM approach. The interFoam solver of OpenFOAM software was developed in the CFD-DEM engine to capture the fluid properties in the two fluid phases. Similar to the experimental results reported in the literature, the specimen with the larger permeability coefficient in the upper layer produced significant heave failure. This study prepared samples with different moisture contents in order to examine the impact of air bubbles on the development of the hydraulic heave phenomenon. The system with increasing content of air bubbles presented a larger particle velocity concentrated around the boundary, which would cause a significant heave phenomenon at the initial stage. Moreover, regarding the stratified soils, the impact of air bubbles in a coarse skeleton on the heave phenomenon is smaller than that contained in fine soil. Special attention was paid to the spatial distribution of the internal structure and drag force to investigate the erosion mechanism. According to the spatial distribution of particles, it is suggested that the fluid flow involving the air phase produces a serious loss of fine particles around the retaining wall and causes rearrangement of the interactions. Based on the contact orientation behavior, the impact of air bubbles produced a higher probability in the vertical direction for the region below the retaining wall, particularly at t = 0.1 s. With the development of seepage erosion, the influence of air bubbles on the migration of particles and the contact orientation was not that significant. The spatial distribution of drag force was further investigated to analyze the occurrence of seepage failure at the initial stages of the seepage process. Compared with the saturated condition, samples with a great number of air bubbles caused a higher ratio of drag force and soil weight adjacent to the wall, increasing the velocity of particles and rapidly

Conclusions
This paper studied the occurrence of heave failure under multiphase fluid flow using the CFD-DEM approach. The interFoam solver of OpenFOAM software was developed in the CFD-DEM engine to capture the fluid properties in the two fluid phases. Similar to the experimental results reported in the literature, the specimen with the larger permeability coefficient in the upper layer produced significant heave failure. This study prepared samples with different moisture contents in order to examine the impact of air bubbles on the development of the hydraulic heave phenomenon. The system with increasing content of air bubbles presented a larger particle velocity concentrated around the boundary, which would cause a significant heave phenomenon at the initial stage. Moreover, regarding the stratified soils, the impact of air bubbles in a coarse skeleton on the heave phenomenon is smaller than that contained in fine soil. Special attention was paid to the spatial distribution of the internal structure and drag force to investigate the erosion mechanism. According to the spatial distribution of particles, it is suggested that the fluid flow involving the air phase produces a serious loss of fine particles around the retaining wall and causes rearrangement of the interactions. Based on the contact orientation behavior, the impact of air bubbles produced a higher probability in the vertical direction for the region below the retaining wall, particularly at t = 0.1 s. With the development of seepage erosion, the influence of air bubbles on the migration of particles and the contact orientation was not that significant. The spatial distribution of drag force was further investigated to analyze the occurrence of seepage failure at the initial stages of the seepage process. Compared with the saturated condition, samples with a great number of air bubbles caused a higher ratio of drag force and soil weight adjacent to the wall, increasing the velocity of particles and rapidly inducing severe seepage failure. The results suggest that the existence of air bubbles may exacerbate the process of heave or boiling failure and affect the infrastructure stability at an early stage.