Three-dimensional seismic slope stability assessment with the application of Scoops3D and GIS: a case study in Atsuma, Hokkaido

Abstract

Background

The 2018 Hokkaido Eastern Iburi earthquake (Mj 6.7, Mw 6.6), occurred on the 6th of September 2018 and triggered thousands of landslides in pyroclastic fall deposits. The triggered landslides caused destructive damages to the structures and resulted in serious causalities. Hundreds of earthquakes persisted after the mainshock and the Iburi region is of high possibility to suffer major earthquakes in the future. An effective method to assess the seismic slope stability is of great importance for the disaster prevention and mitigation in the Iburi region.

Results

Based on the isopachs of different pyroclastic fall deposits mantled in the study area, GIS was employed to process the input soil layers and construct the 3D domain. By applying different horizontal pseudo-acceleration coefficients in the Scoops3D program, the factor-of-safety maps of eight cases were obtained. After validating by the coseismic landslide inventory, the performance of the computed results was evaluated.

Conclusions

A horizontal pseudo-acceleration coefficient between 1/2 and 2/3 of PHGA is suitable for seismic slope stability assessment in pyroclastic fall deposits. The catastrophic Tomisato landslide and Yoshinoya landslide were correctly predicted. Scoops3D proves to be an effective and efficient method for guiding disaster mitigation and management.

Introduction

Slope stability assessment on a regional scale represents a vital aspect of geoenvironmental disaster prevention and mitigation, and has been commonly utilized in slope stability analysis especially during critical rainfall events. Dozens of infinite slope analysis approaches and models (one-dimensional or two-dimensional), such as Shallow Landsliding Stability Model (SHALSTAB, Montgomery and Dietrich 1994), the distributed Shallow Landslide Analysis Model (dSLAM, Wu and Sidle 1995), the Stability Index Mapping (SINMAP, Pack et al. 1998), the Transient Rainfall Infiltration and the Grid-Based Regional Slope-Stability Model (TRIGRS, Baum et al. 2002), have been proposed and applied in previous researches on the basis of the limit equilibrium theory. These models have advances in assessing slope stability under intense rainfall, as they incorporate the variation of groundwater table or soil moisture in response to rainfall. However, one-dimensional or two-dimensional models can not consider the three-dimensional variations of topography and soil conditions in actual slopes and commonly cause conservative computation results (Cavoundis 1987; Duncan 1996).

One challenge in slope stability analysis is how to locate the potential sliding surface. Scoops3D, developed by the U.S. Geological Survey, can regionally evaluate three-dimensional slope stability throughout the digital elevation model (DEM) utilizing 3D method of columns approach (Reid et al. 2015). Scoops3D allows the user to define a series of horizontally and vertically extended points (centers the spheres) and a certain radius increment. Then the spherical surfaces intersected by the spheres and the DEM will serve as the potential sliding surfaces, and the stability of each potential landslides encompassing many DEM cells will be computed. In addition to incorporate complex topography and 3D distributions of subsurface material parameters, Scoops3D can also include the effect of earthquake by applying a horizontal seismic loading to the potential sliding mass in a pseudo-static analysis (Reid et al. 2015). Moreover, the wide application of Geographic Information Systems (GIS) and the availability of Digital Elevation Model (DEM) have significantly facilitated the application of Scoops3D in the assessment of slope stability on a regional scale.

Scoops3D was applied for stability analysis of various aspects and areas in previous studies. It has been employed to evaluate the stability of volcano edifices (Vallance et al. 1998; Reid et al. 2001; Vallance et al. 2004; Reid et al. 2010), coastal bluffs (Brien and Reid 2007) and loess slopes (Xin and Zhang 2018). Tran et al. (2018) utilized Scoops3D and TRIGRS to predict rainfall-induced landslides. Liu et al. (2018) used Scoops3D to evaluate regional slope stability considering variation of water level in reservoir. While Scoops3D has been validated to be an effective way for slope stability analysis and landslide prediction especially in response to rainfall infiltration in previous studies, it has not been applied to slope stability assessment in pyroclastic fall deposits under seismic loading yet.

The 2018 Hokkaido Eastern Iburi earthquake (Mj 6.7, Mw 6.6), which occurred on the 6th of September 2018, triggered thousands of landslides in pyroclastic fall deposits. The triggered landslides caused destructive damages to the structures and resulted in serious causalities. Hundreds of earthquakes persisted after the mainshock and the Iburi region is of high possibility to suffer major earthquake hitherto. An effective method to assess the seismic slope stability is of great importance for the disaster prevention and mitigation in the Iburi region. The aim of this work is to utilize the Scoops3D to conduct the slope stability analysis in Atsuma, Hokkaido, where destructive landslides were resulted during the Iburi earthquake. The high-resolution DEM (5 × 5 m) was used to construct the surface topography and the isopachs of pyroclastic fall deposits were used to construct the subsurface structures. Then a series of horizontal pseudo-acceleration coefficients proposed in previous literatures were selected to compute the slope stability under seismic loading and the results of the calculation were validated based on the landslides triggered by the Iburi earthquake.

Study area and coseismic landslides

The 2018 Hokkaido Eastern Iburi earthquake triggered thousands of landslides in the vicinity of Atsuma, Hokkaido, Japan. In order to evaluate the seismic stability of slopes on a regional scale, four towns in western Atsuma (Fig. 1c, Tomisato, Yoshinoya, Sakuraoka and Horosato), where catastrophic landslides occurred, were selected as this study area. The study area is approximately 12 km away from the epicenter of the Iburi earthquake mainshock (Fig. 1b). The study area is characterized by moderate terrain with the elevation range of 18 m to 244 m (Fig. 2a) and the predominant slope angle below 40° (Fig. 2b). As to the slope aspect of the study area, there is no apparent and preferred inclinations (Fig. 2c). The majority of study area is underlain by Neogene sedimentary rock (sandstone, siltstone, mudstone and conglomerate) and details of the bedrock classification in the study area are shown in Fig. 2d and Table 1.

Fig. 1

Location maps and coseismic landslides in the study area (a) General view of the study area. b Realtime ground-motion monitoring stations (K-NET station and KiK-net station) distributed around the study area. c Coseismic landslides triggered by the Iburi earthquake in the study area. The shaded relief map was generated from a 5 × 5 m DEM

Fig. 2

Topographic and geological information of the study area (a) Elevation. b Slope angle. c Slope aspect. d Geological map (based on the 1:200,000 seamless geological map published by the Geological Survey of Japan, AIST). The detailed explanations of the geological codes (N2sn, N3sn, Hsr, Q2th, N1 sr, Q2sr and Q3tl) are listed in Table 1

Table 1 Classification of geological units in the study area

The Iburi earthquake caused enormous economic losses and resulted in 41 casualties (based on reports by the Ministry of Internal affairs and Communications, Japan). During the Iburi earthquake sequence, the highest ground motion was recorded by a K-NET (Kyoshin network) station (HKD127, Fig. 1b) in Abira and the peak ground acceleration (PGA) of this station was approximately 1.83 g (1,796 gal). The intense ground shaking resulted in severe liquefaction and catastrophic damages to houses, roads and other engineering structures. As is revealed in Fig. 3, road was severely destroyed and jetted sand, uplift of inlet well and surface ruptures were observed in the paddy field in Yoshinoya. Moreover, thousands of slope failures were triggered due to the strong Iburi earthquake and 36 persons were killed by the destructive slope failures (Yamagishi and Yamazaki 2018). The damage of the slope failures is more destructive than that resulted directly from the earthquake. The study area includes 345 landslides triggered by the Iburi earthquake, and all of these landslides were shallow landslides with long run-out distance and high mobility. Figure 4 shows two most destructive landslides, i.e., the Yoshinoya landslide and the Tomisato landslide in the study area. These two landslides destroyed and buried dozens of houses and other engineering structures at the lower footslope.

Fig. 3

Damages and liquefactions resulted from the Iburi earthquake sequence. a Liquefaction and jetted sand in the paddy field. b Road damage. c Uplift of inlet well. d Surface rupture in the paddy field

Fig. 4

Panoramic views of the destructive Yoshinoya landslide and Tomisato landslide (the base maps are from Google Earth). a Yoshinoya landslide before failure. b Yoshinoya landslide after failure. c Tomisato landslide before failure. d Tomisato landslide after failure

Material and methods

The Fortran program Scoops3D, developed by the U.S. Geological Survey has advantages in seismic slope stability calculation on a regional scale, as it can compute the factor of safety of millions of potential landslides based on a predisposed DEM utilizing 3D method of columns approach (Reid et al. 2015). In this work, the high-resolution DEM was served as the ground surface and the isopachs of pyroclastic fall deposits around the study area were used to define the depths of soil layers. Then a series of horizontal pseudo-acceleration coefficients proposed in previous studies were applied in seismic slope stability calculation of the study area.

Soil structure and geotechnical parameters

The study area is mantled with late Pleistocene and Holocene pyroclastic fall deposits (Hirose et al. 2018) originated from the mounts Tarumae and Eniwa, and the Shikotsu caldera. According to the eruption ages of the craters, the surface pumice layers around the study area can be classified as six layers, i.e., Tarumae-a (Ta-a, AD1739), Tarumae-b (Ta-b, AD1667), Tarumae-c (Ta-c, 2.5–3 ka), Tarumae-d (Ta-d, 8–9 ka), Eniwa-a (En-a, 20 ka) and Shikotsu-1 (Spfa-1, 40 ka) (Hirose et al. 2018; Tajika et al. 2016). The thicknesses of these six pyroclastic fall deposits in the vicinity of the study area were illustrated in Fig. 5. The thicknesses of Ta-a, Ta-b, Ta-c, Ta-d, En-a and Spfa-1 pyroclastic fall deposits are 0~0.16 m, 0~0.32 m,0.32 m, 0.3~0.5 m,0.5~1.0 m and 4 m respectively. Considering that the formation age of Ta-a is close to that of Ta-b and the thicknesses of these two layers are very thin compared with other pyroclastic fall deposits, the Ta-a and Ta-b layers are regarded as one layer, Ta-a, b, in this work. During the computation of slope stability analysis, boundary condition is a major concern. To reduce the effect irregular shape of the targeted area on the computation, the rectangle encompassing the study area is chosen as the computation area. Based on the isopachs of pyroclastic fall deposits, the rectangular computation area can be divided into 6 sub-areas. The soil structures of the corresponding sub-areas are depicted in Fig. 6.

Fig. 5

Isopachs of pyroclastic fall deposits (Ta-a, Ta-b, Ta-c, Ta-d, En-a and Spfa-1) in the vicinity of the study area based on Machida and Arai (2003), Furukawa and Nakagawa (2010) and Hirose et al. (2018). a Isopachs of Ta-a layer. b Isopach of Ta-b layer. c Isopach of Ta-c layer. d Isopach of Ta-d layer. e Isopach of En-a laer. f Isopach of Spfa-1 layer

Fig. 6

Soil layers and corresponding depths of the six sub-areas classified based on the isopachs of pyroclastic fall deposits

On the basis of our on-site field reconnaissance from 10 September 2018 and the preliminary report (Hirose et al. 2018), the sliding zone of the majority of the Iburi landslide is located in the Ta-d layer. Geotechnical parameters of the pyroclastic fall deposits in the study area are listed in Table 2. The basic parameters (bulk weight, saturated unit weight, water content, degree of saturation and void ratio) of Ta-a, b, Ta-c, En-a and Spfa-1 layers are based on previous literature (Miura et al. 2003; Miura and Yagi 2005) and the basic parameters of Ta-d layer are based on samples taken from Tomisato during the field work. The effective strength parameters (effective cohesion and effective friction angle) of En-a and Spfa-1 layers are from the consolidated undrained triaxial test (Miura and Yagi 2005) and the strength parameters (cohesion and friction angle) of Ta-a, b, Ta-c and Ta-d layers utilized in the work are the empirical value from the report of Chuo Kaihatsu Corporation.

Table 2 Geotechnical parameters of the pyroclastic fall deposits in the study area

Seismic loading

Earthquake is the predominate factor controlling the occurrence of the Iburi landslides. In order to incorporate the effect of seismic loading in the slope stability computation, the ground acceleration of a distal K-NET station (HKD128) was described in the work. K-NET is a nation-wide strong-motion monitoring network operated by the National Research Institute for Earth Science and Disaster Prevention (NIED) and can send strong-motion data on the Internet. More than 1,000 K-NET stations are distributed uniformly in Japan with a distance of 20 km and the seismograph installed on the ground surface will record the ground accelerations of three directions (east-west (EW), north-south (NS) and up-down (UD)) at each station (Kinoshita 1998). Figure 7 shows the monitored ground accelerations of three orthogonal directions (EW, NS and UD) observed by a distal K-NET station (HKD128) during the mainshock of the Iburi earthquake. The HKD128 station is about 6 km west to the study area (Fig. 1). The maximum ground accelerations in EW, NS and UD directions are 0.681 g, 0.567 g and 0.404 g respectively (Fig. 7). The corresponding time of the maximum ground accelerations in three directions emerges at 26.3 s, 26.5 s and 24.3 s respectively, which indicates the propagation velocity difference of the S wave and P wave.

Fig. 7

Ground accelerations of three orthogonal directions (EW, NS and UD) observed by a K-NET station (HKD128) during the Iburi earthquake. a Ground acceleration in EW direction. b Ground acceleration in NS direction. c Ground acceleration in UD (up and down) direction. The ground acceleration data are derived from the National Research Institute for Earth Science and Disaster Prevention (NIED)

In seismic assessment of slope stability, it is difficult to consider the seismic effect in three directions and the synthetic ground acceleration in horizontal direction is usually employed to simulate the effect of seismic loading. Figure 8 displays the synthetic ground accelerations of two dimensions (EW and NS) and three dimensions (EW, NS and UD). The peak horizontal ground acceleration (PHGA) and the total peak ground acceleration (PGA) are observed at 26.3 s and 26.5 s respectively. The corresponding values of PHGA and PGA are 0.715 g and 0.730 g. As the horizontal seismic loading plays more important role in triggering slope failures and value of PHGA is close to the value of PGA, only horizontal seismic loading was considered in the calculation and the value of PHGA was applied in this work.

Fig. 8

Synthetic ground accelerations of two dimensions (EW and NS) and three dimensions (EW, NS and UD) observed by a K-NET station (HKD128) in Iburi earthquake. a Horizontal ground acceleration (synthesis of EW and NS directions). b Total ground acceleration (synthesis of EW, NS and UD directions). The ground acceleration data are derived from the National Research Institute for Earth Science and Disaster Prevention (NIED)

Three-dimensional slope stability analysis

Scoops3D computes the stability of rigid masses encompassed by the spherical trial surfaces (potential sliding surfaces) and assesses slope stability by extending conventional limit-equilibrium analysis to three dimensions. It uses 3D method-of-columns approach to compute the factor of safety in a 3D domain (Reid et al. 2015). By defining a series of rotational center and gradually increased radii, numerous potential sliding surfaces intersected by the spheres and DEM are generated. Then the factor of safety of each potential sliding mass (grid-based) is calculated by either Bishop’s simplified method or the ordinary method (Reid et al. 2015). As Bishop’s simplified method provides more accurate solutions of factor of safety (F_S) and the value of F_S is similar to some rigorous limit equilibrium methods (Spencer 1967; Hungr 1987; Lam and Fredlund 1993; Reid et al. 2015; Tran et al. 2018), the slope stability calculation of this work is conducted utilizing Bishop’s simplified method.

In moment equilibrium methods, the factor of safety can be expressed as the ratio of shear strength (τ) to the shear stress (s) (Reid et al. 2015). Thus, the factor of safety is expressed in Eq. 1. Figure 9 shows the 3D view of one column intersected by the potential sliding surface and DEM as well as the forces acting on this column. The shear strength (τ) is calculated based on the Coulomb-Mohr failure criterion in Eq. 2.

$$ {F}_s=\frac{\tau }{s} $$

(1)

Fig. 9

Schematic diagram of forces acting on one column of (modified from Reid et al. (2015))

$$ \tau ={c}^{\prime }+\left({\sigma}_n-u\right)\tan {\phi}^{\prime } $$

(2)

where c^′ is the effective cohesion, ϕ ′ is the effective internal friction angle, σ_n is the normal stress, and u is the pore water pressure. Based on the global moment equilibrium, the total resisting moment is equal to the driving moment. Considering all the columns encompassed by the spherical trial surface, the factor of safety can be expressed in Eq. 3.

$$ F=\frac{\sum {R}_{i,j}\left({c}_{i,j}{A}_{i,j}+\left({N}_{i,j}-{u}_{i,j}{A}_{i,j}\right)\mathit{\tan}{\upphi}_{i,j}\right)}{\sum {W}_{i,j}\left({R}_{i,j}\mathit{\sin}{\upalpha}_{i,j}+{k}_{eq}{e}_{i,j}\right)} $$

(3)

where N_i,j is the normal force of column i,j, W_{i, j} is the weight, A_i,j is the area of the trial surface, R_i,j is the distance from the rotation axis to the geometric center of the potential sliding face of column i,j, e_i,j is the horizontal driving force moment arm, α_i,j is the apparent dip angle and k_eq is the horizontal pseudo-acceleration coefficient. The factor of safety expressed in Eq. 3 is suitable for all general moment equilibrium method, and the normal force can be further obtained by the force equilibrium in vertical and horizontal directions in Bishop’s simplified method. Then by substituting the calculated normal force from the force equilibrium into Eq. 3, the factor of safety is finalized as Eq. 4.

$$ {F}_s=\frac{\sum {R}_{i,j}\left[{c}_{i,j}{A}_{h_{i,j}}+\left({W}_{i,j}-{u}_{i,j}{A}_{h_{i,j}}\right)\mathit{\tan}{\upphi}_{i,j}\right]}{\sum {W}_{i,j}\left({R}_{i,j}\mathit{\sin}{\upalpha}_{i,j}+{k}_{eq}{e}_{i,j}\right)\ }\frac{1}{m_{\alpha_{i,j}}} $$

(4)

where ${A}_{h_{i,j}}$ is the projected area of the trial surface in horizontal plane (${A}_{h_{i,j}}={A}_{i,j}\mathit{\cos}{\varepsilon}_{i,j}$) and ${m}_{\alpha_{i,j}}=\cos {\varepsilon}_{i,j}+\left(\mathit{\sin}{\upalpha}_{i,j}\mathit{\tan}{\upphi}_{i,j}\right)/{F}_s$. As is described in Fig. 9, ε_{i, j} is the true dip angle of the trial surface. More details of the theory and operation of Scoops3D can be accessible in the manual book (Reid et al. 2015) published by the U.S. Geological Survey (https://www.usgs.gov/software/scoops3d).

Input parameters and assumptions

In general, the operation of slope stability calculation in Scoops3D consists of three parts, i.e., construction of 3D domain, input of subsurface parameters and seismic loading, and definition of search configuration. In this work, the 3D domain was layer-based and the ascii DEM (5 × 5 m) was inputted as the ground surface. The base (bottom) elevations of the pyroclastic fall deposits below the ground surface were calculated using the raster calculator tool in ArcGIS based on the isopachs of pyroclastic fall deposits. In order to ensure the accuracy of base elevations of the soil layers, 30 points (5 points in each sub-area) distributed randomly in the calculation area were selected. Then the elevation of each layer at each point was extracted and the corresponding depths of the soil layers of the 30 points were calculated (Table 3). The depths of soil layers of the 30 points show high consistency with depths of soil layers depicted in Fig. 6. Another consideration during the construction of 3D domain is the location of groundwater table. Based on the samples from Ta-d layer in Tomisato, the Ta-d layer is highly saturated but not fully saturated (Table 2). Due to the scarcity of the information about the variation of groundwater in the study area, which is being loosely populated, the ground water configuration in this work was assumed to be a piezometric surface located at the bottom of the Ta-d layer.

Table 3 Elevations and depths of the pyroclastic fall deposits in the six sub-areas

The input geotechnical parameters of soil layers are listed in Table 2. A major concern in this part is the selection of a proper horizontal pseudo-acceleration coefficient for the regional stability assessment. Several plausible horizontal seismic coefficients were proposed in previous literatures. Terzaghi (1950) proposed horizontal seismic coefficient of 0.5 is suitable seismic stability assessment for catastrophic earthquakes and Marcuson and Franklin (1983) noted that a value between 1/2 of PHGA and 1/3 of PHGA could be applied as the horizontal seismic coefficient. Other widely applied horizontal seismic coefficients are described in Table 4 and eight cases with different horizontal seismic coefficients were selected to simulate the seismic loading in the calculation.

Table 4 Eight horizontal pseudo-acceleration coefficients applied in the calculation in eight cases

The search configuration includes size criteria of potential slope failures and definition of search extent. 345 coseismic landslides are distributed in the study area and these landslides can be classified as two types, i.e., coherent shallow debris slide and disrupted mobilization of valley fill. 473 landslide source areas were confirmed after a further division of source area and deposition area. More than 90% of the landslide source areas are between 50 m² to 20,000 m², thus this range is set as the size criteria for potential failures. The horizontal search extent of rotational centers is same with the extent of the study area. The vertical extent of rotational centers ranges from 20 m to 500 m.

Results and discussion

In general, the predicted areas with F_S less than 1.0 are considered as unstable and the predicted areas with F_S more than 1.0 are regarded as stable area. In order to evaluate the performance of the calculation results, a more detailed classification of slope stability and instability was introduced (Table 5).The areas with F_S less than 0.75 and between 0.75 and 1.0 indicate the very unstable and unstable areas. In this work, the calculation of case 1 without seismic loading was conducted to simulate the state before the Iburi earthquake occurrence. Case 2 to case 8 with the increasing horizontal seismic coefficients were conducted to check the performance of the calculated results. By validating using the coseismic landslide inventory, a proper horizontal seismic coefficient range suitable for slope stability assessment can be obtained.

Table 5 Classifications of slope stability and instability based on Ray and De Smedt (2009) and Teixeira et al. (2015)

The factor-of-safety map calculated by Scoops3D is grid-based, thus, GIS can be applied to validate the calculation result. Figure 10 shows the factor-of-safety maps calculated with different horizontal pseudo-acceleration coefficients. In case 1, the factor-of-safety map indicates that the whole study area is in stable state without seismic loading and only very limited areas, i.e., the Tomisato sliding area and the eastern boundary of the study area, have high possibility for slope failure (Fig. 10a). Case 2 to case 8 illustrate that slope instability increases with the increase of horizontal pseudo-acceleration coefficients. A horizontal pseudo-acceleration coefficient of 0.5 (case 8) will result in that more than half of the whole study area is in unstable state. In order to quantitively evaluate to calculation results, landslide points distributed in each stability classes were extracted and areas of slope stability classes were calculated for eight cases. Figure 11 displays the cumulative landslide point percentage of the slope stability classes. For case 1, no landslide is distributed in the very unstable class and 99.8% of the coseismic landslides are distribute in the stable classes (quasi stable, stable and very stable). This is in well accordance with the real case. With the horizontal pseudo-acceleration coefficients of between 0.1 and 0.238 (case 2 to case 5), 6.1% to 22.2% of the Iburi landslides are spread in the unstable classes (very unstable and unstable), which is far from the satisfactory results. The cumulative landslide point percentage of unstable classes is 53.5% with the horizontal pseudo-acceleration coefficient of 0.358 (1/2 of PHGA) and the corresponding cumulative class area percentage is 31.5% with same seismic coefficient. In case 7 and 8, the calculation results seem to be perfect as 75.7% and 80.6% of the triggered landslides are located in the unstable classes. However, the cumulative percentage of the unstable classes in case 7 and case 8 are as high 49.1% and 53.4%. Thus, the seismic coefficients applied in case 7 and case 8 (especially in case 8) may result in overestimation of the slope instability.

Fig. 10

Factor-of-safety maps calculated with different horizontal pseudo-acceleration coefficients

Fig. 11

Cumulative landslide percentages of five stability classes for eight cases

The sizes of Iburi landslides in the study area range from 77 m² to 115,214 m², therefore the landslide area should be also considered in the validation of the calculated results. Besides, the aforementioned validation method can not include the false predicted area, i.e., the predicted unstable area in the area where no slope failure occurs during the Iburi earthquake. In view of this, the example contingency table for a two-class problem, also termed as an error matrix (Stehman 1997) or a confusion matrix, were utilized to validate the calculation results. Even though the factor-of-safety map calculated by Scoops3D is for the whole potential sliding area, each grid cell of factor-of-safety map has a unique value representing the factor of safety and all grid cells can be determined as one of the four elements (true positive, false negative, false positive or true negative) in Fig. 12. The grid cell computed to be unstable is counted as true positive if it falls in real landsliding area and it is counted as false positive if it falls outside the real landsliding area; correspondingly, the grid cell calculated as stable is counted as a true negative if it is located outside the triggered landslides and it is counted as a false negative if it is located within the triggered landslides (Godt et al. 2008).

Fig. 12

Schematic diagram of the confusion matrix (modified from Fawcett (2006))

The example contingency table of eight cases in this work is listed in Table 6, and the values in the column TP, FN, TN, P, N and P + N are the grid cell counts of corresponding cases. Similar with Fig. 10a, the TPR for case 1 without seismic loading is as low as 0.002, which means very few landslides are correctly predicted. The TPRs of cases 2 to 8 increase from 0.07 to 0.87 as the seismic coefficients increase from 0.1 to 0.5, but the FPRs also increase from 0.02 to 0.49. The high FPR of case 8 indicates almost half of the area without landslide occurrence during the Iburi earthquake was predicted as unstable area. In general, a perfect result should meet the requirements of maximum TPR and minimum FPR when the TPR/FPR is more than 1 (Fawcett 2006). However, in most models a higher TPR is accompanied by a higher FPR. The perfect predicted results should incorporate high veracity and acceptable error rate. The TPR of case 6 is approximately same to the predicted results using Scoops3D in loess area, while the FPR is much lower than that predicted in loess area (Xin and Zhang 2018). The horizontal seismic coefficient of between 1/2 PHGA to 2/3 PHGA seems to be effective to obtain good results in pyroclastic fall deposits under seismic loading.

Table 6 Calculation results of eight cases with or without seismic loading

The occurrence of coseismic landslides in area with relatively low slope angle is one of the typical characteristics during the Iburi earthquake and the existence of (act) faults also play an important role in the occurrence of the Iburi landslides. More than 50 landslides are distributed in the flat west part of the study area and these landslides are not well predicted even with a horizontal seismic coefficient of 0.5. The occurrence of these landslides in flat topography decreases the predicted accuracy. Even so, a horizontal seismic coefficient between 1/2 PHGA to 2/3 PHGA can correctly predict 53.5% to 75.7% of the landslide points, and correctly predict 63.5% to 82.7% of the total landsliding area (represented by the TPR).

Conclusions

In this work, GIS was applied to process the input the ground surface, soil layers and the output factor-of-safety maps, and Scoops3D was employed to conduct the three-dimensional seismic slope stability calculation. The soil structure of the study area was constructed based on the isopachs of the pyroclastic fall deposits and the groundwater configuration was assumed as a piezometric surface located at the bottom of the Ta-d layer. Regional slope stability assessment of eight cases with different horizontal pseudo-acceleration coefficients were conducted to simulate the states with or without a seismic loading. 345 coseimic landslides in the study occurred during the Iburi earthquake were utilized to validate the accuracy of the calculated results.

The calculated result without seismic loading in case 1 shows high consistency with the real case and more than 99.8% of the study area is in stable state. The horizontal seismic coefficients of 0.1 to 0.238 applied in the slope stability computation result in 6.1% to 22.2% of total landslide number and the corresponding TPR and FPR range from 0.07 to 0.3 and 0.02 to 0.1. Obviously, these calculation results are far from satisfaction. The predicted result with a seismic coefficient of 0.5 manifests 80.6% of the triggered landslides with landslide area percentage of 86.6% are correctly predicted. But this predicted result overestimate the unstable area as the FPR value of this case is as high as 0.49. Based on the modelled results of case 6 and case 7, a horizontal pseudo-acceleration coefficient between 0.358 (1/2 of PHGA) and 0.477 (2/3 of PHGA) may be suitable for slope stability evaluation since 53.5% to 75.7% of the total landslide (landslide point) can be predicted and 63.5% to 82.7% of the total sliding area can be modelled. The computed result a horizontal pseudo-acceleration coefficient between 1/2 and 2/3 of PHGA is satisfactory, especially in hilly regions. Besides, the factor-of-safety maps in Figs. 10f, g indicate that two most destructive landslides, i.e., the Tomisato landslide and the Yoshinoya landslide are correctly predicted.