Exploring the submarine landslide is challenging due to the invisibility nature and the complex soil-water interaction and large deformation throughout its runout process. The purpose of this paper is to investigate the ability of the coupled material point method (MPM) in modeling the soil flows under water. A sand-column collapse experiment is performed fully under water initially, with the results used to benchmark the MPM analysis. Thereafter, the whole failure process of a real submarine landslide in the South Mediterranean sea is simulated using MPM. The results show that MPM can be a reliable tool in capturing the postfailure behaviors of the submarine landslide. The failure mode of the landslide is flow-type, with an initial translational slide moving to a diffusive one eventually.

1. Introduction

Submarine landslides are the most common and frequent geological hazards in the ocean [1, 2]. Both natural and human activities, such as earthquake, gas hydrates dissociation [3, 4], and hydrocarbon production [5], are able to trigger submarine landslides. As opposed to terrestrial landslides, submarine landslides normally affect a wider range of areas [6, 7], not only posing a great risk to submarine pipelines and other seafloor structures [8, 9] but also causing severe casualties due to tsunamis. Compared with forecasting the occurrence probability of a marine landslide, estimating the postfailure behaviors and characteristics of a marine landslide is of equal importance [10, 11].

Evolutions of sliding distance and velocity during a landslide runout process are necessary for disaster prevention and mitigation [12]. The mechanism of submarine landslides has been studied by a couple of scholars. Koning [13] proposed that the marine sediments may turn into flows given that their magnitude of porosity exceeds a critical value. Middleton and Hampton [14] studied the dynamics of flows and concluded that gravity-driven process plays an important role in the transport of sediments. Whitman [15] classified submarine landslides into disintegrative and nondisintegrative failure, making a distinction between flows and less mobile types of landslides. In this vein, the notion of a steady state of deformation was developed by Poulos [16] to quantify the development of disintegrative failure. Besides, many empirical formulations and models have thereby been proposed [17]. However, they are usually of high approximation as well as poor applicability and generality. In contrast, numerical methods permit a better assessment of submarine landslides [18, 19]. Classical mesh-based methods, such as the finite element method (FEM), are normally concerned about when the slide starts, but it is difficult to capture the whole runout process due to the potential severe mesh distortions in large deformations. Discrete element method (DEM) is able to simulate soil flows [20], but its utilization may be restricted due to a relatively prohibitive computational cost.

Material point method (MPM), originating from the fluid mechanics, was first applied to the solid mechanics in 1994 [21], and, since then, increasing attention has been given to it for a wide range of slope failure analyses [22]. By utilizing a hydromechanical coupled MPM formulation, Alonso and Zabala [23] analyzed the progressive failures of Aznalcóllar dam. Bandara and Soga [24] presented a fully coupled MPM formulation based on the mixture theory, where two sets of Lagrangian material points are utilized, representing the soil and water, respectively, and the progressive failure of a river levee was simulated. Yerro et al. [25] further extended the coupled MPM formulation to model the unsaturated soils by integrating three phases (solid, liquid, and gas) within a set of material points, with each occupying a certain fraction. Wang et al. [26] adopted a one-set two-phase coupled MPM formulation, and a rainfall-induced slope failure was analyzed by incorporating a strain-softening model to describe the strength reduction of the soil due to the rainfall infiltration.

Due to the fact that the submarine landslides occur fully under water, their invisibility makes it hard to determine the runout distances as opposed to the landslides. The main objective of this paper is to investigate the ability of MPM to simulate submarine landslides in order to quantify the calamity caused by it. A series of parametric analyses are conducted as well in order to gain a better understanding of the influencing factors to the slide consequences.

2. Overview of MPM

2.1. Basic Principle

MPM is one of the mesh-free methods, which is described by a Lagrange particle system and Euler grid. Although MPM utilizes a background mesh, it is still regarded as a mesh-free method because no polygonal tessellation is involved in the initial discretization of the material. A schematic diagram of MPM is shown in Figure 1, which demonstrates a typical two-dimensional discretization of a continuum. In MPM, all the material information of the medium is carried by the discrete particle system, such as the initial value of position, mass, volume, stress, and velocity. Additional parameters like pore water pressure can be specified at the material points depending on the material being simulated [27]. Background grid covers the whole space particles can reach and is only used for calculating spatial derivatives and solving momentum equations.

At every time step, the material information is first mapped to the background mesh using the shape functions, which is derived from the position of each material point with regard to the mesh. After gradient terms are calculated and governing equations are solved on the background mesh, the updated information will be mapped back to the particles. At last the distorted mesh is abandoned and a new one is initialized.

Because of this feature, MPM avoids difficulty caused by mesh distortion, interface tracking, and nonlinear convection, making it outstand other numerical methods in simulating geotechnical problems involving large displacement. Meanwhile MPM can adopt history-dependent constitutive model. The governing equations and spatial discretization of MPM can be referred to the principle parts of literature [21].

Based on Biot’s theory [28] of porous media, Bandara and Soga [24] developed the principle theory of coupled MPM. The following assumptions are made in the coupled MPM formulation described in the next section:(a)Saturated soil is considered to be made up of solid framework and water phases, both of which are described in a Lagrangian formulation and assumed as a continuous medium.(b)The solid framework is considered as incompressible.(c)An effective stress model is adopted for the solid framework.(d)There is no mass exchange between the solid framework and water phases.(e)Variation of water density is negligible. Variation of porosity and volume of solid framework with time is considered.

2.2. Governing Equations for Coupled MPM

Porosity of soil can be updated from equilibrium equation of solid framework:where and are the mass and volume of the solid framework and n is the porosity. Dividing formula 1 by , based on the hypothesis that solid framework is incompressible, we obtain the equation in which porosity changes with time:

Porosity of soil can be obtained with formula (2). Based on mass balance equation of liquid part, the state equation can be expressed aswhere and are the mass and volume of the liquid part. Divide formula (3) with :where is the bulk modulus of water. Introducing formula (2) to (4), the state equation of water phase is obtained:

Momentum balance equation of soil framework can be expressed as

Momentum balance equation of liquid part can be expressed as

and are acceleration of solid framework and liquid part. is the fluid-coupled momentum exchange term of solid framework, is equal and opposite to , and are partial stress tensor of solid framework and liquid part, and is the body force. and are defined aswhere k is the soil permeability. The second term in formula (8) is called “buoyancy term” in the mixture theory. and can be written aswhere is Terzaghi’s effective stress, is the water pressure, and is the Kronecker delta vector ().

Introducing formulas (8) and (9) into formula (6), the ultimate forms of momentum balance equation of solid framework can be expressed as

Introducing formulas (8) and (10) into formula (7), the ultimate forms of momentum balance equation of liquid part can be expressed as

Combining formulas (11) and (12), the momentum balance equation of mixture part can be expressed as

is total stress of mixture part and it is defined as

is effective water pressure. In saturated soils, .

2.3. Computational Cycle

The calculation procedures in each computational cycle are described by the following steps:(a)The material information is projected to the background mesh using shape functions. Internal forces of the water phase and mixture parts are evaluated in the nodes.(b)The momentum balance equations of solid framework and water phase are solved and nodal accelerations are calculated.(c)With accelerations obtained in step (b), the momentum balance equation of the mixture part is solved.(d)Velocities of particles and nodes are updated with the forward Euler scheme, while the displacements are updated with the backward Euler scheme.(e)The mass balance equations are solved and the material information is updated.(f)The updated information is mapped back to material points. The distorted mesh is abandoned and a new background mesh is initialized for next step.

3. Sand-Column Collapse

The collapse process of a sand column includes some key stages similar to the runout process of a landslide [29, 30] so that it is widely used as a simplified landslide model in geotechnical tests [31]. A submerged sand-column collapse example is presented in this section for an initial demonstration of the ability of the MPM in modeling the large deformation of submarine soils.

3.1. Experiment Setup

A small-scale in-door experiment was conducted to serve as a benchmark, as shown in Figure 2. The model box has a size of 50 cm × 25 cm × 50 cm. A baffle is placed at a certain distance from the left side of the box to keep the sand column steady initially. Supported by the baffle, the sand forms a columnar accumulation with an initial length of L = 7.5 cm and a height of H = 15 cm, respectively. Water is thereafter poured into the box and the water depth is kept at 10 cm to make sure the whole collapse process occurs under water. According to Brok and Spiers [32], the samples need to be immersed in water for 2 hours to ensure that they are fully saturated.

Once the experiment starts, the baffle is quickly pulled up, rendering the sand-column collapses to the right direction, and forms a granular flow in fluids. After the system is stabilized, the sliding distance and the morphology of the sample in its final configuration are measured.

The material used here is silty sand, with the material properties summarized as shown in Table 1.

3.2. MPM Simulations

In order to verify the reliability of the MPM model, the whole failure process is simulated with the dynamic explicit MPM software Anura3D [33], which is capable of dealing with fluid-solid coupling problems [34, 35]. Different from previous researches, where the sand-column collapse is normally carried out in dry conditions, the column is designed to collapse fully under water in the experiment. A one-layer two-phase (solid and fluid) MPM model is adopted for the simulation.

The deformation of the sand column at different times is shown in Figure 3. As can be seen, the upper right part of the column starts to destabilize initially due to the removal of restoring forces and finally collapses under the action of gravity. No displacement occurs in the lower left part of the column. A shear band seems to be formed between the two areas, and the material points above it continue to fall and move along the band until a final equilibrium system is reached at time t = 5.0 s. The accumulation is similar to a triangle with a relatively narrow front end and a depressed upper surface.

3.3. Results and Discussion

Both the sliding distance and the morphology in the final configuration are compared between the numerical solutions and laboratory results. The total displacement of the simulation model is 19.1 cm, which is similar to the experiment, with the sliding distance being measured as 19.2 cm, as shown in Figure 4. As for the comparison of the morphology, Pearson correlation coefficient [36, 37] is chosen to characterize the similarity between the experimental and numerical results. The correlation coefficient . The closer the correlation coefficient is near to 1, the higher the similarity is [38]. By selecting 52 sampling points at the same positions on both the experimental and simulation morphology contours, a linear correlation analysis is carried out. The similarity coefficient between them is 0.945, exhibiting a high similarity. To this end, the ability and adequacy of the MPM in modeling the soil flows underwater are clearly demonstrated.

4. Investigation of the Southern Mediterranean Submarine Landslide

Submarine landslides usually end up in a longer runout distance and affect a wider range of areas than the ashore landslides. In order to have a better understanding of the kinematics of submarine landslides, a real case of submarine escarpment failure in the Southern Mediterranean is selected in this study. The Mediterranean Sea is a land-locked basin that is called a “natural geological laboratory” due to a variety of sedimentary and tectonic environments coexisting in the same place. Submarine landslides are ubiquitous in this region [39]. The margin type of Southern Mediterranean is mainly the transcurrent tectonic margin. Besides, this region is considered as a tectonically active area according to a map of instrumentally recorded seismicity [40].

The detected slope morphology of the submarine landslide is shown in Figure 5. The detailed information of this real case is provided by Dong et al. [41], which was obtained from a routine exploration and the detailed information of the location is difficult to backtrack now. A dotted line and a solid line are utilized here to represent the profile before and after failure, respectively, surmised from the interpretation of the adjacent transept of submarine deposits. Elevations of the idealized transect are shown in Figure 6. The landslide occurred in the Southern Mediterranean at a depth of 630–690 m below the sea level, sliding from a steep slope at an inclination of 13.75° to a gently inclined (at an inclination of 1.1°) horizontal base, finally ending up with a runout distance of 420 m. Soil parameters were obtained through laboratory tests after the soil samples were retrieved from the study site, with the buoyant density of the soil  = 525 kg/m3, cohesion  = 15 kPa, sensitivity (the strength ratio between the intact and fully remoulded state)  = 40, and viscosity coefficient  = 0.0167. The acceleration of gravity is  = 9.81 m/s2. In order to approximate constant volume under the assumed undrained conditions, Poisson’s ratio is taken as 0.49. Young’s modulus of the soils is taken as 3000 kPa.

Several predisposing factors contributed to the destabilization of this slope. The most obvious one is the steep slope itself, which can be a result of a relatively high sedimentation rate [42]. The rapid deposition of sediment will generate excess pore water pressure in the underlying sediments, resulting in an unstable state of the slope. On the other hand, the study area is tectonically active with relatively high earthquake frequency. According to Ai et al. [43], a moderate earthquake (M = 4.0–4.8) can generate enough horizontal acceleration to cause failure in submarine slopes, which makes it highly possible that seismic activities also play a role in the destabilization of this slope.

4.1. Runout Process Simulation

A single-layer two-phase material point model based on the Mohr-Coulomb failure criterion for soil is used. A total stress analysis is utilized due to the undrained situations. Both physical and mechanical parameters of the slide necessary for the simulation are shown in Table 2. The damping ratio here is chosen as 0.05.

After a convergence study, a mesh spacing of 10 m for this simulation is selected. A total of 2613 material points are generated for the numerical model, which are distributed uniformly in the background grid, with 10 material points per cell in the slide body and 4 material points per cell in the bedrock.

The whole simulation process is shown in Figure 7, with the coloring indication representing the displacement relative to the initial configuration. It can be seen that the whole landslide starts as a translational movement along the inclined seabed, as represented by Figures 7(a)7(c), where a uniform displacement field, without perceptible differences, can be found within the sliding body. After entering the gently inclined horizontal area, the relative displacement among different parts within the sliding body becomes larger. The points at the front part move further away as the time goes, while the points at the back gradually stop, rendering a diffusive pattern of the movement, as shown in Figures 7(d)7(f). Eventually, at time t = 57.0 s, the whole motion ceases, as shown in Figure 7(g).

Likewise the sand-column collapse example, the sliding distance and the morphology in its final configuration are compared to verify the simulation results [44, 45]. The sliding distance of the simulation is 420.3 m, similar to the field data, which is recorded as 420.0 m, as shown in Figure 8. The similarity coefficient of the landslide morphology between them is 0.983. Moreover, the sliding time and maximum velocity are 57 s and 14.52 m/s, consistent with the simulation results obtained by Dong et al. [41] who adopted the strain-softening constituive model, which are around 60 s and 14 m/s. The correlation coefficients between them are 0.95 and 0.964, respectively.

4.2. Analysis of the Kinematics

As shown in Figure 6, three monitoring points PA, PB, and PC, located at the slope top, middle, and bottom parts, respectively, are selected to study the motion characteristics of the submarine landslide. Figure 9 depicts the velocity curves of the 3 monitoring points during the whole runout process. It can be seen that the 3 monitoring points accelerated immediately once the gravity is applied. At first, the velocities tend to be uniform and fluctuate in a small range. After approximately 10 s, differences of the velocities begin to be shown, and, gradually, PA obtains the highest velocity among the 3 measuring points, reaching 14.52 m/s; PB, in the middle of the slope, is the second; PC, at the slope shoulder, is the lowest, but with little difference between PB and PC. It is also interesting to find that the material points in the back of the slide body, that is, PB and PC, stop before the material point in the front, that is, PA. This may be due to the fact that the points at the back are sliding on top of the points at the front part rather than on the seabed directly, thereby causing more resistances to the points at the back. Figure 10 shows the displacement curves of the measuring points, where it is shown that the sliding distance of PA at the foot of the slope is the longest, reaching 420.3 m; PB is the second and Pc has the shortest sliding distance. By combining the displacement contours of the sliding body during the whole runout process, it may be concluded that the slope failure is a progressive one. Initially, the slide moves in a uniform way, leading to a translational slide of the whole body, which can be verified by the same displacements seen on the measuring points as depicted in Figure 10. Due to the velocity differences shown afterwards, the whole body movement is shifted to be a diffusive sliding pattern in the horizontal direction, with the front end moving far away, leaving potential gaps within the sliding body.

4.3. Slope Failure Mechanism

The mechanism of slope failure is of great importance for assessing geologic hazards caused by submarine landslide. In this section, the underlying mechanism will be discussed with the above simulation results. So far, knowledge of mechanism involved in submarine landslide generation remains poorly known. One common practice in the analysis is to borrow some established theories from subaerial landslides.

Submarine landslides are initiated when the shear stress downslope exceeds the shear strength of the material forming the slope. A failure criterion that represents the material’s physical characteristics is required to define the shear strength. The most general choice here is the Mohr-Coulomb failure criterion:where is the shear strength, is the friction angle, is the effective cohesion, is the stress that acts normal to the failure surface, and is the pore water pressure. The term is generally summarized as , known as the effective normal stress. All causes of submarine landslides can be concluded as two types: reduction of shear strength and increase of shear stress. According to Mohr-Coulomb failure criterion, the shear strength and effective stress follow a linear relation. Thus, any cause of a reduction of the effective stress can lead to a reduction of the shear strength.

Three types of shear stress downslope are significant with regard to submarine landslides: gravity, storm-wave-induced stress, and seismically induced stress. Detailed site-specific analyses must be performed in order to evaluate the postfailure behaviors of submarine slopes considering three types of stress, which is beyond the scope of our study. As it is a preliminary study, only gravity loading is considered in our simplified model. The seismic loading can generate vertical and normal accelerations and induce normal and shear stress in the slope. The storm-wave loading is much like that caused by an earthquake, which induces shear stress that varies cyclically between the crests and troughs, and may alternate in direction. These two types of loading are complex and cannot be described easily with a few parameters. However, it is certain that the presence of them will cause a reduction in the stability of submarine slopes [46].

From Figures 710, the total failure mechanism can be described in the following stages. At first is the local deformation. In natural slopes, the initial stress field is usually nonuniform due to the complex geological process during the slope formation [47]. Then the softening effect of seawater reduces the shear strength of sediments. The seawater and gravity loading increase the differences in the stress levels and the induced soil deformation is therefore nonuniform too. As a result, plastic deformation first occurs at the foot of the slope. This local deformation lasts for a short time. After the local plastic zone is formed, the deformation further enlarges due to the inability of soil to sustain the stress associated with gravity. The soil at the foot of the slope begins to creep forward as the shear stress exceeds the shear strength.

With the development of deformation, a sliding body is formed in the slope and it begins to slide along the seabed. The sliding body moves forward essentially as a block, with the height reducing very slowly. This stage lasts for about 13 seconds, and the velocity of each part of the landslide increases rapidly until it reaches the maximum velocity in the sliding process. The displacement and velocity curves are consistent with the experiment results by Olivares and Damiano [48], confirming that the soil mass is in a clear flow-like movement. The sliding body elongates and its thickness reduces. This flow-like movement is caused by the sudden increase in excess pore pressure, which reaches such a high value that the shear strength of soil almost vanishes, culminating in development of such a flowing mass with a high degree of fluidity.

Then the kinetic energy of the sliding body dissipates continuously due to the friction resistance of the interface exceeding the driving force caused by gravity. This stage lasts for a relatively long period (about 35 s). Although the velocity is decreasing, there is still considerable displacement in this stage. The displacement curve changes from upper concave to lower concave. Finally the movement of the landslide has almost stopped and the displacement curves of each point tend to a certain value. According to Whitman’s terminology [15], this submarine landslide can be classified as a disintegrative failure as its deformation caused by a transient loading event continues throughout the whole sliding process.

4.4. Parametric Analysis

Orthogonal experiments are used for parametric analysis in this section [49, 50], but note that the parameters variation should be within a reasonable range [51]. Due to the lack of detailed survey data for this submarine landslide, refer to [52, 53] for the range of values taken by previous researches, and it is found that 75%∼125% of the benchmark values for the parametric analysis are considered reasonable for the selection.

Six influential factors, namely, initial porosity, Young’s modulus, buoyant density, cohesion, friction coefficient, and slope angle, are chosen for the sensitivity analysis. Based on the principle of orthogonal experiments, the following orthogonal experimental table with 6 factors and 5 value levels is designed, leading to a total of 25 MPM analyses. The maximum sliding distance and velocity are recorded, with the results summarized and shown in Table 3.

After completing 25 numerical experiments and collecting all the data, the range analysis and variance analysis are conducted to study the influential degree of 6 factors on the target indices. The statistical results of the range analysis and variance analysis are shown in Tables 4 and 5, respectively, which are consistent with each other. According to the values of range in the range analysis and the F values in the variance analysis, the sensitivity of the influential factors on the sliding distance, in a decreasing order, can be concluded as follows: friction coefficient, slope angle, cohesion, buoyant density, Young’s modulus, and initial porosity. Meanwhile the sensitivity of the 6 factors on the maximum velocity is as follows: slope angle, friction coefficient, cohesion, buoyant density, initial porosity, and Young’s modulus.

The result of parametric analysis shows that the postfailure behavior of submarine landslides depends on the combined influences of different factors. The parameters can be classified into two types: the inner and the outer factors. The inner factors are the physical and mechanical properties of a slope. They can affect the initial stability of a submarine slope and the runout mechanism after a landslide is triggered by external factors. A sliding mass with high initial porosity will generate intense excess pore pressure, which can greatly reduce the effective stress as well as the shear strength, since they follow a linear relation. From the parametric analysis, a slope with high initial porosity and low cohesion and buoyant density often ends up in a higher velocity and longer runout distance as its shear strength decreases significantly and can easily change from a block sliding into a flow-like movement.

The outer factors include the slope angle and friction coefficient, which are exactly the most significant two influential factors to the slide. The slope angle represents the geomorphology of the slope when the slide initially happens, and the friction coefficient refers to the pathway resistance to the slide. During the sliding process, the initial gravitational potential would be transferred to the kinetic energy due to the particle movement and the strain energy, and part of dissipation due to the friction occurred in between the slope base and the floor. The slope angle can affect the initial morphology, and thereby a higher slope angle may represent a higher gravitational potential, from which more kinetic energy could be transferred, that is, higher velocity. Likewise, the friction would have very few effects on fluid flow with high Reynolds when its velocity is high. As the velocity decreases, the influence of friction becomes more dominant. Moreover, a higher friction coefficient corresponds to a higher energy dissipation rate, which may influence the time the friction takes to have an effect. Therefore, the friction appears to control the sliding distance the most.

5. Conclusion

By utilizing a hydromechanical coupled MPM formulation, two cases of soil flows fully under water have been analyzed, in which the ability and adequacy of MPM in modeling the submarine landslides are demonstrated. As for the South Mediterranean landslide case, the simulation results match well with the real situation and it can be concluded that the failure mode of the slide is a disintegrative one. The seawater and gravity loading increase the differences in the stress levels of the initial nonuniform strain field and triggered the landslide. Plastic deformation first occurs at the foot of the slope, and the deformation zone becomes larger and extends into the slope body. Due to the sudden increase in excess pore pressure, the block sliding gradually evolves into a flow-like movement with the sliding body elongating and its thickness reducing, ending up with a high maximum velocity of 14.52 m/s and final runout distance of 420.3 m. Parametric analysis results show that the postfailure mechanism of submarine landslides depends on the combined influences of different parameters. The slope angle may influence the maximum sliding velocity in the greatest extent; meanwhile, for the sliding distance, the friction coefficient in between the seabed and sliding materials may have the biggest influence.

Our study proves that MPM is suitable for simulating large deformation and multiphase problems such as submarine landslides by employing both the Lagrangian and Eulerian descriptions. Its ability to reproduce the entire runout process provides an insight into analyzing the postfailure mechanism of submarine slopes, which helps people better assess the risk caused by a potential landslide in offshore or abyssal areas. Our work presents a preliminary application of MPM in studying submarine landslides. However, a fully coupled MPM model still requires further development for refinement. More complex factors should be taken into consideration such as the interaction between the flowing mass and submarine structures to ensure accurate simulation results in geotechnical engineering problems.

Data Availability

The data (including all the figures, tables, and MPM simulation results) used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.


The authors would like to acknowledge the financial supports of the National Key R&D program of China (Grant no. 2018YFC1505104), the CAS Pioneer Hundred Talents Program, the Natural Science Foundation of Jiangsu Province (Grant no. BK20181182), and the Science & Technology Program of Suzhou (Grant no. SYG201613). They also thank Professor Dong Youkou from the China University of Geosciences for kindly providing them with geological information of the submarine landslide in the Southern Mediterranean.