Estimating the Velocity of Landslide Movement Using Visco-Plastic Model in Jeruk Sub-village, Kulon Progo District, Yogyakarta, Indonesia

A ground movement occurred in March and November 2017 on the hills and paddy fields in Jeruk Sub-village, Kulon Progo District, Yogyakarta Special Province. The landslide movement destroyed two houses in the village and the land is still moving especially in the rainy season. The mitigation of landslide hazard requires understanding of landslide triggering factors and its movement mechanism. This paper applies the slope stability analysis and visco-plastic model to predict the movement mechanism and velocity of a translational landslide. The sliding mass is modeled as a low plasticity silt (homogenous soil). The Limit Equilibrium Method is used to estimate the safety factor, whereas the shear strength parameters on the slip surface were determined by using the back analysis approach. The results of the slope stability analysis showed that the shear strength parameters and the fluctuation of groundwater level strongly influence the stability of the landslide. From visco-plastic model simulation, this slope has slow movement velocity with the range of 11.31 to 175.88 mm/day. It is clarified that the velocity of landslide movement is influenced by soil strength parameters, coefficient of dynamic viscosity, and groundwater level fluctuation.


INTRODUCTION
Indonesia is one of the countries affected by the landslide disaster, resulting in casualties and economic losses. Landslide is one of the natural disasters which threatens most areas in Indonesia because of the geological condition of the country which is mainly covered by weathered volcanic rocks in the mountainous and hilly areas intersected by faults and joints. Moreover, the high rainfall intensity is one of triggering factors of landslide with land use change and deforestation causing the increase of landslide hazard frequency. The occurrence of these natural phenomena cannot totally stop, but the potential landslide mechanism should be identified to minimize the losses of human life and economic value.
In March and November 2017, Wilopo and Fathani (2017) reported the land movement occurred which is identified by the cracks in the hills and paddy fields in Jeruk Sub-village, Gerbosari Village, Samigaluh District, Kulon Progo Regency, Yogyakarta Special Region. The landslide is located on the eastern side of Menoreh Mountains where the morphology is steep hillsides (Hadmoko, et al., 2010). The hilly area is generally used by the community as a plantation area and residential area. This can contribute to the high landslide susceptibility. The landslide destroyed two houses in the village and is still moving during the rainy season.
One of the mitigations efforts to reduce risks is by analyzing the causing factors and mechanism of landslide that may be used to develop a warning system. In order to obtain an effective and reliable result, there is a necessity to investigate the geomorphology, geology and geotechnical conditions of the landslide area and to predict the landslide movement and its correlation with the triggering factors such as rainfall and slope hydrological condition.

GEOLOGICAL CONDITION OF THE STUDY AREA
The study area is located in Jeruk Sub-village, Gerbosari Village, Samigaluh Sub-district, Kulon Progo Regency, Yogyakarta Special Region at S07° 41'16.8" E110° 09'47,5" Zone 49S (Figure 1). The morphology of Jeruk Sub-village, Gerbosari Village is steep hillsides (ranging between 19-58°) located on the eastern side of the Menoreh Mountains. The morphology in this area is controlled by lithology and geological structure. The elevation of the study area is between 312 and 620 m from sea level. The hilly area is generally used as a plantation and residential area. The Tinalah River flows permanently in the valley of this area.
According to the regional geology, this area lies in the Kebobutak Formation which is generally composed of andesite breccia, tuff, lapilli tuff, agglomerate and intercalation of andesite lava flow (Rahardjo, et al., 1995). Based on the site investigation, the lithology is dominated by the andesite breccia with fragment size up to gravel. The rocks have been intensively weathered to form a thick layer of soil. Moreover, the rocks in some places were also found to have alteration process with clay mineral composition. This thick layer of soil (around 7 m thick) is one of the controlling factors of the landslide movement. Weathered soil easily absorbs and store water so it can increase groundwater level and potentially reduce the slope stability.
Geological structures such as fractures and tension cracks can be the controlling factors of the landslide movement. The water infiltrates into the discontinuity weak plane in the landslide body, resulting in the progressive movement of landslide at this study area.
Based on the results of geological survey, two major cracks were found on the hill slope and in paddy field as shown in Figure 1. The major crack on the hill with the slope inclination of 55° ( Figure 2) is identified as the crown of landslide with the movement direction to N30°E (relatively northeast). This crack has a length of 137 m with a high difference of land subsidence from 3 to 7.3 m. Based on villagers' information, this crack appeared one month before the mass movement. Due to the heavy rainfall on 28th to 30th November, 2017, landslide movement occurred with new major cracks. It is located over the older major crack with the dimension of approximately 130 m in length and 1-3 m in depth. Minor cracks are also found around the crown of landslide with varying dimensions. Minor cracks have a length of up to 73 m in the same direction. This minor crack has a depth of up to 1 m.  Cracks in the paddy fields ( Figure 3) have a length of up to 143 m in the direction of movement of N310°E (relatively northwest). Cracks were formed in the middle of the paddy fields and resulted in land subsidence in the western part. The cracks have a width of 30 cm with a depth of up to 1 to 5 m. In addition, a spring also found in the paddy fields, where in certain places, the water flows directly into the cracks.
In the vicinity of the valley at the western part of the paddy field, there are small avalanches with rotational sliding where the landslide material is in the form of soil and rock ( Figure 4). This landslide has the same direction as the existing crack in the paddy fields, i.e. N310°E (relatively northwest). The crown length of the landslide reaches 30 m with the slope inclination of 44°. Cracks were also found around residential houses close to paddy fields. This crack has two main directions, namely N54°E (relatively northeast) and N330°E (relatively northwest). This crack has a length of up to 18 m with a fractional opening width of 5 cm. This cracks also resulted in slight damages to the residential houses.

METHOD OF ANALYSIS
3.1 Visco-Plastic Model for Dynamic Simulation Ranalli, et al., (2009) explained slow slope movement is typically associated with "creep" behavior, since the soil can be characterized by a viscous response, if the soil mass starts to move slowly. For continuously moving landslides, a dynamic analysis should be adopted instead of a classical static approach (Ranalli et al., (2009) and Corominas et al., (2005)). The classical static analysis is suitable to determine the slope stability but is not able to model the actual kinematics of the soil mass behavior (Faris & Fathani, 2013).
In the limit equilibrium method, the soil shear strength is usually defined by the Mohr-Coulomb criterion, and the instability condition occurs when the equilibrium is changed by pore water pressure increase and a consequent reduction of the effective stress level. A constant instability force could exist for a given piezometric level and initiate a slope movement with constant acceleration and a corresponding velocity, linearly increasing with time. A possible explanation is the effect of a viscous resisting component of the material. In this case, the mass velocity can be related to the excessive shear stress by different viscous laws, like the Bingham's law, which shows a yield point and a subsequent linear relationship (Ranalli, et al., 2009).
Soils at the critical state are like a visco-plastic fluid, which will flow for applied stresses greater than the critical state shear strength (Angeli, et al., 1996). The shear viscosity () at critical state is the desired parameter for post-failure analysis of soil (Locat and Demer (1998) and Komamura and Huang (1974)). According to Edger and Kalrsrud (1985), the shear viscosity of soil plays an important role in landslide.
In this condition of research area, the movement of landslide mass is controlled by viscous resisting force which depends on the coefficient of dynamic viscosity (C) parameter. This parameter was obtained from the calibration process (Ranalli, et al., 2009) by simulating the velocity of landslide taking into account the groundwater level fluctuation generated from visco-plastic model. The coefficient of dynamic viscosity depends on shear viscosity of the soil and thickness of shear band z which is difficult to determine during the site investigation and laboratory test. Corominas, et al., (2005) adopted the equation for the dynamics velocity of the landslide that has creep behavior which is considered sensitive to water pressure at the slip surface. The equation of the dynamic velocity can be written as: In order to change the equation to become general form (Corominas, et al., (2005); Faris and Fathani (2013)), the Equation (1) needs to be modified so that it can be understood easily by dividing both sides with terms m. Then, the Equation (1) can be written as follows: where, = .
where, C is the coefficient of dynamic viscosity, e is exponential for t (time function), v is the velocity (m/s), z is the shear band thickness (m), η is the dynamic viscosity (N.s.m -2 ), m is mass of the soil, γ is unit weight of soil, α is the inclination of slope.
The pore water pressure was not measured directly, but it was calculated by observing the depth of groundwater level. Assuming a parallel flow to the slope surface: where, γw is the specific weight of water, l is the thickness of the sliding mass and Dw is the depth of groundwater level.

Stability Analysis for Infinite Slope
There are several methods to determine slope stability. One of the most commonly used methods is the Limit Equilibrium Method (LEM). The slope that extends for a relatively long distance and has a consistent subsurface profile may be analyzed as an infinite slope (Naresh & Edward, 2006). The failure plain for this case was considered as parallel to the surface of the slope and LEM can be applied readily.
The factor of safety (FS) is defined as the ratio of resisting shear strength (τr) to driving shear stress (τd). Thus, the factor of safety is given by: If a saturated slope formed by cohesive soil has seepage parallel to the surface of the slope as shown in Figure 5, the effective shear strength parameters are used in the analysis. The pore water pressure (Pw) was also considered in normal force. Thus, Equation (7) is introduced.
where, l is the height of the slice measured in the field.
For an infinite slope analysis, the factor of safety is independent of the slope depth (h) and depends only on the effective strength parameters (c-). Also, in the critical condition of the slope model (FS = 1.0), the shear strength parameters can be determined by using back analysis (Parmar, 2016).

Soil Model
Pastor and Picarelli (2010) generally classified run-out model for landslide hazard and risk mapping into two main models, i.e., empirical models and rational models. In the rational model, there are sub-models such as a continuum model (3D models based on mixture theory, velocity-pressure models, depth integration model) and a discrete model, which is used in this study. The soil model was constructed based on the assumptions that there are two different layers forming the slope. The first one is the upper layer which is an unstable zone acting as sliding mass and the second one is the lower layer which is the stable zone and composed by andesite rock and breccia.
In this study, the upper layer is composed by homogeneous soil type, and the slip surface is assumed as a complex type of translational and earth flow landslide movement, which is suitable to be used in the LEM. The soil is classified by USCS as low plasticity silt obtained from laboratory tests. The elevation models that are through the cross-section line A-B varies between 420 m at the toe and 622 m at the top above the sea level ( Figure 5). After constructing the model, the parameters (bulk weight, cohesion, internal friction angle) are then used in calculating the factor of safety factor using LEM in order to obtain the most critical safety factor. In this case, groundwater level and the depth of slip surface were considered. The depth of slip surface was verified from the site investigation. Figure 5. The detail soil model design with two different layers (section A-B in Figure 2).

Slope Stability Analysis
The input parameters used in this study is based on the existing geotechnical data (see Table 1) and the parameters to determine the safety factor such as the variation of the cohesion and internal friction angle are used to obtain the critical factor of safety in static condition (see Table 2). Based on the field investigation, the groundwater level is 0.36 m from the surface. Using the cross-section in Figure 5, then the movement velocity of landslide can be calculated. The factor of safety of slope model was lower than critical condition when the parameters from the result of laboratory tests were used. The shear strength parameters from the laboratory test did not give a reliable result for slope stability analysis. Thus, the back analysis is used to determine the shear strength parameters at a critical condition (FS = 1.0).   Based on slope cross-section in Figure 5 and the result of calibration of shear strength parameters in Figure 6 and Figure 7, the internal friction in Table 2 is lower than the result in Table 1 and the cohesion in Table 2 is greater than in

Estimating the Velocity of Landslide
Using the cross-section in Figure 5, the velocity of landslide can be simulated by visco-plastic model. Coefficient of dynamic viscosity (C) is needed to simulate the velocity of landslide movement. Coefficient of dynamic viscosity (C) was defined as the ratio of viscosity of soil and thickness of shear band z. It was obtained from calibration process (Ranalli, et al., 2009) by simulating velocity of landslide taking into account daily groundwater level fluctuation (Figure 8) from 11/9/2017 to 12/10/2017. To determine the most reliable viscous parameter through calibration process, slow movement of landslide (Cruden & Varnes, 1996) and shear thickness of slip surface (0.5 m) were considered. The model simulated the velocity of the landslide, on the basis of soil parameters that have been considered as the most appropriate parameter via back analysis previously performed.
The value of coefficient of dynamic viscosity was 4.93×10 8 Nsm -3 , which was based on the calibration process. The value of coefficient of dynamic viscosity provided a good agreement of landslide velocity to the fluctuation of groundwater ( Figure 9). Moreover, the value of coefficient of dynamic viscosity was suitable for the low plasticity silt which formed the slope in the study area. The result of the value of coefficient of dynamic viscosity was verified by movement displacement of landslide which occurred on 28th November 2017. The displacement of movement had been measured manually by head of village. Therefore, the value of coefficient of dynamic viscosity of 4.93×10 8 Nsm -3 was used for further analysis. Moreover, it has been proven that the fluctuation of groundwater gives a great influence on the landslide velocity.   Figure 9 shows the correlation between landslide velocity and groundwater level fluctuation. The uses of these parameters are limited for local conditions as representative for particular area of landslide. The results of this study can be implemented to design the monitoring and mitigation system of the landslideprone area.

CONCLUSION
Based on the site investigation, the significant landslide movement occurred twice (in March and November 2017) with major cracks destroying a major road of the village. The landslide is still moving during the rainy season. The visco-plastic model can predict the velocity of landslide movement and is strongly depending on the fluctuation of groundwater level, the coefficient of dynamic viscosity and engineering properties of the sliding mass. The coefficient of dynamic viscosity is determined by a calibration process. The shear strength parameters are determined from back analysis to reach a critical condition of the slope. The result of visco-plastic model shows the velocity of Jeruk Landslide can be classified as slow movement with the range of 11.31 to 175.88 mm/day. The proportional relationships between the velocity of landslide movement and the increase of groundwater level were caused by heavy rainfall events were found. Moreover, the Jeruk Landslide is defined as complex type of translational and earth flow landslide movement.
The approaches used in this study could not estimate the distribution of landslide material. In the future, it is suggested to use a detail geotechnical analysis and real-time monitoring devices to obtain the actual physical properties of the sliding mass and sliding surface and to determine the warning criteria for slope movement.