Numerical Investigation of Surge Waves Generated by Submarine Debris Flows

Abstract

Submarine debris flows and their generated waves are common disasters in Nature that may destroy offshore infrastructure and cause fatalities. As the propagation of submarine debris flows is complex, involving granular material sliding and wave generation, it is difficult to simulate the process using conventional numerical models. In this study, a numerical model based on the smoothed particle hydrodynamics (SPH) algorithm is proposed to simulate the propagation of submarine debris flow and predict its generated waves. This model contains the Bingham fluid model for granular material, the Newtonian fluid model for the ambient water, and a multiphase granular flow algorithm. Moreover, a boundary treatment technique is applied to consider the repulsive force from the solid boundary. Underwater rigid block slide and underwater sand flow were simulated as numerical examples to verify the proposed SPH model. The computed wave profiles were compared with the observed results recorded in references. The good agreement between the numerical results and experimental data indicates the stability and accuracy of the proposed SPH model.

1. Introduction

Submarine debris flows are widely distributed on continental shelves, continental slopes and in deepwater areas, where they pose a serious threat to offshore infrastructure, such as submarine pipelines and cables, offshore oil and gas platforms, and offshore wind farms [1]. In addition, submarine debris flows can generate huge surge waves that can seriously threaten the safety of infrastructure and people living in coastal areas. For example, the 1958 Lituya Bay debris flow-generated tsunami produced a run-up of over 400 m [2]. Another significant tsunami was generated by a submarine debris flow at the Nice airport in France on 16 October 1979, and swept away 11 people [3]. The 1994 Skagway debris flow-generated tsunami destroyed a railway dock, damaged the harbor, and killed a construction worker in Alaska, USA [4,5]. A submarine debris flow in Papua New Guinea generated a wall of water measuring 15 m high that killed over 2100 people in July 1998 [6]. Furthermore, the 2006 Java tsunami, generated by a submarine debris flow off of Nusa Kambangan Island, produced a run-up in excess of 20 m at Permisan [7]. In 2010, the slope failures of river deltas along the Haitian coastline triggered tsunamis with a height of 3 m, which caused at least three fatalities and some damage to infrastructure [8]. In 2018, the submarine debris flow induced by the Palu earthquake resulted in devastating tsunamis and caused more than 2000 fatalities in Sulawesi, Indonesia [9,10,11]. These debris flow-generated tsunami disasters have therefore attracted a lot of attention from researchers.

Catastrophic tsunami waves generated by submarine debris flows are difficult to observe and record because of their inaccessibility and unpredictability. Current approaches to investigating the surge waves generated by submarine debris flows mainly focus on physical model tests and numerical simulations. For instance, physical model tests were conducted based on the generalized Froude similarity to study the complex wave patterns generated by debris flows at Oregon State University [12]. Later, McFall and Fritz [13] conducted physical modeling using gravel and cobble materials in different topographic conditions, and analyzed offshore tsunami waves. Wang et al. [14] conducted physical model tests to analyze the effect of volume ratio on the maximum tsunami amplitude, and obtained an optimization model to predict the maximum tsunami amplitude. Miller et al. [15] designed a large flume model test to measure the shape and amplitude of the near-field waves and investigate the fraction of debris flow mass that activates the leading wave. Meng [16] conducted a series of model tests to investigate the waves generated by debris flow at laboratory scale, and the effect of the sliding mass on the wave features was analyzed. Takabatake et al. [17] conducted subaerial, partially submerged, and submarine debris flow-generated tsunamis to investigate the different characteristics of these three types of tsunamis. Although some promising results have been obtained, the physical model tests conducted in the abovementioned work require a great deal of manpower and material resources. In addition, the size effect is an inevitable problem in physical modeling.

With the rapid development of computer techniques and numerical algorithms, numerical modeling of such events is critical to better understand the mechanism of surge waves and predict their propagation behavior, and become more prevalent in recent years. For example, Lynett and Liu [18] derived a numerical model to simulate the waves generated by a submarine debris flow and to predict the run-up. Yavari–Ramshe and Ataie–Ashtiani [19] presented a comprehensive review on past numerical investigations of debris flow-generated tsunami waves. In this review study, numerical models based on depth-averaged equations (DAEs) were tabulated according to their conceptual, mathematical, and numerical approaches. Recently, a novel multiphase numerical model based on the moving particle semi-implicit method was proposed and applied to simulate submerged debris flows and tsunami waves. The model was validated and evaluated through comparison with the experimental data and other numerical models [20]. Sun et al. [21] solved the Navier–Stokes equations and incompressible flow continuity equation and developed a numerical model to investigate the propagation of a debris flow-generated tsunami under different initial submergence conditions. Rupali et al. [22] developed a coupled model based on the boundary element method and spectral element method to describe the characteristics of debris flow-generated waves. Baba et al. [23] adopted a two-layer flow model to simulate the submarine mass movement and surged wave propagation. In addition, computational fluid dynamics (CFD) has been widely applied to analyze tsunami wave propagation [24,25]. In the field of CFD, a mesh-free particle method named smoothed particle hydrodynamics (SPH) was proposed as an astrophysics application and has been widely used in various engineering fields [26]. For example, Iryanto [27] established an SPH model to simulate the waves induced by aerial and submarine landslides. Wang et al. [28] proposed a coupled DDA–SPH model to deal with the solid–fluid interaction problem, and applied it to simulate a block sliding along an underwater slope. However, the landslide mass was simulated as a rigid body, and its deformation was neglected in these two models. This work aimed to establish an SPH model to simulate surge waves generated by submarine debris flows, and to show the accuracy and stability of this model in the simulation of multiphase flow problems.

2. Materials and Methods

2.1. SPH Algorithm

The SPH method was proposed in 1977 for astrophysical applications [29]. Recently, this method has been widely applied to a large variety of engineering fields [30,31,32,33]. Compared to the mesh-based methods, the major advantage of this method is that it bypasses the need for numerical meshes and avoids the mesh distortion issue and a great deal of computation for mesh repartition [34].

In the SPH method, the governing equations can be transformed from partial differential equation form into SPH form by two approximation methods. The first one is to produce the functions using integral representations, which is usually called kernel approximation. The second one is to discretize the problem domain into a set of particles, which is called particle approximation. Based on these two approximations, field variables and their derivative in SPH form can be expressed by the following equations [35]:

$$f\left({\mathit{x}}_i\right)\approx {\displaystyle \sum _{j=1}^N\frac{m_{j}f\left({\mathit{x}}_j\right)}{{\rho }_j}W\left(\left|{\mathit{x}}_i-{\mathit{x}}_j\right|,h\right)}$$

(1)

$$\frac{\partial f\left({\mathit{x}}_i\right)}{\partial {\mathit{x}}_i}\approx -{\displaystyle \sum _{j=1}^N\frac{m_{j}f\left({\mathit{x}}_j\right)}{{\rho }_j}\cdot \frac{\partial W\left(\left|{\mathit{x}}_i-{\mathit{x}}_j\right|,h\right)}{\partial {\mathit{x}}_i}}$$

(2)

where f is an arbitrary field function, x is the coordinates of the SPH particles, m is the mass of particles, ρ is the density, W is the smoothing function, h is the smoothing length, and N is the total number of neighboring particles.

2.2. Governing Equations

In the presented model, soil and water in a large deformation situation are assumed as two different fluids that occupy the space by a certain volume fraction. They satisfy the conservation law of mass and momentum, respectively [36]:

$$\frac{1}{{\tilde{\rho }}_s}\frac{D\left({\tilde{\rho }}_s{\varphi }_s\right)}{Dt}+{\varphi }_s\nabla \cdot {\tilde{{v}}}_s=0$$

(3)

$$\frac{1}{{\tilde{\rho }}_w}\frac{D\left({\tilde{\rho }}_w{\varphi }_w\right)}{Dt}+{\varphi }_w\nabla \cdot {\tilde{{v}}}_w=0$$

(4)

$${\varphi }_s\frac{D{\tilde{{v}}}_s}{Dt}=\frac{1}{{\tilde{\rho }}_s}\nabla {\tilde{\mathbf{\sigma }}}_s+{\varphi }_{s}g+\frac{1}{{\tilde{\rho }}_s}{\mathit{f}}_s$$

(5)

$${\varphi }_w\frac{D{\tilde{{v}}}_w}{Dt}=\frac{1}{{\tilde{\rho }}_w}\nabla {\tilde{\mathbf{\sigma }}}_w+{\varphi }_{w}g+\frac{1}{{\tilde{\rho }}_w}{\mathit{f}}_s$$

(6)

where t is the time, g is the gravitational acceleration. v denotes the velocity vector, and σ denotes the stress tensor. The subscripts w and s represent water and soil, respectively. Moreover, φ_s and φ_w are the volume fractions of soil and water, and f_s is the force exerting on the soil by the water, which can be determined by:

$${\mathit{f}}_s=-{\varphi }_s\nabla p+{\mathit{f}}_d$$

(7)

where f_d is the Darcy penetrating force, and p is the pressure term of the fluid, which is calculated by an equation of state [37]:

$$p=\frac{{\rho }_0{c}_s^2}{\gamma }\left[{\left(\frac{\rho }{{\rho }_0}\right)}^{\gamma }-1\right]$$

(8)

where ρ₀ is the reference density, which can be measured through laboratory tests; ρ is the density obtained through the mass conservation equation; c_s is the sound speed, which can be set equal to ten times the maximum velocity [38]; and γ is the exponent of the equation of state, and is usually set to 7.0 [39].

2.3. Material Model

In the presented SPH model, the submarine debris flow is assumed to be a Bingham viscous fluid, and the relationship between the shear stress and the shear strain rate is defined as:

$$\tau =\eta \dot{\gamma }+{\tau }_y$$

(9)

where τ denotes the shear stress, η represents the viscosity coefficient in fluid dynamics, $\dot{\gamma }$ means the shear strain rate, and τ_y is the Bingham yield stress.

The ambient water is assumed as a classical Newtonian fluid, and the deformation behavior can be described as:

$$\tau =\eta \dot{\gamma }$$

(10)

The viscosity of water η is 1.7 × 10⁻³ Pa·s.

2.4. Approach for Multiphase Granular Flow Modeling

In this SPH model, the approach for multiphase granular flow proposed by Tajnesaie et al. [20] is adopted. As shown in Figure 1, the debris flow is represented by granular-type (in brown color) particles, and the water is represented by fluid-type (in blue color) particles. As the granular flow is quite fast, and the permeability is low, it is assumed in this model that the mass exchange between the ambient and pore water is negligible. Each particle has its own physical and mechanical properties. The fluid-type particles are regarded as Newtonian fluid, and the viscosity is constant. In contrast, the granular-type particles are regarded as Bingham viscous fluid, and the viscosity is variable and determined by the rheological model, as described in Section 2.3.

2.5. Treatment on Water–Soil Interface

As shown in Figure 2a, density and mass are discontinuous near the water–soil interface, reducing the computational accuracy and stability when solving the governing equations. Therefore, boundary treatments on the water–soil interface should be made when solving the governing equations.

As shown in Figure 2a, a is the central particle (in red) and b represents the neighboring particles of the other phase (in yellow). In this model, particle b is converted into particle b’ (in blue) of the same phase with particle a, as shown in Figure 2b, which occupies the same position, velocity, and volume as the original particle b. Only the position, volume, velocity, and pressure of these neighboring particles are taken into account in particle approximation.

2.6. Boundary Condition

During the submarine debris flow propagation, the fluid particles face resistance from the local seabed. Therefore, to truly simulate the debris flow motion, it is important to accurately calculate the repulsive force from the solid boundary. In the SPH method, repulsive particles are widely used on the boundary contours to exert a force on the fluid particle in the direction of the line connecting both particles, as shown in Figure 3.

The repulsive force, F_ij, can be determined by Equation (11) [40]:

$$F_{ij}=\{\begin{array}{l}D\left[{\left(\frac{r_0}{\left|{\mathit{r}}_{ij}\right|}\right)}^{n_1}-{\left(\frac{r_0}{\left|{\mathit{r}}_{ij}\right|}\right)}^{n_2}\right]\frac{{\mathit{X}}_{ij}}{{\left|{\mathit{r}}_{ij}\right|}^2},\begin{array}{cc}& if\begin{array}{cc}& \frac{r_0}{\left|{\mathit{r}}_{ij}\right|}\ge 1\end{array}\end{array}\\ 0\begin{array}{cc}& if\begin{array}{cc}& \frac{r_0}{\left|{\mathit{r}}_{ij}\right|}<1\end{array}\end{array}\end{array}$$

(11)

where D is a parameter that depends on the highest flow velocity, r₀ is the cutoff distance, |r_ij| is the distance between a fluid particle and a boundary virtual particle, X_ij is the difference between the position vectors of a fluid particle and a boundary virtual particle, and n₁ and n₂ are parameters that are usually equal to 12 and 4, respectively.

3. Model Validation and Application

3.1. Waves Generated by Underwater Rigid Block Sliding

To validate the developed SPH model, an underwater rigid block sliding experiment performed by Heinrich [41] was simulated, and the numerical results were compared with test data. Figure 4 shows the experimental setup. In the experiment, a block slid down a slope of 45° on the horizontal and generated water waves. The shore was simulated by another incline with a 15° slope. Water surface was level with the intersection of these two inclines, and the water depth in the channel was 1 m. The block was triangular in the cross-section with dimensions of 0.5 m × 0.5 m. The density was 2000 kg/m³. During the experiment, the block slid into the water under the effect of gravity, and then stopped when it ran into a small buffer at the channel bottom. As shown in Figure 4, the computational domain was 5.0 m × 1.3 m in the x and y directions.

Table 1. Parameters used in the SPH simulation of the rigid block sliding experiment.

Rigid block density	ρs (kg/m3)	2000
Water density	ρw (kg/m3)	1000
Water viscosity coefficient	η (Pa·s)	1.7 × 10−3
Acceleration of gravity	g (m/s2)	9.8

lists the calculating parameters used in the SPH simulation.

Figure 5 shows the SPH simulated results for the rigid block sliding experiment. When the block slid down, surge waves were generated because the kinetic energy of the sliding block was transferred to the ambient water body. The surge waves propagated out from the near field to the far field because of the height differences of the water surface. Energy conversion finished after the rigid block stopped, and then both the kinetic and potential energy of the water waves started to decay. The positions of the block and the wave profiles at different times were observed.

To verify the simulation accuracy of the SPH model, the numerical results were compared with the experimental data from Heinrich [41]. Figure 6 compares the numerical and experimental vertical displacement time history of the rigid block. The numerical results coincided with the experiment data well. After an acceleration phase at the beginning, the rigid block reached a constant velocity of 0.6 m/s. At approximately 1.0 s, the block reached the flume bottom in the SPH model as well as in the experiment. Figure 7 shows the comparison between the simulated and observed wave profiles at different times (t = 0.5 s, t = 1.0 s and t = 1.5 s). Although some discrepancy can be observed, there is close overall agreement between the experimental and numerical wave profiles. The maximum run-up of the wave during the simulation was about 0.10 m, which matches the experimental data well.

To evaluate the performance of the presented model furtherly, the numerical results are compared with those calculated by Tajnesaie et al. [20] using a MPS model, and Wang et al. [28] using a coupled DDA-SPH model. Figure 8 compares the water profiles calculated based on the above-mentioned numerical models at t = 0.5 s and t = 1.0 s. It shows that all the numerical models produce similar water profiles which are consistent with that observed in the experiments.

To quantitatively compare the results, the Pearson correlation coefficient of the experimental and numerical water profiles in Figure 8 are calculated. In this case, 40 points with an equal horizontal distance on the water profile are selected as the analysis objects. Assuming that the vector X (x₁, x₂, x₃, ……, x₄₀) represents the calculated elevations of these points, and Y (y₁, y₂, y₃, ……, y₄₀) represents the measured ones, the Pearson correlation coefficient of the two vectors can be calculated as:

$$Corr\left(\mathit{X},\mathit{Y}\right)=\frac{{\displaystyle \sum _{i=1}^n\left(x_i-\overline{x}\right)\left(y_i-\overline{y}\right)}}{\sqrt{{\displaystyle \sum _{i=1}^n{\left(x_i-\overline{x}\right)}^2\times {\displaystyle \sum _{i=1}^n{\left(y_i-\overline{y}\right)}^2}}}}$$

(12)

where $\overline{x}$ and $\overline{y}$ are the average values of the two data sets, respectively.

According to Equation (12), the correlation coefficients between the experimental and numerical results are listed in

Table 2. Correlation coefficients between the experimental and numerical results on the water profiles.

Case	Current Study	Tajnesaie et al. [20]	Wang et al. [28]
t = 0.5	0.9336	0.9389	0.9297
t = 1.0	0.9132	0.9206	0.9040

. It shows that all the correlation coefficients calculated are close to 1, which means that the numerical results have a strong correlation with the experimental profiles. Therefore, the presented SPH model can successfully simulate such complex flows.

3.2. Waves Generated by the Underwater Sand Flow

In Section 3.1, the block was assumed to be rigid, which is not consistent with the actual situation of submarine debris flow. The large deformation of the solid phase should be taken into account in the simulation of the submarine debris flow propagation. Therefore, an experiment of the surge waves generated by underwater sand flow, conducted by Rzadkiewicz et al. [42], was simulated. Figure 9 shows the experimental setup. In this experiment, the maximum water depth was 1.6 m. Under the influence of gravity, a mass of sand (0.65 m × 0.65 m in cross section) slid down along a slope of 45°, and generated waves propagated outward.

In the underwater sand flow experiment, the density of the sand was 1950 kg/m³, and the grain diameter ranged from 2 to 7 mm. In the numerical simulation, the sand flow was assumed to be a Bingham fluid. The parameters are listed in

Table 3. Parameters used in the SPH simulation of underwater sand flow.

Density of sand	ρs (kg/m3)	1950
Viscosity coefficient of sand flow	ηg (Pa·s)	0.15
Bingham yield stress of sand flow	τy (Pa)	750
Density of water	ρw (kg/m3)	1000
Viscosity coefficient of water	ηw (Pa·s)	1.7 × 10−3
Acceleration of gravity	g (m/s2)	9.8

. In the absence of measurements, the viscosity coefficient and yield stress of Bingham model can be determined by trial and error. According to Ataie–Ashtiani and Shobeyri [43], under the conditions of viscosity coefficient η_g = 0.15 Pa·s and Bingham yield stress τ_y = 750 Pa, the numerical results match the experimental data well. Therefore, the same parameter values were used in this study.

Figure 10 shows the SPH simulated results for the underwater sand flow. Unlike the previous numerical example, the solid phase no longer behaved as a rigid block but deformed significantly when sliding down along the slope. The profiles of water and sand slope at different times are presented. The sand slope keeps its initial shape at the beginning of slide motion, as shown in Figure 9a,b, while at t = 0.8 s, most of the sand particles are concentrated at the flow front, as shown in Figure 9d.

Figure 11 compares the simulated slope configurations with experimental data [42] and those calculated by the VOF model proposed in [42] and the ISPH model proposed in [43]. It shows that the simulated slope configurations are coincident with experimental results.

Table 4. Correlation coefficients between the experimental and numerical results on the water profiles.

Case	Current Study	VOF Results [42]	ISPH Results [43]
t = 0.4	0.9092	0.8641	0.9197
t = 0.8	0.8761	0.8318	0.8804

lists the Pearson correlation coefficients of the numerical and experimental slope configurations shown in Figure 11. The comparison shows that the SPH model proposed in this study has the equivalent simulation accuracy with the existing numerical models.

4. Conclusions

Submarine debris flow-generated water waves are natural phenomena that occur under certain conditions and can result in damage of infrastructure and loss of life in coastal areas. In this study, an SPH model with a multiphase granular flow algorithm is presented for numerical simulation of surge waves generated by submarine debris flows. The Bingham fluid model and Newton fluid model are used to describe the motion behavior of submarine debris flow and ambient water, respectively. A simple treatment is proposed to ensure the continuity of the density and pressure near the interface of these two fluids and to avoid numerical instability.

This model was first used to simulate an experiment of underwater rigid block sliding. The simulated wave profiles were compared with the observed results to verify the stability and accuracy of the SPH model.

Then, this model was applied to simulate submarine debris flow. Surged waves generated by the underwater sand flow were computed and compared with the tested results. The good agreement between the numerical and tested results proves that the proposed SPH model is suitable to simulate such complex flows.

In fact, submarine debris flows move across 3D submarine topography and may change direction and split or join in response to the complex topography. Therefore, the 2D model presented in this study cannot truly reproduce the complex dynamic process, and a 3D model with parallel computing techniques is necessary to improve the calculation accuracy and efficiency.