When a reservoir bank landslide body into the water, the generated waves will endanger the safety of passing ships and surrounding buildings, as well as the safety of the dam. Therefore, calculating the waves height is important for evaluating the risks of landslides1.
The methods of calculating the landslide generated wave height can mainly be divided into analytical method2–6, numerical simulation method7–8, and physical modelling method9–10. Analytical method is widely used in engineering applications because of its simple modelling processes, which has few requirements for engineers and high precision.
The analytical method originated from Noda2. Noda proposed an approximate method to find the amplitude of the largest wave in the nonlinear region by utilizing the solutions obtained from linear theory. Since then, many scholars have conducted more in-depth research. For example, Academician Pan Jiazheng divided the landslide body into many two-dimensional (2D) vertical strips and calculated the wave height by considering the horizontal and vertical movement of the landslide. This method is called the Pan Jiazheng method3. Huang et al. improved the Pan Jiazheng method by considering the resistance of water and the change in the friction coefficient4. Miao et al. proposed a sliding block model based on the 2D vertical strip method to predict the maximum wave height5. The American Civil Engineering Society recommends a prediction method of the wave height that assumes the landslide results in the particle motion with a centre of gravity, and Newton's law of motion is used to calculate the wave height6.
The above methods are all 2D analysis methods. In the vertical strip method, the calculation results will differ with the selection of the 2D section. The 2D analysis methods cannot effectively simulate the actual spatial state of three-dimensional (3D) landslide. Hu11 proposed that the value obtained by 2D analysis method is approximately 70% of the value based on 3D analysis method. To date, analytical method based on the 3D landslide body model has not been studied by scholars.
Geographic information systems (GIS) is widely used in geotechnical engineering. The most notable feature of GIS is that they can transform vector data into grid data sets based on a grid column unit model12. Because of the high 3D spatial data processing capability of GIS, many scholars have added geotechnical professional models to their respective systems. For example, Professor Xie's team established a 3D limit equilibrium method based on GIS, and developed a slope stability analysis module called 3Dslope13–15. Jia et al. proposed a slope stability analysis method by coupling a rainfall infiltration model and 3D limit equilibrium method within the GIS environment16. Mergili et al combined GRASS GIS and the 3D Hovland model to implement a 3D slope stability model capable of considering shallow and deep-seated slope failures17. Therefore, to develop a wave height calculation module in GIS, it is necessary to first establish a force analysis model of the 3D landslide in GIS.
Based on the spatial data processing capability of GIS, this paper applies the grid column unit model to establish a 3D landslide model, and proposes a method for calculating the wave height. Compared with 2D analysis methods, the 3D method proposed in this paper better represents the actual spatial state of landslides. Simultaneously, the resistance of the water is considered to improve the accuracy of the calculation result. To make the calculation more convenient, an expansion module is developed to calculate the wave height in GIS, and the feasibility of the module is verified by a case study.
Gis-based Method Of Calculating The Wave Height
Grid column unit model based on GIS. For a slope, the representation of data is mainly in the form of vectors. These data include but are not limited to slip surface, strata, groundwater, fault, slip, and other types of data. These vector data layers can be converted to raster data layers using the spatial analysis capabilities of GIS to form a grid data set. The grid data structure consists of rectangular units. Each rectangular unit has a corresponding row and column number and is assigned an attribute value that represents the grid unit18. Therefore, the slope can be divided into square columns based on the grid units to form a grid column unit model, as shown in Fig. 1.
The spatial relationships between parameters. Figure 2 shows the spatial relationship between parameters. θ is the dip of the grid column at the slip surface; α is the dip direction of the grid column at the slip surface; β is the sliding direction of the landslide; θr is the apparent dip of the main inclination direction of the landslide; αx is the apparent dip of the X-axis; and αy is the apparent dip of the Y-axis. Of these six parameters, α, θ, and β is known and is calculated in a paper published by Xie 's team14. Other three parameters αx, αy and θr can be calculated according to the spatial relationship in Fig. 2(b), which are calculated as follows.
$$\alpha x=\arctan (\cos \alpha \tan \theta ),\alpha y=\arctan (\sin \alpha \tan \theta )$$
1
$${\theta _r}=\arctan \left\{ {\tan \theta } \right.\left. {\left| {\cos \left( {\alpha - \beta } \right)} \right|} \right\}$$
2
The bottom area of one grid column is calculated by
$$A=cellsiz{e^2}\left[ {\frac{{\sqrt {\left( {1 - {{\sin }^2}\alpha x{{\sin }^2}\alpha y} \right)} }}{{\cos \alpha x\cos \alpha y}}} \right]$$
3
where cellsize represents the length of the side of each grid column.
The weight W of the grid column is expressed as
$$W=cellsiz{e^2}\sum\limits_{{m=1}}^{n} {{h_m}{r_m}}$$
4
where m is the number of strata, hm is the height of each stratum, and rm is the unit weight of each stratum. For the grid column units above the water, rm is calculated from the natural unit weight. For grid column units under water, rm is calculated from the buoyant unit weight.
The pore water pressure is obtained as follows19.
u = R/cosθ (5)
where R is the distance from the centre bottom of the grid column to the water surface.
When the sliding body enters the water, the resistance of the water is calculated as follows20.
$$G{\text{=}}\frac{1}{2}{c_w}{\rho _f}{v^2}S$$
6
where G is the resultant force of the resistance of the water to the sliding body (mN); cw is the viscous resistance coefficient, which is 0.15 to 0.18; ρf is the buoyant density (g/m3), taking the average of all stratum; v is the velocity of the landslide (m/s); and S is the surface area of the grid column in the water (m2).
Dynamic equation based on grid column units. As shown in Fig. 3, it is the force analysis of one grid column. In the 3D sliding body, one random grid column ABCDA1B1C1D1 is selected, and the force analysis is explained in the supplementary file.
We assume that all the grid column units move continuously, do not separate in the macroscopic dimension and remain vertical after sliding, as also assumed by Pan Jiazheng3. The force analysis of one grid column and the spatial relationships between parameters at the slip surface are shown in Fig. 2 and Fig. 3, respectively.
We arbitrarily selected a grid column unit (the grid column unit in row i and column j). According to Newton’s laws of motion, dynamic equilibrium equations are established in the sliding direction of the landslide and the vertical direction. The force analyses in the sliding direction of the landslide and vertical direction are shown in Fig. 4.
$${A_{i,j}}{\tau _{i,j}}\cos \theta {r_{i,j}} - {A_{i,j}}{\sigma _{i,j}}\sin {\theta _{i,j}}\cos \left( {{\alpha _{i,j}} - \beta } \right) - k{W_{i,j}}+\Delta {D_{i,j}} - {G_{i,j}}=\frac{{{W_{i,j}}}}{g}{a_x}$$
7
$${A_{i,j}}{\tau _{i,j}}\sin \theta {r_{i,j}}+{A_{i,j}}{\sigma _{i,j}}\cos {\theta _{i,j}} - {W_{i,j}} - {P_{i,j}}+\Delta {V_{i,j}} - \Delta {R_{i,j}}=\frac{{{W_{i,j}}}}{g}a{y_{i,j}}$$
8
$${\tau _{i,j}}={c_{i,j}}+\left( {{\sigma _{i,j}} - {u_{i,j}}} \right)\tan {\varphi _{i,j}}$$
9
where ax and ayi, j are the horizontal acceleration and vertical acceleration of the grid column, respectively; φi, j is the effective friction angle of the grid column at the slip surface; g is gravitational acceleration; ci, j is the effective cohesion of the grid column at the slip surface.
According to this assumption, the horizontal acceleration ax of each grid column unit is the same, and the vertical acceleration ayi, j of each grid column unit varies. Pan Jiazheng3 suggested that there is a certain proportional relationship between ax and ayi, j, that is, ayi, j/ax=tanδi, j. δi, j is the horizontal inclination angle of the line connecting the centre bottom of the grid column to the centre bottom of the next grid column in the sliding direction of the landslide. The effect of vertical tangential forces is ignored, namely, ΔVi, j-ΔRi, j=0; therefore, Eq. (8) can be transformed as follows.
$${A_{i,j}}{\tau _{i,j}}\sin \theta {r_{i,j}}+{A_{i,j}}{\sigma _{i,j}}\cos {\theta _{i,j}} - {W_{i,j}} - {P_{i,j}}=\frac{{{W_{i,j}}}}{g}ax\tan {\delta _{i,j}}$$
10
The simultaneous Eqs. (9) and (10) can be obtained as follows.
$${\sigma _{i,j}}=\frac{{{A_{i,j}}\sin \theta {r_{i,j}}\left( {{u_{i,j}}\tan {\varphi _{i,j}} - {c_{i,j}}} \right)+{W_{i,j}}+{P_{i,j}}+\frac{{{W_{i,j}}}}{g}ax\tan {\delta _{i,j}}}}{{{A_{i,j}}\left( {\sin \theta {r_{i,j}}\tan {\varphi _{i,j}}+\cos {\theta _{i,j}}} \right)}}$$
11
For the entire sliding body, the forces between the grid columns are internal forces, that is, the resultant force is 0, yielding Eq. (12).
$$\sum\limits_{I} {\sum\limits_{J} {\Delta {D_{i,j}}=0} }$$
12
By summing all the grid column units, the horizontal acceleration ax can be determined by Eq. (7).
$${a_x}=\left[ {\sum\limits_{I} {\sum\limits_{J} {\frac{{{A_{i,j}}{\tau _{i,j}}\cos \theta {r_{i,j}} - {A_{i,j}}{\sigma _{i,j}}\sin {\theta _{i,j}}\cos \left( {{\alpha _{i,j}} - \beta } \right) - k{W_{i,j}}{\text{-}}{G_{i,j}}}}{{{W_{i,j}}}}} } } \right]g$$
13
Substituting Eqs. (9) and (11) into Eq. (13) yields the following equation.
$${a_x}=\left[ {\sum\limits_{I} {\sum\limits_{J} {\frac{{{B_{i,j}}+{E_{i,j}} - {F_{i,j}} - ({G_{i,j}}\tan {\delta _{i,j}}{H_{i,j}} - k{W_{i,j}}){L_{i,j}}}}{{{W_{i,j}}\tan {\delta _{i,j}}{H_{i,j}}{Q_{i,j}}}}} } } \right]g$$
14
where
$${B_{i,j}}={A_{i,j}}\cos \theta {r_{i,j}}\left( {{u_{i,j}}c\tan {\varphi _{i,j}} - {c_{i,j}}^{2}} \right)$$
15
$${E_{i,j}}={A_{i,j}}\cos \left( {{\alpha _{i,j}} - \beta } \right)\sin \theta {r_{i,j}}\sin {\theta _{i,j}}\left( {{u_{i,j}}\tan {\varphi _{i,j}} - {c_{i,j}}} \right)$$
16
$${F_{i,j}}=\left[ {\cos \theta {r_{i,j}}\tan {\varphi _{i,j}} - \sin {\theta _{i,j}}\cos \left( {{\alpha _{i,j}} - \beta } \right)} \right]\left( {{W_{i,j}}+{P_{i,j}}} \right)$$
17
$${L_{i,j}}={c_{i,j}}+\sin \theta {r_{i,j}}\tan {\varphi _{i,j}}$$
18
$${H_{i,j}}=\cos \theta {r_{i,j}}\tan {\varphi _{i,j}} - \sin {\theta _{i,j}}\cos \left( {{\alpha _{i,j}} - \beta } \right)$$
19
$${Q_{i,j}}={c_{i,j}}+\sin \theta {r_{i,j}}\tan {\varphi _{i,j}}$$
20
Calculation of the sliding velocity. The steps in calculating the landslide sliding velocity are as follows.
(1) Using the spatial analysis capability of GIS, the landslide body is rasterized, and the size of the grid column unit (i, j) can be set to an arbitrary square. A partitioning line is drawn from the bottom to the top of the landslide every ΔL in the sliding direction of the landslide, and the resulting regions are numbered zone 1, zone 2, zone 3, ..., zone (n-1), zone n. Each partition includes numbers of grid column units, and the length of zone n is less than or equal to ΔL, as shown in Fig. 5. For a grid column unit that is not completely contained within a partition, if the area within the partition is greater than half of the total area, the unit is divided into that partition; otherwise, the unit is divided into the next partition.
(2) For each grid column unit, the parameters required in Eq. (14) are calculated.
(3) t0 is the starting point of when the landslide body begins to slide, and t0 = 0. When the landslide body moves distance ΔL sequentially in the sliding direction of the landslide, the corresponding time is recorded as t1, t2, t3…tn, and the corresponding velocity is expressed as vx1, vx2, vx3…vxn.
(4)The horizontal acceleration at t0 can be calculated by Eq. (14) and is denoted as ax0, and the velocity at time t0 is zero. After sliding distance ΔL is reached, the following equations can be obtained.
$${v_{x1}}=\sqrt {2{a_{x0}}\Delta L}$$
21
$${t_1}={t_0}+\sqrt {\frac{{2\Delta L}}{{{a_{x0}}}}}$$
22
(5) At t = t1, the landslide body has horizontally moved by a distance ΔL in the sliding direction of the landslide, zone 1 has slipped from the sliding surface. The horizontal acceleration ax1 at t1 is still calculated by Eq. (14). Unlike t0, the weight for zone (n-1) changes to the weight for zone n, and the weight for zone (n-2) becomes the weight for zone (n-1), and so on (At this time, there is no grid column for zone n). After ax1 is calculated, the following can be established.
$${v_{x2}}=\sqrt {2{a_{x1}}\Delta L+{v_{x1}}^{2}}$$
23
$${t_2}={t_1}+\frac{{{v_{x2}} - {v_{x1}}}}{{{a_{x1}}}}$$
24
(6) The calculation is continued in turn. When the obtained horizontal acceleration is negative, the maximum velocity can be obtained. Finally, ax and vx in the calculation process can be plotted as respective curves versus the sliding time.
Calculation of the wave height. The China Institute of Water Resources and Hydropower Research proposed an empirical formula for wave height calculation21. In the formula, the main factors that affect the wave height are the sliding velocity and volume of the landslide. The formula for calculating the maximum wave height is as follows.
$${\xi _{\hbox{max} }}=d\frac{{{v_m}^{{1.85}}}}{{2g}}{V^{0.5}}$$
25
where ξmax is the maximum wave height (m); d is the comprehensive influence coefficient, with an average value of 0.12; vm is the maximum sliding velocity (m/s); V is the volume of the landslide body in the water (m3); and g is gravitational acceleration, which equals 9.8 m/s2.
The formula for calculating the wave height at different distances from the landslide body is as follows.
$$\xi ={d_1}\frac{{{v_m}^{n}}}{{2g}}{V^{0.5}}$$
26
where ξ is the wave height at a point from the landslide body L metres (m); n is the calculation coefficient, which is 1.4; and d1 is the influence coefficient related to distance L, which is determined by the following formula.
$${d_1}=\left\{ \begin{gathered} 0.5{\text{ , }}\left( {L \leqslant 35} \right){\text{ }} \hfill \\ 6.1274{L^{ - 0.5945}}{\text{, }}\left( {L>35} \right) \hfill \\ \end{gathered} \right.$$
27
Program Implementation
Combined with the wave height calculation method, an expansion module was developed based on component object model (COM) technology in the ArcGIS environment. Figure 6 illustrates the computational process.