E3S Web of Conferences 415, 02026 (2023)
DFHM8

https://doi.org/10.1051/e3sconf/202341502026

Mixture theory-based SPH model for submerged landslide
Chengwei Zhu1, Chong Peng2, and Wei Wu3*
1Research Center of Coastal and Urban Geotechnical Engineering, Zhejiang University, 310058 Hangzhou, China
2ESS Engineering Software Steyr GmbH, Berggasse 35, 4400 Steyr, Austria
3Institut für Geotechnik, Universität für Bodenkultur, Feistmantelstraße 4, 1180 Vienna, Austria

Abstract. A novel SPH model aimed at solving the coupled water-soil problems is proposed based on themixture 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 ismodeled 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 todemonstrate the capability of the proposed SPH model in solving the dynamic soil – water coupling
problems.

1 Introduction
Submerged landslide, widely observed in thecoastal region, is destructive to underwater facilities
such as submarine cables and pokes critical threat to theresidents 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-3].
Experimental studies have shed significant light on thisproblem but are often economic- and time-consuming
compared to numerical methods. A natural way ofinvestigating the water – soil coupling problem is todirectly consider the flow in the pore network and themotion 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 forengineering problems.
A feasible solution is to study the submergedlandslide dynamics within the unified continuumframework. Two schemes are popular in the previousstudies. One is the multi-phase model, which treats thepure water and sediments as two immiscible materials[6]. This method is straightforward and computationallyefficient, but limited in applications due to the neglectof soil water interactions. The second scheme, based on
the mixture theory, is featured with the spatiallyoverlapped continua representing the fluid and solid
phases [7]. This method is relatively more efficient thanthe 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

hydrodynamics is chosen as the numericalapproximation scheme due to its excellent performanceon dealing with large deformation problems [8].
The remaining part of this paper is outlined asfollows. 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.

2 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 solidphases. This practice is applied throughout this paper
unless otherwise stated. By definition, we can obtain thefollowing expressions between the intrinsic density, 𝜌̃,and the apparent counterpart, 𝜌,𝜌s = 𝜙s 𝜌̃s , 𝜌f = 𝜙f 𝜌̃f .
(2)
2.1 Conservation equations
The Lagrangian form of the mass conservation equation
for the solid phase readsds  s us ,dt(3)

where d(∙)/d𝑡 = ∂(∙)/d𝑡 + 𝐮s ∙ 𝛻(∙) is the material
derivative, the symbol ∇denotes the gradient operator
and · indicates the inner product. As the soil grains areincompressible, the intrinsic density of soil 𝜌̃s isconstant. Therefore, the evolution of the volume fractionof the soil phase is as follows* Corresponding author: wei.wu@boku.ac.at© The Authors, published by EDP Sciences. This is an open access article distributed under the terms of the Creative Commons Attribution License 4.0(https://creativecommons.org/licenses/by/4.0/).

E3S Web of Conferences 415, 02026 (2023)
DFHM8ds s  us .
dt

https://doi.org/10.1051/e3sconf/202341502026

2.3 SPH formulations(4)

Within the framework of SPH, the approximation to thefield function 𝑓(𝐱) and its gradient ∇𝑓(𝐱) at particle iwrite(12)f fWV ,The linear momentum conservation equation for thesoil phase in the Lagrangian form isdu s(5)s   σ s   s g  s  p  f d .
dtHerein, σ′s is the effective stress tensor, g denotes thegravity, p refers to water pressure; and fd represents theinteraction force between solid and fluid phases, owninga quadratic form with Ergun’s coefficients in this study
as [9]ij(13)jwhere 𝑓𝑖 and f j are the abbreviated form of 𝑓(𝐱 𝑖 ) and f(xj), in which xi and x j are the position vectors ofparticles 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 thesize of the support domain. In this work, the Wendland
C2 kernel with a radius of 2h is employed. Vj is thevolume 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 followsda(14) a   u a  ub  j  aWabVb ,dtbin 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 bechosen 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]
(7)

du a1

  a   b   ab I  aWabVb dt
a b(8)a1
  pi  pa  aWaiVi   i f aiWaiVi  ga (15)a baNote that the shear stress of fluid is neglected in Eq. (8)as the interstitial fluid is water and its viscosity ismarginal.

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 sideof Eq. (15) are the inter-phase interactions computed
based on the fluid particles in the support domain of soil
particle a. If there is no fluid particles in the support
domain, these two terms disappear in the momentum
equation. The drag force is computed as

2.2 Constitutive models
The fluids are assumed weakly compressible
dictated by an equation of state (EOS) such that  
p  B  f   1 , f0 ijfi   f j WijV j .

 1  f  1   
uf  us   d f f uf  us uf  us  (6)2f Dc
Dc,df  f   f uf  sus  ,
dt
fduf
f f
 f p  f f g  fd .
dt

jj2

fd   d(9) 1  i  ui  u a 
i Dc2
i 1  i 2f ai  150

where 𝜌̃f0 is the reference intrinsic density and is set as1000 kg/m3 for water, γ is a constant usually taken as 7.
B is a parameter determining the compressibility. It isrelated to the artificial speed of sound cf, which is chosen
to allow approximately 1% of density fluctuation in SPHcomputations.
For soils, an elastoplastic constitutive model with
the Drucker – Prager yield surface is employedspecified as(10)f  I1 , J 2   J 2   I1  kc ,where I1 is the first invariant of the stress tensor σ′s, J2 isthe 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, inwhich the plastic potential function is(11)g  I1 , J 2   J 2  3I1 sin ,where ψ is the dilatancy angle.

ui  u a  ui  u a  . (16)
Dc
Likewise, the fluid phase has the following SPH
discretization,di   i   j u j  i ui  iWijV j dt
i  j1.75

. (17)  1   j  u sj  i WijV j    hc f  Ψij  iWijV jjj
Herein, 𝐮𝑗𝑠 denotes the interpolated soil velocity at thewater particle; 𝛿 = 0.1 and 𝛹ij are the items from the δSPH [12].
dui1    pi  p j   ij  iWijV j
dt
i j1

 fiaWiaVa  gi .

i i(18)aIn 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, andmove independently according to their own governing2

E3S Web of Conferences 415, 02026 (2023)
DFHM8

https://doi.org/10.1051/e3sconf/202341502026

equations. A particle in one layer mostly interacts withthe particles in the same layer but also has connections
with particles in the other layer when computing theinter-phase interactions. The soil volume fraction isobtained as 𝜙𝑠 = 𝜌s /𝜌̃s , where ρs is obtained bynumerically integrating the apparent density using thecontinuity equation. According to the saturation
assumption, the fluid volume fraction can be obtained as(19)i  1  s  xi   1  aWiaVa .

layer SPH method are rather close to the analyticalsolutions.

awhere i indicates a fluid particle and a denotes a soil
particle. The index a loops over all the soil particles inthe support domain of fluid particle i.
2.4 Initial and boundary conditions

Fig. 1. Sketch of the problem of submerged soil.

At the beginning of the simulation, the particle
mass is initialized by(20)m    0  r  ,where α denotes s or f. In each calculation step, theparticle’s volume is updated withmm(21)V     .
3 

In this study, the boundary implementation
proposed by Adami[13] is adopted.

3 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 aparticle density of 𝜌̃𝑠0 = 2750 kg/m3 with a solid
volume fraction of 𝜙𝑠 = 0.5. The mean grain size isassumed to be D50 = 0.05 m. Other material parametersof the soil are: Young’s modulus E = 1.5 × 106 Pa,
Poisson’s ratio ν = 0.2, frictional angle ϕ = 25◦, dilatancy
angle ψ = 0◦ , and cohesion c = 0 Pa. Initially, all thestress and pressure for both the water and soil phases are
set to zero. The gravitational load is then applied in asingle step at the beginning of the simulation. Under thegravitational load, the pressure and stress in the waterand soil phases will build up to balance the gravity, andfinally reach a steady state, under which the exact
pressure of water and the vertical stress of the soil canbe analytically calculated. In this case, the waterpressure and vertical soil effective stress arep   H w  z   f 0 g ,     H  z        g . (22)
zz

w

ss0

Fig. 2. Change of water pressure and vertical effective stressat the point x = 1.0 and z = 0.2 m. Pressure is positive, andeffective stress is negative.

We further show the water pressure and soil verticalstress along the central line in Fig. 3. The numericalresults are in good agreement with the analyticalsolution. Along the central line, the normalized root
mean square error is 0.5% for the water pressure, and1.1% for the soil vertical stress. It is demonstrated againthat the presented SPH method can correctly model thedynamic coupling of water–soil coupled problem.

f0A simulation with particle size ∆r = 0.01 m is performed
with 2.0 s of physical time. After the application of thegravitational load, the pressure and stress in the waterand soil phases accumulate and start
oscillating. This can be observed in Fig. 2, where thetime history of the water pressure and soil vertical stressare plotted. Due to the numerical damping, viscous
effect, and inter-phase drag force, the fluctuations inpressure and stress are quickly damped out, reaching a
steady state after around 0.35 s. The analytical values ofpressure and stress are obtained using Eq. (22). It isfound that the numerical results obtained from the multi-Fig. 3. The water pressure soil vertical effective stress along
the central line at t = 2 s. Pressure is positive, and effective
stress is negative.

3

E3S Web of Conferences 415, 02026 (2023)
DFHM8

https://doi.org/10.1051/e3sconf/202341502026

4 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 thesimulations.
Table 1. Data input in the submerged landslide simulation.

Parameters𝜌̃𝑠 (kg/m3)
𝜌̃𝑓0 (kg/m3)
𝜙𝑠E (MPa)
υ

Value
2550
1000
0.62
1.5
0.3Parameters𝜑 (°)
c (kPa)𝜓 (°)
𝐷𝑐 (mm)Value
35004.5

Fig. 5. Reprinted snapshots of the submerged landslide from[14] (left) and numerical simulation (right). (a) t = 0.4 s, (b) t= 0.8 s.

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 resultscompared with experiment data. On the other hand,drastic fluctuations on the water surfaces are observedin [15], which is probably induced by the overestimationof the mass movement as mentioned above. The good
agreement between the multi-layer SPH andexperimental results demonstrates the feasibility andcapability of the proposed method in studying thesubmerged landslide.

As shown in Fig. 4, the experiment tank consists ofmainly three parts. In the middle is an inclined slope of◦
45 , above which is seated a soil body and a removablegate. After lifting the gate, the soil body can slide down
towards the pure water domain on the left. At the rightend of the tank there is a platform, which is used tomodel the shore line in reality. The dimensions of thecomputational 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 thesimulation.

Fig. 4. Initial configurations of the submerged landslides.

Fig. 5 presents the snapshots from the experimentand the simulation at t = 0.4 s and t = 0.8 s. At the firststage, 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 isobserved 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 andthicker because of the pronounced water resistance.
Concurrently, the water surface, in a “W” shape, hasstrong fluctuations. These phenomena are also observed
from the experiment. By virtue of the theoreticalframework of the mixture theory, the proposed method
is capable of capturing the kinematics of both phases
within the sediment simultaneously, as shown by thesmooth distributions of the fluid pressure and the soil
velocity in Fig. 5.

Fig. 6. Submerged landslide-induced wave surface fromdifferent methods at (a) t = 0.4 s and (b) t = 0.8 s.
The authors would like to acknowledge the funding from
National Natural Science Foundation of China under project
number 51709230. The work in this paper was also partly
funded by the EU Horizon 2020 RISE projects HERCULES(778360) and FRAMED (734485), and the Nazarbayev
University Research Fund (SOE2017001). The first author
wishes to thank the Otto Pregl Foundation for financial supportin Austria.

References
1.
2.
3.

4

M. A. Cabrera, G. Pinzon, W. A. Take, R. P.
Mulligan, J. Geophys. Res-Oceans, 125, 12 (2020)
J. Wang, L. Xiao, S. N. Ward, J. Du, Eng. Geol.
284, 267-277 (2021)
R. Duarte, J. Ribeiro, J. L. Boillat, A. Schleiss, J.
Hydraul. Eng. 137(11), 1385-1393 (2011)

E3S Web of Conferences 415, 02026 (2023)
DFHM8

4.

K. Han, Y. T Feng, D. R. J. Owen, Comput. Struct.
85(11-14), 1080–1088 (2007)

5.

C. Peng, L. Zhan, W. Wu, B. Zhang, Powder
Technol. 387, 509-526 (2021)
T. Capone, A. Panizzo, J. J. Monaghan, J. Hydraul.
Res. 48, sup1, 80–84 (2010)

6.
7.

C. Peng, G. Xu, W. Wu, H. S. Yu, C. Wang,
Comput. Geotech. 81, 239–248 (2017)8.

H. H. Bui, R. Fukagawa, K. Sako, S. Ohno, Int. J.
Numer. Anal. Met. 32(12), 1537–1570 (2008)
S. Ergun, AA. Orning, Ind. Eng. Chem. 41(6),1179–1184 (1949)9.

https://doi.org/10.1051/e3sconf/202341502026

10. C. Zhu, C. Peng, W. Wu, C. Wang, Comput.
Method Appl. M. 396, 115106 (2022)
11. J. J. Monaghan, J. Comput. Phys., 110(2), 399–406 (1994)
12. M. Antuono, A. Colagrossi, S. Marrone, Comput.
Phys. Comm. 183(12), 2570–2580 (2012)
13. S. Adami, X. Y. Hu, N. A. Adams, J. Comput. Phys.
231(21), 7057-7075 (2012)
14. SA. Rzadkiewicz, C. Mariotti, P. Heinrich, J.
Waterw. Port Coast. Ocean Eng. 123(4) 149–157(1997)
15. GF. Ma, JT. Kirby, FY. Shi, Numerical simulationof tsunami waves generated by deformablesubmarine landslides, Ocean Model. 69, 146–165(2013).

5