Evaluating earthquake-induced widespread slope failure hazards using an AHP-GIS combination

Landslides and slope failures are often caused by earthquakes. This study proposes a method to map earthquake-induced slope failure hazards that uses the analytic hierarchy process (AHP) and a geographic information system (GIS) for four districts where many slope failures were induced by earthquakes (the 2018 Hokkaido Eastern Iburi, 2016 Kumamoto, 2008 Iwate-Miyagi Nairiku, and 2004 Mid Niigata Prefecture earthquakes). The assessment system, which was based on the National Research Institute for Earth Science and Disaster Resilience landslide distribution maps, was analyzed using the methods of previously published. We considered the relationships between the earthquake-induced slope failure distributions and landslide hazard factors (elevation, slope angle, slope type, catchment degree, geology, and vegetation). These relationships were utilized for pairwise comparisons of the factors in the AHP analysis. The slope angle, slope type, and catchment degree exerted the highest effects on the slope failure distribution in the four districts. The four earthquake-induced slope failure distributions were highly consistent with the slope failure hazard rank. These results provide a practical method for evaluating earthquake-induced slope-failure hazards.


Introduction
Japan is among the most earthquake-and landslide-prone countries worldwide. The geology and topography of the Pacific Ring of Fire, which includes the Japanese archipelago, is rich in magmatic activity, crustal movements, and weathering; moreover, fragile geology and steep topography are widely distributed. Therefore, many slope disasters triggered by earthquakes have occurred in these regions. For example, we remember the recent 2018 Hokkaido Eastern Iburi earthquake ( Fig. 1A and Table 1; epicenter: central-eastern Iburi region, Hokkaido, earthquake magnitude: 6.7, epicenter The 2018 Hokkaido Eastern Iburi Earthquake *1 Mj 6.7 *1 About 6,000 places *2 Tephra *2 The 2016 Kumamoto Earthquake *3 Mj 7.3 *3 190 or more places *4 Tephra *4 The Iwate-Miyagi Nairiku Earthquake in 2008 *5 Mj 7.2 *5 2000 or more places *6 Tephra *6 Mid Niigata Prefecture Earthquake in 2004 *7 Mj 6.8 *7 About 4,000 places *8 Neogene sediments *8 depth: 37 km, maximum seismic intensity 7) (Headquarters for Earthquake Research Promotion 2018; JMA 2018; AIST 2013), which occurred at 03:07:59 a.m. on September 6, 2018, which caused widely and densely distributed slope failures (GSI 2018;Murakami et al. 2019;Umeda et al. 2019) ( Fig. 2A and Table 1). Immediately after the occurrence of this earthquake, academic research teams from various fields visited the site and performed surveys and field observations of the distribution and damage caused by these slope failures (Yamagishi and Yamazaki 2018;Hirose et al. 2018;Osanai et al. 2019;JGS 2019;JSCE 2019;Tajika et al. 2020;Ishikawa et al. 2021). In addition, many researchers and engineers have conducted surveys and studies Umeda et al. 2019;Kawamura et al. 2019;Kohno et al. 2020a;Li et al. 2020;Wang et al. 2021) to investigate the causes of slope failures and elucidate the mechanisms of their occurrence. In such cases, it is important to clarify the mechanism of the slope failure and thereby establish effective methods, such as risk assessment for slope failure, prediction of the reach of collapsed sediment using numerical analysis, priority evaluation of disaster prevention measures, and others. As shown by the slope failure cases from the 2016 Kumamoto ( Fig. 1B   . The location of these areas is shown in Fig. 5 the 2008 Iwate -Miyagi Nairiku, and 2004 Mid Niigata Prefecture earthquakes, most slope failures occurred in areas with a Japan Meteorological Agency (JMA) seismic intensity between 6+ and 6−. In the 2016 Kumamoto earthquake, most of the slope failures occurred in areas with a JMA seismic intensity between 6− and 5+ . Although the damage caused by slope failures due to earthquakes is enormous, it is financially difficult to take preventive measures in all slope-failure risk areas. Therefore, it is necessary to identify slopes with a high risk of slope failure during an earthquake and considering this priority, conduct surveys, implement countermeasures, and perform slope stability analysis. In such cases, one of the first items to consider is topography. In recent years, aerial video surveys using drones (UAVs) have attracted attention and demonstrated their benefits in terms of time and cost when conducting topographical surveys over a wide area (Niethammer et al. 2012;Lucieer et al. 2014;Tanaka et al. 2016;Yamamura et al. 2016;Peternel et al. 2017;Yamasaki 2017;Cui et al. 2020;Ma et al. 2020). However, there are many hurdles to overcome, such as aviation law, radio law, and ordinances. Therefore, as an alternative approach, identifying areas of slope failure hazards from earthquakes objectively using easily available and versatile data can provide useful basic or auxiliary data.
In this study, we first applied the landslide hazard mapping method proposed by Kohno et al. (2020b), which uses an AHP analytic hierarchy process (AHP) analysis and a geographic information system (GIS), to the slope failure distributions (Fig. 2) caused by the above four earthquakes (the 2018Hokkaido Eastern Iburi, 2016Kumamoto, 2008Iwate-Miyagi Nairiku, 2004. Subsequently, the consistency between the landslide hazard rank and slope failure distribution due to the earthquakes was examined. The slope failure due to earthquakes, which was targeted in this study, has different causes and processes compared to landslides. However, a landslide hazard map was first created using the landslide topography distribution, which was then modified by changing the evaluation factors concerning slope failure to map earthquake-induced slope fail hazards. In addition, after confirming the degree of influence of topography, geology, vegetation, etc., on the distribution of earthquake-induced slope failure, we explored the common slope failure evaluation factors among the four earthquakes. Based on the above results, we developed a method to evaluate earthquake-induced widespread slope failure hazards caused by earthquakes using an AHP-GIS combination analysis.

Overview of the landslide hazard mapping method by using AHP analysis and GIS proposed by Kohno et al. (2020b)
This section briefly describes the landslide hazard mapping method, using AHP analysis and GIS, proposed by Kohno et al. (2020b). AHP analysis is a decision-making method developed by Saaty (1980) that assigns weights to evaluation factors based on paired comparisons. In the AHP method, the importance weight W i of each element is generally calculated using the n-th order square matrix A (= [a ij ] i,j = 1, 2,…,n ; paired comparison matrix) [Eq.
(1)] of n elements in the evaluation item. (1) where a ij > 0, a ii = 1, a ji = 1/a ij . The calculation of the weight of each element includes the application of the main eigenvector method, the geometric mean method, and the harmonic mean method. Of these, the geometric mean method is the simplest and the most commonly used method. Furthermore, Kohno et al. (2020b) implemented the geometric mean method to calculate the weights of the elements. The importance weight W i of each element on the i-th row using the geometric mean method can be calculated by Eq. (2).
where R.I. is the random index that represents the average of random consistency index.
Many previous studies on landslide hazard assessment have utilized AHP analyses (e.g., Miyagi et al. 2004;Yoshimatsu and Abe 2006;Yagi et al. 2009;Hasekioğulları and Ercanoglu 2012;Pourghasemi et al. 2012; Geological Survey of Hokkaido 2013; Mansouri Daneshvar 2014; Roodposhti et al. 2014;Myronidis et al. 2016;Bera et al. 2020). However, the evaluation targets of each study (Miyagi et al. 2004;Yoshimatsu and Abe 2006;Yagi et al. 2009; Geological Survey of Hokkaido 2013) were limited to landslide topography. Thus, these methods primarily focused on reactive landslides. They were also established through brainstorming with several skilled engineers who had conducted landslide surveys for many years. Therefore, they depended largely on the subjective and empirical judgments of the evaluators. Furthermore, the evaluation factors were subdivided, thus requiring advanced technology and additional time. Hence, Kohno et al. (2020b) proposed a new landslide hazard mapping method using AHP analysis. A Frow chart of this method is shown in Fig. 3. In this method, data that are easy to obtain and highly versatile are set as the evaluation factors. In the paired comparison of evaluation factors related to the landslide hazard by AHP analysis, the relationship between the landslide topography distribution (NIED 2014) and the evaluation factors was quantified (ratio between elements in the evaluation factor).
This method first investigates the relationship between the landslide topography and evaluation factors/elements (i.e., the area ratio of both) based on the landslide topography distribution map (NIED 2014) in the target area. The six evaluation factors (hierarchy level 1) related to the landslide hazard by the AHP analysis were elevation, slope angle, slope type, catchment degree, geology, and vegetation, with an emphasis on easy data acquisition and high versatility. A landslide hazard hierarchy system was constructed based on these evaluation factors (Fig. 3A). From this hierarchy system, the importance weight between each element (hierarchy level 2) was first calculated, and then the importance weights among the higher evaluation factors (hierarchy level 1) were calculated. In general, in AHP analysis, the standard scale for paired comparisons when calculating the importance weight of an evaluation factor (or element) is as follows: "equal importance" is 1; "moderate importance" is 3; "strong importance" is 5; "very strong importance" is 7; and "absolute importance" is 9 (2, 4, 6 and 8 are used complementarily). If the importance is low, then the reciprocal is used. Thus, this standard scale inevitably incorporates the subjective and empirical judgments of the engineer. Therefore, in the method by Kohno et al. (2020b), the importance weight was calculated by introducing the quantified relationship between the landslide topography distribution and the elements in the evaluation factor into a paired comparison (Fig. 3B). The landslide hazard level for each evaluation factor was scored based on the obtained importance weights of each layer (Fig. 3C). Landslide hazard mapping was then performed by totaling the corresponding landslide hazard scores for each evaluation factor using GIS (Fig. 3D, E and Fig. 4). This method finally investigates the relationship between the landslide topography and landslide hazard rank based on an obtained landslide hazard map (Fig. 3F). According to the landslide hazard map obtained using this method, the landslide topography distribution ratio tended to increase as the landslide hazard increased. Hence, this method can be expected to extract unstable slopes and can contribute to the development of efficient measures, assigning priorities, and slope disaster prevention technology (Kohno et al. 2020b). Therefore, in this study, we proceeded with the study of the earthquake-induced slope failure hazard evaluation method based on the method by Kohno et al. (2020b). That is, the method by Kohno et al. (2020b) was applied to the slope failure area caused by four earthquakes, and subsequent landslide hazard mapping was conducted (

Choosing the research area in this study
In this study, the study area ( Fig. 5) was set as the northern Kamuikotan metamorphic, South Kitakami, Joetsu, and Ryoke belt districts, with reference to the geotectonic subdivisions of the Japanese archipelago (AIST 2001;Isozaki et al. 2010). Each of these study areas includes a slope failure distribution area due to one of the above earthquakes. The Kamuikotan metamorphic belt district runs from north to south in central Hokkaido and is a metamorphic belt mainly composed of crystalline schist (Kato et al. 1990). It also Fig. 4 A schematic of the thematic layer overlay analysis using a geographic information system (GIS) to derive landslide slope hazard map 1 3 has the largest area of serpentinite distribution in Japan, and many landslides and deepseated landslides have resulted from this serpentinite. The South Kitakami belt district is distributed from the southwestern part of Iwate Prefecture to the northeastern part of Miyagi Prefecture and consists of Pre-Silurian basement rocks (mainly granites) and shallow marine sedimentary rocks (Oide et al. 1989). The Joetsu belt district is distributed from the northwestern part of Gunma Prefecture to the southeastern part of Niigata Prefecture and consists of ultramafic rocks, crystalline schists, hornfelsized schists, and granites (Omori et al. 1986). As shown in Fig. 2D, this area is one of the most landslideprone regions in Japan. The Ryoke belt district is in the central-eastern part of Kyushu and consists mainly of Ryoke granite and localized Ryoke metamorphic rocks (Karakida et al. 1992).

Overview of data used in this study
The GIS software used was ArcGIS 10.2.2, from ESRI, Japan. In this study, raster data were used to input and superimpose the landslide hazard scores of each evaluation factor calculated by AHP analysis into each cell of the GIS data model (Fig. 4). The cell size of the raster data was set to 20 m, 30 m, and 40 m, and the degree of influence of the cell size on the analysis results was also examined. The landslide topography distribution map used the digital archive for landslide distribution maps (SHP format) published by the NIED (2014). This map is a drawing wherein 1:40,000 monochrome aerial photographs taken in the 1970s are topographically interpreted using a simplified actual condition mirror at four magnifications. Furthermore, its outline and basic structure (scarp and moving body) are mapped regarding "landslide topography," which is a topographical trace formed by landslide fluctuation, and its distribution is shown in a 1:50,000 topographic map. However, we did not determine whether the mapped landslide topography occurred because of an earthquake, heavy rain, or other factors. This drawing only presents relatively large landslide areas with a width of ≥ 150 m. Therefore, very small slope fluctuations, such as surface collapse, debris flow, and falling rock, and landslide topography measuring less than 150 m in width, are not subject to interpretation, and no field survey was conducted on the interpreted topography. Therefore, the absence of landslide topography on the map does not indicate the absence of a landslide (NIED 2014). GIS data published by the GSI (2018) were used for the slope failure distribution map for the 2018 Hokkaido Eastern Iburi Earthquake, and the GIS data published by the NIED (2016) were used for the slope failure distribution map for the 2016 Kumamoto Earthquakes. GIS data were obtained from the GSI Maps (GSI 2016) for the slope failure distribution map for the Iwate-Miyagi Nairiku Earthquake and Niigata Chuetsu Earthquake. The outline of the slope failures was directly traced from an aerial photograph [GSI Maps (GSI 2016)] and made into polygons (Fig. 6). Table 2 lists the number of confirmed slope failures caused by the earthquakes in each district. The scale of these confirmed landslide topographies and slope failures due varies in size. Therefore, in this study, the landslide or slope failure topography (vector data) was not treated as a single cell; however, the moving body part (polygon data) of the topography was represented by raster data (Figs. 4 and 6). Naturally, jaggies appear in the outline of the terrain as expressed by the raster. The scalp of the landslide (line data) was not included in the analysis because it did not contain area information.  For topographical data, a digital elevation model (10 m mesh, Digital Elevation Model: DEM10B) of the basic map information published by the GSI (2019) was used. Table 2 lists the number of DEM10B files that comprised each target area. Based on DEM10B, a GIS surface analysis was performed, and elevation, slope angle, slope type, and catchment degree maps of each target area were created. The slope type is an undulating form, Suzuki (1997) in which the unit ground surface is three-dimensionally divided by a combination of the tilt angle of the cell and the change state in the tilt direction. The catchment degree was calculated using GIS surface analysis and Terrain Analysis Using DEM (TauDEM), a topographic analysis tool developed by Tarboton (2006).
For the geological data, we used the Seamless Digital Geological Map of Japan (1:200,000) (Basic version, data update date: 29 May 2015) published by the AIST (2015).
For the vegetation data, we used 1/50,000 Natural Environmental Information GIS (maintained from 1979 to 1998) published by the MOE (1998).

Setting the elements (hierarchy level 2)
The above-mentioned six evaluation factors (hierarchy level 1) were classified into each element (hierarchy level 2). The set evaluation factors (elevation, slope angle, slope type, catchment degree, geology, and vegetation) are traditionally considered to be more or less related to each other. However, it is extremely difficult to clarify in detail whether or not these factors are independent of each other when conducting landslide hazard assessments. These factors have often been utilized in previous studies, are closely related to landslides/ slope failures, and represent versatile data. Therefore, these evaluation factors cannot be excluded during landslide hazard mapping. The AHP method assumes that the evaluation factors/elements to be incorporated into the same hierarchy level have high independence. Therefore, in this study, the set evaluation factors were assumed to be independent of each other.
The evaluation factor of elevation was divided into elements from 0 m to 2600 m (the highest altitude in the four areas was 2575.3 m) at intervals of 200 m. The slope angle was divided elements from 0° to 80° at intervals of 10°. The slope type was divided into nine types (convex ridge, rectilinear ridge, concave ridge, convex straight, rectilinear straight, concave straight, convex valley, rectilinear valley, concave valley type slopes) based on the basic classification scheme of Suzuki (1997) (Fig. 7). The catchment degree was divided into five equal classifications. The geological evaluation factors were divided based on the geological age and rock type using the legend of the Seamless Digital Geological Map of Japan Based on these classifications, the proportion of the landslide topography within the divisions of each factor was investigated. Figure 8 shows an example using the relationship 1 3 between the landslide topography distribution and slope angle at a cell size of 20 m. Landslide areas were most prevalent among all areas at a slope inclination of 10.0001°-20°. When the inclination angle was 20.0001° or greater, the landslide area tended to decrease as the slope angle increased. Therefore, a close relationship exists between the landslide distribution and slope angle. The same tendency was observed for cell sizes of 30 m and 40 m.

Paired comparisons between evaluation factors (hierarchy level 1) or elements (hierarchy level 2)
First, the paired comparisons between the elements (hierarchical level 2) related to the evaluation factors is explained here using the slope angle (Table 3) with a cell size of 20 m in the Kamuikotan metamorphic belt district as an example. As mentioned above, the numerical value of the standard scale for paired comparison uses the ratio of the landslide occupancy between the elements. For example, the paired comparison values a 12 in Table 3 for "0°-10°" (4.66%, Fig. 8) and "10.0001°-20°" (13.95%, Fig. 8) are 4.66%/13.95% = 0.33. In contrast, the paired comparison between the "10.0001-20°" and "0°-10°" values is the reciprocal (a 21 = 13.95%/4.66% = 2.99 in Table 3). Based on this idea, paired comparisons between elements were performed in a similar manner to complete the paired comparison matrix (8-by-8 square matrix; Table 3). The importance weights W 1 to W 8 were calculated by the geometric mean method in Eq.
(2). However, since the geometric mean method was used to calculate the importance weights, if no landslide topography was identified ("70.0001°-80°"), the landslide area occupancy in the element was calculated as 1.0 × 10 -10 % instead of 0% in the ratio between elements (e.g., a 18 = 4.66%/1.0 × 10 -10 % = 4.7 × 10 -10 in Table 3). The order of the importance weights (Table 3) (Saaty 1980). In this study, the ratio between the elements to be paired was used; thus, C.I. = 0 and C.R. = 0. Second, a paired comparison between the evaluation factors (hierarchy level 1) was performed (Table 4), and the importance weight was calculated. The coefficient of variation of the importance weight of hierarchy level 2 was used as the standard scale for the paired comparisons between evaluation factors at hierarchy level 1. An example of this is shown in Table 4 for a cell size of 20 m in the Kamuikotan metamorphic belt district. For details on setting the standard scale for hierarchy level 1, please refer to Kohno et al. (2020b). The importance weight (0.3602 in Table 3) with a slope angle of 10.0001°-20° was significantly larger than the other elements (0°-10°, 20.0001°-30° to 70.0001°-80°). This means that a slope angle of 10.0001°-20° is closely related to the distribution of landslides. Therefore, the variation in the importance weights of the elements of the slope angle becomes large owing to the influence of the values of the elements in 10.0001°-20°. In other words, if there is a large variation among the importance weights of the elements at hierarchy level 2, the landslide topography is concentrated and distributed based on a specific element. Therefore, the relationship between the evaluation factor and landslide distribution is strong. AHP analysis is an expert system, and hierarchy level 1 paired comparisons must rely on empirical judgments. However, as described above, it may be possible to adequately express the importance of evaluation factors at hierarchy level 1 by using a variation index, such as the coefficient of variation for the relationship of the importance between elements. Therefore, the proposed method is effective. Hence, Table 3 An example of the pairwise comparisons of evaluation factors, using slope angle, of hierarchy level 2 in the Kamuikotan metamorphic belt district derived using raster data with a cell size of 20 m 0°-10°10.0001°-20°20.0001°-30°30.0001°-40°40.0001°-50°50.0001°-60°60.0001°-70°70.0001°-80°Weight 0°-10°1 among the indexes showing the variation of data, this study uses the coefficient of variation C v to set the numerical values of the standard scale for the paired comparisons between the evaluation factors at hierarchy level 1. In addition, if the importance weights of all elements at hierarchy level 2 are the same, the relationship with the landslide distribution is low, and the coefficient of variation C v is 0. Here, the numerical value (a 21 = 1.50 in Table 4) of the paired comparison between the slope angle and the elevation in Table 4 is explained as an example. First, the coefficient of variation C v of the eight importance weights W 1 -W 8 (0°-10° to 70.0001°-80°) of the slope angle can be calculated as 88.80% from Table 3. Similarly, the coefficient of variation C v of the elevation was 59.30%. The coefficient of variation for the slope angle is larger than that for the elevation; the variation is also large. Therefore, the landslide hazard evaluation is very important. Similar to the paired comparisons at hierarchy level 2, the paired comparison value was set to a 21 = 1.50 in Table 4, using the ratio of the two values (80.80%/59.30%). Based on this idea, paired comparisons between evaluation factors (hierarchy level 1) were performed in a similar manner to complete the paired comparison matrix (6-by-6 square matrix; Table 4). Similar to hierarchy level 2, C.I. = 0 and C.R. = 0.

Relationship between landslide hazard rank and earthquake-induced slope failure distribution
From the importance weights of each level (evaluation factors and elements), the landslide hazard score p element for each element related to each evaluation factor was calculated using the following formula: where W 1 is the importance weight of each evaluation factor at hierarchy level 1, W 2 is the importance weight of each element at hierarchy level 2, and W 2MAX is the maximum value of the importance weight of each element at hierarchy level 2. The landslide hazard score pelement for each element is shown in Table 5. The total of all the landslide hazard scores (overlaid using GIS) for each evaluation factor on a target slope (one cell) was considered the total landslide hazard score P [Eq. (6)].
(5)  Since the score is calculated from the relationship between the element and the landslide topography, using the AHP method, the total landslide hazard score P naturally depends on the accumulated score of all elements. That is, a high total score indicated a high risk of landslides, and a low total score indicated a low risk. Landslide hazard mapping was performed based on the hazard ranks, and as an example, the total landslide hazard score was classified from I to V (out of P = 100 points divided into 5 ranks). Figure 9 shows a landslide hazard map for each district with a cell size of 20 m. The characteristics of these maps are summarized below. In the Kamuikotan metamorphic belt district, hazard ranks III (41.2%) and IV (49.9%) occupied most of the area. Hazard rank V (3.6%) was distributed harmoniously with the Yubari and Teshio Mountains. Hazard ranks II (4.9%) and I (0.4%) were distributed in localized areas throughout the region. In the South Kitakami belt district, hazard rank IV (20.3%) was harmoniously distributed in the Mahiru Mountains, located in the northwestern and southwestern regions of the Kitakami Mountains. Hazard rank V (0.4%) was distributed within a narrow subsection of the hazard rank IV distribution. Hazard ranks I (9.8%) and II (34.6%) were distributed in the Sendai Plain, which occupies the eastern part of the Mahiru Mountains and southern half of the region. In the Joetsu belt district, the hazard rank IV (72.6%) distribution occupied most of Fig. 9 Representative examples of landslide hazard maps created by the AHP model for the Kamuikotan metamorphic, South Kitakami, Joetsu, and Ryoke belt districts the area. The hazard ranks III (8.5%) and V (18.8%) were locally distributed. The hazard rank II area (0.1%) was extremely narrow, and hazard rank I was not present. In the Ryoke belt district, hazard ranks IV (61.8%) and III (31.2%) occupied most of the area. Hazard ranks V (2.0%) and II (5.0%) were distributed within a narrow range. Hazard rank I was not present. Figure 10 shows the relationships between the landslide hazard rank and landslide distribution. In all districts, landslide areas tended to increase as the hazard rank increased. Therefore, a good correspondence was found between the landslide hazard map and the landslide distribution. This was consistent with the tendency shown by Kohno et al. (2020b). The landslide hazard map shown here was verified using the ROC curve and AUC values. The ROC curve (red symbols) demonstrating the relationship between the landslide hazard map (Fig. 9) and landslide distribution is shown in Fig. 11. In addition, the AUC values obtained from the ROC curve are shown in Fig. 12. The ROC curve is a diagram wherein the false positive rate (i.e., 1-specificity) is plotted on the horizontal axis and the true positive rate (sensitivity) is plotted on the vertical axis. It is typically used to evaluate the performance of multiple methods. Values closer to the upper left (the point where a low false positive rate and a high true positive rate are achieved at the same time) signify a better model. The AUC value indicates the quality of the method. The value ranges from 0 to 1, with 1 indicating that complete classification is possible, and 0.5 indicating random classification. The AUC was from 0.76 to 0.83 (Fig. 11), being high in all regions. Further, it can be inferred that it is as high as that in previous studies (Hasekioğulları and Ercanoglu 2012;Pourghasemi et al. 2012;Roodposhti et al. 2014;Myronidis et al. 2016;Pereira et al. 2017;Bahrami et al. 2021;Salehpour Jam et al. 2021) on landslide hazard mapping. Therefore, this result is considered to have the ability to classify landslide hazards. In each district, the characteristics of the landslide hazard maps with a 20 m cell size were approximately the same as those of the 30 m and 40 m cell sizes. Furthermore, the relationship with the landslide topography distribution showed the same tendency, as shown in Fig. 10. In addition, there was no significant difference in the AUC values ( Fig. 12) with a change in the cell size. Therefore, when this method is used with a raster data cell size of 20 to 40 m, the actual cell size has almost no effect on the landslide hazard map. Figure 13 shows the relationships between the landslide hazard rank and slope failure distribution due to the four earthquakes considered in this study. The following trends Fig. 10 Percentage of landslide areas in areas with different levels of hazard rank in the landslide hazard maps shown in Fig. 9. The results were derived using raster data with a cell size of 20 m, 30 m, and 40 m, respectively were observed in three of the regions, excluding the Joetsu belt district. In the Kamuikotan metamorphic belt district, the slope failure area due to the earthquake peaked in hazard rank III areas and decreased in hazard rank IV and V areas. In the South Kitakami and Ryoke belt districts, the slope failure peaked in the hazard rank IV areas and decreased in the hazard rank V areas. In contrast, in the Joetsu belt district, good correspondence was found between the landslide hazard rank and slope failure distribution. As mentioned above, this may be because the hazard rank IV areas occupied most of this area. Furthermore, the slope failures caused by the Mid-Niigata Prefecture Earthquake in 2004 were concentrated in the hazard rank V areas. Hence, in the created landslide hazard map, slope failures due to earthquakes did not necessarily occur in the areas where the landslide hazard was high. The ROC in this case is shown by the blue symbols in Fig. 12. All AUC values were lower (AUC = 0.46 to 0.66) than the red symbols except in the case of the Joetsu belt district (AUC = 0.86). In addition, the AUC value was less than 0.5 in the case of the Kamuikotan metamorphic belt district. Therefore, the landslide hazard map created by this method does not properly represent the slope failure distributions caused by the earthquakes.

Extracting the evaluation factors that most affected the earthquake-induced slope failure distribution
Among the evaluation factors for the landslide hazard map created in the previous section, we investigated which factor had the greatest effect on the slope failure distributions due to the earthquakes. In this method, the landslide hazard score was calculated with the importance weight of one evaluation factor set at 99.9, and the importance weight of the other five evaluation factors set at 0.02 at hierarchy level 1 in the AHP analysis. Because there were six evaluation factors, six corresponding landslide hazard maps were also created. We thereby confirmed the relationships between the landslide hazard rank in the landslide hazard maps and the slope failure distributions due to the four earthquakes. This makes it possible to understand which evaluation factors affected the earthquake-induced slope failure distributions. In addition, we examined whether the evaluation factors have common characteristics among the four regions. Based on this knowledge, the evaluation factors were narrowed down and set again, and the landslide hazard mapping was performed again, allowing us to examine the relationships between the factors and slope failure distributions due to the earthquakes.  Figure 14 shows the relationships between the landslide hazard ranks and the slope failure distributions for the landslide hazard maps that were obtained when the importance weight of each evaluation factor (hierarchy level 1) was set to 99.9. In addition, Fig. 15 displays the AUC values obtained from these ROC curves. In all areas, when the importance weight of the slope type or catchment degree was 99.9, most slope failures to occurred in the hazard rank IV or V areas. Therefore, these two evaluation factors are closely related to the slope failure distributions due to the earthquakes. As per these factors, the AUC values were high in all four districts. That is, when the importance of the slope type is 99.9%, AUC = 0.68-0.80, and when the importance of the catchment degree is 99.9%, AUC = 0.60-0.80 (Fig. 15). However, in the Kamuikotan metamorphic and Ryoke belt district, for altitude, geology, and vegetation, 80% or more of slope failures occurred in the relatively low-risk hazard rank I and II areas. As per these factors, the AUC values were not high in all four districts, with several districts presenting AUC values less than 0.5 (Fig. 15). Therefore, the involvement of these three factors in the slope failure distribution is not significant. Generally, there is a close relationship between slope failure and geology/vegetation. However, it is probable that their influence was not discernable because the slope failures occurred in areas with the same geology and vegetation in these cases. Regarding the slope angle, approximately 70% of the slope failure occurred in the landslide hazard rank III area in the Joetsu belt district. However, as shown in Fig. 8, the landslide topography was mostly distributed at a specific slope angle (10.0001°-20°) in all areas. In addition, the AUC values (= 0.63 to 0.74; Fig. 15) were highly similar, in all four districts. Furthermore, examining many studies (e.g., Osanai et al. 2007;Wakai et al. 2008;Sakai et al. 2013;Hamasaki et al. 2015) on slope failure hazard evaluation due to earthquakes reveals many cases in which the slope angle is mentioned as an evaluation factor. Therefore, we determined that it was not appropriate to evaluate the slope failure hazards due to earthquakes without the slope angle and decided to add it to the evaluation factors.
Hence, landslide hazard mapping was performed again using only three evaluation factors (hierarchy level 1), which were the sum of the slope type, catchment degree, and slope angle, as discussed in the following section. Moreover, all three factors (slope angle, slope type, and catchment degree) are topographic primary causes. As mentioned above, these evaluation factors had similarly high AUC values (Fig. 15) in all four districts. Then, the relationship between this landslide hazard map and slope failure distribution due to the earthquakes was investigated.

Landslide hazard mapping based on slope angle, slope type, and catchment degree
Based on the methods described in Sects. 2 and 3, landslide hazard mapping was performed using the three evaluation factors. Because all three evaluation factors listed here are common topographic primary causes, the importance weights of the evaluation factors (hierarchy level 1) were all the same. Figure 16 shows an example of these hazard maps (slope failure distribution due to the 2018 Hokkaido Eastern Iburi Earthquake, cell size: 20 m). Hazard ranks II and III were found on the west side of the Atsuma Fault extending in the NNW-SSE direction, and landslide hazard ranks IV and V were found on the east side. As shown in the figure, most of the slope failures occurred in the landslide hazard rank IV and V areas. Figure 17 shows the relationships between the landslide hazard ranks and slope failure distributions caused by the four earthquakes. In all areas, as the landslide hazard rank increased, the slope failure area tended to increase. The ROC in this case is shown by the green symbols in Fig. 11. Except for the Joetsu belt district, the ROC curve was closer to the upper left than the blue symbols. Furthermore, the AUC value was also higher than the blue symbols. Thus, the landslide hazard maps corresponded well to the slope failure distributions. This is because landslide hazard mapping was performed by collecting and setting only the factors that were closely related to the slope failure distribution. Therefore, evaluating the slope failure hazard due to an earthquake using the slope angle, slope type, and catchment degree is effective. Moreover, this method can be effective for the advance determination of areas with a high risk of slope failure due to sudden phenomena such as earthquakes. In all areas, there was no significant differences in the relationship between the landslide hazard rank and the slope failure distribution due to the earthquake for cell Fig. 16 An example of a landslide hazard map based on slope angle, slope type, and catchment degree, i.e., a part of Kamuikotan metamorphic belt district, which is also shown in Figs. 5 and 9 Fig. 17 Percentage of areas with earthquake-induced slope failures in areas with different levels of hazard rank in the landslide hazard maps created by the AHP model, considering slope angle, slope type, and catchment degree. The results were derived using raster data with a cell size of 20 m, 30 m, and 40 m, respectively sizes from 20 to 40 m. Therefore, setting the cell size to 40 m is better for the computational load times of GIS analysis. In developing countries, the accuracy of the topographical data may be low or it may be difficult to obtain data, such as geology and vegetation. Even in such a case, this method, which only requires DEM data of a relatively large 40 m mesh to perform the evaluation, is both advantageous and effective.

Conclusions
In this study, we applied the landslide hazard mapping method of Kohno et al. (2020b) to widespread and high-density landslide distribution areas caused by four earthquakes (the 2018Hokkaido Eastern Iburi, 2016Kumamoto, 2008Iwate-Miyagi Nairiku, and 2004 Mid Niigata Prefecture Earthquakes) based on several hundreds to tens of thousands of landslide topographies. Then, the relationship between the landslide hazard ranks and slope failure distribution due to the earthquakes was investigated. Furthermore, the degree of influence of topography, geology, vegetation, etc., on the slope failure distribution was confirmed. Subsequently, we explored the commonalities among the evaluation items for the slope failure distribution and proposed a new earthquake-induced slope failure hazard evaluation method based on the results.
When the method of Kohno et al. (2020b) was used, it was unable to properly express the state of the slope failure distributions during an earthquake. However, the primary topographic causes played major roles in the slope failure distribution during an earthquake. Therefore, it is effective to evaluate the slope failure hazard due to an earthquake using three factors: slope angle, slope type, and catchment degree. We concluded that this method is effective the advance determination of areas with a high risk of slope failure due to sudden phenomena such as earthquakes. In addition, this method has the advantage that the analysis can be performed using only DEM data without the need for other data, such as geology and vegetation. Therefore, it can effectively be applied to areas where data are scarce or inaccurate. Hence, this method can be utilized to locate slopes that will be unstable during an earthquake and require a field survey. Because this method enables efficient priority measurements, we think that it can greatly contribute to the development of slope disaster prevention technology. In future, the relationship between seismic motion and slope failure occurrence mechanism, the spatial correlation between earthquake events and slope failure, and time lag of landslide occurrence after earthquake occurrence will be investigated in detail via time series analysis, and our data can be used for these analyses. By incorporating this information, we aim to improve the accuracy of the method in this study. Furthermore, we plan to further improve this method by exploring its applications and potential in areas subject to various restrictions, such as developing countries.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.