Seepage failure prediction of breakwater using an unresolved ISPH-DEM coupling method enriched with Terzaghi’s critical hydraulic gradient

This study develops a new numerical simulation model for rubble mound failure prediction caused by piping destruction under seepage flows. The piping has been pointed out as a significant cause of breakwater failure during tsunamis. Once boiling and heaving occur on the mound surface, the piping suddenly propagates in the opposite direction of seepage flow. For the seepage failure prediction, a coupled fluid-soil-structure simulator is developed by combining the ISPH for fluid and the DEM for rubble mounds and caisson blocks. The ISPH, a Lagrangian particle method for incompressible fluids, can simulate seepage and violent flows such as tsunamis. The DEM has been applied for discrete particle and rigid body simulations that include discontinuous deformation, as in the rubble mounds failure and large displacement of the caisson block. ISPH-DEM coupling simulations have already been proposed as a technique for multi-phase flows. Still, the technique cannot reproduce the sudden onset of piping from a stable mound. Two simple assumptions are applied to reduce the numerical cost for the fluid-soil-structure simulators of a breakwater structure composed of a rubble mound and the caisson block. Firstly, each rubble is modeled as an idealized spherical DEM particle with the mean diameter of the rubble. The ISPH particle size is assumed to be the same size as the DEM particle. Under these assumptions, the unresolved coupling model between rubble mound particles and fluid, which obtains the interaction through empirical drag force, should be applied. At the same time, the interaction between the fluid and the caisson block is fully resolved with the spatial resolution with the ISPH and DEM particle size. Our new contribution in this paper is how to model the interaction as an unresolved coupling between seepage flow simulated by ISPH and rubble mound particle modeled with DEM. Our original seepage failure experiment is simulated using the proposed ISPH-DEM coupling simulator. We identified the conventional drag force models as the unresolved coupling model are insufficient to initiate the boiling and piping observed in the experiment. It may be due in one part to excessive averaging of flow velocities caused by unresolved coupling. Therefore, Terzaghi’s critical hydraulic gradient is introduced to initiate the boiling and heaving. Unstable DEM particles, judged by Terzaghi’s critical hydraulic gradient, gradually lose their mass to represent unresolved suspended fine rubble mound particles. Our models qualitatively reproduce the sand boiling and backward erosion in the opposite direction of the seepage flow, as shown in the experiment.


Introduction
On March 11, 2011, The huge tsunami induced by the Tohoku Earthquake(Mw 9.0) caused severe damage to many port structures, such as breakwaters in the coastal area in Japan. A proper understanding of the failure mechanism of caisson-type breakwaters is one of the significant tasks to reduce the damages that following millennium tsunami disasters may bring. According to Arikawa et al. [1], Breakwater failure, in general, has multiple causes; (i) Horizontal force due to the water-level difference between the front and backside of a caisson block, (ii) Scour induced by tsunami over-topping behind the caisson block [2], (iii) Piping destruction induced by seepage flow beneath the caisson block [3]. Developing a numerical simulation model, which can simultaneously cover the above three causes, is desired for the practical design of resilient breakwater systems. The numerical prediction of the failure of soil structures in the case of multiple factors is being developed as a general-purpose technique that can be applied not only to tsunami disasters but also to clarify the causes of embankment failure in recent torrential rain disasters. However, it is challenging to develop such multi-physics simulation tools due to the complexity of the failure mechanism of a caisson-type breakwater on the rubble mound. Therefore, numerical simulations are needed to investigate the above three failure factors individually. As Sassa et al. [4] have shown experimentally, seepage flow in rubble mounds reduces the bearing capacity of the mound. It contributes to scour when a tsunami flows over the caisson. For these reasons, in this paper, the third cause, (iii) seepage-induced piping caused by the bearing capacity degradation of a rubble mound, is solely considered.
In preparation for the catastrophic tsunamis that are expected to occur in the future, there is a need to construct a resilient breakwater system that is tenacious enough to prevent total collapse even if hit by tsunamis and still function in a disaster-prevention ability. Numerical analyses have been conducted by many researchers in the field of breakwaters.Takahashi et al. [3] showed through centrifuge experiments that seepage flows through a mound caused by tsunamis reduce the bearing capacity of rubble mounds and reproduce the characteristics of the reduction in bearing capacity by FEM analysis. Ding et al. [5] studied the stability of the back of a caisson breakwater. They performed a numerical analysis using the DDA: Discontinuous Deformation Analysis [6] to account for the sliding of the caisson under hydrodynamic loading. They also simulated the effect of the armor units' geometry that covers and protects the breakwater mound using DDA. However, it is desirable to consider the displacement of the caisson and armor units and the deformation of the rubble mound that serves as the foundation of the caisson and armor unit to study the damage caused by the tsunami in detail. As Sassa et al. [4] have shown, seepage flow through rubble mounds reduces the bearing capacity of the mound and contributes to scour. Caisson blocks placed on a mound with reduced bearing capacity can also be expected to experience more significant displacement. Therefore, it is essential to select an accurate simulation method that can track the deformation of the rubble mound to design breakwaters. For this purpose, numerical simulation should simultaneously support the discontinuous deformation analysis of rubble mounds composed of gravel and other granular materials and the violent fluid flow analysis. The Lagrange-type particle methods, which directly track the behavior of the particles, are effective for both analyses. A typical method for analyzing discontinuous materials composed of granular particles, such as rubble mounds, is DEM: Discrete Element Method proposed by Cundall and Strack [7]. The DEM can represent discontinuities' behavior by placing spring and damping elements between rigid particles and has been used in numerical simulations in powder engineering and geotechnical engineering. The Lagrange-type particle method can be applied to continuum mechanics and fluid simulation. One of the famous particle methods for continuum mechanics is the Smoothed Particle Hydrodynamics (SPH) [8,9], which was developed for compressible fluids and was already modified for incompressible fluid simulation as the Incompressible SPH (ISPH) [10,11]. The alternative method is the Moving Particle Semi-implicit method [12] for incompressible fluids. Due to their Lagrangian nature, these particle methods are well suited for calculating free surface flows with complex interfacial geometries, such as tsunamis. The recent research has expanded the range of applications to include geomechanical applications [13][14][15], thermal melting problems in solids [16,17], and medical problems [18].
Recent advances in computing technology have made CFD-DEM coupled approaches to particle-fluid mixed flows more practical and efficient. In CFD-DEM, the particle trajectories are usually simulated by the Discrete Element Method (DEM). Then, the fluid part is solved by some CFD methods using the Finite Volume Method(FVM) [19], Finite Difference Method(FDM) [20], and meshfree particle methods [21][22][23] or Lattice-Boltzmann method(LBM) [24][25][26]. Whether the fluid around the particles is resolved, CFD-DEM can be classified into two types: resolved coupling model and unresolved coupling model, as shown in Fig. 1.
The resolved CFD-DEM is an extension of the direct Navier-Stokes solver to the Fluidparticle Interaction without empirical law for the interaction force between fluid and particles. Since the conventional CFD methods are based on the Eulerian description with the fixed meshes or points and moving particle boundary of DEM should be modeled by immersed boundary method (IBM) [19,27] and fictitious boundary method (FBM) [28]. Although the resolved CFD-DEM approach has a clear advantage in its accuracy without additional empirical laws related to the interaction forces, the fixed meshes or spacial points should be smaller than 1/10 of the particle diameters, as discussed [29]. For example, Fukumoto et al. [26] demonstrated the applicability of direct simulations using the Lattice Boltzmann method, an alternative CFD method to the Navier-Stokes equation solver, without macroscopic empirical law to seepage failure of saturated granular soils. This simulation is, however, limited to two-dimensional calculations using soil particles with an average size of 550 µm in a 200 × 100 mm 2 area.
In contrast, the unresolved CFD-DEM [30,31] can reduce the numerical cost in particulate flow simulations with a macroscopic empirical drag force model. In this approach, the spatial resolution of CFD simulation can be expanded to at least three times larger than the particle diameters [32]. In addition, the fluid flow passes through these particles as a permeable flow. The drag force converts from the average permeable flow velocity to the average force acting on the structure. Unresolved models use empirical models to exchange interaction forces, so the fluid around the particles does not need to be resolved finely. Therefore, the unresolved coupled model has an advantage for coupled fluid-particle flow simulations due to the lower resolution requirements on the CFD side.
In most existing unresolved CFD-DEM, DEM is solved as a Lagrangian description, and CFD is mostly solved as an Eulerian description, such as FDM or FVM. In this case, it is necessary to check whether the overall mass conservation is satisfied for the transi-tion from free surface to seepage flow. For this reason, it seems essential to have special treatment of the advection term and the boundary between the seepage flow and the free surface flow, according to Fujisawa et al. [33]. However, SPH-DEM or MPS-DEM, which perform both in a Lagrangian description, track the movement of water particles with individually representative volumes, automatically ensuring that mass conservation is satisfied in the calculation process. This is one of the advantages of using the particle method.
The Lagrange-type particle methods have been applied to research related to breakwater problems. Matsuda et al. [34] evaluated the performance of a caisson breakwater against earthquake and tsunami damage using the SPH. The failure region propagates from a localized area near the caisson in this simulation. It is, however, impossible to represent seepage failure directly because the rubble mound is modeled as an impermeable continuum domain. Iwamoto et al. [35] proposed a WCSPH-DEM coupled simulation for scouring prediction after the overtopping on the top of the caisson block. Using the DEM, they modeled a rubble mound, and fluid flow was simulated by a Weakly Compressible SPH (WCSPH). This simulation deals with scouring failure caused by tsunami currents overtopping the caisson, and it does not care about seepage flow, which is the leading cause of the local piping failure.
In this study, we perform a fluid-soil-structure coupling simulation in which the fluid, rubble mound, and rigid body are all solved using the particle method, the ISPH method for the fluid, and DEM for the solid. In the fluid simulation, among the ISPH methods, we employed the stabilized ISPH method proposed by Asai et al. [36], one of the authors of this paper, which enables stable pressure calculations. The entire rubble mound of the breakwater is modeled using three-dimensional DEM with resolved coupled approach for interactions between fluid and caisson block but with an unresolved coupled approach for rubble mound particles. After we check the performance of the conventional drag force model in the heaving and boiling behavior during the piping destruction, some modifications will be included based on the traditional critical hydraulic gradient criteria proposed by Terzaghi [37]. Figure 2 shows the schematic diagram of the seepage failure simulation of the breakwater mound. Here, the domain to be simulated is defined as f for the fluid-only domain outside the porous rubble mound, m for the domain inside the porous rubble mound, and c for the caisson block. Also, the boundary of the caisson block surface is defined as c . Two simple assumptions are applied to reduce the numerical cost for the fluid-soilstructure simulators of a breakwater structure composed of rubble mound and caisson block. Firstly, each rubble is modeled as an idealized spherical DEM particle with the average particle size of the rubble. The ISPH water particle size is assumed to be the same size as the DEM particle. Under these two assumptions, the unresolved coupling model between rubble mound particles and fluid, which obtains the interaction force through empirical drag force, should be applied. On the other hand, the interaction between the fluid and the caisson block is fully resolved with the spatial resolution with the ISPH and DEM particle size. Therefore, the Darcy-Brinkman-type governing equations, which describe the seepage and surface flows continuously, are solved by the stabilized ISPH method [36] to represent the failure process of the entire breakwater by applying both unresolved coupled and resolved coupling techniques. Our new contribution in this paper is how to model the interaction forces as an unresolved coupling between seepage flow Coupling types for fluid-particle two-phase flows simulated by ISPH and rubble mound particle modeled with DEM. A seepage collapse experiment conducted by one of the authors [38] on the Kamaishi breakwater is selected as the validation test.
We identified the conventional drag force models as the unresolved coupling model are insufficient to initiate the boiling and piping destruction observed in the experiment. It may be due in one part to excessive averaging of flow velocities caused by unresolved coupling. Therefore, Terzaghi's critical hydraulic gradient is introduced to initiate the boiling and heaving. In addition, unstable DEM particles, judged by Terzaghi's critical hydrodynamic gradient, gradually lose their mass to represent unresolved floating rubble mound particles. Our proposed models qualitatively reproduce the sand boiling and backward erosion in the opposite direction of the seepage flow, as shown in the experiment. This paper is organized as follows. "Porous media and particulate flows simulations with the Darcy velocity" section and "Numerical method for rubble mound and caisson block motion based on DEM" section presents the basic mathematical equations for calculating fluid flow using ISPH and the calculation of ground deformation and caisson behavior using DEM. "Preliminary validation test for the unified ISPH through a fixed porous media" section presents a dam-break simulation of flow through a porous media using a coupled unresolved ISPH-DEM coupling model. Finally, a simulation for a seepage failure experiment of a rubble mound of a caisson breakwater is performed using an unresolved ISPH-DEM coupling model in "ISPH-DEM coupled simulation for estimating the seepage fail-ure of rubble mound" section. Simulations were also conducted using the proposed method, in which Terzaghi's critical hydraulic gradient provides a quantitative starting point for failure, changes the mass of soil particles in the failed area, and shows the usefulness of the proposed method.

Porous media and particulate flows simulations with the Darcy velocity
In this study, the entire rubble mound of the breakwater is modeled with DEM particles. Since the caisson block, larger than the soil particles, is not permeable to fluids, a resolved coupled model is applied to treat it as a rigid body and a moving wall boundary for the fluid. On the other hand, the fluid problem inside the porous media should be governed by the Navier-Stokes equations in a microscopic sense. This requires a high computational cost if a coupled resolving model is used, which resolves the details of the inside of the porous media. Therefore, this study applies an unresolved coupling model to the interaction between soil particles and fluid. To reduce the computational cost, a macroscopic fluid analysis is carried out to obtain the Darcy flow velocity, which is the volume-averaged flow velocity. Note that in the unresolved coupled model, the drag forces acting on the particles are calculated from a semi-empirical model. This drag force model for simple shape structures is already established, but the application to the porous media flow composed of fine particles is still discussing. The drag force model for the porous media should be intensely dependent on the porosity and its micro-structures. The drag force models in the geomechanics are nicely summarized [39]. The selection of the drag force models will be discussed again in the next section.
After confirming the performance of the conventional drag model for heaving and boiling behavior during piping destruction, some modifications will be made based on the traditional critical hydraulic gradient criterion proposed by Terzaghi [37].

An unified governing equation for Navier-Stokes flow and permeable porous media flows
Before we discuss the drag force model in the unresolved CFD-DEM coupled problem, we summarize the governing equations for a unified equation for Navier-Stokes flow and permeable porous media flows governed by the extended Darcy law. In the seepageinduced failure simulation of the rubble mound, water modeled as seepage flow goes through the rubble mound and becomes general surface flows outside the mound. The Darcy-Brinkman type equations are applied in this research to unify these equations. According to Akbari et al. [40,41], a set of unified governing equations modeling fluid flows in the saturated porous domain m and outside porous domain f are derived as: where ρ f and g represent the fluid density and the gravitational acceleration vector, respectively(ρ f = 1.0 g/cm 3 , |g| = 981 cm/s 2 ). p is the fluid pressure, and ε f is the porosity. v D is the Darcy velocity, defined as a locally averaged velocity, defined as the intrinsic velocity v f multiplied by the porosity i.e., v D = ε f v f . Here,ρ denotes the appar-ent density, which is given byρ f = ε f ρ f . This relation regarding the apparent density must be employed to guarantee the volume conservation of fluid inside the porous domain with the SPH. In this study, the fluid is assumed to be incompressible, and the porosity is also assumed to be unchanged at each time step. Therefore, the equation of continuity becomes a divergence-free condition of the Darcy velocity, as follows: The inertial coefficient C r (ε f ) to add resistance force related to the virtual mass and the effective viscosity ν E (ε f ) including the kinematic viscosity ν f and an eddy viscosity ν T modeled by the Smagorinsky turbulent model are given as: When the porosity ε f = 1.0, C r = 1.0 and F r = 0 , the unified Eq. (1) is consistent with the original incompressible Navier-Stokes equation. Then, the last term in Eq. (1), resistance force F r (v D , ε f ) for the 'fixed' porous with the porosity, is assumed by: b where a(ε f ) and b(ε f ) are the linear and non-linear coefficients in the extended Darcy law, respectively. These two coefficients are determined from the average diameter D 50 of rubble mound particle and the empirically derived parameters α and β.
The above equations are assumed that the rubble mound domain is fixed with an initial porosity distribution. Most of them can be applied to the deformable and movable rubble mound situation. However, the extended Darcy law is not applicable in such situations, and the resistance model must be modified. In this research, deformation and motion in the rubble mound are represented by DEM spherical particles with a constant diameter equal to the average diameter D 50 of the rubble mound particle. The resistance force should be determined by considering the porosity change referring to the DEM particle motion and the relative velocity of fluid flow and DEM particle velocity. We selected one of the most fundamental and the most cited drag force models proposed by Wen and Yu [42] in the particulate multi-phase flow simulation as the modified resistance force. Wen and Yu proposed a drag force model with two different low and high porosity states. In this study, the resistance force is assumed by Eq. (6) in the low porosity seepage state, and the drag force in the high porosity suspended state is given as the equation proposed by Wen and Yu. The resistance force F r (v D , ε f ) in the Eq. (1), defined by Eq. (6), is replaced by this modified resistance force F * r (v r , ε f ): Re where C d represents the drag coefficient depending on the Reynolds number Re. μ f represents the fluid viscosity. Note here that the first equation in Eq. (9) has the same function shape as Eq. (6), but the variable was changed from the Darcy velocity v D to a relative velocity v r := v f − v s . The modified resistance force can treat a more comprehensive condition ranging from flow through the mound as a seepage flow to flow with gravels suspended in the mound.
For the seepage flow problem on the porous structure, the coefficients of the resistance model are often determined as α = 150 and β = 1.75 based on the Ergun model [43], for example, by Larese et al. [44] and Peng et al. [45]. However, These empirical coefficients α and β have varieties in previous studies, and the determination should be carefully considered, as discussed by Losada et al. [39]. In this study, a preliminary study to determine this coefficient is carried out in the dam-break simulation through the porous media in "Preliminary validation test for the unified ISPH through a fixed porous media" section.

Spatial discretization of fluid using the SPH based on the Lagrangian description
The SPH is based on the Lagrangian description and divides a continuum into a set of discrete particles. These particles have a spatial distance, known as the smoothing length, over which a kernel function smoothes their properties. In a general SPH analysis, the physical quantity of the target particle i can be obtained by summing all the neighboring particles j, which exist within the range of the kernel area. Then, the contribution of each neighboring particle is weighted according to its distance from the target particle using a smoothing function. Finally, the motion of the particles can be described by the interpolated physical quantity based on a simple algorithm. Here, the basic concept of the SPH method is described by referring to Asai et al. [36].
A spatial discretization using scattered particles is summarized. First, The scalar function φ(x i , t) of the field represented by the domain can be expressed as a volume integration using the Dirac delta function δ at any point x in the domain . The Dirac delta functions δ is a sharp function that is ∞ at the origin and 0 elsewhere. It has the property of becoming 1 when integrated over a domain with φ(x) = 1 as the function to be integrated.
Instead of this sharp delta function δ, the SPH method deals with approximations using a smooth kernel function: where W is a weight function called the smoothing kernel function. In the smoothing kernel function, r ij and h are the particle distance between neighbor particles and the smoothing length, respectively. The distance r ij is simply the length of the relative coordinate vector r ij = x j − x i (the relative distance r ij = |r ij |). In SPH discretization, the differential operators (e.g., the gradient ∇φ, the divergence ∇ · φ and the Laplacian ∇ 2 φ) acting on a target particle i are evaluated using the neighboring particles j within an effective radius r e . Let S i be defined as the set of neighbor particles j of the target particle i, as follows: where N SPH is the number of SPH particles (serial ID number). For an SPH simulation, the volume integration in Eq. (15), ∇ · φ, ∇φ and ∇ 2 φ can be approximated as: The subscripts i and j indicate the positions of a labeled particle; for example, m j and ρ j mean the representative mass and density of a neighbor particle j, respectively. Note that the triangle bracket · indicates the SPH approximation of a particular function. Furthermore, note that the two sets of expressions of the gradient models have various properties that can be converted to each other analytically. η in Eq. (21) is the parameter to avoid division by zero and is defined by the following expression, i.e., η 2 = 0.0001(h/2) 2 . In this study, the kernel function adopts the cubic spline curve with the smoothing length h = 1.2r 0 (r 0 being the initial distance between SPH particles) and an effective radius r e = 2h.

Discretization of the unified governing equations with the ISPH method
In the ISPH: Incompressible SPH method, the unified governing equations, Eqs. (1) and (2), are first discretized in time by following the projection method based on the predictor and corrector scheme. Then, the spatial Discretization is implemented for the time-discretized equations.
To begin with the time discretization, v D at n + 1 step is written as: where v * D is the predictor term calculated explicitly from the physical quantities at n step, while v * D is the corrector term which is implicitly given from the physical quantities at n + 1 step to correct the predictor term. Based on the projection method, Eq. (22) can be separated as: The pressure p n+1 in Eq. (24) Next, the PPE is discretized into the particle quantities using Eqs. (18) and (21) as: Here,ρ f denotes the apparent fluid density that is defined by the following equation, i.e., ρ f = ε f ρ f . This concept must be employed to guarantee the volume conservation of fluid inside the rubble mound. The representative volume changes as the fluid density changes depending on the porosity. Consequently, in an entire analysis domain, including outside and inside of the rubble mound, the total amount of the fluid volume can be theoretically conserved.
Using the above SPH approximate operators, the PPE in Eq. (25) is described as: This derivation of the PPE is well-known as the formulation under the divergence-free condition. Next, the velocity at n + 1 step v n+1 D is determined from the obtained pressure calculated through Eq. (28) by using Eq. (22) and Eq. (24). The pressure gradient term needs to be determined to implement this updating procedure. By referring to Eq. (20), it can be written as: Finally, the position of a particle is updated from the updated velocity: A particle's intrinsic Lagrange velocity v f is the Darcy velocity v D divided by the porosity ε f . In the next step, this intrinsic velocity v f is used to calculate the particle position r n+1 i .

Stabilized ISPH method for the unified equation
The heterogeneity of the particle distribution directly affects the accuracy of the SPH discretization. Therefore, the stabilized ISPH method proposed by Asai et al. [36] is arranged for the unified equation. In the incompressible Navier-Stokes equation, numerical density ρ f i should be preserved with the true fluid density value ρ f . On the other hand, the numerical density in the pressure Poisson equation using the Darcy-Brinkman type unified equation must beρ fi = ε n fi ρ f , consistent with the apparent density following the surrounding porosity distribution. The discrete PPE (28) is slightly modified by adding a stabilization term, which tries to maintain the reasonable density as: where γ (0 ≤ γ ≤ 1) is called the relaxation coefficient and is generally set to be much less than 1.

Porosity and apparent density/volume
In this study, the ISPH and DEM are coupled through drag forces in an unresolved coupling scheme. Therefore, the porosity ε f , which is the volume fraction of voids in a unit fluid volume, is necessary to calculate the space-averaged fluid densityρ f and Darcy velocity v D . Based on the concept of interpolation approximation of the SPH method, the porosity is obtained from the total volume of wall particles and DEM particles contained within the influence zone of the target particles as follows: where V s represents the volume of DEM particle. The sets S c i and D i , which are subject to the sum of the numerator-denominator, are the sets of all DEM particles and the sidewall/caisson particles of the SPH contained within the influence area of the target particle i, respectively. Figure 3 shows that the fluid deforms to avoid solid particles in the gap. On the other hand, the unresolved coupling allows for the overlap of solids and fluid, so the apparent increase in fluid volume and an apparent decrease in fluid density must be considered. Therefore, under the constant mass of fluid particles, spatial averaging is performed using porosity ε f , and the apparent densityρ f , volumeV f , and Darcy flow velocity v D are defined as follows:

Numerical method for rubble mound and caisson block motion based on DEM
The rubble mound's deformation and a caisson block's motion are represented by rigidbody calculations using DEM. Although the caisson block is cracked and damaged by the tsunami's impact, the caisson block is assumed to be a rigid body without deformation in this study. The rubble mound is represented by rigid spherical particles frequently used in general DEM simulation. The equations of motion of the rigid elements consider the

Defomation analysis of the rubble mound
The equations of motion for the translational and rotational directions of a rubble mound particle in the fluid can be expressed as follows: where v s and ω s represent the velocity and the angular velocity. Furthermore, m s , V s and I s represent one DEM particle's mass, volume and moment of inertia. Note that the contact forces f c are summed only for the set D con i of particles j in contact with the particles i. As in the conventional DEM, the contact is judged to have occurred if there is a positive amount of overlap δ n = d s −|x ij |. The general spring-dashpot model calculates the contact forces f c between the rubble mound particles or the caisson block. The contact forces are divided into normal force f cn and tangential contact forces f ct . The frictional forces are obtained according to Coulomb's friction law in the tangential direction.
where k, η, μ s , e and δ represent the spring stiffness, the damping coefficient, the friction coefficient, the restitution coefficient and the contact overlap, respectively. The subscripts n and t denote the normal and tangential components. e n and e t represent the unit normal and tangential vector. v rs is the relative velocity between soil particles. The fluid forces acting on the soil particles include the buoyancy force ∇pV s caused by the pressure gradient and the drag force f d , which depends on the relative velocity v rs between the neighboring fluid particle and soil particle. This drag force is the interaction force in coupling with the fluid and is calculated as the reaction force acting on a single particle in the opposite direction of the fluid resistance force described in Eq. (6).
Rubble mounds are packed with spherical particles, and the torque m c that produces rotational motion is generated only by tangential contact forces f ct . Spherical particles cannot represent the interlocking effects resulting from the actual uneven shape of the particles. Therefore, rolling friction m r is introduced to suppress rotation as follows artificially: where μ r andω represent the coefficient of rolling friction and the axis vector of rotation.

Tracking of the caisson block motion
The caisson block motions as the rigid body give moving wall boundary conditions for the fluid. Particles are placed on the surface of the caisson block that acts as both DEM contact spheres and ISPH wall particles. The contact forces with the rubble mound particles and the hydrodynamic forces acting from the fluid are calculated simultaneously. The equations of motion for the translational and rotational directions of the caisson block are shown below: where M, I, V and represent the mass of the caisson block, the tensor of inertia, the velocity and the angular velocity. F c , M c , F f and M f represent the contact forces, contact torques acting on the caisson block, the hydrodynamic forces and its torques, respectively. These are then summed up to the forces acting on the surface SPH/DEM particles of the where x j , X g and x cj represent the position vector of the particle j, the center of gravity of the caisson block and the position vector of the contact particle, respectively. a j and n j are required to calculate the hydrodynamic forces representing the area represented by the particle j and the normal outward vector. The rotational motion is calculated according to the quaternion formulation. Due to the explicit time integration with a high penalty parameter, DEM requires sufficiently tiny time increments to maintain the stability condition. In the DEM and SPH coupled simulations, the time increment is generally governed by the DEM stability condition. To overcome this difficulty, Asai et al. [46] successfully apply the ETI: Energy Tracking Impulse method, which is an unconditioned stable scheme for multiple bodies contact simulation, to the fluid-structure interaction simulation. However, the ETI in DEM may require many iterations to simultaneously satisfy all the constraint conditions at each DEM particle. Therefore, in this study, the multiple time-step algorithms are adopted, in which each time increment ( t SPH > t DEM ) is set for ISPH and DEM, and only the DEM calculation is solved with fine time increments. Furthermore, it was assumed that the interaction force with the fluid is constant within the small loop of the DEM.

Dambreak flow passing through a fixed porous layer
First, as a preliminary validation experiment for the ISPH of the unified equation, we select the water column collapse through a porous media experiment performed by Liu et al. [47] as shown in Fig. 4. In this example, we can assume that no DEM simulation was performed and that the porous media was fixed during the simulation. In this experiment, a 25 cm-high water column collapsed due to gravity and passed through a porous layer.
The porosity and mean diameter are 0.49 and 1.59 cm, respectively, placed in the center of the tank to the right-hand side. The conditions for the fluid analysis based on the ISPH method were initial particle spacing of 0.5 cm and time increment of t SPH = 10 −4 s. To show the advantages of the stabilized ISPH method described in Eq. (31), we performed the same simulations for relaxation coefficients γ of 0.00 and 0.01. In addition, as a preliminary study of the drag model coefficients α in Eq. (7) and β in Eq. (8), simulations were carried out here with the combination of α = 150, β = 1.75 and α = 150, β = 0.50. Furthermore, in order to indicate the superiority of the stabilized ISPH method, as shown in equation (31), similar simulations were carried out for stabilization parameters of γ = 0.00 (normal ISPH) and γ = 0.01 (stabilized ISPH). Figure 5 a) and b) shows the pressure distribution for (α = 150, β = 1.75, γ = 0.01) and (α = 150, β = 0.50, γ = 0.01), respectively. The experimentally observed free surface is plotted with red circles in Fig. 5). The smoothed pressure distribution can be evaluated using stabilized ISPH with γ = 0.01. However, compared to the experimental results, the fluid motion is relatively slow, especially after t = 0.4 s. This is because the coefficients α and β in the drag model do not match this experimental result, and the appropriate resisting force is not applied. By calibrating these two parameters, the nonlinear term in the drag force is reduced by decreasing β to obtain the optimal simulation results shown in Fig. 5b). Finally, we conducted a simulation with the same set of optimal parameters alpha and beta, excluding our stabilization term γ = 0. Figure 5c shows this simulation result. The stabilization term plays an important role in these simulations to satisfy the total water volume and to evaluate the smoothed pressure distribution.
To confirm the effect of this stabilization term, the apparent fluid density distribution expressed by Eq. (35) in the simulation results with each parameter set is shown in Fig. 6. For γ = 0.0 in the Fig. 6c), the apparent density cannot keep the value corresponding to the porosity and oscillates significantly. As a result, not only does it decrease the total water volume over time, but it also causes fluctuations in the pressure distribution as shown in Fig. 5. On the other hand, as shown in Fig. 6a, b, the simulation using the stabilized ISPH method with γ = 0.01 shows that the apparent fluid density in the porous media is consistent with the porosity, resulting in no volume reduction. These results show that the stabilized ISPH method can achieve excellent volume conservation and stable pressure calculations even in fluid analysis where there is a transition between the free surface and seepage flow.
The last part of this section discusses the role of drag and parameter settings. As summarised in the review paper by Losada et al. [39], many forms of drag model equations have been proposed, and their coefficients are set in a wide range (α = 150 ∼ 1000, β = 0.2 ∼ 2.0) depending on the porous structure under consideration. In a similar problem using the same coefficients (α = 150, β = 1.75) as in Fig. 5, some papers have shown better results using the high-accuracy MPS method, such as Harada et al. [48]. Coefficient forms that do not treat the drag coefficient as a constant and use the porosity as a variable have also been proposed [49,50], which could be used to set parameters for a wider range of application problems automatically. Incidentally, the coefficients of this drag model are close to those obtained by Korzani et al. in their coupled simulations of sheet pile piping problems (α = 150, β = 0.40). The same coefficients (α = 150, β = 0.50) were used in the following numerical simulation of the seepage collapse of a breakwater, although it is difficult to determine a drag model that can be commonly applied to all cases, and further study is needed.

ISPH-DEM coupled simulation for estimating the seepage failure of rubble mound
The main validation for the seepage failure prediction is given here. Kasama et al. [38] conducted a series of hydraulic model tests on a scale of 1/100 for the Kamaishi Harbor Mouth Breakwater. In this experiment, seepage flow is generated in the rubble mound due to the water head difference in water level between the seaside and inner harbor side. As shown in Fig.7, sand erosion occurred in the rubble mound on the harbor side when the water level difference was 145 mm. The caisson block's sliding, rotation and settlement occurred due to the horizontal hydrodynamic forces caused by the water level difference.
In this study, the ISPH-DEM coupling method, given in the previous section, is first applied to represent this experimental test [38]. Then, the additional failure criteria based on Terzaghi's critical hydraulic gradient will be introduced to improve the accuracy of the seepage failure prediction. The particle simulation model is shown in Fig. 8, and simulation conditions are summarized in Table 1. The initial conditions are created by randomly packing rubble mound particles under gravity. It is desirable to model the rubble mound according to the particle size distribution of the actual mound. However, the number of particles would be enormous, and the computational cost would be unfeasible if even the finest particles were modeled. Therefore, this rubble mound is formed by DEM particles with average size (single particle size) to represent the overall collapse behavior. The boundary condition of the sidewall is set as the slip condition in the fluid analysis and zero friction in the DEM.

ISPH simulation of fluid flow passing through a fixed porous mound
We assumed the fixed porous media for the rubble mound like the above example to investigate pressure and the piezo water head distribution before we introduced the DEM incorporated with the ISPH for the movable mound example. Figure 9 shows the pressure distribution at a hydraulic head difference h of 145 mm and the distributions of the where z is the height from the initial water level at the harbor side. Moreover, the hydraulic gradient of each water particle is defined as the gradient of the piezo head H As shown in Fig. 9, the stabilized ISPH method [36] produces a transition of smooth pressure and piezo head distributions throughout the interior of the rubble mound. Fig. 10 also shows the piezo head distribution obtained on the rubble mound's DEM particles at the hydraulic head difference of 145 mm, and this distribution is normalized by the hydraulic head difference of 145 mm, i.e., h t = H i / h. In this section, the DEM particles are fixed at their initial positions. The porosity values on the DEM particle are evaluated as in Fig. 11a. The range of evaluated porosity is from 0.3 to 0.6, and the value is lower than the average value (0.493) in the experimental setting. Particle packing with spherical particles of average size may give a lower porosity than with random particle shapes. Therefore, we put a limitation of the evaluated porosity value on the DEM particle with the average porosity value (0.493) in the experimental setting. Applying the limitation, the DEM porosity distribution becomes almost homogeneous, as shown in Fig. 11b. The other parameters for the Darcy-Brinkman type unified equation are given only from physical properties, such as porosity and the mean particle size, without any calibrated parameters. Under these assumptions, the evaluated piezo head agrees with the experimental data shown in Fig. 10.

ISPH-DEM simulation of seepage failure in caisson-type breakwater
Next, all DEM particles constituting the rubble mound and caisson block were set to be movable, and the ISPH-DEM coupling analysis was performed to simulate the seepage failure of the rubble mound. The simulation conditions are the same as in the fixed mound simulation case listed in Table 1. The drag forces in the rubble mound give the deformation of the mound, and the caisson block is also moved by receiving the fluid force acting on its surface. Deformation of the mound and the caisson block and water The rubble mound shows a small deformation by receiving the drag force, but we cannot represent the caisson block's soil heaving behavior and rotation as shown in Fig. 7. The rubble mound is modeled in this simulation with single-size DEM particles to reduce computational cost. Therefore, smaller particles than the DEM particle cannot directly affect this simulation. In addition, the DEM particle movements are modeled as an unresolved coupling framework with an empirical drag force model. Our simulation results gave a good piezo water head compared with the experimental measurement in the preliminary validation with the fixed mound. However, the selected drag force due to the averaged Darcy velocity may need to be modified, or the contribution of particles smaller than the average size may need to be included. Cheng et al. [51] also pointed out the limitation of the unresolved coupling model. They conclude that unresolved coupled simulations with the conventional drag force may be relatively smaller than the actual averaged force, thus overestimating the critical hydraulic gradient of the subject ground. In order to represent in detail localized failures, such as the experimentally confirmed sand boiling, we propose a simple modification of the drag force modeling to consider the contribution of the unresolved fine particle effects as a reduction of mass in the DEM particles.
In geomechanics, the critical hydraulic gradient is often applied to predict the initiation of piping failures. Piping is considered to occur when the hydraulic gradient of the seepage flow through the mound exceeds the critical hydraulic gradient. In this study, piping failure at each soil DEM particle is judged by comparing their hydraulic gradient i s with the critical hydraulic gradient i c proposed by Terzaghi [37]. The hydraulic gradient at DEM particle is given as: Fig. 10 The comparison of piezo head distribution between proposed simulation and the experimental data of Kasama et al. [38] Here, G s is the specific gravity of the mound material, which is given by G s = ρ s /ρ f . Moreover, e is the void ratio defined as e = ε s /1 − ε s . This study's G s is set to 2.63. When the ground is horizontal, the vertical component of the hydraulic gradient vector i s is compared with the critical hydraulic gradient, and the judgment is made as follows: Breakwater mounds are generally trapezoidal with a slope, and their shape changes over time due to piping failure. Therefore, the piping criteria for horizontal surfaces should be generalized to any shape with a slope. In this study, the normal outward vector orthogonal to the rubble mound surface is calculated for each DEM particle, and the inner product of the normal vector and the hydraulic gradient vector is used to determine whether the .
The unit normal vector n s is updated every step according to the following equation using a set D stable i that includes only stable DEM particles j in its neighborhood (j ∈ stable m ): After judgment with a generalized piping failure criteria defined in Eq. (61), we assume that the mass of the DEM particles may decrease as the unresolved fine particles flow out. In our experimental setting, the inside of the mound was given the experimentally obtained density ρ s 0 = 1.86 g/cm 3 and mean porosity ε s,ave = 0.493. The DEM density value after the piping failure was corrected to linearly decrease to ρ s,min = 1.1 g/cm 3 following their porosity ε s , according to the following formula: The criteria for the minimum density ρ s,min is the same value (ε si = 0.80) as the criteria of the drag force defined in Eq. (44) for suspended particles. The density of DEM particles, which once decreased, is not allowed to recover.
The same numerical example is solved here with the modified drag force model with Terzaghi's critical hydraulic gradient. Figure 13 shows simulation results with the proposed model. In Fig. 13a, red and blue particles indicate unstable and stable particles. Figure 13b shows the reduced density of unstable particles. As in the experiment, the updated results show that piping failure begins beneath the caisson block at the harbor side. The localized sand blowing from the failure point is qualitatively expressed by decreasing the mass of fine particles near the mound surface, assuming they are lost due to the piping failure caused by seepage flow. Furthermore, backward erosion, in which the failure proceeds in the opposite direction to seepage flow toward the seaside, is expressed qualitatively from where the sand boiling occurs. The clockwise rotation of the caisson block and deformation of the mound is similar to those observed in the experiment. Notes that the mass and volume of the rubble mound do not conserve because the fine suspended particles are never traced in this simulation.
Through numerical tests with the conventional drag force model and the proposed one, we confirmed that the conventional unresolved simulations could not reproduce the internal destabilization of the ground caused by the detachment and movement of fine particles in the soil skeleton due to seepage flow. We proposed a simple modification in the drag force model to reduce mass after the piping failure judgment with Terzaghi's critical hydraulic gradient. Although the proposed method still has room for improvement, such as satisfying the mass and volume conservation for suspended fine particles, it can demonstrate the backward erosion, sand boiling, and caisson motion observed in our piping failure experimental tests.

Conclusion
This study simulates the piping failure of a breakwater rubble mound caused by a seepage flow using an unresolved ISPH-DEM coupling method. In the proposed method, Darcy-Brinkman-type governing equations were employed to describe surface and seepage flows in a unified manner with porosity as a parameter. The interaction force between the fluid and soil particles is assumed to be a drag force that depends on the relative velocity of the two, and the unresolved coupling method is used to reduce the computational cost. On the other hand, for caisson blocks that do not allow water to pass through to the inside, a resolved coupling is implemented in which ISPH wall particles are placed on the surface and treated as a moving wall boundary. We performed a numerical analysis corresponding to the seepage collapse experiment by Kasama et al. [38]. The simulation qualitatively reproduced the overall failure morphology of the breakwater, including deformation of the rubble mound due to seepage flow and sliding and rotating of the caisson block due to horizontal forces. However, it wasn't easy to represent the localized sand boiling observed in the experiments. This problem is that the unresolved coupling model implemented in this study excessively averaged out the locally large velocity flows that are the starting point for blowing sand. In addition, we considered that the modeling of rubble mounds using particles more significant than the average particle size from the viewpoint of the computational cost was a problem because it could not represent the instability of the ground due to the detachment of fine particles. Therefore, a piping judgment based on the critical hydraulic gradient was performed, and particles judged to be unstable were assumed to have been washed out by seepage flow from the surrounding fine particles. A mass reduction process was applied to the particles. This method represented backward erosion propagating beneath the caisson block in the opposite direction to the seepage flow and reproduced the caisson settlement caused by the piping.
As fundamental issues in our SPH simulation, the method for setting drag coefficients and improving the SPH method need to be discussed in-depth in the future. In particular, to solve detailed flow inside the ground, it is necessary to improve the method to represent complex flows, including negative pressure regions. For these problems, it is necessary to improve the accuracy of the gradient model [52] and the Laplacian model [53,54], which are also problems of the conventional SPH method.
As this study shows, internal erosion that occurs in microscopic areas of the ground, such as soil particle loss due to seepage flow, significantly influences the macroscopic stability of the ground. Therefore, analyzing the actual particle size distribution, including even the finest soil particles, is desirable, using DEM, and calculating detailed pore flow using resolved coupling simulation. In addition, we aim to develop a simulator that can efficiently solve internal erosion in the ground by implementing a semi-resolved coupling [51,55] that combines a resolved coupling and an unresolved coupling. If this semi-resolved coupling method becomes available, it should be able to efficiently and in more detail determine the interaction between flow and soil particles in the ground. It should apply to the numerical analysis of large-scale structures such as breakwaters.