Simulating landslide-induced tsunamis in the Yangtze River at the Three Gorges in China

Abstract

Landslide-induced tsunamis may cause fatalities, damages and financial losses. In the Three Gorges Reservoir Area of China, several large landslides are still unstable and persistently creeping toward the Yangtze River. In this paper, we investigate the impacts of landslide-induced tsunamis in the Three Gorges Reservoir by using a hybrid numerical approach. One of the largest unstable mass in this area, the Huangtupo landslide, is chosen as the study object. First, the landslide deformation and initiating velocities are obtained by using the finite-discrete element method. The landslide-induced tsunamis and their impacts on shipping on the Yangtze River are then investigated through smooth particle hydrodynamics modelling. Our results reveal that an approximately 80% reduction in shear strength of the tip in the landslide will lead to catastrophic failure of the landslide, with sliding velocities of up to 8 m/s. Subsequently, such a collapse may initiate a river tsunami, propagating up to 9 m on the nearby reservoir banks within 3 km. The impacts on surrounding floating objects, such as surges and sways, heaves and rolls, are up to 110 m, 8 m and 6°, respectively. The simulations indicate that although the likelihood of a catastrophic failure of the whole landslide is low, the partial sliding still poses severe threat to the nearby reservoir banks and shipping on the Yangtze River. Thus, we recommend continuous monitoring as well as landslide early warning systems at this and also other hazardous sites in this area.

1 Introduction

With the world’s largest hydropower station, the Three Gorges Dam is the key project to manage and develop the Yangtze River in China. However, the narrow and long shape reservoir behind the dam with approximately 660 km length reactivated severe geohazards such as landslides, ground subsidence and earthquakes [23]. Landslides, rock falls and debris flows occur most frequently in the Three Gorges Reservoir Area (TGRA). It is reported that 90% of the 4000 documented geohazards are landslides and rock falls with total volume of up to 4 billion $\hbox {m}^3$ [43]. The failure of landslides results in not only damage such as deformation, collapse, and destruction of buildings and facilities but also impulse waves (or tsunamis), causing great harm to human live, property and important water infrastructure [18, 53]. The failure of the Vajont landslide in 1963, for instance, caused an impulse wave with a maximum surge height of 175 m, overflowing the dam and killing more than 2600 people in the Italian Alps [34].

Similar catastrophic events have also been witnessed in the TGRA of China. Since the impoundment of the reservoir in 2003, there have been two major catastrophic landslide events: the Qianjiangping and Gongjiafang landslides [45, 52]. The Qianjiangping landslide was located close to the estuary of Qinggan River, a main tributary of the Yangtze River. The landslide occurred on 13 July 2003 after the reservoir water level rose from 95 to 135 m above sea level (a.s.l.). A rock mass of 24 million $\hbox {m}^3$ and more than 250 m wide slid into the river, killing 24 people and destroying 346 houses and four factories. The recorded data indicated that the highest induced wave was up to 30 m [45]. The Gongjiafang slope, located on the north bank of a main stream of the Yangtze River, collapsed on 23 November 2008. The fatal sliding process of $380{,}000\,\hbox {m}^3$ of rock was recorded by the smartphone of a passerby. According to the video, the maximum velocity of the sliding mass reached 12 m/s, and the highest wave amplitude was 32 m with a propagation velocity of approximately 18 m/s. The subsequent investigation revealed that the wave height remained 1.1 m at 4 km away. This incident fortunately caused no casualties, but economic losses of approximately 5 million RMB ($\sim 650{,}000 \,\hbox {euro}$) were incurred [18].

Currently, many large landslides in the TGRA remain unstable and exhibit a creep movement pattern with episodic acceleration and deceleration [43]. Although, some of these landslides were partially reinforced to slow their movement, the potential risk of catastrophic collapse is considerably high, owing to the periodic variations in water pressure [48], sliding materials deterioration [10] and/or extreme events such as earthquakes and rainstorms [54]. Nevertheless, analyses, simulations and predictions of the risk of impulsive waves caused by large landslides in the TGRA remains major challenges.

The threats to nearby water areas and banks are mainly caused by landslide-induced tsunamis [18, 53]. Since the properties of such impulsive waves are closely related to landslide geometry, volume, sliding velocity, topography and water level, physical models are often employed [15, 19]. On the other hand, owing to the ability to simulate large deformations and free fluid surfaces, discrete and point-based numerical methods, such as the finite-discrete element method (FDEM), particle finite element method (PFEM), material point method (MPM) and smoothed particle hydrodynamics (SPH), are also widely used for modelling landslide processes and impulsive waves [3,4,5,6].

Existing studies on landslide-induced tsunamis are mainly focused on the landslide sliding velocity, wave propagation, wave height, wave run-up, and related influencing factors. Most of these studies are based on failed landslides, and the simulation results can be verified with recorded data [17, 18, 52, 53]. However, few studies predict the induced tsunamis of high-risk landslides that have yet to occur. Moreover, the effects of landslide-induced tsunamis on the motion of floating objects such as ships in nearby water, are rarely reported, although this aspect remains an important issue for the shipping lanes in the TGRA.

In this paper, we study landslide failure, the induced tsunami and the subsequent impacts on shipping on the Yangtze River through a hybrid numerical approach. A large landslide in the TGRA, Huangtupo landslide, is taken as an example for this simulation. A 2D numerical model using the finite-discrete element method (FDEM) is adopted to predict the landslide kinetics at failure, which are then used as the initiation for smoothed particle hydrodynamics (SPH) modelling. The wave propagation and the dynamic motions of floating objects (i.e. ships) are numerically investigated by 2D and 3D simulations.

2 Background

The Huangtupo landslide is one of the largest and the most harmful landslides in the TGRA (Fig. 1). It has attracted worldwide attention due to its large size and complex structure, as well as the story behind it (e.g. [11, 14, 42, 47]). In the early 1980s, the Huangtupo area was selected as the location of Badong County, which was below the designed reservoir water level in the construction of the Three Gorges Project. However, the Huangtupo landslide was reactivated after the County seat was nearly established, resulting in a second resettlement of approximately 18,000 residents.

Fig. 1

Location (a) and images (b) of the study area. The photographs of Huangtupo area show the original view (in 2001), the view after reservoir impoundment and revetment construction (in 2014), and the view after evacuation (in 2019). The photographs of the Leijiaping area show the residential area and the Earth surface deformation monitoring system (IBIS-FL) located on the opposite bank

To minimize the landsliding risk, an evacuation of the residents on the Huangtupo was started in 2008. As shown by the insets in Fig. 1b, the residents were completely evacuated during the last decade and all buildings were cleared by the end of 2018. Nevertheless, the landslide still poses a significant threat to neighbouring residential areas. Once failure occurs, the landslide-induced tsunami may severely threaten the Lejiaping district, a small town with a population of 1200, located opposite to the Huangtupo landslide. In addition, the induced tsunami may wreck the vessels on the Yangtze River. Notably, there are hundreds of inland water vessels passing through the Huangtupo water area every day, and the number of large ships has continuously increased since the opening of the Three Gorges ship lock [24]. Hence, the evacuation may partially eliminate the threats of landslides.

Fig. 2

Plan (a) and representative cross section M–M’ (b) of the Huangtupo landslide

The Huangtupo landslide is composed of several sliding masses: No. 1 riverside sliding mass, No. 2 riverside sliding mass, No. 3 substation landslide and No. 4 garden spot landslide (Fig. 2a). The No. 1 sliding mass can be further divided into two partially overlapping secondary sliding masses, namely No. 1-1 and No. 1-2 sliding masses (Figs. 2 and 3) with volumes of approximately 13 and 5 million $\hbox {m}^3$ and a maximum thicknesses of 90 and 65 m, respectively [46]. The bedrock at Huangtupo area consists of clastic and carbonate rocks of the Badong group, middle Triassic formation ($T_2b$), which can be divided into five segments. The general dip direction of bedrock is $\hbox {N25}^\circ \hbox {E}$ with the dip angle varying between $30^\circ$ and $50^\circ$. The front part of the riverside sliding masses (No. 1) is composed of loose rocks and soil debris, originating from the third segment of the Badong formation ($T_2b^3$). The rear and deep regions are disturbed rock mixed with gravels and soils (Fig. 2b). Developed between the disturbed rock and bedrock, the sliding zone of the riverside sliding masses is composed of greyish, yellow, silty clays with gravels. Other than the riverside sliding masses, the materials of the upper and lower parts of the No. 3 and No. 4 landslides are mainly fractured rocks originated from amaranth mud and siltstones of $T_2b^2$ and grey limestone of $T_2b^3$, respectively. The sliding zones are developed in the fractured rocks with different sources, and the material is mainly brownish, red, silty clay with gravels.

Fig. 3

a Earth surface deformation data from January 2016 to December 2017, b location of deep cracks and c cracks deformation data from 2018 to 2019

The deformation of the Earth’s surface caused by landslide movement was monitored by the IBIS-FL system, a ground-based interferometric synthetic aperture radar developed by IDS GeoRadar s.r.l. The monitoring system was installed in the Leijiaping district to cover the entire landslide (Fig. 1b). The surface displacement accumulated from January 2016 to December 2017 is shown in Fig. 3a. It reveals that the center of the No. 1 riverside sliding mass contributed the largest deformation, up to 9 cm. More evidence was observed beneath the landslide. The underground tunnel, intersecting the landslide shear zone, was constructed beneath the landslide in December 2012. Since then, cracks in the investigation tunnel have started to appear and exhibit continuous deformation of up to 1 mm per year (Fig. 3b, c; monitored cracks in the Main tunnel, No. 3 Branch Tunnel, No. 5 Branch Tunnel). The development of these cracks offers direct an evidence of the landslide movement along the basal shear zones. All evidence indicates that the No. 1-1 riverside sliding mass exhibits the largest deformation, revealing the most fragile part of the landslide [47, 50].

Currently, the Huangtupo landslide has become a key study object to better understand the mechanisms of similar slides in the TGRA. During the last two decades, there has been ongoing research on the Huangtupo landslide. A three-stage evolution model was proposed to illustrate the involved rock mass creep, primary landsliding and partial reactivation of this landslide in 2000 after Badong County experienced the first relocation [14]. In 2012 the Badong field test site was established. The test site includes a large-scale tunnel with auxiliary adits and an extensive in-situ monitoring system [43]. Since then, research on the internal structure, failure mechanisms and kinematic features has intensified. Through a uranium–thorium disequilibrium dating method, it was found that the riverside sliding masses underwent at least two periods of movement approximately 100 and 40 ka ago [41]. The structure and shear surface distribution were further investigated according to the tunnelling excavation [46]. The strength and rheological behaviours of the sliding zone materials were evaluated by different testing methods both in the laboratory and field [13, 40, 49, 50]. The results indicate that the low shear strength with continuous reduction of sliding zone materials is one of the fundamental mechanisms for landsliding. Apart from the gravitational effect and the defect of weak sliding zones, initial reservoir impoundment, annual water-level fluctuation, and rainfall precipitation may also play important roles in reactivating the landslide [48]. Numerical investigations have proven that the safety factor of the landslide obviously decreases during the reservoir water level drawdown and intense rainfall infiltration [11, 38]. The systematic monitoring data including Earth surface deformation, underground deformation, groundwater table and rainfall also verified these conclusions [28, 42, 44]. A close inspection of the existing literature indicate that there is no research on the potential sliding process and its future consequences. Moreover, the effects of landslide-induced tsunami on the riverside communities and infrastructure, and the vessels sailing along the Yangtze River waterway have to be discussed yet. The insufficiency in the research of the aforementioned aspects brings about the necessity of performing numerical simulations to investigate the impacts of induced tsunami.

3 Numerical modelling approaches

3.1 Modelling concept

Typically, continuum, discrete and hybrid methods are available to simulate landslides [5]. Discrete element method (DEM) treat the material as separate blocks or particles, which therefore enables the displacements and rotations of discrete elements. Although the DEM can easily deal with the large displacement of landslides, it cannot simulate the deformation of elements during sliding and cracking phenomena once the landslide has initiated. Combining the advantage of both continuum and discrete methods, the finite-discrete element method (FDEM) is used to simulate elastic deformation, fracturing, translation, rotation and velocity of the Huangtupo landslide.

Fig. 4

Conceptual modelling approach including the two numerical methods (FDEM and SPH) and the main governing equations

Furthermore, for free-surface flows simulations and induced-tsunami analysis, the SPH method generally presents some remarkable advantages over mesh-based methods: (1) no special treatment to detect the free surface; (2) straightforward modelling of moving complex boundaries and interfaces; (3) no need for special variables to detect different phases in the space since each individual particle possesses the material properties of its phase; (4) natural incorporation of coefficient discontinuities and singular forces into the numerical scheme [2]. The meshless nature of an SPH-based model allows for capture of the violent hydrodynamics of the tsunami waves acting on river bank structures and ships on the river. Hence, the SPH-based model is adopted to study the landslide-induced tsunami and its impacts on ships.

The conceptual modelling approach of the present study is shown in Fig. 4. Generally, the FDEM modelling is performed to capture the kinematic features of the landslide during failure. After failure, the intrusion of the landslide into the river, and the subsequent induced tsunamis are simulated by 2D and 3D SPH models. The kinematic information obtained from the FDEM model is used as the initial condition to simulate the landslide motion in the SPH modelling. Since the focus of this study lies on the latter aspect, the FDEM simulation is simplified by a 2D model. Section M–M’ (see Fig. 2), crossing the central No. 1 sliding mass, was chosen as the representative profile of the FDEM modelling. It is worth noting that the sliding mass is treated as a rigid body in the SHP modelling. Thus, the motion of the landslide can be determined by the time serials movement data of the 2D cross section in the FDEM simulation. Please note that this 2D section represents the worst scenario of landslide failure. It can provide a reasonable estimation of the kinematics for the intrusion, so that a conservative prediction of the induced tsunami is achieved.

3.2 FDEM with a nonlinear fracture model

In a 2D FDEM simulation, the intact body is discretized with a mesh comprising three-noded triangular elements. Four-noded crack elements are embedded between the edges of all adjacent triangle pairs. An explicit time integration scheme is applied to solve the equations of motion for the discretized system and to update the nodal coordinates [25, 26, 35, 36].

3.2.1 Governing equation

The generalized governing equations of motion can be expressed as:

$$\begin{aligned} \varvec{M}\frac{\partial \varvec{x}}{\partial t^2}+\varvec{C}\frac{\partial \varvec{x}}{\partial t}+\varvec{F}_{\rm{int}}(\varvec{x})-\varvec{F}_{\rm{exp}}(\varvec{x})-\varvec{F}_c(\varvec{x})=0 \end{aligned}$$

(1)

where $\varvec{M}$ is the system mass matrix; $\varvec{C}$ is the damping matrix; $\varvec{x}$ is the vector of nodal displacement; $\varvec{F}_{\rm{int}}$, $\varvec{F}_{\rm{ext}}$, and $\varvec{F}_c$ are the vectors of internal resisting forces, applied external loads, and contact forces, respectively. The internal resisting forces $\varvec{F}_{\rm{int}}$ include the contributions from the elastic reaction forces and the inter-element bonding forces, which are used to simulate elastic deformation and fracture, respectively. The contact forces $\varvec{F}_c$ are calculated either between contacting discrete bodies or along internal discontinuities (i.e. pre-existing or newly created fractures).

Numerical damping is introduced into the model to account for energy dissipation due to nonlinear material behaviour or to model quasi-static phenomena by dynamic relaxation. The matrix $\varvec{C}$ is equal to:

$$\begin{aligned} \varvec{C} =\mu \varvec{I} \end{aligned}$$

(2)

where $\mu$ and $\varvec{I}$ are a constant viscous damping coefficient and the identity matrix, respectively. The theoretical critical damping coefficient can be obtained from [36].

3.2.2 Fracture model

The progressive failure of materials is simulated using a cohesive-zone approach [35]. During elastic loading, the deformation of the material is solved by the model based on constant strain triangular finite elements by applying linear elasticity continuum theory. On the other hand, the initiation and propagation of fractures are simulated by using the concepts of nonlinear fracture mechanics. This approach assumes that the plastic strains are localized within a narrow fracture process zone when the stress exceeds the peak strength of the material. The mechanical response of the fracture process zone is captured by a nonlinear interdependence between stress and crack displacement implemented at the crack element level [25, 26].

The softening response of a crack element after fracturing is defined in terms of a variation in the bonding stresses between the edges of the triangular element pair as a function of the crack relative displacements in the normal and tangential directions, respectively. In tension, the response of each crack element depends on the cohesive tensile strength, $\sigma _t$, and the opening mode energy. In shear, the behaviour is governed by the peak shear strength, $\sigma _s$, and the sliding mode energy [25]. The peak shear strength is defined as:

$$\begin{aligned} \sigma _s=c + \sigma _n \text {tan}\phi \end{aligned}$$

(3)

where c is the cohesion, $\phi$ is the internal friction angle, and $\sigma _n$ is the normal stress acting across the crack elements. Material interfaces, such as the sliding surfaces between sliding mass and bedrock (Fig. 4), are the typical preexisting fractures, on which frictional forces are applied.

3.3 SPH with rigid body dynamics

SPH is a fully Lagrangian meshless method that discretizes the continuum material into a set of particles. For the simulation of fluid dynamics, the set of neighbouring particles is defined by a function with an associated characteristic distance, called smoothing length (h). The contributions of neighbouring particles for the computation of position, velocity, density and fluid pressure are dependent on the distance between particles and the corresponding parameters are calculated by using a weighted kernel function (W) in the area defined by smoothing length. Please note that particles are initially created with an interparticle distance, $d_p$, which is used as a reference value to define the smoothing length using $h =2d_p$.

3.3.1 Governing equations in SPH

In the classical SPH formulation, the Navier–Stokes equations can be written in a discrete SPH formalism, as follows:

$$\begin{aligned} \frac{\rm{d}\varvec{r}_i}{\rm{d}t}= & {} \varvec{v}_i \end{aligned}$$

(4)

$$\begin{aligned} \frac{\rm{d}\varvec{v}_i}{\rm{d}t}= & {} -\sum _{j=1}^{n}m_j\left ( \frac{p_j}{\rho _j^2}+ \frac{p_i}{\rho _i^2} + \Pi _{ij}\right )\nabla _i W_{ij}+\varvec{g}\end{aligned}$$

(5)

$$\begin{aligned} \frac{\rm{d}\rho _i}{\rm{d}t}= & {} \sum _{j=1}^{n} m_j\varvec{v}_{ij}\nabla _i W_{ij}\nonumber \\&+\,2\delta h c_0 \sum _{j=1}^n m_j \left( \frac{\rho _i}{\rho _j}-1\right) \frac{\varvec{r}_{ij}\cdot \nabla _i W_{ij}}{\vert \varvec{r}_{ij}\vert ^2+\eta ^2} \end{aligned}$$

(6)

where $d(\cdot )/\rm{d}t$ represents the material time derivative, $\nabla$ is the gradient operator. $\varvec{r}$ is the position, $\varvec{v}$ is the velocity, p is pressure, and $\varvec{g}$ denotes acceleration induced by body force, e.g. gravity; $\rho$ is the density, m is the mass, $c_0$ is the numerical speed of sound.

The kernel function $W_{ij}$ depends on the normalized distance between particle i and neighbouring particle j, i.e., $\varvec{r}_{ij} = \vert \varvec{r}_{i} - \varvec{r}_{j} \vert$. The Quintic kernel proposed by Wendland [51] is adopted for the present study. This weighting function vanishes for interparticle distances greater than 2h.

The artificial viscosity ($\Pi _{ij}$) proposed by Monaghan [32] is used here.

$$\begin{aligned} \Pi _{ij} ={\left\{ \begin{array}{ll} \alpha {\bar{c}}_{ij}\mu _{ij}/{\bar{\rho }}_{ij} &{} \varvec{v}_{ij}\cdot \varvec{r}_{ij<0}\\ 0 &{} \varvec{v}_{ij}\cdot \varvec{r}_{ij>0} \end{array}\right. } \end{aligned}$$

(7)

where $\alpha$ is a parameter set as 0.01, which has proven to yield the best results in the validation of wave flumes to study wave propagation and wave loadings exerted on river bank structures [1]; ${\bar{c}}_{ij}$ and ${\bar{\rho }}_{ij}$ are the mean speed of sound and mean density of the particles i and j, respectively; $\mu _{ij} = h\varvec{v}_{ij}\cdot \varvec{r}_{ij}/(\vert \varvec{r}_{ij}\vert ^2+\eta ^2)$ with $\eta ^2=0.01h^2$ adopted to avoid singularities in the fraction caused by very close particles.

The second term of the right-hand side of Eq. (6) is a diffusion formulation [31] designed to increase the accuracy of impact force computation. This so-called $\delta$-SPH approach is applied to smooth out the oscillations for better pressure results, with $\delta = 0.1$ recommended for most applications. The latest $\delta$-SPH formulation free from calibration parameters can be found in [30].

The fluid is considered as weakly compressible. This assumption allows us to compute the pressure directly from density using an equation of state as:

$$p = \frac{{c_{0}^{2} \rho _{0} }}{\gamma }\left[ {\left( {\frac{\rho }{{\rho _{0} }}} \right)^{\gamma } - 1} \right]$$

(8)

where $\gamma = 7$ for a fluid like water and $\rho _0 = 1000\, \hbox {kg/m}^3$ is the reference density of the fluid when there is no pressure. According to Eq. (8), the compressibility of the fluid depends on $c_0$, such that the fluid is virtually incompressible for a sufficiently high sound celerity [9].

3.3.2 Boundary conditions

The dynamic boundary condition (DBC) is used in this study. This method treats the boundary as particles that fulfill the same equations as the fluid particles. The boundary particles, however, do not move individually. Instead, they remain either fixed in position or are moved according to prescribed velocities. Once the distance between a fluid particle and a boundary particle decreases beyond the kernel range, the density of the boundary particles increases, leading to an increase in pressure. This leads to a repulsive force being exerted on the fluid particle, as calculated from the momentum equation (5). Hence, these fixed and moving DBCs are respectively utilized to simulate the stable Earth surface and the moving landslide. In this way, the moving path and velocities determined by landslide simulations with the FDEM are therefore assigned to the DBC to simulate the moving sliding mass.

3.3.3 Motion of floating structures

One of the specific capabilities of SPH models is to simulate the motion of floating objects, e.g. ships, represented by a set of boundary particles that are equi-spaced around the boundary [33]. This is done by considering their interaction with fluid particles and using the sum of fluid pressure for the entire floating body to derive its movement [12]. The floating objects are assumed rigid, and the net force on each object particle is computed according to the sum of the contributions of all surrounding fluid particles according to the designated kernel function and smoothing length. The force per unit mass $\varvec{f}_k$ on the boundary particle k is given by:

$$\begin{aligned} \varvec{f}_k= \frac{\rm{d} \varvec{v}_k}{\rm{d}t}=\sum _{i=1}^m \varvec{f}_{ki} \end{aligned}$$

(9)

where $\varvec{f}_{ki}$ denotes the force per unit mass on boundary particle k exerted by the fluid particle i on the boundary particle k and the summation is over the surrounding fluid particles (particle number m).

For the motion of the floating object, Newton’s equations for rigid body dynamics in the domain reference frame are used, and the discretization consists of summing the contributions from each SPH node, as follows:

$$\begin{aligned} M_J\frac{\rm{d} \varvec{V}_J}{\rm{d}t}= & {} \sum _{k\in J}m_k\varvec{f}_k, \end{aligned}$$

(10)

$$\begin{aligned} I_J\frac{\rm{d}\varvec{\Omega }_J}{\rm{d}t}= & {} \sum _{k\in J}m_k(\varvec{r}_{k}-\varvec{R}_J) \times \varvec{f}_k \end{aligned}$$

(11)

where rigid body J possess a mass $M_J$, velocity $I_J$, moment of inertia $\varvec{V}_J$, angular velocity $\varvec{\Omega }_J$, and center of gravity $R_J$. Equations (10) and (11) are integrated in time to predict $\varvec{\Omega }_J$ and $\varvec{V}_J$ for the beginning of the next time step. Each particle within the body J then has a velocity $\varvec{v}_k$ through:

$$\begin{aligned} \varvec{v}_k=\varvec{V}_J+\varvec{\Omega }_J \times (\varvec{r}_k-\varvec{R}_J) \end{aligned}$$

(12)

Finally, the floating object particles within the rigid body are moved by integrating Eq. (12) in time. This technique conserves both linear and angular momentum [33], and is successfully validated by other numerical methods and experimental measurements [7, 8].

4 Numerical simulation

In this section, the failure initiation of the landslide and the subsequently induced-tsunamis are simulated. In particular, the kinematics of the landslide are simulated with FDEM using Irazu2D (Version 4.00) commercial software developed by Geomechanica Inc. Furthermore, the analysis of the resulting impulse waves and behaviour of floating objects is carried out by using DualSPHysics (Version 4.4.12), an open-source code developed to study free-surface flows based on SPH [12].

4.1 Simulation of failure initiation

A 2D numerical simulation model is set up according to the representative cross section M–M’ (Fig. 5). The model ranges 2500 m along the sliding direction, with the elevation from 0 to 600 m above sea level. By simulating this section, it is possible not only to study the deformation and failure characteristics of the most dangerous area of Huangtupo landslide but also to determine the behaviour of the No. 4 landslide, which is located behind and overlaps with the No. 1-1 sliding mass. Considering that the failure of the No. 1-1 sliding mass could result in the loss of support at the toe area of the No. 4 landslide, this may increase the risk of a subsequent slide.

Fig. 5

2D finite-discrete element model of the Huangtupo landslide

According to the lithology and geological structure of the Huangtupo landslide, the numerical model is subdivided into the following geotechnical materials: (1) debris, (2) disintegrated rock, (3) disturbed rock, (4) fractured rock, and (5) bedrock. In addition, the disintegrated and disturbed rocks are distinguished between, above and below the water level. The contact surface between the sliding body and bedrock is the main sliding surface. The input parameters for the numerical model are presented in Table 1. The initial shear strength parameters (c and $\phi$) of the main sliding surface material are obtained from in situ direct shear tests during tunnelling [50]. The residual friction angle of the sliding surface is only $14 \pm 1$ degree, as determined by averaging the results of the in situ reversal shear and laboratory ring shear tests. The properties of the bedrock are tested in the laboratory by triaxial compression ($c, \phi , E$ and $\nu$) and uniaxial tensile tests ($\sigma _t$). The properties of the sliding masses, including the debris, disintegrated, disturbed and fractured rocks, are empirically estimated by reducing the values of the original intact rock data according to their weathering and the state of disintegration.

Table 1 Material parameters for the FEDM simulation

The distribution of the sliding surface of the landslide is clearly exposed based on previous geological survey. In the FDEM simulation, the largest deformation is expected to occur along the sliding surface, where the mesh size is relatively finer. The sliding mass is meshed with 10 m, which is sufficient to reveal the deformation phenomena such as geometry and cracks illustrated in Fig. 7. Since the bedrock remains stable, the maximum mesh size at the boundaries is 100 m (Fig. 5). Ultimately, the model is meshed into 5459 elements with 2857 nodes. The bottom boundary of the model is assigned as pin condition (i.e. zero displacement in both the x- and y-directions), while the lateral boundaries are assigned as vertical roller (i.e. zero displacement in the x-direction).

According to recent studies [10, 20], mechanical properties of the rocks and soils close to the reservoir banks gradually degrade due to the periodic water level fluctuations, resulting in drying and wetting circles, as well as continuous creep. To study the evolution of the Huangtupo landslide under such conditions, the mechanical strength parameters such as cohesion and friction angle were systematically reduced from 95 to 80% of the initial values. Note that the failure criterion used in FEM does not apply here, because FDEM simulation can always continue at any reduction factor. We proceed to define the initiation of landslide failure by examining the maximum displacement during a specified time duration. Consequently, at this time the shear strength parameters of the sliding surface are set to residual strength to simulate the deformation and motion characteristics of the landslide after the failure.

4.2 Simulation of induced-tsunamis and impacts

In the current study, a digital elevation model with a resolution of 30 m is used to build the 3D SPH model. The model covers a water area of 3 km around the Huangtupo area. Figure 6 shows the SPH model with the potential sliding mass, i.e. the No. 1-1 sliding mass. In addition, two reservoir conditions at low and high water levels (145 m and 175 m above sea level) are considered to reflect normal operation of the Three Gorges Reservoir.

Fig. 6

3D model for impulsive wave and floating objects simulation showing a the virtual reality scene, and b the solid 3D model

The floating objects on the water surface in the model are represented by rectangular solids with length, width and height of 60 m, 20 m and 20 m, respectively, referring to a typical size of the ships that navigate on the Yangtze River (Fig. 6). A regular spacing of 180 m is set between the 4 floating objects in front of the landslide, whereas all others are set as 500 m. To restrict the reservoir, water cannot spill from the model during the simulation, and hence virtual dams with a top elevation of 175 m are set as model boundaries at the ends of the river on the east, west and north sides. The specific geometries of the Earth surface, the simulated landslide, the water level and ships are imported into DualSPHysics using CAD files. The simulated objects can be generally divided into two types according to their properties: fluid or rigid bodies. In this model, water is set as fluids, while other objects are set as rigid bodies.

The main variables that must be defined in DualSPHysics cases are interparticle distance, liquid viscosity, liquid density, density of float objects and object motion. The interparticle distance ($d_p$) is the initial distance between each particle generated in the domain of objects. A smaller value of $d_p$ increases the resolution of the model and shows more details but implies a much higher computational cost. The ratio of wave height to interparticle distance ($H/d_p$) is used as the factor to evaluate the resolution of the model. Convergence studies can be used to analyse the influence of the model resolution (different $H/d_p$ values) on the simulation results. According to the works of Altomare et al. and Roselli et al. [2, 37], a relatively coarse resolution ($H/d_p=3.0$) causes approximately 10% error in comparison with the theoretical value of wave height calculating results. It is suggested that $H/d_p$ should be larger than 10 to achieve a more accurate modelling, which remains an error less than 3% comparing with the theoretical value [2, 37]. Additionally, an integer value of $d/d_p$ (ratio of water level to $d_p$) is recommended to ensure that the top layer of particles exactly matches the initial still water surface [2]. Because the water levels of the reservoir in our simulated scenarios are 175 m and 145 m, the recommended $d_p$ is 1 m or 5 m; considering the limitation of computing capacity, here we used $d_p = 5$ m for the 3D model and created 5,207,283 and 6,783,774 particles for the two reservoir levels, respectively. It is worth noting that the maximum wave height simulated by our model is approximately 20 m (measured from wave crest to trough), resulting in $H/d_p=4.0$. Hence, the error of the 3D simulation should be restricted within 10%.

Generally, floating objects are sensitive to the model shapes. To more accurately simulate the roll motions of the floating objects, a 2D SPH model describing more detailed shapes of the ships ($d_p = 1\,\hbox {m}$, $H/d_p = 20$) is established based on the typical section that crosses from ship N1 to ship N4 (Fig. 10a). During the actual navigation, ships N1 to N4 are passing perpendicular to the sliding direction of the landslide and to the propagation direction of the induced impulsive wave, implying that they will be more vulnerable to significant rolling motions (Fig. 10). The 2D model profile is consistent with the model shown in Fig. 5, while the positions, densities and overall sizes of the ships N1 to N4 are consistent with the same floating objects in the 3D model in Fig. 6.

The density of the fluid is set to $1000 \,\hbox {kg/m}^3$. The bounded objects of a simulation can be still, floating or in motion. The modelled reservoir bank and rock mass are fixed on the initial location during the simulation process. The ship models are configured in a floating state, and the density of these floating objects is set to be 25% that of water, i.e., $250 \,\hbox {kg/m}^3$. Thus, the flotation depth of the ships is 5 m, which is a common value for most ships on the Yangtze River.

The motion of the sliding mass is modelled using the function “set a motion from a file”, which contains a time series of coordinates describing the dynamics of the object. The input for the velocities (movements) is obtained from the FDEM simulation. As illustrated in Fig. 7d, the serial time movement data of the sliding object are decomposed into 3 components along the x- (east), y- (north), z- (elevation) directions.

5 Results and discussion

5.1 Landslide simulation

Figure 7 shows the results of the FDEM simulations. The final maximum displacements dependent on the reduced strength parameters of the sliding surface are illustrated in Fig. 7b. When the strength parameters decrease from 100 to 90%, the landslide deformations are less than 2.6 m. During this stage, the landslide becomes unstable, with continues creep. This movement pattern agrees with the current monitoring data (Fig. 3). With a further decrease in the strength from 82 to 81%, however, the final displacement of the landslide increases dramatically to approximately 40 m, which can be treated as a catastrophic failure. After failure, the shear strength of the sliding surface of No. 1-1 sliding mass is set to residual strength. The final deposition of the landslide with a runout distance approximately 230 m is shown in Fig. 7a.

Fig. 7

Landslide simulations: a final displacement and geometry of the landslide after failure, b the maximum displacement under different strength parameters, c sliding velocity over time for the No. 1-1 sliding mass after failure and d the time serials movement data of the sliding mass

In contrast, the No. 4 landslide does not fail, however, the stress state is redistributed resulting in compressive shear cracks distributed at the toe of this landslide followed by vertical shallow tensile cracks (Fig. 7a). The latter might act as preferential rainfall infiltration passages and therefore could destabilize this landslide. Two representative points, A and A’ shown in Fig. 7a, are considered to show the velocities of the failed landslide. The velocities gradually increase over time, giving rise to a maximum velocity of approximately 8 m/s at 27 s. Then, the sliding slows and halts at 51 s.

5.2 Tsunami simulation

Figure 8 shows the propagation and wave height data over time for both reservoir water levels at 145 m and 175 m above sea level. The resulting water levels due to the landslide-induced tsunami are also shown. The black dots with labels at different locations on the reservoir bank show the maximum wave run-up elevations above sea level.

Fig. 8

Wave propagation and run-up simulation results in the 145 m and 175 m reservoir water level scenarios due to the landslide-induced tsunami. (The red contours show the wave height, while black dots and numbers indicate the wave run-up at different locations on the river banks. All values are expressed as the elevations above sea level)

The maximum wave height occurs at the adjacent position between the landslide and reservoir. The maximum wave run-up is approximately 12 m indicating the possible magnitude of the tsunami. Only after approximately 40 s does the tsunami reach the river banks of the Leijiaping area on the opposite side of the Huangtupo landslide. Here, the maximum wave run-up still remains approximately 9 m. At 80 s, the impulsive waves propagate to 1.5 km away from the landslide with a velocity of approximately 19 m/s. The wave height and wave run-up in the north side tributary of Dongrang River are obviously higher than those in the Yangtze River. In particular, the wave height and run-up increase to 3 m in the turning and narrowing area of the river. If consider the interaction between the riverbank slope and tsunami, much higher local wave run-ups could be achieved.

Apart from the landslide adjacency, the highest wave run-ups are found in the Leijiaping area. Under the high reservoir condition (i.e. 175 m above sea level), the wave run-up may not seriously threaten the buildings and roads with elevations of above 200 m above sea level. However, areas where local residents live and work close to the reservoir are still under risk, such as the docks on the water surface as well as fruit trees and vegetable fields, located between 175 and 185 m above sea level. For the other parts of the reservoir bank, the wave run-up might reach locally to higher values due to the influences of local terrain turning and narrowing. According to our worst-case simulations, humans and properties within 6 m above the water levels of the Yangtze River and 8 m above the water levels of its tributaries must pay special attention to the landslide-induced tsunamis.

5.3 Impacts on shipping

The anti-wave capacities of inland river vessels are usually designed to be lower than those for marine vessels [39]. However, in narrow river channels such as those in the Three Gorge Reservoir, landslide-induced tsunamis might exert larger impacts on the vessels. In general, the most common risks of vessels facing in water waves include collision, sinking and capsizing caused by different motions. Typically, the motion behaviours of floating objects are described by translational (surge, sway and heave) and rotational motions (pitch, roll and yaw). Sway and surge are the linear transverse (side-to-side motion) and the linear longitudinal (front/back) motions. Heave is the linear vertical (up/down) motion [22]. As the surge and sway are both horizontal movements that travel in two orthogonal directions, here we treat them as a combination value by surge and sway, indicating the total horizontal displacement. Pitch, roll and yaw are the up/down rotation about the transverse (side-to-side), tilting rotation about the longitudinal (front-back), and the turning rotation about the vertical axis (z-axis), respectively. Here, only the roll is analysed,as it is the largest threat to the safety of a ship, considering that the roll can overturn a ship

Fig. 9

Simulation results of the surge and sway and heave of floating objects. (The red and orange curves show the surge and sway and heave values of the floating objects on different locations of the river in the 145 m and 175 m water level scenarios, respectively)

Generally, the potential damages caused by surge and sway are the following: the collision of ships with nearby banks, submerged rocks, fixed structures such as bridges and docks or other anchored vessels. The harm caused by heave is mainly the water inflow into cabins, especially for fully loaded ships, which have large draft depths and whereby the boards are very close to the water’s surface. Massive water inflows can destroy the cargo, and might also cause the ships to sink.

With a low reservoir water level (145 m above sea level), the maximum surge and sway of 110 m is observed for ship N1 closest to the landslide (Figs. 6b and 9a). With the wave propagation, the values on the east and west sides gradually decrease to approximately 24 m and 34 m, respectively. Notably, the north side shows different phenomena: the surge and sway values of ships gradually increase after entering the tributary. The horizontal movements of ship N5 close to the northern model boundary also increase to 110 m. Similar phenomena are observed with high reservoir water levels (175 m above sea level); however, the simulated surge and sway at the ship N1 (approximately 61 m) is much lower. Therefore, the ships navigating in the main stream within 500 m surrounding the landslide and in the tributary have higher risk of collision with the nearby river banks, submerged rocks, bridges, docks or other vessels because of larger surge and sway values when exposed to the landslide-induced tsunami.

Depending on the reservoir water level, the maximum heave of the ships may reach between 7 and 8 m (Fig. 9b). Towards the east and west, the heaves gradually reduce to approximately 2 m. In general, the surge and sway and the heaves of floating objects decrease away from the wave source, and exhibit more rapidly decreasing velocity in the close vicinity of the landslide ($\leqslant 500 \,\hbox {m}$). However, the water depths of the river may also influence the heave of the floating objects, which can be observed in the positions of ships N4 and N5. Here, the relatively shallow water depths result in local increases in heaves. Hence, close to the landslide ($\leqslant 500 \,\hbox {m}$), there are greater risk of water inflow and the sinking of ships due to the larger heaves. This information is important with respect to installation of an early warning system for passing ships, which for example, could be achieved by sirens installed along the river bank.

Fig. 10

a The 2D model for roll simulation and b simulation results of rolls for floating objects

Figure 10b shows the time series of the roll data of each floating object during both high and low reservoir water levels. The rolls of ships N1 to N3 are less than $3^\circ$. A maximum roll of $6^\circ$ is observed for ship N4, which is close to the opposite river bank approximately 1 minute after the landslide event, therefore providing sufficient time for an early warning system. This is caused by the shallow water depth at this location. Generally, a roll of $6^\circ$ is less than the anti-capsize capacity of a normal ship. Hence, in the case of a sudden landslide failure, the possibility of ships capsizing near the landslide is small. However, the simulated roll of $3^\circ$ to $6^\circ$ may still cause the movement of cargo and the falling of passengers.

Although our predictions cannot be validated, the results provide an initial estimation, as the deployed modelling methods have already demonstrated their reliability already in previous studies [16, 21, 27, 29, 55]. Furthermore, the simulation results of the sliding velocity of the landslide and the wave propagation parameters are within similar ranges as actual recorded data of previous cases in nearby areas [18, 45], as well as the inverse estimation based on tsunami sediments [20].

6 Conclusions

In this study, we predict and evaluate the possible sliding and risk of the impulsive wave (tsunami) induced by the Huangtupo landslide at the Three Gorges Reservoir. 2D numerical simulationy with strength reductions are performed to predict possible deformations and sliding velocities. Furthermore, an SPH approach is deployed to simulate the landslide-induced tsunami and its impact on shipping on the Yangtze River. From our simulations, the following conclusions can be drawn:

• 1. From the entire Huangtupo landslide, the No. 1-1 riverside sliding mass exhibits the largest ongoing deformations and therefore imposes the highest risk of triggering a landslide-induced tsunami. If the shear strength of the sliding zones is reduced to about 80% of initial values, immediate failure could occur, resulting in maximum displacement of approximately 230 m.

• 2. The velocities, heights and run-ups of the tsunami are similar for both studied typical reservoir conditions 145 m and 175 m above sea level. The speed of the tsunami is approximately 19 m/s, and the run-up on the nearby river banks ranges between 3 and 9 m depending on the distance and the local topography.

• 3. For the floating objects, the surge and sway range from 18 to 110 m, and the heaves range from 2 to 8 m. The ships within 30–60 m of reservoir banks or any other water structures in the main stream and within 50–100 m in the Dongrang River face the risk of collision, while the ships within 500 m of the landslide have the risk of massive water inflow and even sinking. Finally, the roll of the floating objects close to the landslide is generally low, at less than 6 degrees, which cannot cause the rollover or the capsizing of ships. However, the roll may still present risks of loss and injury due to the movement of cargo and the falling of passengers and crew.

Hence, based on the outcome of our study, we strongly recommend continuous monitoring of hazardous landslides in the Three Gorges Reservoir Area (TGRA), the implementations of early warning systems for the landslides and the passing ships, along with the establishment of emergency plans at the studied and other possible hazardous sites.