Scenario-Based Risk Assessment of Earthquake Disaster Using Slope Displacement , PGA , and Population Density in the Guyuan Region , China

Mega-earthquakes that occur in mountainous areas of densely populated cities are particularly catastrophic, triggering large landslides, destroying more buildings, and usually resulting in significant death tolls. In this paper, earthquake scenarios in the Guyuan Region of China are used as an example to study earthquake disaster risk assessment and a method of assessment is proposed that uses the peak ground acceleration (PGA), landslides triggered by the earthquake, and the effects on the population. The method is used to develop scenarios for earthquake disaster risk assessment along the Haiyuan and Liupanshan Faults for earthquake magnitudes of Ms 7.0, 7.5, 8.0, and 8.5 triggered by one of the two faults. The quantitative earthquake disaster risk maps in the study area were developed by integrating the values of the at-risk elements for the earthquake factor, population, and landslide hazard. According to the model results, the high-hazard zone was mainly located in the severely affected areas along the faults and on the western side of the faults. These results can be useful for emergency preparation planning, response plans, and resource assessment.


Introduction
Earthquakes are considered to be one of the most serious natural hazards worldwide, due to their potential for causing huge loss of life, extensive property damage, and significant economic loss [1][2][3][4].Particularly mega-earthquakes, that occur in the mountainous areas and densely populated cities, are more dangerous because of their potential to trigger large landslides, destroy numerous buildings, and cause more deaths [5,6].A number of mega-earthquakes have occurred in mountainous and densely populated areas in recent decades, including the Wenchuan Earthquake (Ms 8.0) in China, the Nepal Earthquake (Ms 7.5), the Haiyuan Earthquake (Ms 8.5), and the Kashmir Earthquake (Ms 7.6) in Pakistan.Strong earthquakes have occurred in the past and may re-occur in areas such as Vancouver [7] and Montreal [8] in Canada and in the Ningxia Province [9] and the Guanzhong Plain [10] in China.
Since the reoccurrence intervals between mega-earthquakes in most areas are more than 100 years, local populations tend to forget the previous major earthquakes and think it is a safe place, paying no attention to earthquake mitigation.Unfortunately, it is still not possible to control or accurately forecast earthquakes, and so they remain a disaster that cannot be avoided.In China, the primary method used to improve safety and minimize loss and injury is hazard zoning, which provides information on the intensity of previous earthquakes for a given location.However, this method is not sufficient for informed decision-making, especially in mountainous and densely populated areas [11].In addition, most earthquake loss assessment is only conducted several days after the occurrence of an earthquake, ISPRS Int.J. Geo-Inf.2019, 8, 85 2 of 16 which can significantly delay relief supplies.For example, the risk map for the Wenchuan Earthquake was developed 10 days after the event [12].The scenario methods are the acknowledged basis of emergency planning.
One of the most effective ways to reduce the impact of earthquake disaster on people and property is through scenario-based earthquake disaster risk assessment [13,14].Scenario-based earthquake disaster risk assessment can provide good support to local governments for budget planning, determining appropriate levels of relief supply reserves, raising public awareness, allocating human resources for mitigation and disaster management operations, educating the public and professionals on preparation and mitigation, and the prioritization of retrofit applications [15].Thus, scenario-based seismic disaster risk assessment is an important tool for reducing earthquake-induced losses [16].
Because scenario-based seismic disaster risk assessment is an effective method for improving earthquake relief, a variety of analytical tools, including both qualitative and quantitative risk analyses, are now available for the risk assessment of seismic loss [17].Qualitative earthquake disaster risk analyses primarily study the zoning of earthquakes at a regional scale, based on the earthquake magnitude and geological characteristics of the land use [16].Quantitative methods require data of extremely high quality, which make their general use impractical, and so are typically used only for studies of small, specific areas.The more widely used quantitative risk analysis methods benefit from simplified models, using a Hazards US (HAZUS) analysis under earthquake conditions, which can be used for the evaluation of specific sites [18].HAZUS was built by the US Federal Emergency Management Agency (FEMA) [19], based on numerical modeling with geographic information system (GIS) data for assessing the potential losses from earthquakes.While the HAZUS model is widely used in the USA and Canada due to the advantages of quantifying risks [4,18], it requires very detailed data about the potential earthquake sites.Thus, the model is not necessarily useful worldwide and is not promoted in other counties because the required data cannot be obtained easily or is too expensive to acquire [17].The other widely used earthquake hazard assessment method is PSHA (probabilistic seismic hazard assessment), which mainly focuses on the seismic ground motion assessment and gives the scenario seismic intensity map [20].
Scenario-based seismic disaster risk assessment is widely used to quantify the potential social and economic losses from earthquakes and has become a very important tool for local people and the governments [15,16].However, seismic disaster risk assessment is a complex process because it is necessary to integrate several spatial parameters from different locations.These parameters include (1) geology, (2) geomorphology, (3) engineering geology, (4) structural engineering, (5) human behavior, and (6) population density [2].Despite their complexity, seismic risk assessment studies have proven to be a very useful tool for developing emergency preparation plans and for promoting seismic risk mitigation [15,16].

Region of Study
The study area of the "Guyuan Region", located in the south of the Ningxia Province in China, is a tectonically active plate boundary, separating the Qing-Tibet Plate and the Erodos Sub-plate.There are two famous active faults in the area; the Haiyuan Fault and the Liupanshan Fault (Figure 1).Both faults have been generating intense earthquakes that have affected the Guyuan Region for centuries, including the destructive Ms 8.5 earthquake in 1920 [21] and the major Ms 7.0 earthquakes in 1219, 1306, and 1622 [22].The ~200 km long, active, left-lateral Haiyuan Fault is a major geological structure in northern Tibet and is a well-known left-lateral zone at the eastern edge of the Tibetan Plateau [21].On December 16, 1920, an earthquake of Ms 8.5 occurred in Haiyuan County, in the Ningxia Hui Autonomous Region.The earthquake caused 234,117 deaths and many more were injured [23].Several paleoseismic investigations have been conducted on various points along the fault [24], with three major earthquake events extrapolated to have occurred during the periods A.D. 350-986, B.C. 1325-2031, and B.C. 1985-2351 [24].Thus, the time intervals between these events are 1252 ± 318 years, 2346 ± 475 years, and 490 ± 397 years, respectively, indicating considerable variations.Even on the less active central and eastern sections of the Haiyuan Fault, the mean estimated Quaternary slip rate is 12 ± 4 mm/yr [25], while the modern fault slip rate or strain accumulation rate, as derived from GPS and interferometric synthetic aperture radar, suggest a lower rate of ~5-8 mm/yr [26,27].Taking into account these past events and the Haiyuan Fault's slip rate, an analysis by Song et al. [28] suggested that a 1920-type earthquake would likely occur again within 500 years, with a mean interval time of about 420 years.500 years, with a mean interval time of about 420 years.
The Liupanshan Fault is a thrust fault exposed along the eastern foot of the Liupan Mountain and is characterized by the eastward displacement of the Cretaceous Liupan mountain group over the Tertiary red-beds.The Liupanshan Fault is approximately 130 km long, with the strikes at 330°-335°, and dips about 45 o W (Figure 1).The mean estimated left and right slip rate is 1-3 mm/yr, with 0.9 mm/yr in vertical displacement [29].Three major earthquake events have been identified in the most recent 800 years, in 1219 (Ms 7.0), 1306 (Ms 7.0), and 1622 (Ms 7.0) [30].Thus, the time intervals between these events are 100 years and 300 years, respectively.Considering these past events and the Liupanshan Fault's slip rate, the other analysts have suggested that the accumulated energy is enough to trigger an earthquake with a magnitude of 7.0 or higher [29,30].
Thus, the Haiyuan and Liupanshan Faults pose a major seismic threat to the populations on both sides of the faults.Scenario-based seismic risk assessment is essential for the region, as it can provide a reference for emergency response training exercises and the determination of appropriate relief supply reserves.

Methods and Data Collection
It is generally accepted that qualitative risk analysis for earthquakes is preferred over hazard analysis because it allows for a more explicitly objective output for hazard relief [31].The method is based on the widely used definitions of the risk assessment method proposed by Varnes [32], Fell [33], and Lee and Jones [34], where: ( ) ( ) Here, H is the earthquake hazard.The value of (VA) is for all of the elements at risk, which is an index derived from the casualties and from the buildings and infrastructure destroyed in the earthquake [35].Seismic risk assessment can provide estimates of the probability losses induced by an earthquake, including to housing infrastructure, the economy, casualties and societal losses [36,37].The Liupanshan Fault is a thrust fault exposed along the eastern foot of the Liupan Mountain and is characterized by the eastward displacement of the Cretaceous Liupan mountain group over the Tertiary red-beds.The Liupanshan Fault is approximately 130 km long, with the strikes at 330 • -335 • , and dips about 45 • W (Figure 1).The mean estimated left and right slip rate is 1-3 mm/yr, with 0.9 mm/yr in vertical displacement [29].Three major earthquake events have been identified in the most recent 800 years, in 1219 (Ms 7.0), 1306 (Ms 7.0), and 1622 (Ms 7.0) [30].Thus, the time intervals between these events are 100 years and 300 years, respectively.Considering these past events and the Liupanshan Fault's slip rate, the other analysts have suggested that the accumulated energy is enough to trigger an earthquake with a magnitude of 7.0 or higher [29,30].
Thus, the Haiyuan and Liupanshan Faults pose a major seismic threat to the populations on both sides of the faults.Scenario-based seismic risk assessment is essential for the region, as it can provide a reference for emergency response training exercises and the determination of appropriate relief supply reserves.

Methods and Data Collection
It is generally accepted that qualitative risk analysis for earthquakes is preferred over hazard analysis because it allows for a more explicitly objective output for hazard relief [31].The method is based on the widely used definitions of the risk assessment method proposed by Varnes [32], Fell [33], and Lee and Jones [34], where: Here, H is the earthquake hazard.The value of (VA) is for all of the elements at risk, which is an index derived from the casualties and from the buildings and infrastructure destroyed in the earthquake [35].Seismic risk assessment can provide estimates of the probability losses induced by an earthquake, including to housing infrastructure, the economy, casualties and societal losses [36,37].
There are three important factors used to assess an earthquake's impact.The peak ground acceleration (PGA) is a very important factor in engineering seismology, as it influences how much earthquake-affected infrastructure is destroyed and its displacement, and it is also a main factor in determining the magnitude of the earthquake.The PGA is mainly affected by the earthquake's magnitude and the distance to the active fault [38].The second factor, the "population", is a direct index related to the social economics of the region (e.g., the county in the west of China for this work).This variable is also a factor in the relationship between the economy and the amount of infrastructure and the number of buildings in the region [39].Both factors reflect the earthquake-associated hazard and all elements at risk.Finally, earthquake-induced landslides in mountain areas can cause half of all of the casualties in an earthquake, as well as destroying buildings and infrastructure, as seen in recent mega-earthquakes which occurred in mountainous areas [5].Thus, any risk assessment of an earthquake should also consider the potential effects of landslides triggered by this earthquake in a mountainous area.
Therefore, in this paper, we propose a method that uses the PGA of an earthquake and the landslides triggered by an earthquake as the primary factors, along with the population as the element at risk.These three factors can be used to quantify the potential earthquake effects and the social economic elements at risk.More importantly, all three quantitative factors can be obtained easily.
The approach used in this study is shown schematically in Figure 2. The spatial and magnitude probabilities of landslides are estimated using the scenario-based earthquake hazard analysis, according to the stress-strain curve based on triaxial testing.The consequence analysis uses the peak ground acceleration (PGA) for each element at risk, along with the population density in the study unit.The last step is the risk analysis, which involves estimating the values of the elements at risk and then combining these values with the vulnerability and hazard values to produce risk maps for multiple scenarios.
ISPRS Int.J. Geo-Inf.2019, 8, x FOR PEER REVIEW 4 of 16 There are three important factors used to assess an earthquake's impact.The peak ground acceleration (PGA) is a very important factor in engineering seismology, as it influences how much earthquake-affected infrastructure is destroyed and its displacement, and it is also a main factor in determining the magnitude of the earthquake.The PGA is mainly affected by the earthquake's magnitude and the distance to the active fault [38].The second factor, the "population", is a direct index related to the social economics of the region (e.g., the county in the west of China for this work).This variable is also a factor in the relationship between the economy and the amount of infrastructure and the number of buildings in the region [39].Both factors reflect the earthquakeassociated hazard and all elements at risk.Finally, earthquake-induced landslides in mountain areas can cause half of all of the casualties in an earthquake, as well as destroying buildings and infrastructure, as seen in recent mega-earthquakes which occurred in mountainous areas [5].Thus, any risk assessment of an earthquake should also consider the potential effects of landslides triggered by this earthquake in a mountainous area.
Therefore, in this paper, we propose a method that uses the PGA of an earthquake and the landslides triggered by an earthquake as the primary factors, along with the population as the element at risk.These three factors can be used to quantify the potential earthquake effects and the social economic elements at risk.More importantly, all three quantitative factors can be obtained easily.
The approach used in this study is shown schematically in Figure 2. The spatial and magnitude probabilities of landslides are estimated using the scenario-based earthquake hazard analysis, according to the stress-strain curve based on triaxial testing.The consequence analysis uses the peak ground acceleration (PGA) for each element at risk, along with the population density in the study unit.The last step is the risk analysis, which involves estimating the values of the elements at risk and then combining these values with the vulnerability and hazard values to produce risk maps for multiple scenarios.Flow chart showing the steps involved in producing a seismic risk map (c is the cohesive strength, φ is the internal friction angle, r is the material unit weigh, a is the slope angle, Fs is the static factor of safety, and DN is the landslide movement distance in cm.DEM is Digital Elevation Model).

Landslide Assessment
Among the simplified models published for slope stability analysis, the permanentdisplacement analysis developed by Newmark [40] can be used to estimate co-seismic landslide displacements on regional scales for high-magnitude earthquakes and is widely used in co-seismic landslide assessment [41][42][43].Many studies have used the Newmark method and movement records to estimate the earthquake-induced slope displacement in active fault areas by using GIS, such as along the Longmenshan Fault in the Lushan Earthquake [43].These studies have shown that the Newmark model is one of the most accurate models for assessing earthquake-induced landslides and is a suitable method for studying the triggering factors in different areas.

Landslide Assessment
Among the simplified models published for slope stability analysis, the permanent-displacement analysis developed by Newmark [40] can be used to estimate co-seismic landslide displacements on regional scales for high-magnitude earthquakes and is widely used in co-seismic landslide assessment [41][42][43].Many studies have used the Newmark method and movement records to estimate the earthquake-induced slope displacement in active fault areas by using GIS, such as along the Longmenshan Fault in the Lushan Earthquake [43].These studies have shown that the Newmark model is one of the most accurate models for assessing earthquake-induced landslides and is a suitable method for studying the triggering factors in different areas.
Newmark's approach models a landslide as a rigid friction block sliding on an inclined plane.The block has a critical acceleration, which is the threshold acceleration required to overcome the shear resistance; when the critical acceleration (a c ) exceeds the threshold acceleration, downslope movement occurs and initiates a landslide [41].The model calculates the cumulative permanent displacement of the block relative to its base when it is subjected to the effects of earthquake acceleration over time.The critical acceleration is expressed as: where a c is the critical acceleration (g), g is the acceleration due to Earth's gravity.
The displacement calculation is obtained by the quadratic integration of the instances of seismic acceleration over the critical acceleration.Accelerations below this level cause no permanent displacement of the block.Those accelerations exceeding the critical acceleration are integrated to obtain the velocity profile of the block; a second integration is performed to obtain the cumulative displacement of the block.
Researchers have developed several simplified models for estimating the Newmark displacement (D N ) based on different earthquake parameters, such as Newmark displacement models based on the Arias intensity [41], the critical acceleration ratio (i.e., the ratio of the critical acceleration to the maximum acceleration) [41,44], and multivariate parameter synthesis [45,46].Unfortunately, most models developed to date have only used magnitudes of 6.0-7.3, and so are not suitable for accurate displacement prediction when the magnitude is above 7.4.This is particularly true in China, as the data used to build the models does not include Chinese earthquake data.Instead, some researchers built a new method based on data from the Wenchuan and Yushu, Ya'an earthquakes that is more suitable for studying mega-earthquakes in China.The most widely used model in China, suitable for mega-magnitude earthquakes, was developed by Xu et al. [47] using 1253 data points from 420 monitoring stations all over China that collected information during a recent mega-earthquake in China.Using this data in the following equation improves the R 2 value to 91.4%, while the standard deviation decreases slightly.The equation is as follows: Here, a max is the peak ground acceleration, and a c is the critical acceleration.This empirical equation was derived from an extensive dataset of seismic records collected in different regions.The equation is used in the present study as a reasonable and simple computational tool for deriving information about the slope displacement for a given set of geotechnical criteria and assuming a specific earthquake intensity.

Peak Ground Acceleration (PGA)
The value of a max is a very important factor in engineering seismology, as it influences the earthquake-induced landslide displacement.The relationship of a max to the earthquake magnitude and distance was first reported by Gutenberg and Richter in 1942 and 1956.Since 1965, numerous other correlations of the peak acceleration with the earthquake magnitude and distance have been presented, and the regression equation is generally presented in the following form [48,49]: where R is the distance to the fault, M S is the earthquake magnitude, T is 0.35 (for a horizontal earthquake vector), and c 0 , c 1 , c 2 , and c 3 are the regression coefficients for various parameters.Many values of these regression coefficients have been presented for different earthquake areas due to variations in geology and landform.In a loess plateau, Yang et al. [50] obtained the relationship between the peak acceleration and the earthquake magnitude and distance by using 1221 datasets from 49 regional earthquakes.For such regions, the best-fit equation for the regression equation is as follows: lga max = 0.71178 + 0.25898M s + 0.62321lgR + 1.985T (5)

Critical Acceleration
The critical acceleration of a potential landslide block is a simple function of the static factor of safety (F S ), which can be calculated using the following equation proposed by Van Westen and Terlien [51]: where r w is the unit weight of water, t is the slope-normal thickness of the failing slab, and m is the proportion of the slab thickness that is saturated.Since the study region for this work has a semi-arid climate, a negligible pore-water pressure (m = 0) was assumed for this model.For simplicity, the parameters of c, ϕ, and r were derived from the rock hardness category (Table 1).The sliding surface thickness was 6 m, based on the available statistics for the landslides triggered by the Haiyuan Earthquake in the loess plateau.All parameters were obtained from a site category analysis and a DEM and then converted to a grid with the same resolution.The F S value was calculated using Equation (6).The subsequent production of the critical acceleration grid used Equation ( 2) to combine the slope angle with the calculated factors of safety.

Population Data
The population data for the region was obtained from the sixth census data in 2015, which used the township as the smallest administrative unit and also the smallest unit in terms of natural hazard mitigation and assessment.Hence, this study uses the township as the assessment unit for scenario-based seismic risk assessment.

Geological Data and DEM
The site category was identified using 1:5000-scale digital geological maps.The weathering degree determines the hardness of the rock and can be estimated by testing this hardness.The hardness of the rock can be classified into five categories: (1) Very hard, (2) hard, (3) soft-hard, (4) soft, and (5) very soft, according to the Code for the Investigation of Geotechnical Engineering and Testing [52].Representative values for the frictional and cohesive components of the shear strength were assigned to each geological unit.The rock hardness affects the strength of the lithology and the values of the cohesive strength (c) and the internal friction angle (ϕ).We compiled the results from a series of shear tests on samples from a variety of geological units in the study area and the surrounding region.In addition, we queried several papers and reports describing the strength of sites with the same lithology that were similar to those of the study area.The values of c and ϕ for different lithologies are presented in Table 1.For this study, we used the mean values of c, ϕ, and r derived from the strengths of the different reported lithologies.
The DEM, with the grid size of the raster resolution of 25 m, was produced using 1:50,000-scale Quadrangle maps.The slope map and distance to the rupture plane were derived from the DEM, using spatial analysis tools in GIS (Figure 3).The DEM, with the grid size of the raster resolution of 25 m, was produced using 1:50,000-scale Quadrangle maps.The slope map and distance to the rupture plane were derived from the DEM, using spatial analysis tools in GIS (Figure 3).

Results
The Haiyuan Fault and the Liupanshan Fault have previously generated intense earthquake activity, ranging from Ms 5.0 to 8.5.In this study, the seismic magnitude scenarios studied are Ms 7.0, 7.5, 8.0, and 8.5 to provide a risk assessment according to the respective active feature.

Earthquake-Induced Landslides
All of the grid maps were constructed with the same resolution, after which an empirical value of the landslide movement distance (DN) was calculated and used to estimate the regional landslide displacement.This analysis was conducted using the grid calculation tool in ArcGIS to link the slope critical acceleration map and the peak acceleration, using Equation (3).Landslide displacement classes were identified and mapped as follows: 0-1 cm, 1-5 cm, 5-10 cm, 10-15 cm, 15-20 cm, 20-30 cm, and >30 cm (Figure 4).
The maps exhibit consistent levels of DN > 20 cm along the faults and on the western sides of the faults.This is similar to the geological distribution, where loess-covered areas are more susceptible to landslide activity [53][54][55].This behavior is probably due to the geotechnical features of loessmacro-pores, vertical joints, loose texture, and a silt content greater than 50%, which makes it sensitive to earthquakes and prone to landslides [53][54][55].Furthermore, the model predicts that the displacement decreases with an increase in the distance from the fault.The percent of the slope displacement above 20 cm which is triggered by the Haiyuan Fault and the Liupanshan Fault, respectively, increases with the magnitude of the earthquake, from 0.208% and 0.199% at Ms 7.0 to 0.813% and 0.783% at Ms 8.5 (Figure 5).

Results
The Haiyuan Fault and the Liupanshan Fault have previously generated intense earthquake activity, ranging from Ms 5.0 to 8.5.In this study, the seismic magnitude scenarios studied are Ms 7.0, 7.5, 8.0, and 8.5 to provide a risk assessment according to the respective active feature.

Earthquake-Induced Landslides
All of the grid maps were constructed with the same resolution, after which an empirical value of the landslide movement distance (D N ) was calculated and used to estimate the regional landslide displacement.This analysis was conducted using the grid calculation tool in ArcGIS to link the slope critical acceleration map and the peak acceleration, using Equation (3).Landslide displacement classes were identified and mapped as follows: 0-1 cm, 1-5 cm, 5-10 cm, 10-15 cm, 15-20 cm, 20-30 cm, and >30 cm (Figure 4).
The maps exhibit consistent levels of D N > 20 cm along the faults and on the western sides of the faults.This is similar to the geological distribution, where loess-covered areas are more susceptible to landslide activity [53][54][55].This behavior is probably due to the geotechnical features of loess-macro-pores, vertical joints, loose texture, and a silt content greater than 50%, which makes it sensitive to earthquakes and prone to landslides [53][54][55].Furthermore, the model predicts that the displacement decreases with an increase in the distance from the fault.The percent of the slope displacement above 20 cm which is triggered by the Haiyuan Fault and the Liupanshan Fault, respectively, increases with the magnitude of the earthquake, from 0.208% and 0.199% at Ms 7.0 to 0.813% and 0.783% at Ms 8.5 (Figure 5).Stress-strain tests in loess areas have shown that most failures in loess occur after 3% of the strain threshold has been attained [56][57][58].In our study, the slope failure depth was an average of 6 m due to landslides induced by Haiyuan Fault earthquakes, with the landslides' occurrence at a displacement of 18 cm.Model results indicate that a landslide displacement of 18 cm results in a landslide density of 0.803% (i.e., 0.00803 landslides per km 2 ) in the study area, assuming a seismic magnitude of Ms 8.5 triggered by the Haiyuan Fault (Figure 6).For comparison, the landslide densities triggered by the Haiyuan and Wenchuan Earthquakes were 0.844% (0.00844/km 2 ) and 1.41% (0.0141/km 2 ), respectively [59].The landslide density in the present study was, therefore, similar to that of the Haiyuan Earthquake.This suggests that the proposed model, as well as using a displacement above 18 cm as the slope failure, provides an accurate co-seismic landslide assessment and can be used for risk assessment.Stress-strain tests in loess areas have shown that most failures in loess occur after 3% of the strain threshold has been attained [56][57][58].In our study, the slope failure depth was an average of 6 m due to landslides induced by Haiyuan Fault earthquakes, with the landslides' occurrence at a displacement of 18 cm.Model results indicate that a landslide displacement of 18 cm results in a landslide density of 0.803% (i.e., 0.00803 landslides per km 2 ) in the study area, assuming a seismic magnitude of Ms 8.5 triggered by the Haiyuan Fault (Figure 6).For comparison, the landslide densities triggered by the Haiyuan and Wenchuan Earthquakes were 0.844% (0.00844/km 2 ) and 1.41% (0.0141/km 2 ), respectively [59].The landslide density in the present study was, therefore, similar to that of the Haiyuan Earthquake.This suggests that the proposed model, as well as using a displacement above 18 cm as the slope failure, provides an accurate co-seismic landslide assessment and can be used for risk assessment.Stress-strain tests in loess areas have shown that most failures in loess occur after 3% of the strain threshold has been attained [56][57][58].In our study, the slope failure depth was an average of 6 m due to landslides induced by Haiyuan Fault earthquakes, with the landslides' occurrence at a displacement of 18 cm.Model results indicate that a landslide displacement of 18 cm results in a landslide density of 0.803% (i.e., 0.00803 landslides per km 2 ) in the study area, assuming a seismic magnitude of Ms 8.5 triggered by the Haiyuan Fault (Figure 6).For comparison, the landslide densities triggered by the Haiyuan and Wenchuan Earthquakes were 0.844% (0.00844/km 2 ) and 1.41% (0.0141/km 2 ), respectively [59].The landslide density in the present study was, therefore, similar to that of the Haiyuan Earthquake.This suggests that the proposed model, as well as using a displacement above 18 cm as the slope failure, provides an accurate co-seismic landslide assessment and can be used for risk assessment.The slope failure density for each township was also calculated as a factor for the risk assessment.The slope failure density was classified into five degrees, according to the landslide densities induced by the Haiyuan and Wenchuan Earthquakes (Table 2).

Peak Ground Acceleration and Population
The grid map of the peak ground acceleration (amax) was generated by combining the distance maps for the scenario, according to Equation (5) (Figure 7).The slope failure density for each township was also calculated as a factor for the risk assessment.The slope failure density was classified into five degrees, according to the landslide densities induced by the Haiyuan and Wenchuan Earthquakes (Table 2).

Peak Ground Acceleration and Population
The grid map of the peak ground acceleration (a max ) was generated by combining the distance maps for the scenario, according to Equation (5) (Figure 7).The slope failure density for each township was also calculated as a factor for the risk assessment.The slope failure density was classified into five degrees, according to the landslide densities induced by the Haiyuan and Wenchuan Earthquakes (Table 2).

Peak Ground Acceleration and Population
The grid map of the peak ground acceleration (amax) was generated by combining the distance maps for the scenario, according to Equation (5) (Figure 7).The population and infrastructure are mainly concentrated in towns (as identified by the government census) in China.Thus, the peak ground acceleration in these locations was selected as the input factor for the risk assessment.The study area belongs to an eight-intensity fortification area, where the government requires buildings and infrastructures to be designed to resist an intensity eight earthquake (equal to a peak ground acceleration of 0.2 g [60]).This means that an earthquake with an intensity level of less than eight would not damage any buildings or infrastructure in this earthquake zone [61].The hazards for the peak ground acceleration can be divided into five degrees of intensity, as shown in Table 3.The population data includes both the number of people and the population density in each township unit.The average population densities used for this work were 143 people/km 2 in China and 134 people/km 2 in the study area.In order to determine a quantitative risk assessment, we divided the population densities into five levels according to the average population density in China and the study area.This classification is shown in Table 4.These three factors (slope failure density, peak ground acceleration, and population density) were normalized and then the risk value was calculated in each township.The risk classes were identified and mapped according to the three classification factors.The population and infrastructure are mainly concentrated in towns (as identified by the government census) in China.Thus, the peak ground acceleration in these locations was selected as the input factor for the risk assessment.The study area belongs to an eight-intensity fortification area, where the government requires buildings and infrastructures to be designed to resist an intensity eight earthquake (equal to a peak ground acceleration of 0.2 g [60]).This means that an earthquake with an intensity level of less than eight would not damage any buildings or infrastructure in this earthquake zone [61].The hazards for the peak ground acceleration can be divided into five degrees of intensity, as shown in Table 3.The population data includes both the number of people and the population density in each township unit.The average population densities used for this work were 143 people/km 2 in China and 134 people/km 2 in the study area.In order to determine a quantitative risk assessment, we divided the population densities into five levels according to the average population density in China and the study area.This classification is shown in Table 4.These three factors (slope failure density, peak ground acceleration, and population density) were normalized and then the risk value was calculated in each township.The risk classes were identified and mapped according to the three classification factors.

Scenario-Based Risk Assessment
The maps exhibit consistent levels of seriously affected townships, mainly along the fault and on the western side of both the Liupanshan and Haiyuan Faults.This is consistent with the geological features and the landslide distribution, where loess-covered areas are more susceptible to slope failure.Furthermore, the model predicted the decreasing displacement with an increase in the distance from the fault.
The study results in Figure 8 show that the number of high-risk townships increases with the magnitude of the earthquake.The number of high-risk townships is much larger for earthquakes occurring from the Liupanshan Fault than from the Haiyuan Fault for a given earthquake magnitude.The high-risk boundary for an Ms 8.5 magnitude earthquake is in agreement with the area covered by the level IX intensity of the Ms 8.5 Haiyuan earthquake in 1920.The severely-affected counties of the Ms 8.0 Wenchuan Earthquake in 2008 are also in agreement with the area covered by the level IX intensity of the earthquake according to the risk assessment conducted after the Wenchuan Earthquake [12].These results show that the factors selected and the proposed quantitative classification of the assessment results are in practical agreement with known events and are appropriate for creating earthquake risk scenarios.Therefore, the results can be used to develop response plans for earthquake hazard relief.

Scenario-Based Risk Assessment
The maps exhibit consistent levels of seriously affected townships, mainly along the fault and on the western side of both the Liupanshan and Haiyuan Faults.This is consistent with the geological features and the landslide distribution, where loess-covered areas are more susceptible to slope failure.Furthermore, the model predicted the decreasing displacement with an increase in the distance from the fault.
The study results in Figure 8 show that the number of high-risk townships increases with the magnitude of the earthquake.The number of high-risk townships is much larger for earthquakes occurring from the Liupanshan Fault than from the Haiyuan Fault for a given earthquake magnitude.The high-risk boundary for an Ms 8.5 magnitude earthquake is in agreement with the area covered by the level Ⅸ intensity of the Ms 8.5 Haiyuan earthquake in 1920.The severely-affected counties of the Ms 8.0 Wenchuan Earthquake in 2008 are also in agreement with the area covered by the level Ⅸ intensity of the earthquake according to the risk assessment conducted after the Wenchuan Earthquake [12].These results show that the factors selected and the proposed quantitative classification of the assessment results are in practical agreement with known events and are appropriate for creating earthquake risk scenarios.Therefore, the results can be used to develop response plans for earthquake hazard relief.

Haiyuan Risk Zone
The very high-risk area covers between 0 and 25 counties, depending on the earthquake magnitude, and accounts for 0-37.76% of the severely affected area.The high-hazard zone is mainly located in the severely affected areas along the Haiyuan Fault and the western side of the fault.The areas where the landslide concentration is the highest are covered with loess and lie within a high density of population.The area of the very high-risk zone increases as the earthquake magnitude increases and sharply increases when the magnitude is higher than 8.0.

Liupanshan Risk Zone
The high-risk area covers 20 counties with an area of 7691.64 km 2 , accounting for 10.5% of the severely affected areas.The high-hazard zone is mainly located in the severely affected areas along the Liupanshan Fault.The percentage of high-risk areas in the Guyuan Region associated with the Liupanshan Fault is more than the percentage affected by the Haiyuan Fault because the total area affected by the Liupanshan Fault is larger than that affected by the Haiyuan Fault.The very high-risk zone is also mainly located on the western side of the fault, due to the concentration of landslides in this area.

Discussion and Conclusions
Scenario-based earthquake risk assessment is a very important tool for reducing earthquake-induced losses in active fault areas [15,16].However, the methodology of seismic risk assessment, from seismic hazard assessment to the evaluation of the potential losses, encompasses numerous uncertainties [15].Many researchers have been attempting to develop a new method that is simple, accurate, and easy to quantify factors to enable scenario-based earthquake risk assessment, but to date, no single method can be universally used to quickly conduct a risk assessment for proper budgetary planning and earthquake mitigation [17,18].The HAZUS model, which is widely used and provides an accurate assessment, needs a considerable amount of data and uses many uncertain methods that prevent the process from being used for all cases or creating an assessment quickly [14,17].In contrast, the method proposed in this study needs only three factors to quickly develop a scenario-based earthquake risk assessment, all of which can be quantified and can be easily obtained.
The displacement criteria used to assess slope failure is based on the stress-strain curve derived from triaxial tests using loess soil.There is no specific value used to assess the slope failure and most studies take 5 cm as the displacement criteria to determine slope failure [62].Displacement, in terms of slope failure, is primarily affected by the geology of the landslide sliding surface and the slope volume.In mountain areas, most landslides triggered by earthquakes and precipitation slide along the interface between the regolith and rock [63,64].In these situations, the displacement can be over 2 m prior to the occurrence of the slope failure, while some slopes fail when the displacement is only 1 mm [65].Thus, displacement cannot easily determine the slope failure, due to the significant spatial variation found in slope geological characteristics, which are also affected by many other factors [65].In the present study, the proposed method assumes that the slope depth and the criteria displacement are consistent in all areas.However, this assumption is not necessarily true and may result in prediction errors.
The value of a max is different at the plate and at the footwall, given the same distance to the fault and the same geological characteristics.The relationship between a max at different plates remains unclear [66,67].Meanwhile, the typical movement direction of strike-slip faults, the propagation direction of the seismic wave, and the upper and lower plates also affect a max and the slope failure distribution.This was proven in the Wenchuan and Chi-Chi Earthquakes [68][69][70].Therefore, future studies should ideally revise the prediction model with larger datasets.
In this paper, a method was proposed to provide scenario-based earthquake risk assessment for earthquakes along the Haiyuan and Liupanshan Faults, using peak ground acceleration (PGA), the density of the landslides triggered by an earthquake, and the population for earthquake scenarios of magnitude Ms 8.5, 8.0, 7.5, and 7.0 triggered by one of the two faults.We used the Guyuan Region as an example to develop a scenario-based earthquake risk assessment for the area.The township units at very high-risk ranged from 0 to 25 townships and the affected area ranged from 0-37.76% for scenarios of earthquake magnitudes Ms 7.0-8.5,respectively, for an earthquake triggered by the Haiyuan Fault.For an earthquake triggered by the Liupanshan Fault, 1-45 townships and an area of 0.09-58.31% of the total area were at very high-risk for earthquakes of magnitude Ms 7.0-8.5,respectively.The high-hazard zone was mainly located in the severely affected areas along the faults and along the western side of the fault.This is consistent with the geology of the region and the landslide distribution, where loess-covered areas are more susceptible to slope failure.Maps of the estimated risk assessment for the earthquake scenarios of magnitudes Ms 8.5, 8.0, 7.5, and 7.0, generated by the proposed method, will be useful for emergency preparation planning, response plans, and resource assessment.

Figure 1 .
Figure 1.The engineering geology and faults around the study region: (A) Hard rock; (B) loess; (C) moderate rock; (D) very soft rock; and (E) soft rock.(The earthquake data is from the China Earthquake Administration).

Figure 1 .
Figure 1.The engineering geology and faults around the study region: (A) Hard rock; (B) loess; (C) moderate rock; (D) very soft rock; and (E) soft rock.(The earthquake data is from the China Earthquake Administration).

Figure 2 .
Figure 2.Flow chart showing the steps involved in producing a seismic risk map (c is the cohesive strength, φ is the internal friction angle, r is the material unit weigh, a is the slope angle, Fs is the static factor of safety, and DN is the landslide movement distance in cm.DEM is Digital Elevation Model).

Figure 2 .
Figure 2.Flow chart showing the steps involved in producing a seismic risk map (c is the cohesive strength, ϕ is the internal friction angle, r is the material unit weigh, a is the slope angle, F s is the static factor of safety, and D N is the landslide movement distance in cm.DEM is Digital Elevation Model).

Figure 3 .
Figure 3.The slope, population density, and distance to faults; (a) the slope of the study area, (b) the population density of the study area, (c) the distance to the Haiyuan Fault, and (d) the distance to the Liupanshan Fault.

Figure 3 .
Figure 3.The slope, population density, and distance to faults; (a) the slope of the study area, (b) the population density of the study area, (c) the distance to the Haiyuan Fault, and (d) the distance to the Liupanshan Fault.

Figure 4 .
Figure 4.The distribution of the landslide movement distance (DN) along the Haiyuan and Liupanshan Faults: The seismic magnitude of Ms 7.5 along the Haiyuan Fault (a), Ms 8.5 along the Haiyuan Fault (b), Ms 7.5 along the Liupanshan Fault (c), and Ms 8.5 along the Liupanshan Fault (d).

Figure 5 .
Figure 5.The area percentage of the displacement triggered by the Haiyuan/Liupanshan Fault at different magnitudes; the grey gradual change bar is the area percentage of the displacement associated with the yaxis, triggered by the Haiyuan Fault at the seismic scenario, and the blue gradual change bar is the area percentage of the displacement associated with the y-axis, triggered by the Liupanshan Fault at the seismic scenario).

Figure 4 .
Figure 4.The distribution of the landslide movement distance (D N ) along the Haiyuan and Liupanshan Faults: The seismic magnitude of Ms 7.5 along the Haiyuan Fault (a), Ms 8.5 along the Haiyuan Fault (b), Ms 7.5 along the Liupanshan Fault (c), and Ms 8.5 along the Liupanshan Fault (d).

Figure 4 .
Figure 4.The distribution of the landslide movement distance (DN) along the Haiyuan and Liupanshan Faults: The seismic magnitude of Ms 7.5 along the Haiyuan Fault (a), Ms 8.5 along the Haiyuan Fault (b), Ms 7.5 along the Liupanshan Fault (c), and Ms 8.5 along the Liupanshan Fault (d).

Figure 5 .
Figure 5.The area percentage of the displacement triggered by the Haiyuan/Liupanshan Fault at different magnitudes; the grey gradual change bar is the area percentage of the displacement associated with the yaxis, triggered by the Haiyuan Fault at the seismic scenario, and the blue gradual change bar is the area percentage of the displacement associated with the y-axis, triggered by the Liupanshan Fault at the seismic scenario).

Figure 5 .
Figure 5.The area percentage of the displacement triggered by the Haiyuan/Liupanshan Fault at different magnitudes; the grey gradual change bar is the area percentage of the displacement associated with the y-axis, triggered by the Haiyuan Fault at the seismic scenario, and the blue gradual change bar is the area percentage of the displacement associated with the y-axis, triggered by the Liupanshan Fault at the seismic scenario).

Figure 6 .
Figure 6.The slope failure distribution with a displacement above 18 cm; (a) seismic magnitude of Ms 7.5 along the Haiyuan Fault, (b) seismic magnitude of Ms 8.5 along the Haiyuan Fault, (c) seismic magnitude of Ms 7.5 along the Liupanshan Fault, and (d) seismic magnitude of Ms 8.5 along the Liupanshan fault.

Figure 6 .
Figure 6.The slope failure distribution with a displacement above 18 cm; (a) seismic magnitude of Ms 7.5 along the Haiyuan Fault, (b) seismic magnitude of Ms 8.5 along the Haiyuan Fault, (c) seismic magnitude of Ms 7.5 along the Liupanshan Fault, and (d) seismic magnitude of Ms 8.5 along the Liupanshan fault.

16 Figure 6 .
Figure 6.The slope failure distribution with a displacement above 18 cm; (a) seismic magnitude of Ms 7.5 along the Haiyuan Fault, (b) seismic magnitude of Ms 8.5 along the Haiyuan Fault, (c) seismic magnitude of Ms 7.5 along the Liupanshan Fault, and (d) seismic magnitude of Ms 8.5 along the Liupanshan fault.

Figure 7 .
Figure 7.The distribution of the peak ground acceleration (amax) along the Haiyuan Fault and the Liupanshan Fault: (a-d): The seismic magnitude is Ms 7.0-8.5 along the Haiyuan Fault, and (e-h): The seismic magnitude is Ms 7.0-8.5 along the Liupanshan Fault.

Figure 7 .
Figure 7.The distribution of the peak ground acceleration (a max ) along the Haiyuan Fault and the Liupanshan Fault: (a-d): The seismic magnitude is Ms 7.0-8.5 along the Haiyuan Fault, and (e-h): The seismic magnitude is Ms 7.0-8.5 along the Liupanshan Fault.

Figure 8 .
Figure 8. Risk assessment zoning maps in the Guyuan Region, at the township level, for a seismic scenario of Ms 7.0-8.5 along the Haiyuan Fault and the Liupanshan Fault: (a-d): A seismic scenario of Ms 7.0-8.5 along the Haiyuan Fault; and (e-h): A seismic scenario of Ms 7.0-8.5 along the Liupanshan Fault).

Figure 8 .
Figure 8. Risk assessment zoning maps in the Guyuan Region, at the township level, for a seismic scenario of Ms 7.0-8.5 along the Haiyuan Fault and the Liupanshan Fault: (a-d): A seismic scenario of Ms 7.0-8.5 along the Haiyuan Fault; and (e-h): A seismic scenario of Ms 7.0-8.5 along the Liupanshan Fault).

Table 1 .
Values of the cohesive strength (c) and the internal friction angle (ϕ) and the site categories for different lithologies classified by rock hardness.

Table 2 .
The classified landslide density categories.

Table 2 .
The classified landslide density categories.

Table 2 .
The classified landslide density categories.

Table 3 .
Classification of the peak ground acceleration (amax) categories.

Table 4 .
Classification of the population density categories.

Table 3 .
Classification of the peak ground acceleration (a max ) categories.

Table 4 .
Classification of the population density categories.