Landslide hazard probability and risk assessment at the community level: a case of western Hubei, China

Small communities living in mountainous terrain in Hubei province are often affected by landslides. Previous studies by the China Geological Survey focused on the 1:100 000 scale. Therefore, a more detailed assessment, especially at the community level, is urgently required by local governments for risk management. In this study, we conducted a more detailed semiquantitative landslide and risk assessment at the community level using a scale of 1:10 000. We applied the probabilistic method to assess landslide spatial, temporal, and size probabilities, while the hazard and risk assessment were considered for four return periods (5, 10, 20, and 50 years) and two size scenarios (landslide volume). The spatial probability from susceptibility mapping with an accuracy of 84 % indicates that the major controlling factors are Quaternary deposits and weathered eluvium from Ordovician limestones. This study revealed that most building areas in hazard maps are at the foot of major slopes with very high hazard probabilities, and therefore we computed the potential loss of life and property for each slope. The results reveal that 1530 people and USD 18 million worth of property were at risk of landslides within a 50-year return period and a landslide volume of 50 000 m³. The longer the return period is, the higher the hazard probability is. Compared with the classic inverse gamma and power law distribution of landslide magnitude and frequency, the function by the ordinary least squares method is more suitable for landslide size probability analysis of the study area. According to these methods, the proposed procedure of landslide risk assessment proves more useful than the existing data from the 1:100 000 scale in western Hubei, China.

1 Introduction

Risk analysis and assessment in China is an effective means of reducing casualties and economic losses induced by landslides. Although theory and techniques applied worldwide are available (Van Westen et al., 2005; Lee et al., 2007; Neuhäuser and Terhorst, 2007; Erener et al., 2016; Huang et al., 2017; Jiménez-Perálvarez et al., 2017; Van Westen and Greiving, 2017), these are not yet well utilized in western Hubei, China, where the current research was undertaken. To date, very little scientific work is documented at the community level for the study area. Mountainous communities in the area are exposed to landslides because of high rainfall and urbanization. Annually, road construction and anthropogenic modifications at the community level (e.g., excavations in search of building materials), the degree of urbanization, and subsequent population growth have promoted the frequency of landslides, with immense risk to the communities, causing death and unaccountable property loss. According to the Chinese geological disaster notification report (2017), for example, around 850 disasters occurred in 2017 in Hubei province, China, causing 23 deaths and about CNY 254 million (equal to USD ∼36.29 million by November 2019 conversion) of economic losses. Quantifying landslide risk and developing a reduction strategy remains a challenging issue.

Currently, guidelines for landslide risk zoning and land use planning with the framework, definitions, and recommendations are available for clearly defined scales (Fell et al., 2008a). Also, research work highlighting landslide risks at a community level has recently been conducted, and some of the results reached are available to the public (Erener and Düzgün, 2012; Abdulwahid and Pradhan, 2016; H. X. Chen et al., 2016; Liu et al., 2016; McAdoo et al., 2018; Paliaga et al., 2019). Paliaga et al. (2019), for example, used a spatial multicriteria analysis technique to propose geohydrological risk mitigation measures in a small but densely populated catchment, with descriptive parameters involving the extent of urban development and elements at risk. In Spain, a quantitative assessment of landslide risk for the road network of the Basque Country was used for calculating hazard probability and expected consequences (Mavrouli et al., 2019).

Our objective in this study is to conduct a community level landslide hazard probability and risk assessment of an area with limited landslide data and damage records. Despite these limitations, in this study, we attempted to quantify the landslide risk for the Yuyangguan community, in Hubei province, China. The probabilistic method was used to assess landslide spatial, temporal, and size probabilities. Landslide hazard and risk assessment are considered for four return periods and two magnitude scenarios. From these, we propose risk reduction strategies that stakeholders can use for risk management and control.

2 Study area and data

The study area, Yuyangguan community (30^∘07^′–30^∘12^′ N, 111^∘01^′–111^∘07^′ E), is in Wufeng county (29^∘56^′–30^∘25^′ N, 110^∘14^′–111^∘20^′ E), western Hubei province, China (Fig. 1). It covers an area of about 34 km² with 15 villages, including Sanfangping, Dafangping, and Caojiaping. This area has been inhabited since 750 years ago, but intense urbanization development involving the construction of national-class roads began only in 2012. The area was selected for this study because of frequent landslides responsible for tremendous damage in recent years. The residential area is surrounded by steep slopes, with elevations ranging from 180 to 680 m a.s.l. The climate is typical monsoonal with annual average precipitation of about 1500 mm.

Figure 1 (a) Insert representing Wufeng county hosting the Yuyangguan and neighboring Chengguan community; (b) the distribution of lithological units and recorded landslides in the case community. (The lithology codes in the legend are described in Table 1.)

The main lithological units outcropping in the study area comprise Ordovician (O₂g) limestones and shales, Silurian (S₂s) sandstones and shales, Devonian sandstones with coal layer, Permian sandstones, and recent deposits. The S₂s sandstones dominate the rock types along the community slopes, and intense weathering explains the low mass strength of the bedrock. Consequently, weathered S₂s and O₂g rocks are the primary sources of landslides, with the records and bedrocks presented in Table 2.

The Xiannvshan fault (Fault I in Fig. 1) and Yuyangguan fault (Fault II in Fig. 1) of the Changleping tectonic belt constrain the bulk architecture of the area. The Xiannvshan fault, striking approximately 340^∘ N–345^∘ W and dipping 60–70^∘ NE, is a transpressional fault terminating at the center of the community. The fault is a seismically active belt, exemplified by the Panjiawan earthquake (M_s=4.9) of 1961 and the Zigui earthquake (M_s=3.3) of 1972. Conversely, the Yuyangguan fault, striking E–W and dipping 60–70^∘ S, is a transtensional fault. Rocks associated with the Yuyangguan fault comprise a 20–50 m wide area of cataclasite and brecciated mylonite, with several secondary faults merging into it, and, together, go through the center of the Yuyangguan community. According to the China Earthquake Administration, the studied area is in a weak-seismic-activity region with a basic earthquake intensity of VI, and the 50-year 10 % probability exceedance of the peak ground acceleration is 0.05 g. No historical record is available for earthquake-induced landslides in the area.

Rainfall and anthropogenic activities contribute significantly to triggering mass movements, mainly landslides. In the rainy season of 2013, a slope along the main road collapsed, breaking and causing lengthy traffic jams and transportation problems (Fig. 2). A landslide occurred in a new residential quarter of the community on 3 January 2013, due to an unstable foundation pit after a 1 d excavation (Fig. 3). These two landslides are examples, with more landslide records prepared from aerial photograph interpretations and validated by field investigations.

Table 1 Characteristics of lithology distributed among Yuyangguan.

Figure 2 Landslide YYG01 occurred in the rainy season of 2013, causing damage on the national road (G318) in Yuyangguan. (Travel distance = 57 m.) (a) An overview of the landslide YYG01, (b) left side of the loose landslide body, (c) river blockage caused by the landslide YYG01, and (d) the retaining wall.

Figure 3 Landslide YYG03 caused by slope incision in a residential quarter of Yuyangguan.

Table 2 shows all historical landslide data from 1976 to 2013 in areas dominated by Quaternary deposits. Besides the landslide inventory database, other datasets collected for landslide risk mapping include the following.

• A 10 m × 10 m resolution digital elevation model (DEM) generated from a topographic map obtained using an unmanned aerial vehicle (UAV) was used. The DEM allowed the extraction of the slope, elevation, aspect, and curvature data using the surface analysis tool in ArcGIS (a geographic information system for working with maps and geographic information, http://www.esri.com/software/arcgis/arcgisonline, last access: September 2017; Fig. 4a–d).

• A Geological map at a 1:50 000 scale (Fig. 1) was used to extract datasets, including lithology, faults, and slope structure. The slope structure map (Fig. 4e) was generated using the standard and stratigraphic altitude advocated by Cruden (1991). The land use map provided the distribution of rivers and roads (Fig. 4h and i).

• Landslide inventory databases of the case community (named Yuyangguan) and a neighboring community (named Chengguan) were utilized in this study to analyze the size probability in hazard assessment (see Table 2). We used the hazard database of the neighboring community for the following reasons: (1) the landslides in the case community are limited for probability analysis; and (2) the two communities are similar in geomorphology, geology, climate, and landslide types. The location of the Chengguan community is displayed in Fig. 1a, and administratively the two communities are in Wufeng county. Both communities belong to the same structural belt named the Changleping anticline fold that extends in a nearly E–W direction and comprises the Silurian and Ordovician rocks. Previous landslides in these communities involve shallow and soil slope movement of weathered bedrock. The factors triggering the landslides are rainfall and slope cutting, associated with urbanization development.

• A building footprint map (Fig. 5) was interpreted and checked in the field by the authors, with most buildings on or at the toe of the first slope zones, with an elevation of up to 350 m a.s.l. The entire built area is 757 000 m², with data on the economic value of buildings obtained from the Department of Natural Resources of Hubei Province (see Table 3).

• Census data were obtained by integrating the information derived from the China population data (2010) (http://www.stats.gov.cn/tjsj/pcsj/rkpc/6rp/indexch.htm, last access: September 2016) and sampling survey (Fig. 6), amounting to a population of 45 914 in the area.

Figure 4 Thematic maps for landslide susceptibility mapping: (a) slope, (b) elevation, (c) slope aspect, (d) curvature, (e) slope structure, (f) lithology, (g) distance to fault, (h) distance to river, and (i) distance to road.

Figure 5 Building footprint map with economic values for the Yuyangguan community.

Figure 6 Population distribution map of the Yuyangguan community.

Table 2 Historical landslides investigated in the field by the authors.

R – rainfall; SC – slope cut. The code of bedrock is listed in Table 1. The projected coordinate system is Xian_1980_3_Degree_GK_Zone_37 with projection of Gauss_Kruger.

3 Methodology

This section summarizes the methods employed for semiquantitative risk analysis for landslides at a community level. Initially, the slope-unit-based hazard probability was calculated by integrating spatial, temporal, and size probabilities. We generated eight hazard probability maps involving four return periods (5, 10, 20, and 50 years) and two size scenarios. For each map, we determined the potential landslide-influencing area semiquantitatively by calculating the traveling distance of the slope unit. Subsequently, the element-at-risk map was interpreted from the image data, population census data, and sampling survey. The vulnerability map was created by assimilating assets and landslide-influencing areas. Values in the vulnerability map were semiquantitatively determined for buildings and people in the buildings. The combination of hazard maps, vulnerability maps, and element-at-risk maps produced the risk value for each slope unit, contributing to the final risk maps of the population and buildings in the study area. The flowchart of the methodology is depicted in Fig. 7.

Figure 7 Flowchart of the methodology for landslide hazard risk assessment.

3.1 Landslide hazard from spatial, temporal, and size probability analysis

Hazard assessment is an essential step in landslide risk assessment. For the community level, this is achievable through the deterministic model (Gokceoglu and Aksoy, 1996; Qiao et al., 2019), given enough engineering geology data. For the study area, soil or rock strength parameters were unavailable for use in the deterministic method. We therefore used three probabilities (spatial, temporal, and size) to resolve questions on where and how potential landslides will occur with absolute magnitude in a given time (Guzzetti et al., 2005) using the following expression:

$$\begin{array}{}\text{(1)}& H=P\left(S\right)\times P\left(N_{\mathrm{L}}\right)\times P\left(A_{\mathrm{L}}\right),\end{array}$$

where H represents the hazard probability, P(S) is the spatial probability, P(N_L) is a temporal probability, and P(A_L) is size probability.

3.1.1 Spatial probability

Spatial probability, based on the concept of susceptibility, assesses the locations where a mass movement exists or may potentially occur. Landslide susceptibility mapping (LSM) is now widely used by researchers (Ayalew et al., 2004; Fell et al., 2008b; Van Westen et al., 2008; Guzzetti et al., 2012). In this study, morphometric and geoenvironmental factors, including altitude, slope, aspect, curvature, slope structure, distance to rivers, and proximity to roads, were chosen as variables. The morphometric factor maps were derived from the DEM with a 10 m × 10 m resolution obtained by a UAV. Geoenvironmental factors, such as lithology and faults, were prepared and transformed from shapefiles to grid-based maps using a 1:50 000 scale geological field map from the China Geological Survey (http://www.cgs.gov.cn/, last access: September 2016). Detailed processing steps in ArcGIS are provided in Catani, Casagli, Ermini, Righini, and Menduni (2005). The commonly applied weight of evidence (WoE) method was used to assess landslide susceptibility in this study. This is a probabilistic model considering evidence factors of landslides, based on the conditional independence hypothesis (Hong et al., 2017). In the WoE method, W⁺ and W⁻ were used as the weights where the evidence was present or absent, respectively, and the contrasts (differences between W⁺ and W⁻) were used as the weight for the morphometric and geoenvironmental factors.

The effectiveness of the LSM was tested using the receiver operating characteristics (ROC) curve (Metz, 1978; Zezere et al., 2017), with the area under the ROC curve used to assess the success rate.

The above grid-based LSM data were then converted into a slope-unit-based susceptibility map. The susceptibility value of each slope was calculated from the average value of the susceptibility of the grids within the slope. The slope unit was subdivided using the hydrology analysis method in the ArcGIS platform. The slope-unit-based spatial probability map was classified into five classes, including very high, high, moderate, low, and very low.

3.1.2 Temporal probability

Landslide temporal probability P(N_L) is evaluated by assuming that slope failures are independent random point events in the time domain (Crovelli, 2000; Guzzetti et al., 2006). In this study, the Poisson model (Crovelli, 2000) was adopted for constructing temporal probability. It is the exceedance probability of landslide occurrence during a given period (see Eq. 2), meaning the probability of experiencing one or more landslides during a given time, and derived as

$$\begin{array}{}\text{(2)}& P\left(N_{\mathrm{L}}\right)=\mathrm{1}-e^{-T/\mathrm{RI}},\mathrm{RI}=t/N,\end{array}$$

where T is the return period, e.g., 1, 10, 20, and 50 years; the recurrence interval (RI) is the historical mean recurrence interval for each slope unit; t is the temporal interval of the landslide database; and N is the number of landslides recorded in each slope. If the historical landslide database is incomplete, N on a slope unit with a very high susceptibility class in the LSM is assigned as 1.

3.1.3 Size probability

Landslide size probability is calculated based on the relationship between landslide volume and cumulative frequency. Guzzetti et al. (2005) used the probability density function of the landslide area to predict the probability of a specific landslide area in each slope unit. Stark and Hovius (2001) found that landslides in New Zealand and Taiwan fitted with a double Pareto probability distribution, and two differences in our study will improve the application. The first is that the present study seeks the most suitable distribution for the case study area. Therefore, we compared the distributions stated previously and introduced another type using the ordinary least squares (OLS) method in the MATLAB software, with the best fit used for probability calculation. The second is that landslide volume is an acceptable indicator for risk control practice in the study area. Therefore, we converted the landslide size probability distribution from area to volume using the volume–area relationship simulated by the OLS method. For this, we used the historical landslides in the database of the study area in Table 2. The hazard database of the Chengguan community was implemented because of the landslide records limitation, but the geoenvironment and hazards in both communities are similar. Meanwhile, the two size scenarios were determined from the distribution of landslide volume in the case study area.

3.2 Data preparation for elements at risk

The next step was determining the elements at landslide risk. This study focuses on residential buildings and people within. The building footprint map (see Fig. 5) was interpreted from the 2013 UAV image data. The building structure (reinforced, reinforced concrete, masonry, and wooden) and numbers of floors are involved in the building map database. To express the risk in monetary values, we used the economic value of the buildings (see Table 3) obtained from the Department of Natural Resources of Hubei Province. These data were converted to building values by multiplying the unit economic values, footprint areas, and the number of floors. Data for the population in buildings were obtained by integrating the information from the China population census data (2010) and the sampling survey. The average number of people per building was calculated by applying a dasymetric mapping approach (a methodology for generating a surface-based representation of the population; Mennis, 2003), which contributes to the data of the population in each building. To assess the element at risk, the building footprint map was then combined with the potential landslide influence area at the community level. The influence area of each slope was semiquantitatively determined by calculating the travel distances using the following formula (Hungr et al., 2005):

$$\begin{array}{}\text{(3)}& \log (H/L)=A+B\times \log V,\end{array}$$

where L is the travel distance, H is the slope height, V is slope volume, and A and B are constants. These constants are referred to in Corominas (1996).

Table 3 The economic value of buildings in the Yuyangguan community (the Department of Natural Resources of Hubei Province, 2016).

3.3 Vulnerability analysis and risk assessment

Quantitative vulnerability analysis is still a challenge in landslide risk assessment (Chen et al., 2011; Peduto et al., 2017). Physical vulnerability assessment is performed in a large- or local-scale area (Fell et al., 2008b; Li et al., 2010; Quan Luna et al., 2011). In this study, the physical vulnerability was semiquantitatively determined for buildings, using two indicators. The first indicator relates to building structures, such as reinforced, reinforced concrete, masonry, and wooden, while the second indicator is the landslide travel distance. The assumption is that vulnerability is 1.0 for the buildings on the slope, and it decreases from the toe of the surface rupture to the farthest travel distance. We therefore propose vulnerability values for different types of buildings (see Table 4).

Table 4 Vulnerability value of buildings impacted by landslides (proposed by authors).

L is travel distance; H is slope height.

The vulnerability of the population in buildings follows a power–law relationship with building vulnerability (Li et al., 2010) as follows (see Eq. 4):

$$\begin{array}{}\text{(4)}& V_{\mathrm{p}}=\mathrm{0.0014}\times e^{\mathrm{6.07}\times V_{\mathrm{b}}},\end{array}$$

where V_p is the vulnerability of the population in buildings and V_b is the vulnerability of a building.

A landslide risk map was then generated in ArcGIS based on the concept defined by the IAEG (International Association of Engineering Geology) and Varnes (1984) as “the expected number of lives lost, persons injured, damage to property, and disruption of economic activity due to a particularly damaging phenomenon for a given area and reference period”. The conceptual equation for risk is

$$\begin{array}{}\text{(5)}& R=H\times V\times E,\end{array}$$

where R is the expected loss for some return period, and H is the landslide probability of some return period with a given size scenario. In the present study, V is the physical vulnerability of buildings or the population in the buildings, and E is the quantification of the exposed elements at risk. Using Eq. (5), the risk curve is fitted by plotting the probability versus potential loss, with the annual risk calculated from the area under the risk curve (Van Westen et al., 2002).

4 Results

This section provides the results of the case study to illustrate the application of the proposed framework and methodology in Sect. 3.

4.1 Landslide susceptibility assessment

In assessing landslide susceptibility, we investigated the elevation, slope, aspect, curvature, lithology, slope structure, distance to fault, rivers, and roads. The weights and contrasts values from the WoE method for LSM of Yuyangguan are presented in Table 5. According to the contrast values explained in Sect. 3.1 for lithology, the Ordovician limestones and shales (O₂g) and Quaternary eluvium ($Q_{\mathrm{4}}^{{\mathrm{d}}_l+{\mathrm{e}}_l}$) are the top two units, implying that these units are susceptible to erosion. Under this background, heavy rainfall in the area accelerates erosion, thus triggering landslides.

General road construction is of secondary importance, with a contrast value of 0.95. Data in Table 5 also reveal generally significant morphometric factors in the study area. For aspect and elevation, the contrast values are elevated on the north-facing slopes with an elevation from 0 to 260 m a.s.l., but low on the south-facing slopes with an elevation above 350 m a.s.l. For the slope, the steeper the slope is, the higher the landslide probability is, with the contrast value for slopes ranging from 10 to 30^∘ being 0.19, indicating a relatively high landslide probability.

The grid-based susceptibility map was converted to a slope-unit-based map with 701 slope units in total. The slope-unit-based susceptibility map was ordered into five classes, ranging from very low to very high (Fig. 8). The performance from the map reveals an accuracy of 84 %, using the ROC curve. The landslide susceptibility is very high for the north-facing slopes along the main road, especially where $Q_{\mathrm{4}}^{{\mathrm{d}}_l+{\mathrm{e}}_l}$ and O₂g rocks are present. These results correspond well with the contrast values presented in Table 5.

Figure 8 Slope-unit-based LSM at the community level, converted from the grid-based LSM using the weight-of-evidence method for Yuyangguan, western Hubei, China.

4.2 Landslide hazard probability

The landslide hazard involves spatial, temporal, and size probabilities. The landslide data for Yuyangguan presented in Table 2 cover the 33 years from 1981 to 2013. For each slope unit, the historical mean RI is calculated using Eq. (2). Assuming that the past is an indicator of the future, landslides in the study area may be modeled over the next 50 years, based on the past 50 years.

Four landslide temporal probability maps are displayed in Fig. 9 for four return periods (5, 10, 20, and 50 years). The map for 50 years, for example, shows the highest probability for the slope units experiencing landslide events among the four maps. Slope units with high and very high probability values (>0.5) cluster on the first slope zones around the community. The probability values on the slopes increase from the return period of 5 years. For example, the enlarged windows in maps of Fig. 9 demonstrate that the slope in Huanglongzhai village experienced a very low probability in the 5 years to a high class in the 50-year return period.

Figure 9 Landslide temporal probability maps based on the Poisson model, showing the exceedance probability of landslide occurrence in each slope unit for four return periods (5, 10, 20, 50 years).

Table 5 The weight and contrast values by the weight-of-evidence model for landslide susceptibility mapping of Yuyangguan, western Hubei, China.

The total area of the community is 29.14 km². Total landslide area is 3740 pixels. ^a Variance of W⁺. ^b Variance of W⁻. ^c Standard deviation of contrast. ^d Studentized value of contrast.

The landslide probability distribution curves are created using three different fitting functions, as displayed in Fig. 10. In comparison to the inverse gamma and the power law distributions, the function by the OLS method shows the best fit, with the lowest root mean square error of 0.05. This indicates that the landslide frequency distribution function by the OLS method is the most appropriate technique to apply in Yuyangguan. The volume–area relationship is analyzed in Fig. 11a with an R² value of 0.915. This indicates that converting the size probability distribution from landslide area data to volume is feasible in the study area. Meanwhile, two size scenarios are determined from the cumulative frequency curve (Fig. 11b) based on a landslide volume of 50 000 and 100 000 m³. The number of landslides begins to increase rapidly in Fig. 11b when the volume is greater than 50 000 m³, representing a probable threshold value of landslide hazard volume, while a landside volume greater than 100 000 m³ represents the maximum for the landslides, as shown in Fig. 11b. The value of 100 000 m³ is a standard threshold for landslide classification provided by the China Geological Survey and a widely accepted value in landslide hazard risk control in China.

Figure 10 Magnitude–frequency relationships simulated by the ordinary least squares (OLS) method, inverse gamma, and power law for landslide in western Hubei. The samples are presented in Table 2 and involve the Chengguan database except for the study area.

Figure 11 (a) Relationship between landslide volume and area (R²=0.915) used to convert magnitude–frequency relationships for the area to landslide volume. (b) Frequency distribution of landslide volume. The samples are presented in Table 2, involving the Chengguan database, except for the study area.

Integrating the spatial probability in Fig. 9, temporal probability in Fig. 10, and the magnitude–frequency relationships in Fig. 11, eight hazard maps were produced using Eq. (1). Through these maps, we demonstrate the hazard probability for each slope unit for four return periods (5, 10, 20, and 50 years) and landslide sizes equal to or greater than 50 000 m³ and equal to or greater than 100 000 m³. For example, Fig. 12 shows the four landslide hazard maps for the four return periods and the landslide size scenarios equal to or greater than 50 000 m³. In these maps, the landslide hazard probability values comprise five categories from very low (0.0–0.1) to very high (0.4–0.5). The maps for the 50-year return period show very high hazard probabilities for most building areas in the community at the bottom or near slopes. Conversely, buildings located in the community center are in areas from low to very low hazard probability class.

Figure 12 Landslide hazard maps for four return periods (5, 10, 20, and 50 years) and the landslide sizes scenario equal to or greater than 50 000 m³. The maps were generated by integrating the spatial probability in Fig. 8, temporal probability in Fig. 9, and the magnitude–frequency relationships in Fig. 10.

4.3 Element-at-risk and vulnerability assessment

We assume that the buildings and population in the buildings are exposed to the slopes with high and very high class probability in the hazard maps. Therefore, the element-at-risk maps for buildings and the population generated from the data in Sect. 2 and methodology in Sect. 3.2 are consistent with the four return periods and two size scenarios of the landslide hazard probability maps.

Data in Table 6 show no exposure for the return periods of 5 years and 10 years, while for the 50-year return period and 50 000 m³ size scenario, 570 000 m² of hosing areas and 14 257 persons are exposed to landslide risk. The data also highlight the potential for damaged building areas and the number of persons in the buildings for return periods of 20 to 50 years. The vulnerability value is indirectly calculated by Eq. (4) in Sect. 3.3. In the size scenario of landslides for the volume equal to or greater than 50 000 m³, a sharp increase in exposure exists for built areas and their population, with the exposure more than doubling from 20 to 50 years. Similar results emerge for the size scenario of landslide volume equal to or greater than 100 000 m³. Comparing the two size scenarios, we find that the number exposed for the 100 000 m³ volume is lower than that for 50 000 m³. The percentage of building areas exposed to landslides in the 20-year return period are 8.76 % for the 50 000 m³ volume, whereas it is 5.86 % for the 100 000 m³ volume, probably due to the lower hazard probability in the latter scenario. A similar tendency is observed for the population in buildings, with about 30 % exposed to landslides for the 50-year return period.

Table 6 Exposure and risk of buildings and population for landslides under four return periods (5, 10, 20, and 50 years) and two size scenarios (landslide volume equal to or greater than 50 000 or 100 000 m³) by using the methodology in Sect. 3.2 and 3.3. (Number in brackets is in percentage.)

Accordingly, eight vulnerability maps for buildings and eight vulnerability maps for the population in the buildings are created for the four return periods and two size scenarios. For example, Figs. 13 and 14 show the resulting vulnerability maps for buildings and the population in buildings for the four return periods and when the landslide volume is equal to or greater than 50 000 m³. Most exposed buildings are located on slopes with high or very high hazard probability, including the Caolinwan, Dafangping, and Sanfangping villages. The buildings and population in buildings in these central communities are not exposed to slope hazard because of the very low class (0–0.1) of landslide probability.

Figure 13 Buildings exposed to landslides, vulnerability distribution map for four return periods (5, 10, 20, and 50 years), and size scenario of landslide volume equal to or greater than 50 000 m³.

Figure 14 Population exposed to landslides and vulnerability distribution map for four return periods (5, 10, 20, and 50 years) and size scenarios of landslide volume equal to or greater than 50 000 m³.

As presented in Table 4 in Sect. 3.3, we assigned a vulnerability value of 1.0 to buildings on slopes. This causes buildings outside the flat areas in the community to display very high vulnerability values in Fig. 13, involving five categories ranging from very low to very high. For the 50-year return period and 50 000 m³ size scenario, 18 % of built areas exhibit very high class vulnerability. The vulnerability value for the population was then assigned according to Eq. (4) in Sect. 3.3 from the building vulnerability result. The very high class population vulnerability shown in Fig. 14 was over 0.5, with the area representing about 10 % of the people exposed, concentrated in the slope influence areas of Dafangping, Caolinwan, and Sanfangping or on the slopes with the very high class hazard probabilities displayed in Fig. 12.

4.4 Risk assessment

Landslide risk maps are then generated from the eight hazard maps, eight vulnerability maps, and the values of the elements at risk by Eq. (5). An example of a risk map for the population loss is shown in Fig. 15. The map for the 50-year return period, for example (see Fig. 15d), shows that potential loss is concentrated in the urban and densely populated areas along the Yuyangguan River, especially at the toe of the slope of the Caolinwan, Dafangping, and Sanfangping villages.

Figure 15 Example of a risk map for population loss for four return periods (5, 10, 20, and 50 years) and the landslide size scenario equal to or greater than 50 000 m³.

The final population and building risks for four return periods and two size scenarios are presented in Table 6. Accordingly, no potential losses emerge for the 5- and 10-year return periods. The absence of risk for the 5-year return period is proven by the fact that no casualties or economic losses have been reported in the Yuyangguan community since 2014. For the volume equal to or greater than 50 000 m³, however, the potential casualties are 771 persons, representing 1.68 % of the total population in the community. The economic losses are estimated at USD 8.57 million for the 20-year return period, amounting to 3.7 % of economic values associated with buildings. As presented in Table 6, the risk is expected to double in the next 50-year return period. For the volume equal to or greater than 100 000 m³, the potential loss is lower, with 1.27 % of the population exposed to landslides, amounting to a decrease of 0.41 % relative to the size scenario of 50 000 m³.

5 Discussion

5.1 Discussion on landslide susceptibility map

The landslide susceptibility results demonstrate that lithology is the most important controlling factor. $Q_{\mathrm{4}}^{{\mathrm{d}}_l+{\mathrm{e}}_l}$, weathered Guniutan O₂g limestone, and weathered S₂s sandstone exhibit the top three weight values in lithology contribution in the susceptibility map. This is consistent with findings inform the fieldwork. The residual deposit and eluvium are composed of clay and gravels that are characterized by low strength and cover the main areas with intensive human activities in the Yuyangguan community. For the O₂g bedrock, its high strength makes it less vulnerable to a landslide. However, in the field, we found two groups of surface joints in the O₂g, and this helps to explain why the YYG07 landslide occurred (see Table 2). Field observations also reveal that most landslides occurred in S₂s sandstones and shales. However, the weight value (0.3) of S is lower than those of $Q_{\mathrm{4}}^{{\mathrm{d}}_l+{\mathrm{e}}_l}$ and O₂g. This is because the area covered by the S₂s formation is much larger than the area of $Q_{\mathrm{4}}^{{\mathrm{d}}_l+{\mathrm{e}}_l}$ and O₂g (see Fig. 1). Therefore, high landslide probability for slopes in O₂g limestone must be considered, and attention should be devoted to slopes where historical landslides are not as frequent now.

5.2 Discussion on landslide hazard assessment

Landslide hazard maps are generated for four return periods (5, 10, 20, and 50 years) and two size scenarios. Theoretically, the definition of hazard scenarios, as stated by L. Chen et al. (2016), should be based on the analysis of landslide occurrences and triggering events. Due to incomplete information on landslide dates in the community, it is challenging to establish the relationship between the landslide return period and triggering factor (rainfall in this case study). However, we observed that no landslide has occurred in the Yuyangguan community since 2014, which is consistent with the hazard result for the 5-year return period (see Fig. 12a). This means that the temporal probability approach using the Poisson model is suitable for landslide hazard assessment when the landslide database is missing the occurrence date and triggering event (e.g., rainfall) data.

In the size probability analysis, the landslide probability distribution is the key for quantifications. We found that the classical distribution model (Stark and Hovius, 2001; Malamud et al., 2004) failed to produce an excellent fitting performance in this study. The difference in the landslide size between Malamud's and our landslide database accounts for this inconsistency. No small landslides (<1000 m³ in Malamud's research) are present or recorded in the Yuyangguan community, and the simulation equation (in Fig. 10) is suitable for landslide risk assessment in this study. In the future, however, comparison with classical models should be undertaken, and other factors, such as triggers and landslide types, should be considered using a complete database or one with more landslide events.

5.3 Other limitations in risk results

Uncertainties exist in the final risk maps due to some other factors, such as element-at-risk data and their vulnerability or resilience. In this study, the data for buildings and population in buildings at risk are derived from empirical calculations for the landslide influence area. Further studies considering landslide material, pore water pressure, and ground surface characteristics besides slope height and volume used in this paper will be conducted by numerical modeling.

Meanwhile, the risk results cannot be tested because of a lack of historical damage data in the area, which is a common difficulty in China and other areas (Ghosh et al., 2011; L. Chen et al., 2016). Damage data are also crucial for vulnerability analysis. We considered building typology and location from the landslide to assess the physical vulnerability based on local experts' opinions. However, the resilience of the element at risk contributes to a reasonable decrease in vulnerability. The mobility of persons and their characteristics (e.g., age, education, and physical disability) and disaster prevention capability of the government are not considered in the community. In the future, physical vulnerability curves for buildings and the population will be constructed for the area.

From susceptibility to risk assessment, in general, we assumed that landslides will occur under the same condition as historical landslides. However, morphometric and geoenvironmental conditions will change with time. For example, the slope degree, elements at risk, and land use cover may change because of new infrastructures, or the number of historical landslides may increase due to heavier rainfall. Future studies will be necessary, considering the changing conditions for more accurate results and more practical applications. So far, the series of maps represent a basis for landslide risk control and land use planning in the Yuyangguan community. Risk control measurements can be planned on each slope unit using a matrix combining landslide hazard probability and risk maps. For example, risk management on slope units with very high class hazard probability and very high risk can be suggested as requiring relocation or engineering works. The precondition is that a more detailed geotechnical investigation and comprehensive analysis be complemented. This achievement can be used by professionals on engineering geology. The slopes with high class hazard probability and midclass risk are suggested for monitoring, considering the importance of safety for the surrounding population, roads, or other elements at risk. Meanwhile, the results of the four return periods and two size scenarios are useful for multitemporal land use planning, including short term (5 years), mid-term (10–20 years), and long term (50 years) strategies. The annual risk value on each slope is useful for the cost–benefit analysis of risk decisions. This kind of achievements can be applied for government decision makers.

6 Conclusion

We conducted a semiquantitative risk assessment for landslides at the community level based on the definition of landslide risk provided by Varnes (1984) and the IAEG. In our case study, we focused on the potential damage of buildings and loss of life for the population in the buildings. We generated a susceptibility map, eight hazard maps, eight vulnerability maps, and eight risk maps for four return periods (5, 10, 20, and 50 years) and two size scenarios (equal to or greater than 50 000 m³ and equal to or greater than 100 000 m³). The landslide susceptibility result was tested, and it yielded a success rate of 0.84, highlighting the important contribution of $Q_{\mathrm{4}}^{{\mathrm{d}}_l+{\mathrm{e}}_l}$ and the Guniutan O₂g limestone. The approach for generating hazard maps, which involved integrating three probabilities (spatial, temporal, and size probabilities), proved applicable in the case study area. In the size probability calculation, the use of a normal distribution function for landslides requires caution, with a better fitting function suggested when small landslide data are scarce in an area. Also, the landslide influence area was empirically determined at a community level using simple data (slope height and volume) in the absence of geotechnical parameters. However, for a more accurate vulnerability assessment, numerical modeling on landslide travel distance is suggested because the resulting intensity parameters, such as velocity and depth, are essential input data for vulnerability quantification. Besides the presented limitations, we believe that the proposed risk maps can help local stakeholders in establishing periods of risk planning for the community, including short-term (in 5 years), mid-term (10–20 years), and long-term (50 years) strategies, or provide a reference for the cost–benefit analysis for each slope unit from the quantified annual risk values.