1 Introduction

A submarine landslide is a common marine geological disaster that occurs worldwide. It threatens the safe operation of marine infrastructure such as offshore wind turbines, drilling platforms, submarine pipelines, and cables due to its characteristics of large volume, high speed, and long runout distance (Rzadkiewicz et al., 1997; Perez-Gruszkiewicz, 2012; Dong et al., 2022). For example, the oil pipeline of Texaco Company in the United States was damaged by a submarine landslide, causing the leakage of 2100 gallons of crude oil into the sea (Zhao et al., 2021). In 2006, many submarine cables were damaged by a submarine landslide induced by the Pingtung earthquake in Taiwan, resulting in a major failure of international communication and substantial financial losses (Hsu et al., 2008). In 2009, two deepwater pipelines were destroyed by a submarine landslide triggered by the Suruga Bay earthquake in Japan (Matsumoto et al., 2010). In addition, submarine landslides may generate huge waves and lead to secondary disasters such as tsunamis. For example, the 1998 Papua New Guinea submarine landslide-induced tsunami resulted in over 2000 fatalities, destroyed numerous villages, and left more than 12000 people homeless (Tappin et al., 2001). In 2011, an earthquake-induced submarine landslide occurred in northeastern Japan and triggered a huge tsunami, killing tens of thousands of people and causing extensive damage to coastal buildings (Strasser et al., 2013). Moreover, in 2018, a submarine landslide-induced tsunami was triggered by a 7.5 magnitude earthquake in Sulawesi, Indonesia, which killed and injured more than 2000 people and caused extensive property damage (Pakoksung et al., 2019). Therefore, the safety of marine infrastructure can be ensured by obtaining a comprehensive understanding of the kinetic characteristics of submarine landslides and their impact on surrounding structures.

The extremely complicated submarine environment increases the difficulty of directly observing and monitoring the occurrence and movement processes of submarine landslides. Hence, physical model tests and numerical simulations are still the leading methods for studying the kinetic characteristics of submarine landslides. In terms of model tests, Watts (1998, 2000) established a series of flume model experiments to qualitatively and quantitatively investigate the wave characteristics generated by underwater landslides in solid blocks. Wiegel (1995) studied the influences of different initial submerged states, sliding angles, and water depths on the height and period of tsunami waves triggered by underwater landslides. Ataie-Ashtiani and Najafi-Jilani (2008) conducted 120 laboratory tests to analyze the effect of seabed slope angle, initial inundation, sliding geometry, shape, and deformation on shock wave properties, such as amplitude and period. Pilvar et al. (2019) studied the sliding forms of three different slopes, material types, and bed roughness and established a non-constant theoretical model to forecast the average velocity of underwater granular landslides. Zakeri et al. (2008) investigated the effect of clay-rich submarine mudflows on pipelines with two different placement forms (suspended and laid on the seafloor) by performing extensive flume model tests. They also proposed an equation to estimate the impact force of submarine landslides on the pipeline. Sahdi et al. (2014) explored the influence of soil strength, density, and the relative pipe-soil velocity on the impact force based on centrifuge model tests. They evaluated the impact forces of the submarine landslide by introducing a method that comprehensively considers geomechanics and hydrodynamics. Dai et al. (2017) discussed the relationship between the impact force and the sliding velocity, landslide volume, and cable diameter. Some achievements have been obtained through the above research. However, constraints such as scale effect, cost, and space limitation may be encountered by experimental studies.

Numerical simulation has recently been rapidly developed as an important alternative approach to investigating the motion and impact behavior of submarine landslides. For example, Yavari-Ramshe and Ataie-Ashtiani (2017) used the finite volume method (FVM) with nonlinear wet/dry transition treatment to investigate the submarine landslidetriggered wave. Zakeri et al. (2009) simulated the impact of submarine landslides on the suspended and laid-on-seafloor pipelines by applying the FVM method combined with the Herschel-Bulkley rheological model and validated the numerical results by comparing them with the experimental data. Considering three-phase media, Zhao et al. (2021) developed a coupled immersed boundary model to investigate the hydrodynamic load and vibration of pipelines exerted by submarine landslides. Fu and Jin (2015) employed the moving particle semi-implicit method to simulate the submarine landslide surge phenomenon and improved the model through a smoothing term and a pressure oscillation reduction method. Zhang et al. (2019) simulated the submarine landslide impact on pipelines at various impact angles using the computational fluid dynamics (CFD) model. Guo et al. (2022) introduced a CFD-based approach to quantitatively analyze the composition of drag forces on suspended pipelines exerted by submarine landslides. Dong et al. (2017) simulated the effect of the transient process of submarine landslides on the pipeline embedded in the seabed using the material point method combined with the Herschel-Bulkley rheological model. The smooth particle hydrodynamics (SPH) method has recently been widely used in multiphase or multiphysics problems due to its unique advantages in handling large deformations, free surfaces, and moving interfaces. For example, Ritchie and Thomas (2001) adapted the SPH technique to enable a multiphase fluid, which facilitates the free intermixing of SPH particles with widely different densities. Zhu et al. (2018) presented an SPH model for simulating multiphase fluid flows with large density ratios, which demonstrated the effective capture of inherent multiphase flow physics. Shimizu et al. (2020) introduced an enhanced multiphase ISPH-based method to simulate oil spill problems. Small-scale turbulence and the oil-water mixing effects were reproduced and captured in the simulations. Khayyer et al. (2023) improved the resolution of the continuity equation in the explicit weakly compressible SPH (WCSPH) model for the simulation of incompressible free-surface fluid flows. Enhanced satisfaction with incompressibility conditions and invariant density conditions were achieved. In addition, Luo et al. (2021) provided a state-of-the-art review of the applications of the SPH method in multiphase and multiphysics problems and emphasized the future perspectives of extending this method to a wide range of ocean and coastal engineering applications. The SPH method has been extensively applied in the simulation of submarine landslides in the past decade. For instance, Capone et al. (2010) applied the SPH method and established the underwater landslide as a non-Newtonian Bingham fluid to simulate the process of submarine landslide motion as well as wave generation. Qi et al. (2022) introduced a new multiphase SPH method for simulating various types of submarine landslides in surge experiments. Yeylaghi et al. (2017) simulated subaerial and submarine landslides through the development of an incompressible SPH (ISPH) model. Zhang et al. (2022) simulated the production and propagation of tsunami waves from submarine and subaerial landslides by developing a WCSPH model based on two-phase mixing theory and a large eddy simulation (LES) turbulence model. SPH has recently been coupled with other numerical algorithms for submarine landslide simulation. For instance, Qiu et al. (2017) simulated tsunamis generated by submerged, rigid, and deformed landslides using a combination of the discrete element method (DEM) and SPH. They employed a combination of the Papanastasiou and Herschel-Bulkley rheological models to capture the viscoplastic behavior of non-Newtonian flows in submarine-deformed landslides. Peng et al. (2020) developed an SPH model coupled with (discontinuous deformational analysis) DDA to address the fluid-solid interaction involved in submarine landslide propagation. The above studies reveal the unique advantages of the SPH method in simulating submarine landslide movement due to its meshfree and self-adaptability characteristics. Mean-while, SPH simulations of the impact behavior of fluids on structures have also produced some research results (Marrone et al., 2011; Pan et al., 2016; Li et al., 2020). However, investigations into the impact behavior of submarine landslides are still limited.

An SPH-based multiphase flow model is established in this paper to simulate the motion and impact behavior of submarine landslides. This model considers the landslide body and the ambient water as non-Newtonian and Newtonian fluids, respectively. The influence of particle resolution and rheological parameters on the motion behavior of submarine landslides is investigated, and the accuracy of the model in predicting the submarine landslide impact force on cable is validated.

2 Multiphase Flow SPH Model 2.1 SPH Formulation

The SPH method was originally applied to the simulation of astrophysical problems (Gingold and Monaghan, 1977; Lucy, 1997). Owing to its continuous improvement (Tafuni et al., 2018; Khayyer et al., 2023), the SPH method is increasingly used in the simulation of problems in the multiphase flow field, such as submarine landslides, wave impacts, and oil spills (Shimizu et al., 2020; Luo et al., 2021; Mahallem et al., 2022; Su et al., 2022). In the SPH method, the continuous fluid (or solid) is replaced by finitely interacting particles. Each particle contains field variables such as mass, pressure, and velocity, and the approximate values of the field variables are calculated using the weighted average of the interactions with other particles in the support domain (Gomez-Gesteira et al., 2010; Snelling et al., 2020). The SPH discretization of the governing equation is performed in two separate steps: 1) kernel approximation, which applies an integral representation to an arbitrary function, and 2) particle approximation, which transforms the integral form into a discretized one. Therefore, the approximate field variable functions f (x_i) and their derivatives ∇f (x_i) in the SPH method are shown as follows (Liu and Liu, 2010):

$ f\left({{{\boldsymbol{x}}_i}} \right) \approx \sum\limits_{j = 1}^N {\frac{{{m_j}}}{{{\rho _j}}}} f\left({{{\boldsymbol{x}}_j}} \right){W_{ij}}, $

(1)

$ \nabla f\left({{{\boldsymbol{x}}_i}} \right) \approx \sum\limits_{j = 1}^N {\frac{{{m_j}}}{{{\rho _j}}}} f\left({{{\boldsymbol{x}}_i}} \right){\nabla _i}{W_{ij}}, $

(2)

$ {\nabla _i}{W_{ij}} = \frac{{{{\boldsymbol{x}}_i} - {{\boldsymbol{x}}_j}}}{{\left| {{{\boldsymbol{r}}_{ij}}} \right|}}\frac{{\partial W\left({\left| {{{\boldsymbol{x}}_i} - {{\boldsymbol{x}}_j}} \right|, h} \right)}}{{\partial {{\boldsymbol{r}}_{ij}}}}, $

(3)

where x is the coordinate of the SPH particle, m is the mass of the particle, ρ is the density, h is the smooth length, r_ij = x_i – x_j is the distance between particles, N is the total number of neighboring particles in the support domain, and W is the smooth kernel function. The quintic spline kernel function proposed by Morris (1996) is applied in this study, as shown in Eq. (4):

$ W(q, h)=\alpha_d \times\left\{\begin{array}{l} (3-q)^5-6(2-q)^5+15(1-q)^5, 0 \leq q <1 \\ (3-q)^5-6(2-q)^5, 1 \leq q <2 \\ (3-q)^5, 2 \leq q <3 \\ 0, q \geq 3 \end{array}\right., $

(4)

where q = |r_ij|/h, α_d = 7/(478πh²).

2.2 Governing Equations

The underwater landslide is treated as a weakly compressible viscous flow, which satisfies the law of mass and momentum conservation, and the corresponding governing equations are shown below:

$ \frac{\mathrm{d} \rho}{\mathrm{~d} t}=-\rho \nabla \cdot \boldsymbol{v}, $

(5)

$ \frac{\mathrm{d} v}{\mathrm{~d} t}=-\frac{1}{\rho} \nabla p+\frac{1}{\rho} \boldsymbol{F}_{\mathrm{vis}}+g, $

(6)

where t is the time, v is the velocity vector, g is the gravitational acceleration, F_vis is the viscous force of the fluid, and p is the pressure term of the fluid, which is calculated using the following equation of state (Monaghan, 2012):

$ p = {\rho _0}c_s^2\left[ {\left({\frac{\rho }{{{\rho _0}}} - 1} \right)} \right] + \chi, $

(7)

where ρ₀ is the reference density, c_s is the sound velocity, and χ is the background pressure. c_s should be greater than 10 times the maximum velocity to ensure that the density fluctuations are within 1% and substantially less than the true velocity, thereby avoiding extremely small time steps. c_s is typically limited by the following equation (Meringolo et al., 2018):

$ {c_s} \geqslant 10\max \left({\sqrt {\frac{{{p_{\max}}}}{{{\rho _0}}}}, {v_{\max}}} \right), $

(8)

where p_max and v_max are the estimated maximum pressure and velocity in the flow field, respectively. As shown below, the two phases are required to satisfy the ρ₀c_s² conservation (Colagrossi and Landrini, 2003) and enhance the stability and accuracy of the simulations:

$ {\rho _{{0_1}}}c_{{s_1}}^2 = {\rho _{{0_2}}}c_{{s_2}}^2, $

(9)

where subscripts 1 and 2 indicate the landslide particle and water phases, respectively.

During calculation in the pressure field, numerical oscillations in the traditional WCSPH method may be attributed to the density approximation. Therefore, in this study, an artificial diffusion term D_i is added to the continuity equation to effectively eliminate the high-frequency numerical noise.

$ {D_i} = 2\delta {h_i}{c_s}\sum\limits_{j = 1}^N {{\boldsymbol{\psi}_{ij}}} {V_j}{\nabla _i}{W_{ij}}, $

(10)

where V denotes the volume of the particle, and δ is the dimensionless coefficient, which is usually taken as δ = 0.1. The term ψ_ij uses the modified form proposed by Molteni and Colagrossi (2009) and Krimi et al. (2018):

$ {\boldsymbol{\psi}_{ij}} = \left({{\beta _{{\rho _i}}} - {\beta _{{\rho _j}}}} \right){\rho _{0i}}\frac{{{\boldsymbol{r}_{ij}}}}{{{{\left| {{\boldsymbol{r}_{ij}}} \right|}^2}}}, $

(11)

where β_ρ is the density ratio of a particle in terms of its initial reference density, i.e., β_ρ = ρ/ρ₀.

2.3 Viscous Item

The acceleration due to viscous forces is shown as follows (Adami et al., 2010):

$ \frac{1}{\rho} \boldsymbol{F}_{\mathrm{vis}}=\frac{1}{m_i} \sum\limits_{j=1}^N \frac{2 \mu_i \mu_j}{\mu_i+\mu_j}\left(V_i^2+V_j^2\right) \times \frac{\boldsymbol{r}_{i j} \cdot \nabla_i W_{i j}}{\left|\boldsymbol{r}_{i j}\right|^2} \boldsymbol{v}_{i j}, $

(12)

where μ is the dynamic viscosity of the particle. The viscosity for water particles is calculated based on the Newtonian fluid model, while μ for landslide particles is derived using the non-Newtonian Herschel-Bulkley rheological model. The expression is as follows:

$ \mu = \eta {\left| {\dot \gamma } \right|^{K - 1}} + \frac{{{\tau _y}}}{{\left| {\dot \gamma } \right|}}, $

(13)

where τ_y is the yield stress, η is the viscosity coefficient, K is the power-law exponent, and their magnitude is determined by the specific material. $\dot \gamma $ is the shear strain rate given by:

$ \dot{\gamma}=\nabla \boldsymbol{v}+(\nabla \boldsymbol{v})^{\mathrm{T}},, $

(14)

where ∇ν and (∇ν)^T represent the velocity gradient and transpose tensors, respectively.

2.4 Boundary Conditions

In the simulation of the submarine landslides, the support domain of fluid particles near the solid wall is truncated by the wall boundary, leading to a reduction in computational accuracy. Adami et al. (2012) proposed a solidwall virtual particle boundary condition that discretizes the solid-wall boundary into a series of virtual particles with physical quantities, as shown in Fig.1. This method not only prevents the penetration of the fluid particles near the boundary but also maintains the continuity of the pressure, density, and other physical fields of the fluid particles near the boundary to improve the computational accuracy of the simulation. Therefore, this paper applied the aforementioned boundary condition in the simulation of the multiphase flow problem of submarine landslides. The physical quantity of ghost particles is obtained by interpolating the neighboring fluid particles, and the pressure can be expressed by Eq. (14):

$ {p_b} = \frac{{\sum\limits_f {{p_f}{W_{bf}} + \left({g - {{a}_w}} \right)\sum\limits_f {{p_f}\left({{{x}_b} - {{x}_f}} \right){W_{bf}}} } }}{{\sum\limits_f {{W_{bf}}} }}, $

(15)

where the subscript b denotes solid-wall ghost particles, the subscript f denotes fluid particles, and a_w is the acceleration at the wall boundary.

A free slip boundary condition is implemented by the model by ignoring the viscous force interaction between fluid particles and neighboring ghost particles. Ghost particles for the no-slip boundary condition are imposed via interpolation with the interacting fluid particles to obtain the velocity. The expression for the ghost particle velocity is given as follows:

$ \boldsymbol{v}_b=2 \boldsymbol{v}_w-\frac{\sum\limits_f \boldsymbol{v}_f W_{b f}}{\sum\limits_f W_{b f}}, $

(16)

where v_w is the velocity of the solid wall.

Solid boundary particles are used to build the structure to simulate the impact force of submarine landslides on the structure, and the virtual acceleration of the structure can be obtained by solving Eq. (5). The total impact force can then be calculated by multiplying the acceleration with the particle mass (Crespo et al., 2015; Klapp et al., 2020). The specific expression is shown below:

$\frac{\mathrm{d} \boldsymbol{v}_s}{\mathrm{~d} t}=-\frac{1}{m_s} \sum\limits_f\left(p_s V_s^2+p_f V_f^2\right) \nabla_s W_{s f}+\frac{1}{m_s} \sum\limits_f \frac{2 \mu_s \mu_f}{\mu_s+\mu_f}\left(V_s^2+V_f^2\right) \times \frac{\boldsymbol{r}_{s f} \cdot \nabla_s W_{s f}}{\left|\boldsymbol{r}_{s f}\right|^2} \boldsymbol{v}_{s f}, $

(17)

$ \boldsymbol{F}=m \sum\limits_s \frac{\mathrm{~d} \boldsymbol{v}_s}{\mathrm{~d} t}, $

(18)

where the subscript s denotes the structure particle, the subscript f denotes the fluid particle, and F denotes the impact force on the structure.

2.5 Time-Step Scheme

Restricting the size of the time step by imposing an artificial term that is related to the sound velocity, viscous diffusion term, and fluid maximum acceleration conditions is necessary to ensure the stability of the simulation (Monaghan, 1992, 2005; Morris et al., 1997). The time-step constraints are shown below:

$ \Delta t \leq \min \left(1.25 \frac{h}{c_s}, 0.125 \frac{h^2}{\max \left(\frac{\mu_i}{\rho_i}\right)}, 0.25 \mathrm{~min} \sqrt{\frac{h}{\left|\boldsymbol{a}_i\right|}}\right), $

(19)

where a is the acceleration of the particle.

3 Model Validation and Application 3.1 Two-Phase Hydrostatic Tank

The stability and accuracy of multiphase flow models are validated by widely simulating the two-phase hydrostatic tank as a benchmark problem (Krimi et al., 2018; Rezavand et al., 2020). As shown in Fig.2, a two-dimensional rectangular tank with a width of 0.4 m is filled with two stratified fluids. The upper and lower layers are light and heavy fluids with densities of 1000 and 2000 kg m⁻³, respectively. Both fluids are subjected to a vertical gravitational acceleration of 9.81 m s⁻². The upper fluid behaves as a Newtonian fluid with a constant viscosity of η₁ = 0.02 Pa s, while the lower fluid behaves as a Herschel-Bulkley rheological fluid with a viscosity coefficient η₂ = 0.001 Pa s, a power-law exponent N = 2, and a yield strength τ_y = 250 Pa.

The two-phase hydrostatic tank model is replaced by SPH particles with three different sizes (dp = 0.001, 0.002, and 0.004 m) to investigate the influence of SPH particle size on the simulation accuracy. Fig.3 shows the simulated pressure field of the model with three particle sizes. The pressure distributions are smooth, and pressure discontinuity and oscillation are absent. The model with fine particle resolutions can obtain highly stable results.

The numerical results are compared with the analytical solution to validate the simulation accuracy of the multiphase flow SPH model. As shown in Fig.4, the vertical pressure distribution obtained using SPH simulation is close to the analytical solution.

The L₂ error in pressure is calculated using the following equation to quantify the error between the numerical and analytical results:

$ {L_2}{\kern 1pt} {\rm error} = \frac{{\left\| {p - \hat p} \right\|}}{{\left\| {\hat p} \right\|}}, $

(20)

where p is the calculated pressure in the SPH simulation, and${\text{\hat p}}$is the analytical result. The error paradigm ||p|| can be calculated as follows:

$ \|p\|=\left(\int\limits_{\Omega_0} p \cdot p \mathrm{~d} \Omega_0\right)^{1 / 2}. $

(21)

Fig.5 shows the convergence curves of the multiphase flow SPH model. These curves indicate the small L₂ error of the model, which decreases as the number of particles increases.

The numerical experiment in this section shows consistency between the simulated pressures of the multiphase flow SPH model and those of the analytical solution, and the SPH model with fine particle resolution can obtain highly accurate results.

3.2 Waves Generated by the Submarine Landslide

A physical model test designed by Rzadkiewicz et al. (1997) is simulated as a case study in this paper to validate the computational accuracy of the multiphase flow SPH model. The test was conducted in a 4 m × 2 m water tank with a maximum water depth of 1.6 m. The slide body with a cross-sectional area of 0.65 m × 0.65 m was placed on the slope with an inclined angle of 45˚. Fig.6 shows the schematic of the model test. The densities of the landslide body and ambient water were 1950 and 1000 kg m⁻³, respectively. The Herschel-Bulkley rheological model is introduced to describe the landslide behavior, and the Newtonian fluid model is used for the simulation of the surrounding water.

The rheological parameters used in this paper are selected in accordance with the simulations performed by Ataie-Ashtiani and Shobeyri (2008), τ_y = 750 Pa and η = 0.15 Pa s. Fig.7 shows the pressure field of the water at different moments. The figure shows the absence of pressure oscillation and a relatively continuous and smooth pressure field. Fig.8 shows the velocity vector of the fluid phase, which describes the interaction between a landslide and water. When the landslide body slides downward under the effect of gravity, its kinetic energy is transferred to the ambient water, thus inducing surge waves at the water surface and creating a vortex above the front of the slide.

The landslide profiles at the moments of 0.4 and 0.8 s are compared with the experimental results recorded in Rzadkiewicz et al. (1997) and the ISPH results presented in Ataie-Ashtiani and Shobeyri (2008) to validate the accuracy of the multiphase flow SPH model. Fig.9 shows some discrepancies between measured and simulated landslide profiles. When t = 0.8 s, the sliding velocity of the landslide body is relatively high, resulting in a strong interaction between soil particles and the ambient water. Therefore, water erodes the soil particles in the landslide front, leading to a thick landslide front. However, the SPH particles in the numerical model are substantially larger than the real soil particles and are thus difficult to erode by the water. Therefore, the simulated landslide profile near the front cannot accurately match the experimental results.

The simulated surface wave profiles are compared with the experimental observations in Rzadkiewicz et al. (1997) and the numerical results in Capone et al. (2010) and Bu et al. (2022) (as shown in Fig.10) to further validate the simulation accuracy of the SPH model. The goodness of fit R² is applied for the quantitative analysis of the results, and the calculated expression is shown in Eq. (21). In this case, the y-direction coordinates (y₁, y₂, y₃, ···, y₃₅) of 35 points with equal horizontal distances on the wave profile are selected for the analysis.

$ {R^2} = 1 - \frac{{\sum\limits_{i = 1}^n {{{\left({{y_i} - {{\hat y}_i}} \right)}^2}} }}{{\sum\limits_{i = 1}^n {{{\left({{y_i} - {{\bar y}_i}} \right)}^2}} }} . $

(22)

Table 1 shows the calculated goodness of fit. The calculated fitness is greater than 0.8, which indicates that the numerical results agree well with those from the experimental and existing numerical models. Therefore, the landslide-induced wave phenomenon can be effectively reproduced by the SPH model established in this paper.

Sensitivity analysis of η and τ_y of the Herschel-Bulkley rheological model is conducted based on the above results to analyze the influence of rheological parameters of landslides on the characteristics of landslide motion. The maximum frontal velocity v_f-max of the granular landslide and the landslide body profile at 0.8 s are simulated, and the results are shown in Fig.11. The v_f-max of the submarine landslide gradually increases with the decrease in η and τ_y, and τ_y has a notable influence on the v_f-max, as shown in Figs. 11a and 11b. Fig.11c shows that the landslide configuration at 0.8 s slightly changes with the variation of η when the η value is small. The landslide movement distance considerably shortens only when the η value is large. Fig.11d shows the acceleration of landslide flows and the extension of the sliding distance at 0.8 s with the reduction in τ_y. Therefore, the rheological parameter τ_y has a more notable influence on the mobility of the landslide body compared with the viscosity coefficient η.

3.3 Impact of an Underwater Landslide on a Cable

A rotating flume experiment conducted by Wang et al. (2018) to replicate the impact of an underwater landslide on a cable is simulated to validate the SPH impact model established in this work. Fig.12 shows the schematic of the model test. The main part of the apparatus is an annular trough with an axle in the center. The outer radius, inner radius, and width of the trough are 0.9, 0.6, and 0.4 m, respectively. The axle of the apparatus is connected to a motor to facilitate rotation in a vertical plane at a controlled velocity. A certain amount of water and soil material is placed in the apparatus to simulate the submarine landslide. The fine silica sand used in the experiment as the landslide body is assumed to be a Herschel-Bulkley fluid in the simulation. A cable is installed in the bottom of the annular trough parallel to the axel. The cable comprises stainless steel with a length of 370 mm. The height of the cable from the frame bottom is 0.2 m.

Two load cells are set at two sides of the cable to measure the impact force from the water and landslide mass. The measurement range of the load cell is 20 N, and the measuring precision is 0.1 N. A water pressure transducer with a measurement range of 10 kPa and a precision of 0.1 kPa is fixed inside the frame bed to measure the water pressure during the experiment.

The water and the landslide body remain at the bottom during the rotation due to the effect of gravity. The cable moves at a certain speed and impacts the water and the landslide, simulating the process of a submarine landslide impacting the cable. Fig.13 shows the SPH-simulated velocity fields at different times with the flume rotation at 0.38 m s⁻¹. The dry weight of the landslide body is 30 kg, and the cable diameter is 20 mm. The cable has a certain effect on the velocity field of the surrounding fluid, and a vortex appears near the cable as it moves. The initial state of the soilwater mixture is disturbed when the cable starts to move, and the surrounding water agitates the sand particles. Turbid fluid appears along the cable trace. Fig.14 shows that the initial smooth pressure field of the soil-water mixture is also disturbed due to the impact of the cable.

The impact force and water pressure are compared with the experimental data to validate the simulation accuracy of the SPH model. As shown in Fig.15, the evolution of the simulated impact force and water pressure curves agrees well with the experimental data, and the peak values of the calculated impact force and water pressure are close to the measured data. The goodness of fit R² is calculated and listed in Table 2. The R² of impact force and water pressure are both greater than 0.8, indicating that the SPH simulation results fit well with the measured data in the experiments. Convergence analysis of impact force and pressure time history curves is conducted, and the L₂ error formula (Eqs. (20) and (21)) is used to calculate the error of different particle spacing (dp = 0.001, 0.002, and 0.004 m), as shown in Fig.16. The simulated results of impact force and pressure are closer to the test data as the number of particles increases.

Fig.17 shows a comparison between the measured peak impact force F and the results obtained from the SPH simulation under different silica sand mass m. The diameter of the pipe is 20 mm, and the velocity of the landslide body is 0.38 m s⁻¹. The simulated peak impact forces are closer to the experimental data. Therefore, the submarine landslide-cable interaction can be reproduced by the SPH model with certain accuracy.

4 Conclusions

Submarine landslides seriously threaten coastal marine infrastructure as well as human life due to their immense impact force and the generation of huge waves. An SPH-based multiphase flow model is presented in this paper to simulate submarine landslide motion and impact behavior. The Herschel-Bulkley and Newton fluid models are introduced to simulate the submarine landslide body and the ambient water, respectively. The model is improved by applying artificial diffusion terms and multiphase flow corrections to the equation of state to ensure numerical accuracy and stability.

The accuracy and stability of the multiphase flow SPH model are validated by simulating and analyzing three benchmark cases. The computational accuracy of the SPH model increases with particle resolution. Sensitivity analysis shows that the rheological parameter τ_y has a more notable influence on the mobility of the landslide body compared with the viscosity coefficient η. The numerical results are compared with the analytical solution, experimental data, and simulated results obtained from other numerical models. Therefore, the SPH model is suitable for predicting the impact force of submarine landslides exerted on cables.

The multiphase model presented in this paper demonstrates good performance in simulating the motion and impact behavior of submarine landslides in two dimensions. However, the real submarine landslide occurs in a three-dimensional space with complex submarine topography. Therefore, a 3D model with high-efficiency parallel computing technology is necessary to reproduce the real motion and impact behavior of submarine landslides.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (Nos. 42472332, 42102318 and 42006143), the Open Research Fund Program of Zhoushan Field Scientific Observation and Research Station for Marine Geo-Hazards, China Geological Survey (No. ZSORS-22-07), the Program for Professor of Special Appointment (Eastern Scholar) at Shanghai Institutions of Higher Learning (No. TP2019037), and the Open Research Fund Program of Marine Ecological Restoration and Smart Ocean Engineering Research Center of Hebei Province (No. HB MESO2312).

Author Contributions

All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Baisen Lan. The draft of the manuscript was written by Zili Dai, and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.

Data Availability

All data generated and analyzed during this study are included in this published article and its additional files.

Declarations

Ethics Approval and Consent to Participate

This article does not contain any studies with human participants or animals performed by any of the authors.

Consent for Publication

Informed consent for publication was obtained from all participants.

Conflict of Interests

The authors declare that they have no conflict of interests.