Simulation of the Initiation and Motion of Seismically Induced Aso-Ohashi Landslide During the 2016 Kumamoto Earthquake

To verify the mechanism of Aso-Ohashi landslide which is the seismically induced slope failure triggered by 2016 Kumamoto earthquake, field investigations and numerical simulations with terrain data measured by airborne laser scanners are carried out. A finite difference method (FDM) model is set up by the true terrain data. The stability of landslide region under shear strength reduction method and earthquake acceleration wave input were discussed. The effect of foreshock which caused strength reduction indicated that the critical value of strength was c = 16 kPa and φ = 33.87 degree. The earthquake wave made only steeply dipping region unstable and moved, but landslide hardly occurred in gently dipping region. In order to verify that only the steeply dipping region was moved due to earthquake, a GIS (Geographic Information System)-based depth average Bingham model which estimates the movement of slide mass is presented.


INTRODUCTION
Landslides are geological processes that often occur due to slope failure, and they can happen in almost every conceivable manner, either steeply dipping or gently dipping. Landslides can be induced by human activities such as excavation or cutting into existing slopes, or they can result from strong earthquakes or heavy rainfall. Seismically induced landslides can cause big disasters such as the Hattian Balla landslide triggered in 2004 in Kashmir, Pakistan [1], the 2006 Leyte landslide triggered by a small nearby earthquake after a heavy rainfall in the Philippines [2], and many big landslides triggered by the 2008 Wenchuan earthquake in Sichuan, China which caused 15000 geohazards [3]. Aso-Ohashi landslide was a major slope failure in the 2016 Kumamoto earthquake. The landslide mass travelled a distance of 800 m and destroyed an important bridge connecting Minamiaso village to the city of Kumamoto.
Much study on the mechanism of landslides was conducted. Research results showed that one of the main contributing factors of Tsaoling landslide triggered by rain or earthquakes was the weathering of shale [4]. The Rufi slide was largely due to dipping of interbedded conglomerates and marls. Primary sliding could be attributed to strength degradation in the marls due to weathering processes and saturation during periods of heavy precipitation [5]. The LantaKhola landslide was also caused by weathering and heavy precipitation. The schist was weathered to fine sand and silt with lower shear strength and permeability according to field observation. Numerical modelling confirms that blocks were collapsing and eventually rolling down the slope after cloud bursts [6]. Sliding surfaces of catastrophic landslides triggered by the Mid Niigta prefecture earthquake were formed at the boundary between the overlying permeable sandstone and underlying siltstone or along the bedding planes of alternating beds of sandstone and siltstone [7]. Observations, data and analyses show that the cause of distress to a single-family residence located adjacent to a major highway cut slope was landslide which was induced by a deep-seated translational excavation [8]. The rapid motion of Qianjiangping landslide was caused by the rate effect of the black silt layer during the motion phase after the creep failure [9].
Landslide simulations are studied by some numerical methods, such as the discrete element method (DEM), discontinuous deformation analysis (DDA), the finite element method (FEM) and FDM. FEM and FDM are mesh-based methods in the analysis of initiation (slope failure) and motion of landslide. They have been extensively used as an alternative to and in conjunction with the limit equilibrium method in the analysis of landslide. For example, a finite element analysis of the Senise landslide of 26 July 1986 was performed using a nonlocal elasto-viscoplastic model [10], prediction of slope failure was simulated by using FDM [11]. A number of new and sophisticated numerical models have recently been developed, and modelling of the propagation stage has been largely carried out. Discrete element methods, such as DEM and DDA, have been widely used to analyze the kinematic and run-out behaviour of discontinuous material such as blocky rock masses. For example, dynamic DDA simulations of the Vaiont landslide have been conducted, and the influence of the geometric discontinuity on the kinematic behaviour of landslides has been analyzed [12]. More recently, Utili and Crosta [13] used a DEM model to simulate the evolution of natural cliffs subjected to weathering. The most commonly adopted numerical methods are the continuum approaches, such as FDM and FEM. While FDM uses the differential form of the governing partial differential equations, FEM is based on their integral form and requires solving a global equation system. They are suitable tools for simulating the stress and deformation fields around underground excavations. However, due to the lack of an internal length scale, standard strength-based, strain-softening constitutive relationships cannot reproduce the localization of failure, as the underlying mathematical problem becomes ill-conditioned. In discrete (or discontinuous) modelling techniques, commonly known as the discrete element method (DEM), the material is treated as an assembly of independent, rigid or deformable blocks or particles. Unique features of the DEM are the abilities to: (i) capture finite displacements and rotations of discrete bodies, including complete detachment, and (ii) automatically recognize new contacts as the simulation progresses. Unlike continuum methods, which are based on constitutive laws, DEM relies on interaction laws. DDA is one category of DEM which uses an implicit (and thus unconditionally stable) time integration scheme. The applications of these methods were mainly in the field of granular materials and jointed structures. Further developments made DEMs also capable of explicitly simulating failures through intact rock material.
In this study, the initiation and motion of seismically induced Aso-Ohashi landslide during the 2016 Kumamoto earthquake were simulated which reproduced the destruction process by FDM. It is convenient for the simulation of initiation and motion to shear the same finite difference mesh of FDM. Most areas of this landslide were gentle slopes. Generally, slopes with large inclination angles are extremely sensitive to the acting seismic wave, whereas gentle slopes are less likely to be affected. Therefore, this region could hardly slide as a whole. Failure mechanism of a different region was explored which showed that steeply dipping region was unstable and moved by earthquake wave, and gently dipping region might move by the shear force of the sliding mass. A movement simulation model was presented and this model verified that, due to earthquake, the steeply dipping region was moved, but the gently dipping region was not. The paper is organized as follows. In Sect. 2, we summarized the 2016 Kumamoto earthquake. Sect. 3 is the outline of Aso-Ohashi landslide. Sect. 4 briefly introduces the terrain data measured by airborne laser and compares the terrain before and after earthquake. Sect. 5 presents a FDM model to reproduce the initiation of landslide and a Bingham model to simulate the movement which was induced by landslide. The results are summarized and discussed in Sect. 6. The earthquakes generated many of sediment disasters including debris flows, collapses and landslides (Tab. 1). Most of the disaster happened in Kumamoto. At least 49 people were killed and about 3000 people injured in the earthquakes. The earthquakes seriously damaged the infrastructures in Kumamoto and the circumjacent area. One of the largest seismically induced landslides was located on the National Road 57 and destroyed an important bridge named Aso-Ohashi Bridge.

OUTLINE OF THE ASO-OHASHI LANDSLIDE 3.1 Landslide
The landslide is located at the west of the caldera of Mount Aso where a river named Kurokawa flows through. The gorge of the river is about 80 m deep and a 200 m long Aso-Ohashi Bridge spanned the disaster region and the other side across the river. The Aso-Ohashi Bridge also collapsed in the landslide.
Based on field investigation, there were three regions from above to below: steeply dipping region, erosion region and gently dipping region. In a general way, the acceleration caused by earthquake may expand in convexity of steep slopes. So the steeply dipping region might become unstable first due to the earthquake. There were weathered strata in the bottom and rubble in the side of sliding surface. This showed that the surface of landslide might be a weathered layer. Because there was no evidence to suggest that there was groundwater and there was no rainfall before earthquake, groundwater and rainfall were not taken into account as the cause of landslide. The gently dipping region appeared near the river bank. There was a scarp which was several meters deep and some stratifications outcropped (Fig. 2). This region gently dipped, and the width was continuous with erosion region and steeply dipping region. It was estimated that this region was under shear force induced by the mass movement of upper unstable region. The landslide mass almost vertically flew into the Kurokawa river (the angle was about 70 degree) but did not result in blocking it. The mass was found on the opposite bank which was about 30 to 60 meters height from the riverbed.

Geologic Setting
The rock type of landslide region is pyroxene andesite lava, which belongs to Late Pliocene-pleistocene volcanic rocks, Tertiary to Quatemary Period [14]. Fig. 3 indicates the geology of Aso-Ohashi landslide. In general, the densities of andesite (as base rock, not collapsed) range from 2.1 g/cm 3 to 2.8 g/cm 3 , the cohesions range from 10 MPa to 40 MPa and the friction angles range from 45 degree to 50 degree. Yamashita et al. [15] measured cohesion and friction angle of collapsed andesite rock in Mount Aso with some experiments. Cohesion of collapsed layer is from 0 to 17 kPa, and friction angle is from 15. 5 degree to 19 degree. Dang et al. [16] tested some samples from Aso-Ohashi landslide. Friction angle at peak was 41.8 degrees and friction angle during motion was 36.1 degrees, and cohesion was estimated as 30 kPa.

TERRAIN DATA FROM AIRBORNE LASER SCANNERS SYSTEM 4.1 Airborne Laser Scanners System
The airborne laser scanning system is a compact laserbased system designed for acquisition of topographic and return signal intensity data from a variety of airborne platforms. The data is computed using range and return signal intensity measurements recorded in flight along with position and attitude data derived from airborne Global Navigation Satellite System (GNSS) and inertial measurement unit (IMU). Laser distance measuring device, IMU and GPS receiver antenna are installed in aircraft. GPS antenna measures the position of the aircraft and IMU measures the attitude [17]. Laser is emitting ten thousand or hundred thousand short bursts of light every second, which will measure range and reflectance of objects on the earth surface.
Airborne laser scanners for recording topographic data have been used in various applications [18]. In contrast to microwave radar techniques, lasers are advantageous for wider range measurements because high energy pulses can be realized in short intervals and their comparatively short wavelengths can be highly collimated using small apertures [19]. Laser scanning is not capable of any direct pointing to particular objects or object features. The resulting co-ordinates refer to the footprints of the laser scan as they happen. Laser scanning is characterised by high accuracy, high sampling densities, and a high degree of automation.

Comparison of Terrain Before and After Earthquake
The terrain data used in this paper is provided by the Ministry of Land, Infrastructure, Transport and Tourism of Japan (MLIT) and PASCO Corporation. The data is surveyed by airborne laser scanners system and the accuracy is 1 meter. Fig. 4 shows the shadowgraphs before (obtained on October 27, 2012) and after (obtained on April 16, 2014) the Aso-Ohashi landslide happened.

FDM MODELLING AND REPRODUCTION OF DESTRUCTION PROCESS
The reproduction of collapses was simulated by threedimensional FDM modelling using the Itasca FLAC code. The length of model in x -axis (East-West direction) and in y -axis (North-South direction) was separately 1000 m and 650 m (Fig. 8). The top surface was based on the terrain measured on October 27, 2012 (before landslide happened). The surface terrain data was the same data introduced in Section 4. The top boundary was free. The side boundaries were restrained horizontally and merely allowed to displace vertically. The bottom boundary was fixed in the initial condition, and then imposed by the acceleration of time history to reproduce earthquake loading. In addition, a free field boundary was used to avoid the influence of side boundaries. The model was divided into two groups respectively which were weathered layer and base layer. Weathered layer represented weathered rock or soil in the surface. Weathered layer was of equal thickness in z -axis. The depth was 20 m because the deepest point of collapsed region was about 20 m. In this model, the weathered layer was assumed to be 20 m and the parameters of shear strength referred to collapsed layer mentioned in Section 3.2: cohesion was 20 kPa and friction angle was 40 degree. Base layer was under weathered layer. The deepest was 456.35 m, and the shallowest was 13.57 m. The parameters of shear strength referred to base pyroxene andesite: cohesion was 10 MPa and friction angle was 45 degrees based on the reference range of andesite (Tab. 2). The model had 63360 zones and 71500 grid-points. The constitutive model was Mohr Coulomb model. Mesh size was 10 m (Fig. 9). The more grid-points and smaller mesh produced more accurate solutions, but the computation time was too long in dynamic computation. The strength parameters were reduced when the strength reduction method was used to analyse the slope stability caused by foreshock. After that, earthquake loading was projected on the model.  The earthquake loading used the strong motion records at a station (Nakamatsu station) which was 7 km away from Aso-Ohashi landslide (data from Japan Meteorological Agency). The acceleration time histories of EW, NS directions were projected on the landslide (Fig. 10) in which the maximum acceleration of EW direction is 606.646 cm/s 2 and minimum is −592.832 cm/s 2 , the maximum acceleration of NS direction is 794.265 cm/s 2 and minimum is −733.508 cm/s 2 .
For a slope,

 
If factor of safety resistant force/driving force 1 FS   (1) the slopes will be instable [20], and there are 2 alternatives: 1) Resistant force decreases, when friction coefficient decreases (For example, liquefaction) and underground water table rises by rainfall, 2) driving force increases, when underground water table rises by rainfall and earthquake which generates inertia force. As introduced above, there was no rainfall before landslide which means landslide was not caused by rainfall. In this paper, rainfall was not considered.

RESULTS AND DISCUSSIONS 6.1 Stability Analysis of Foreshock by Using Shear Strength Reduction Method
During earthquake, the shaking of the ground may cause rocks and soils to lose their resistance and behave like liquids. Because landslide mass almost vertically flew into river and landslide mass did not result in dam, it was judged that earthquake wave reduced strength of rock and made collapsed rock liquefaction. Liquefaction is defined as the transformation of a granular material from a solid to a liquefied state as a consequence of increased pore-water pressure and reduced effective stress [21]. This phenomenon will cause landslides and other failures [22]. Because the liquefaction resistance can be consistently related to the strength parameters of friction angle and cohesion in terms of effective stress, shear strength reduction method (SSR) which is an evaluation method of slope stability by reducing shear strength was used [23].

Figure 11
Strength reduction of group 1. Strength of groups 1 and 2 were reduced at the same time. When reduction rate was 80%, the model began to be unstable Weathered layer was given values of cohesion (c) ranging from 24 kPa (120% of initial value) to 12 kPa (60% of initial value) and friction angle (ϕ) from 45.20 degree to 26.72 degree. Cohesion and friction angle were reduced in the same percent. Strength of base layer was not reduced (Fig. 11). When reduction rate was 80% (c was 16 kPa and ϕ was 33.87 degree), factor of safety approached to 1 and the model was unstable. Plastic zones were found in the steeply dipping region but there was almost no plastic zone in gently dipping region (Fig. 12). Fig. 13 showed vector of velocity of the model when reduction rate was 80%. Vector of velocity was shown in green arrow. Fig. 13c, Fig.  13d and Fig. 13e represented vertical sections of steeply dipping region, erosion region and gently dipping region, respectively. Unstable steeply dipping region were obvious because plastic zones were formed and there was vector of velocity in the same zones. The other two regions remained almost stable. There were some plastic zones in the edge of gently dipping region, but there was almost no velocity according to the simulation results. Therefore, we estimated that the gently dipping region could be stable in this condition. Part of the river bank was also unstable, but it was not correlated with landslide. This process could estimate the effect of foreshock which caused strength reduction, and the 80% reduction rate was the critical value. But actually after the foreshock, these two regions remained stable, therefore true strength should be higher than the critical value. Figure 12 Block state of the model (a) and vertical section AA' (b) when reduction rate was 80%

Stability Analysis of Mainshock by Seismic Inputs
The horizontal earthquake acceleration wave of 10 seconds on the bottom boundary was imposed on the basis of result of SSR. The input wave and acquiring wave on the bottom were almost the same (Fig. 14). The amplitude of wave began to get big at about 3.5 sec. At 5 sec and 6.7 sec, acceleration wave was at its two maximum amplitude. From block state in section AA' shown in Fig. 15, the plastic zones appeared at 3.7 sec in active collapsed region and generated in large amounts at 5.3 sec. Through several seconds shock, the plastic zones spread all over active collapsed zone at 9.5 sec. The wave velocity was calculated as about 1.2 × 10 3 m/s, and the travel time from bottom to surface of the model was about 0.27 sec ~ 0.35 sec. Therefore plastic zones at 3.7 sec appeared caused by 3.5 sec input wave, and plastic zones at 5.3 sec were caused by 5 sec input wave. Fig. 16 shows the state and vector of velocity after seismic inputs of 10 seconds. Plastic zones were formed in steeply dipping region and all the river bank, and there was vector of velocity in plastic zones. This indicated that steeply dipping region and all the river bank were unstable after the earthquake action, but gently dipping region remained stable. The gently dipping region was a gently dipping layer. A study of Tang [24] indicated that the primary triggering factor of slide on gently dipping layer was the earthquake. But in this study, earthquake could not make gently dipping region slide based on the result of our simulation. Polemio and Sdao [25] indicated that rainfall was the main triggering factor of landslide.
There was no rainfall before Aso-Ohashi landslide, so rainfall was not the main triggering factor. Stratifications outcropping was found in the scarp of gently dipping region through field investigation (Fig. 2). The earthquake might reduce the shear strength of the layers. The mass of steeply dipping region flew past and produced shear force on the gently dipping layers which made landslide on gently dipping region.

Motion Simulation by a Bingham Model
As mentioned in Section 3.1, the slide mass was found on the opposite bank. Velocity before slide mass flew into river v 0 can be calculated by: where L is the distance between both sides of the river. Using the Eq. (2) and Eq. (3) velocity v 0 can be calculated by: According to the terrain data, 30 m ≤ L ≤ 60 m, 25 m ≤ h ≤ 30m, therefore v 0max = 26.83 m/s, v 0min = 12.24 m/s. For the slide mass to rush to the opposite bank, its velocity needed to be 12.24 m/s. To test if the velocity could reach the 12.24 m/s, we assumed all steeply dipping region was unstable and flew down, and a numerical simulations based on a depth-averaged Bingham model was applied. This model can simulate the propagation of mass movement across complex three-dimensional terrain. A raster grid network from a digital elevation model in GIS can be used as the finite difference mesh, and the governing equations, including the continuity equation and the momentum equations, are depth-averaged and simplified by shallow water theory. [26] The continuity equation of the flow is obtained as: where h is the thickness of flow (m), t is time (s), x and the momentum equation along the y -axis is where α is an ad-hoc coefficient that provides information about the deviation of the vertical velocity profile from uniformity; the coefficient β represents the ratio of the vertical normal stress to the horizontal one; H is the height of the free surface (m); μ is the viscosity (Pa s); θ x and θ y are the angles of inclination of the bed along the x and y directions (°), respectively; tanφ is the dynamic friction coefficient; ρ d is equivalent density of the flow (kg/m 2 ); g is acceleration of gravity (m/s 2 ). The Bingham fluid does not exhibit any shear rate (no flow and thus no velocity) until a certain stress (τ 0 ) is achieved. The viscosity of Bingham fluid is: The fluid flows when: 0 τ τ > (9) For the laminar flow, on the basis of Newton's Law of viscous fluid, the shear stress of horizontal direction is: The Eq. (10) is integrated through the flow depth: The flow is considered as a steady, uniform Bingham flow on a plane at an angle of θ with respect to the horizontal. h(x, y, z) is the total depth of the flow.The yield stress is: The Eq. (14) is integrated through the flow depth, where h = h(x, y, t) =  −  b is the flow depth. When the Bingham fluid starts moving from rest, the determined statements are: In this paper, the DEM (Digital Elevation Model) resolution is 1 m ×1 m. The input parameters are shown in Tab. 3. The depth-averaged velocities are considered as blunt, therefore, α = 1, β = 1. The slide region is steeply dipping region. The effects of seismic waves are not considered, initial velocity of slide mass is assumed to be 0. In 20 seconds, the slide mass will arrive at the bank of river (Fig. 17). At this time, the velocity v 1 ≈ 10 m/s. The velocity v 0 is greater than the final velocity v 1 in the simulation: v 0 > v 1 . We assumed the initial velocity of slide mass was 0. But, if steeply dipping region was unstable and flew down because of the earthquake wave, the initial velocity should not be 0. This verified that the slide mass had an initial velocity. Therefore, the steeply dipping region was unstable when earthquake happened.

CONCLUSIONS
In this study, terrain data was measured by airborne laser scanning. The Aso-Ohashi landslide during 2016 Kumamoto earthquake was reproduced by using numerical simulation. The SSR method was adopted and earthquake wave was input to explore the initiation mechanism of the Aso-Ohashi landslide. From the result of simulation, we could estimate that 1) earthquake wave reduced strength of rock and made collapsed rock liquefaction. 2) The foreshock simulated by SSR indicated that critical values of strength were c = 16 kPa and ϕ = 33.87 degree. 3) Seismic inputs could make steeply dipping region unstable. Through simulation of movement caused by landslide, it was verified that the slide mass was only from the steeply dipping region and had an initial velocity due to the earthquake. Therefore, the steeply dipping region was unstable and flew past. Instability hardly occurred in gently dipping region where the primary triggering factor of slide might produce shear force when mass of steeply dipping region flew past.
Earthquake-triggered landslide is a major secondary hazard following the strong ground motion, especially in hilly regions. The study on the failure mechanism and forming reason of earthquake triggered landslides is of great importance for disaster prevention and mitigation in similar scenarios of massive large landslides of gentle slopes.