Mixture theory-based SPH model for submerged landslide

. A novel SPH model aimed at solving the coupled water-soil problems is proposed based on the mixture theory. This method is featured with the spatially overlapped dual continua for both fluid and solid phases. The water phase is modeled as a weakly-compressible Newtonian fluid, and the soil phase is modeled using an elastoplastic constitutive model. The benchmark problem, fully submerged soil subjected to gravity , is examined to validate this SPH model. Finally, a submerged landslide is simulated to demonstrate the capability of the proposed SPH model in solving the dynamic soil – water coupling problems.


Introduction
Submerged landslide, widely observed in the coastal region, is destructive to underwater facilities such as submarine cables and pokes critical threat to the residents and properties due to the landslide-generated tsunami. As a typical example of coupled soil water dynamic problems, the submerged landslide has drawn wide attention from the academic community [1][2][3]. Experimental studies have shed significant light on this problem but are often economic-and time-consuming compared to numerical methods. A natural way of investigating the watersoil coupling problem is to directly consider the flow in the pore network and the motion of the soil grains at microscale, as in the fully resolved coupled lattice Boltzmann and discrete element method (LBM-DEM) [4], or coupled smoothed particle hydrodynamics and DEM (SPH-DEM) [5]. However, these promising tools are powerful for microscopic scenarios, but highly computationally expensive for engineering problems.
A feasible solution is to study the submerged landslide dynamics within the unified continuum framework. Two schemes are popular in the previous studies. One is the multi-phase model, which treats the pure water and sediments as two immiscible materials [6]. This method is straightforward and computationally efficient, but limited in applications due to the neglect of soil water interactions. The second scheme, based on the mixture theory, is featured with the spatially overlapped continua representing the fluid and solid phases [7]. This method is relatively more efficient than the coupled CFD-DEM procedure because of its continuum mechanics nature. Besides, it can consider the interaction between the fluid and solid phases, thus adopted in this study.
Smoothed particle * Corresponding author: wei.wu@boku.ac.at hydrodynamics is chosen as the numerical approximation scheme due to its excellent performance on dealing with large deformation problems [8].
The remaining part of this paper is outlined as follows. The theoretical framework will be presented in Section 2. Subsequently, in Section 3, the newly proposed method will be examined with a benchmark problem. Section 4 presents the application of the SPH model on the dynamics of a submerged landslide and its induced wave.

Methodology
Based on the mixture theory, the following identify holds everywhere for a saturated soil, s + f = 1, (1) where, denotes the volume fractions; the subscripts, f and s, in the physical quantities refer to fluid and solid phases. This practice is applied throughout this paper unless otherwise stated. By definition, we can obtain the following expressions between the intrinsic density, ̃, and the apparent counterpart, , s = s̃s , f = f̃f . (2)

Conservation equations
The Lagrangian form of the mass conservation equation for the solid phase reads where d(•)/d = ∂(•)/d + s • (•) is the material derivative, the symbol ∇denotes the gradient operator and · indicates the inner product. As the soil grains are incompressible, the intrinsic density of soil ̃s is constant. Therefore, the evolution of the volume fraction of the soil phase is as follows (6) in which is the dynamic viscosity of water, and Dc refers to the characteristic length of the soil grains. For coarse-grained soils, the characteristic length can be chosen as D50, which is the grain size at 50% in the grain size distribution.
=150 and =1.75 are adopted in this work.
The mass and linear momentum conservation equations for the fluid phase are [10]  f Note that the shear stress of fluid is neglected in Eq. (8) as the interstitial fluid is water and its viscosity is marginal.

Constitutive models
The fluids are assumed weakly compressible dictated by an equation of state (EOS) such that where ̃f 0 is the reference intrinsic density and is set as 1000 kg/m 3 for water, γ is a constant usually taken as 7.
B is a parameter determining the compressibility. It is related to the artificial speed of sound cf, which is chosen to allow approximately 1% of density fluctuation in SPH computations. For soils, an elastoplastic constitutive model with the Drucker -Prager yield surface is employed specified as where I1 is the first invariant of the stress tensor σ′s, J2 is the second invariant of the deviatoric stress tensor s′s = σ′s − (I1/3)I. αϕ and kc are Drucker-Prager's parameters, which can be related to the cohesion c and friction angle ϕ by matching the yield surface of Mohr-Coulomb model. A non-associated flow rule is employed, in which the plastic potential function is where ψ is the dilatancy angle.

SPH formulations
Within the framework of SPH, the approximation to the field function ( ) and its gradient ∇ ( ) at particle i write where and f j are the abbreviated form of ( ) and f (xj), in which xi and x j are the position vectors of particles i and j, respectively. The index j loops over all the particles in the support domain of the kernel function centered at particle i. The kernel function is a bellshaped weighting function defined as = ( − , ℎ), where h is the smoothing length determining the size of the support domain. In this work, the Wendland C2 kernel with a radius of 2h is employed. Vj is the volume of particle j. ∇ iWij is the gradient of the kernel function evaluated at the position of particle i.
Accordingly, the SPH discretized form of the soil governing equations, i.e., Eqs. (4)(5) are as follows (15) where pa is the interpolated water pressure at the soil particle and πab is the widely used artificial viscosity in SPH [11]. The second and third terms on the right side of Eq.
Herein, denotes the interpolated soil velocity at the water particle; = 0.1 and ij are the items from the δ-SPH [12].
In this work, two sets of SPH particles are generated to represent the water and soil phases, respectively. These two layers of SPH particles can overlap, and move independently according to their own governing https://doi.org/10.1051/e3sconf/202341502026 , 02026 (2023) E3S Web of Conferences 415 DFHM8 2 equations. A particle in one layer mostly interacts with the particles in the same layer but also has connections with particles in the other layer when computing the inter-phase interactions. The soil volume fraction is obtained as = s /̃s , where ρs is obtained by numerically integrating the apparent density using the continuity equation. According to the saturation assumption, the fluid volume fraction can be obtained as where i indicates a fluid particle and a denotes a soil particle. The index a loops over all the soil particles in the support domain of fluid particle i.

Initial and boundary conditions
At the beginning of the simulation, the particle mass is initialized by where α denotes s or f. In each calculation step, the particle's volume is updated with In this study, the boundary implementation proposed by Adami [13] is adopted.

Validation
In this section a dynamic and fully coupled problem is validated against analytical solution. As shown in Fig.  1, a soil body of height Hs = 0.5 m is fully submerged into a water body of depth Hw = 1.0 m. The soil has a particle density of ̃0 = 2750 kg/m 3 with a solid volume fraction of = 0.5. The mean grain size is assumed to be D50 = 0.05 m. Other material parameters of the soil are: Young's modulus E = 1.5 × 10 6 Pa, Poisson's ratio ν = 0.2, frictional angle ϕ = 25 • , dilatancy angle ψ = 0 • , and cohesion c = 0 Pa. Initially, all the stress and pressure for both the water and soil phases are set to zero. The gravitational load is then applied in a single step at the beginning of the simulation. Under the gravitational load, the pressure and stress in the water and soil phases will build up to balance the gravity, and finally reach a steady state, under which the exact pressure of water and the vertical stress of the soil can be analytically calculated. In this case, the water pressure and vertical soil effective stress are A simulation with particle size ∆r = 0.01 m is performed with 2.0 s of physical time. After the application of the gravitational load, the pressure and stress in the water and soil phases accumulate and start oscillating. This can be observed in Fig. 2, where the time history of the water pressure and soil vertical stress are plotted. Due to the numerical damping, viscous effect, and inter-phase drag force, the fluctuations in pressure and stress are quickly damped out, reaching a steady state after around 0.35 s. The analytical values of pressure and stress are obtained using Eq. (22). It is found that the numerical results obtained from the multi-layer SPH method are rather close to the analytical solutions.  We further show the water pressure and soil vertical stress along the central line in Fig. 3. The numerical results are in good agreement with the analytical solution. Along the central line, the normalized root mean square error is 0.5% for the water pressure, and 1.1% for the soil vertical stress. It is demonstrated again that the presented SPH method can correctly model the dynamic coupling of water-soil coupled problem.

Application
In this section, the proposed multi-layer SPH method is applied to investigate the submerged landslide and its induced waves, which was experimentally studied in [14]. Table 1 lists the material constants used in the simulations. Table 1. Data input in the submerged landslide simulation.

Parameters
Value As shown in Fig. 4, the experiment tank consists of mainly three parts. In the middle is an inclined slope of 45 • , above which is seated a soil body and a removable gate. After lifting the gate, the soil body can slide down towards the pure water domain on the left. At the right end of the tank there is a platform, which is used to model the shore line in reality. The dimensions of the computational domain is 4 m, 0.2 m and 1.7 m in the x, y and z directions. An initial particle spacing ∆r = 10 mm and a smoothing length h = 1.5∆r are employed in the simulation.  At the first stage, the soil mass slides downward slightly and still remains the original triangular shape. Besides, the water surface rises in front of the slide, whereas a drop is observed above the soil mass, resulting in a downward concave shape in the free water surface. Next, the slide tail lengthens and becomes sharp due to the accelerating movement, while the front is comparatively shorter and thicker because of the pronounced water resistance. Concurrently, the water surface, in a "W" shape, has strong fluctuations. These phenomena are also observed from the experiment. By virtue of the theoretical framework of the mixture theory, the proposed method is capable of capturing the kinematics of both phases within the sediment simultaneously, as shown by the smooth distributions of the fluid pressure and the soil velocity in Fig. 5. The prediction of impulsive wave is another key point in this test. As seen in Fig. 6, our simulation and [6] both give reasonable and satisfactory results compared with experiment data. On the other hand, drastic fluctuations on the water surfaces are observed in [15], which is probably induced by the overestimation of the mass movement as mentioned above. The good agreement between the multi-layer SPH and experimental results demonstrates the feasibility and capability of the proposed method in studying the submerged landslide.