A Numerical Analysis of Landslide Movements Considering the Erosion and Deposition along the Flow Path

Landslides are one of the most frequent disasters which occur widespread in Indonesia. This disaster often causes damages and fatalities. One of the mitigations efforts to reduce potential loss is by predicting the area affected by landslide movement. This research developed a numerical model of landslide movement by incorporating the erosion and deposition laws along the flow path. This model improves the accuracy of the previous models which assume that landslide volume is constant without any consideration for the erosion and deposition. The governing equation of this newly developed model uses the Eulerian numerical approach based on the finite difference scheme. The erosion-deposition laws applied in this research are from Egashira et al. (2001), McDougall and Hungr (2005), and Blanc (2008). The simulation program applies Python programming language and examines an imaginary slope with ellipsoid-shaped of source area. The simulation result shows that the additional erosion-deposition formula can enlarge the volume and the affected area of landslide movement. It is clarified that the erosion rate is a determinant factor affecting the results of calculation.


INTRODUCTION
The damage potential of landslides is determined by the velocity and the affected area. Therefore, predictions of landslide velocity, travelling path, depositional area, and flow depth are important keys in landslide risk and hazard assessments. This prediction can be used to estimate the damage potential, to design protective measures, or to identify the possibility of secondary effects such as landslide-generated wave and landslide dam that triggers flood (McDougall and Hungr, 2004;Hungr, 2007). Miyamoto (2010) suggested a 2D numerical model simulate landslide movement in Unzen-Mayuyama, Japan. This model gave the description of the landslide velocity and the affected area with finite difference numerical scheme. Fathani, Legono and Alfath (2017) further developed a numerical model by adding the earthquake factor and rheology parameter from Coulomb and Voellmy into the simulation model. Both numerical models assume that landslide volume is constant. In reality, landslide volume can increase because of erosion and deposition process along the landslide path. The landslide occurring in Pasir Panjang Village, the calculation was the hyper-concentrated solid-liquid mixture (Egashira et al., 1989;Egashira, 1997;Egashira, Miyamoto and Ito, 1997) and the Mohr-Coulomb. The erosion law was added into the numerical model and then was further tested on an imaginary slope. The purpose of this research is to develop landslide movement simulation program by incorporating the erosion laws suggested by Egashira, Itoh and Takeuchi (2001), McDougall and Hungr (2005) and Blanc (2008). This research also discusses the effect of erosion and deposition on the volume and covered area of landslide movement.

Landslide Movement Governing Equations
Landslide movement governing equations suggested by Miyamoto (2010) and Fathani, Legono and Karnawati (2017) is based on the momentum conservation law. The mass movement equation is shown in Eq. (1) and the continuity equation is shown in Eq. (2). The continuity equation is applied when there is no addition or reduction of landslide volume caused by erosion and deposition along the flow path.  In which M is the flux vector;  is the coefficient of momentum; u is the depthaveraged velocity; u t is the transverse vector of u; gz is the gravitational acceleration; h is the thickness of landslide mass; H is the slope height; T is the shear stress acting on the sliding surface; and ρm is the average density of landslide.  is explained in Eq. (3) as follows: Whereas i, j are the vector units on the direction x, y on Cartesian coordinate. Average density (ρm) is described in detail by Miyamoto (2010) and Fathani, Legono and Karnawati (2017).
Shear stress (T) in the hyperconcentrated solidliquid mixture model in Eq. (1) is also described in detail by Miyamoto (2010) and Fathani, Legono and Alfath (2017): where Ts , Td, Tf are the shear stress on a solid phase, colliding particles and supported by the interstitial liquid phase, respectively; ru is the ratio of pore water pressure at a sliding surface to the total pressure above the sliding surface; cs is the concentration of the solid phase of volume in the flow; θ is the eroded slope degree angle; s  is the internal friction angle along the sliding surface; ρs is the solid mass density; ρl is the fluid mass density; d is the diameter of particle in the flow; e is the coefficient of restitution; kg and kf are the empirical constant, kg= 0.0828 and kf = 0.25.
Shear stress (T) for the Mohr-Coulomb model in Eq. (1) was described in detail by Fathani, Legono and Alfath (2017) as follows: ; c is the cohesion of soil along the sliding surface.

Condition of Landslide Movement
Eq. (5) to Eq. (9) can only be used if the landslide moves because the friction on dense phase, which is a part of the shear stress (T) should be in balance with the external force. In the condition where the landslide stops, the shear stress (T) described by Miyamoto (2010) and Fathani, Legono and Alfath (2017) turns into: Shear stress (T) in the hyper-concentrated solidliquid mixture model is equal to the solid friction stress and less than the value on the right side of Eq. (5) as follows (Miyamoto, 2010): As for the Mohr-Coulomb model, shear stress (T) is equal to the solid friction stress and less than the value on the right side of the Eq. (8) as follows (Fathani, Legono and Alfath, 2017): (1 ) cos tan Furthermore, Miyamoto (2010) also described the application of the finite difference method to model the landslide mass when stopping at Dt.

Erosion and Deposition
Depth integrated models usually apply a simple erosion law i.e., the Hungr's erosion law (Hungr and Evans, 2004), the modified erosion law (Egashira, 1993;Egashira, Itoh and Takeuchi, 2001), path-controlled erosion (Chen, Crosta and Lee, 2006), and the erosion law suggested by Blanc (2008). These erosion laws do not give any accurate result and do not have good consistency; however, they are still accepted to be implemented for various simple problems.

Egashira Erosion Law
According to Egashira (1993), the sketch that described the mechanism of debris flow to erode a sloped surface is shown in Figure 1. The sediment deposition rate on the slope eroded by the debris flow is formulated as follows (Pastor et al., 2014): In which er is the erosion rate; v is the average velocity of debris flow; h is the depth of landslide flow; Ds is the distance the debris flow travels in a time period of Dt. In which θe is the equilibrium angle of the eroded slope.
The Egashira erosion law is added with the empirical factor of K to improve its accuracy. This addition is based on the debris flow research in Tsing Shan in 1990 and 2000; therefore, the Egashira erosion law is modified into (Blanc, 2008):

McDougall and Hungr Erosion Law
According to Hungr (1995), the erosion process that occurs during the debris flow path is formulated into: where Es is the rate of erosion.
McDougall and Hungr (2005) assumed that the estimated total volume exiting the zone (Vf) is the sum of the initial volume (V0) and the volume taken from erosion (Ve). Therefore, the Es value can also be obtained from the estimated total volume entering the zone (V0), the estimated total volume exiting the zone (Vf =V0 +Ve), and the approximate average path length of the zone (l), with the following formula, Es is defined as the average addition rate which comes from the natural addition of erosion process with displacement (McDougall and Hungr, 2005).

Blanc Erosion Law
The Blanc erosion law (2008) is a new equation which is a combination of the Egashira erosion law (1993) and Hungr (1995) erosion law. According to Blanc (2008), the erosion process occurring during the landslide flow path is formulated as follows,

Deposition Law
The Egashira law can also be applied to explain the process of deposition of debris flow during its travel. If θe > θ, the erosion rate (er) is negative. This means deposition exists; this then causes the debris flow to decrease during the travel. If er < 0, then < 0. In the calculation, the Egashira law considers that the amount of deposition is not based on the reduction of the particle height (reduction of debris flow depth), but based on the particle velocity (debris flow velocity) equal to 0 m/s. Egashira law is only used to calculate the erosion process. Therefore, if the erosion rate is negative, there is no change in the debris flow volume. The erosion process is more dominant compared to the deposition when debris flow occurs. Therefore, the longer debris flow travels result in a larger volume when it stops.

Erosion Law Adaptation to Numerical Model
For the sake of the simulation accuracy, an erosion law was implemented to the numerical model. A modification was made on the flow continuity equation, which resulted in the following equation:

Simulation Model
The simulation model is based on continuum mechanics in the form of Depth-Integrated Models. This simulation model used Eq.
(1) to calculate flux and Eq.
(2) to obtain the thickness of landslide mass for each time period. Shear stress (T) have resulted from two constitutive equations of landslide material, which are the Egashira and Coulomb constitutive equation. Erosion and deposition occurred along the landslide path were calculated with erosion laws of Egashira, Itoh and Takeuchi (2001), McDougall and Hungr (2005), and Blanc (2008).

Momentum Equation
If the Eq. (1) is described on x -y coordinates, the results are Eq. (20) and Eq. (21). Figure 2 shows the description of finite-difference model meshing (modified by Fathani, Legono and Karnawati, 2017). Flux on x-axis can be rewritten, as shown in Eq. (22).
Whereas n shows a time of n, i and j is the grid on x and y-axis, respectively.
Eq. (22) can be rewritten on Eq. (23) as follows: The PMX value is shown in Table 1, and PMGZ is described in Eq. (25) as follows: In which,   Value of u2,u1,u3 and v2,v1,v3 is shown in Table 2 as follows: The shear resistance value of hyperconcentrated solid-liquid mixture model is obtained from Eq. (4) up to Eq. (7) at axis x and y as follows: If Eq. (44) is fulfilled, the mesh stops, but if it is not, then the flux on the next step can be calculated with the equation as follows: If (0t 't) condition is not fulfilled, then the flux is not equal to zero on the interval of t; therefore, the flux on the next step can be calculated with Eq. (46) as follows:   1/2 1/2 1/2, 1/2, The flux on Eq. (39), (40) and (41) are applied on moving condition, therefore Eq. (46) is semiimplicit, as shown in the equation as follows: 1/ 2 1/ 2, 1/ 2 1/ 2, In which, The flux on y-axis is obtained with a similar method:

Boundary Condition
As described in the continuity equation, erosion (er) is a function of landslide mass height. When er is negative, it is considered that volume reduction does not occur; therefore, there is no reduction of h caused by erosion. Apart from erosion, the change in h value arises from the condition in which the landslide height at one point decreases due to the movement to the next point according to the exiting flux.
Large t can cause a negative h value. A correction is needed to avoid a negative momentum value on the next step. The boundary condition was corrected to generate an accurate value. The correction was carried out by replacing the h < 0 value into h = 0 in the adjacent mesh according to the exiting flux.

Numerical Model Flow Chart
Calculation steps on the simulation program are shown in Figure 3. Input parameters can be determined directly on the program, while topography can be exported with Excel data. The calculation steps were started with an initial condition, where the initial data was used as the input for calculation on 1 st step. The initial condition was in static, while the 1 st step was for starting of the movement. The thickness of new landslide mass was calculated by implementing the three erosion laws. The calculation kept working until stopping term was fulfilled. If the stopping term is not fulfilled, the calculation procedure will redefine the thickness of the new landslide mass on the next step (n = n+1).
Flux calculation flow chart is shown in Figure 4. If the total of the occurring flux is equal to zero, it means there is no movement; thus, the landslide mass stops. The flux with the value of zero results in the flux value on the next step to be equal to zero. If the flux calculation result does not equal to zero, then the flux on the next step is defined based on previous equations.

Manual Calculation
In the simulation program, the solution from the partial differential equation solution was approached with a finite difference scheme. To prove that the simulation program is stable, a manual calculation was also conducted in parallel in one of the observed points. A Calculation was conducted in the first step up to the third step as the description of calculation for static condition and moving condition.

Numerical Simulations with Additional Erosion Formula
In this research, the newly developed model was examined on an imaginary slope. The imaginary slope was inclined plane by 35° angle and combined with a flat plane at the bottom. The landslide mass was located on the top of the ellipsoid-shaped simple slip plane. 2D image of landslide mass and 3D visualization on the imaginary slope are shown in Figure 5. This imaginary slope was also used in the research of Fathani, Legono and Alfath (2017) to model the landslide movement with constant landslide volume.
The erosion formula of Egashira, Itoh and Takeuchi (2001) needs an input of the solid phase concentration of volume in a packed state (c*), angle of the eroded slope (θ), equilibrium angle of the eroded slope (θ e ), and K factor. The numerical simulation result of the landslide movement by the erosion formula of Egashira, Itoh and Takeuchi (2001) showed that the landslide mass started to move in relatively low velocity in the first second and had not reached the flat plane yet. By that time, the maximum landslide mass thickness was 8.612 m. After 3 seconds, the landslide mass moved significantly and spread on the flat plane. The maximum thickness of the landslide movement was 6.607 m. The increasing volume still occurred in the 5 th second with a large increasing rate; the maximum landslide mass thickness was 3.028 m.
This happened because the velocity kept increasing until it reached the maximum velocity. After 7 seconds, the landslide coverage still occurred with a wide range of movement on the flat plane. In the 12 th second, the certain mesh was still moving with a relatively low velocity. The iteration on the simulation program was limited to 2000 th step (in the 20 th second) because of very low velocity (less than 0.1 m/s); therefore it did not change the landslide coverage area. (V0), the estimated total volume exiting the zone (Vf), and the approximate average path length of the zone (l).
The result of the numerical simulation on the landslide movement by McDougall and Hungr (2005) showed that the landslide material started to move in relatively low speed at the time of 1 second. The erosion formula from McDougall and Hungr (2005) needs the input parameters in the form of the estimated total volume entering the zone By that time, the deposit thickness was 8.603 m and still had not reached the flat plane. In the 3 rd second, the landslide mass moved significantly and spread on the flat plane. The maximum thickness of the landslide movement was 6.218 m. The increasing volume still occurred in the 5 th second with a large increasing rate; the maximum landslide mass thickness was 3.027 m. This happened because the velocity kept increasing until it reached maximum velocity. In the 7 th second, the landslide coverage still occurred with a wide range of movement on the flat plane. In the 12 th second, several of certain mesh parts were still moving with a relatively low velocity. The iteration in the simulation program was also limited to 2000 th step (in the 20 th second).
The erosion formula from Blanc (2008) Blanc (2008) showed that the landslide material started to move in relatively low velocity in the 1 st second.
In the 3 rd second, the landslide mass moved significantly and spread on the flat plane with increasing velocity. The volume was still increasing rapidly in the 5 th second. After 7 seconds, the landslide deposit was spreading with la ow velocity. In the 12th second, some of mesh was still moving with a relatively low velocity. The simulation result of the threeerosion formula is shown in Figure 6. Egashira, Itoh and Takeuchi (2001) McDougall and Hungr (2005) Blanc (

DISCUSSION
The result of the simulation suggested that the incorporated erosion law could increase the precision, particularly on the landslide areal extent. The volume addition can be controlled by assuming the erosion rate value where Egashira, Itoh and Takeuchi (2001) and Blanc (2008) use K factor. However, the erosion formula from Blanc (2008) is simpler and does not require any input of the solid phase concentration of the volume in a packed state (c*), and equilibrium angle of the eroded slope (θe). The erosion formula from Blanc (2008)  In the newly developed simulation program, back analysis can be conducted with trial and error in the soil and the erosion parameters to get a landslide coverage area coinciding with the actual condition. Soil parameters can be obtained from laboratory tests or calibration procedures based on field actual conditions. The erosion parameter (Es, K) is obtained by trial and error with pre-determined soil parameter value assumption. In this case, the erosion law from McDougall and Hungr (2005) provides a better solution as the coordinate and the flow path of the eroded area can be determined. Therefore, it is possible to determine which part will be eroded if being passed by landslide flow whereas the erosion model of (Egashira, Itoh and Takeuchi (2001) and Blanc (2008) depended on the empirical factor value (K) with the assumption that the erosion occurring along the eroded flow path are homogeneous.

CONCLUSIONS
This newly developed landslide-movement numerical model that considers erosion and deposition can be used to predict the velocity, affected area, flow depth, and depth of landslide deposit. Using the finite difference method, this model can give a more accurate prediction by considering the volume change in the equation.
The simulation result shows that the landslide affected area is enlarged by adding the erosion formula from Egashira et al. (2001), McDougall and Hungr (2005), and Blanc (2008). The application of McDougall and Hungr (2005) erosion formula is suitable for landslides that have previous event records, or to model landslide movement with back analysis by comparing landslide deposited volume and source volume. The erosion formula from Egashira, Itoh and Takeuchi (2001) can be applied using the rheology model from (Egashira, Miyamoto and Ito (1997) by inputting the concentration of the solid phase of volume in a packed state (c*), and the equilibrium angle of the eroded slope (θe). The erosion formula from Blanc (2008) can be used on various rheology models, as it only needs K factor value as the erosion rate determinant.
To develop the simulation program, calibration on the actual cases is advised. Furthermore, the Lagrangian model or mesh-free method can also be used, particularly on the landslide material threshold which is usually displaced very far with thin deposit thickness.