Risk-based models for potential large-scale landslide monitoring and management in Taiwan

ABSTRACT Several related pre-event environment indices coupled with a 200-year event of Typhoon Morakot were used in this study to analyse landslide risk/scale for mapping and verifying the danger grade of watershed landslides. Risk analysis was modelled using maximum daily rainfall of the tested event as the hazard factor, and a complex indicator derived from indices of return period, road development and green deterioration as the vulnerability factor. The vulnerability factor was then modified using the slope to adjust the weight for excluding the flat areas and increasing impact of steep areas. A comprehensive indicator integrated from indices of river concave, headward erosion and dip slope was provided as the hotspot of potential large-scale landslides which were modified by vegetation index to correct sediment yield due to antecedent effects. Finally, sites with high potential large-scale landslide could be interpreted by the danger grade calculated from collapse risk/scale analysis. The results show that the higher the danger grade, the greater the collapse ratio; and the coefficient of determination for the relationship is greater than 0.9. The model developed in this study for accurately estimating the risk/scale of rainfall-induced landslide can be fulfilled.


Introduction
A landslide could be defined as a certain soil volume located on the weak slope which fails when triggered by external forces and affected by gravity (Cruden 1991). Taiwan is located on the convergent boundary between the Eurasian plate and the Philippine Sea plate, and the tectonic activities result in fragile geological feature. In addition, the climate is influenced by the East Asian Monsoon. Several torrential and/or convectional rainfalls cause frequent landslide occurrence in typhoon seasons. In order to enhance landslide assessment and monitoring in Taiwan effectively, hotspots of landslide risk, where the risks of landslide disasters are particularly high (Dilley et al. 2005;Jaedicke et al. 2014), should be identified and located. Risk analysis has been frequently applied in the fields of economy, health, insurance and disaster prediction (Blanciforti and Luster 2006;Tsai and Chen 2010;Hsu et al. 2011;Marhavilasa et al. 2011;Zahran et al. 2011;Braman et al. 2013;Dong et al. 2015). According to the report of the Office of the United Nations Disaster Relief Coordinator, risk is the potential loss of life and property and/or economic activity at a specific natural disaster, and it generally consists of two components, namely, hazard and vulnerability (UNDRO 1979(UNDRO , 1991UNDP 2004;Rausand 2011;Chen and Wu 2016).
Hazard indicates external factors which specify the type and intensity of the hazard. Vulnerability could be regarded as internal factors which refer to the capacity to anticipate, cope with, resist and recover from the impact (Bohle 2001;Taubenb€ ock et al. 2008). Since the complexity of vulnerability has been stated in site-specific and scale-dependent dimensions (Thywissen 2006), many aspects of vulnerability, developing from various physical, social, economic and environmental factors, were integrated to the holistic approaches (Taubenb€ ock et al. 2008;). Exposure was also discussed in several studies relating to risk analysis (Kron 2005;Kannami 2008), and it means any protected objects such as population, buildings and economic activities located at the specific site that could be susceptible to the risk event (UNDRO 1979;Tobin and Montz 1997). Previous studies reveal that earthquake and rainfall are the main external forces which could result in landslides in Taiwan Feng et al. 2012). Hence, the maximum daily rainfall caused by Typhoon Morakot, which brought recordÀbreaking torrential rainfall in 2009 (Chien and Kuo 2011), was adopted in this study as the index of hazard. Several studies derived landslide vulnerability through the various factors, e.g. topography, geology, hydrology, soil, land use, etc. (van Westen et al. 2008;Liao et al. 2011;Khan and ur Rehman 2014). In this study, a complex indicator derived from indices of return period, road development and green deterioration was used as vulnerability for evaluating the spatial distribution of environmental weak sites.
In the previous studies, the most common descriptions of a large-scale landslide are large-area landslide (Liu et al. 2011), huge debris failure (Liu et al. 2011;Zerathe et al. 2014) and/or the bedrock collapse. Many landslide risk models show the probability of landslide occurrence (Remondo et al. 2005;Lee and Pradhan 2006;Zêzere et al. 2008;Akgun et al. 2012), but are unable to predict the landslide scale. Traditionally, aerial photo interpretation (Chigira and Kiho 1994;Fernandez and Whitworth 2016) and LiDAR (Chen et al. 2014) have been used in large-scale landslide delineation; however, the related indices were still seldom provided for quantitative assessment. Therefore, the other indices should be introduced to portray the potential landslide scale. Generally, areas susceptible to large-scale landslide usually have unique characteristics in topography, geology and vegetation; hence, a comprehensive indicator integrated from indices of river concave, headward erosion and dip slope is also introduced to effectively reflect the hotspots of potential large-scale landslide. Finally, risk-based models will be established and verified in the study, and the sites with high potential large-scale landslide could be interpreted by the danger grade calculated from the analysis of collapse risk and failure scale (Figure 1).

Study area
Typhoon Morakot struck Taiwan in 2009, and brought extreme rainfall causing large-scale landslides in Hsiao-Lin village (Tsou et al. 2011) and Erwanping (Huang 2009). The most serious debris disaster event had been triggered and caused several billions of sediments yield in the upper channels ( Figure 2). Hence, the government has gradually paid more attention to large-scale landslides.  Taiwan was regarded as the study area, and the potential risk and/or large-scale landslides were delineated using related environmental indices extracted from existing big data.

Study materials
Related information used in the models was adopted from big data established by the public sectors of government authorities (Table 1). Historic rainfall data were obtained from the Water Resources Agency, Ministry of Economic Affairs to estimate the spatial distribution of 200-year event rainfall. Multi-temporal SPOT imageries were collected and corrected to extract the vegetation index. A 5 m £ 5 m resolution digital elevation model (DEM) derived from the Department of Land Administration was used to carry out terrain analysis, topographic index extraction and watershed-subdivision delineation. Road maps including national, provincial, county and township roads were acquired from the Institute of Transportation, Ministry of Transportation and Communications; the farm and forest roads were gained from the Soil and Water Conservation Bureau and the Forestry Bureau, respectively. Road density and spatial dispersion index (SDI) could be calculated using road maps coupled with watershed subdivision. River network obtained from the Institute of Transportation could be used to extract areas of river concave and headward watershed. The information of geologically sensitive area announced by the Central Geological Survey was applied to extract the spatial distribution of dip slope.

Analysis unit
A DEM with resolution of 5 m £ 5 m was employed as grid resolution for terrain analysis in the models, while the analysis unit for the indices of road development and headward watershed are based on watershed division. Generally, a larger number of subunit delineated will result in a higher environmentally homogeneous subdivision, and is useful for treatment/management works in watershed landslide restoration. However, the more subdivisions a watershed has, the less accurate in descriptions of landslide characteristics it will have, and an excessive number of subdivisions could cause heavy loads in treatment/management works due to limited topographic factors and/or resources, respectively (Lin et al. 2016). Basically, the size of an analysed watershed should be defined firstly. An optimal watershed scale with the threshold of 1000 ha was specified first because the rational formula can only be applied to estimate surface run-off in a small watershed. The spatial distribution of the specific watershed was delineated in Figure 3(A), and then the specific watershed was further subdivided into right bank, left bank and headward (if existing) areas for the sake of homogeneity (Figure 3(B)). Area of subdivision of the specific watershed could be regarded as the analysis unit to evaluate the potential large-scale landslide.

Risk analysis
Risk analysis is the probability of the occurrence of a certain hazard event (Cao et al. 2010). A quantitative risk model was provided for disaster analysis of the risk issues by the UNDRO (1979). The model could be expressed as follows: The indices used in the study were all normalized individually at the range of 0-1 to display the trend of impacts and/or scale for better comparison of the indices. The normalization equation could be written as follows: where nA is the normalized value of the index, A represents original value, A min and A max indicate the minimum and maximum value of the index, respectively. Equation (3) was employed as a function of linear transformation to derive the reversal index, and iA is the normalized value of the index after linear reversal transformation:

Index of hazard
Hazard degree used in this study implicitly means the intensity of external force that could cause watershed landslide. Maximum daily rainfall of Typhoon Morakot (Figure 4(A)), an extreme event with about 200-year return period, resulted in severe landslides in southern Taiwan and was used as the index of hazard (Figure 4(B)) to verify this conceptual model.

Index of vulnerability
There must be some assumptions and a combined indicator to effectively delineate environmental weakness due to heterogeneity and dynamic changes. In order to derive a composite indicator to service as vulnerability, energy release theory that depicts accidents in terms of energy transference is the main consideration in this study. A synthesis indicator integrating indices of return period, road development and green deterioration is used as the vulnerability index which finally was modified by the slope coverage to adjust the weight for excluding the flat areas and increasing impact of steep areas. Conceptually, heavy rainfall could displace soil particles and cause soil erosion process. Thus, a site located at a high-intensity rainfall region will have a higher resistance to torrential-induced disasters due to its ability to withstand heavier rainfall. In addition, an event with higher occurrence probability will lead to a smaller landslide scale and vice versa (Moon et al 2005). Since concepts of return period do reflect the natural pulsation, it could be adopted as the environmental vulnerability index for evaluating the disaster probability of an area susceptible to external forces. Road development, the second vulnerability index, could intercept overland flow to form channel flow and cause slope failure due to run-off concentration under extreme rainfall ( Figure 5). In Taiwan, the design return period of road side ditch on the slopeland is usually not more than 25 years. During high return period rainfall event, the excessive drainage water could overflow the road side ditches at the sharp curve (valley or ridge), which is likely to erode the slope toe and cause slope failure Lin et al. 2013;Lin and Lan 2015). In addition to road density, the exposure of constructed roads to internal forces should also be concerned for assessing the area with man-made vulnerability. For deriving the third vulnerability index, areas of vegetation deterioration are usually vulnerable to erosion and could cause slope failure (Collison and Anderson 1996;G€ okceoglu and Aksoy 1996). The index of green deterioration, which is calculated from satellite-based vegetation index, could be applied to show the degree of vulnerability.
2.4.2.1. Index of return period. An event with certain return period could be directly calculated from the frequency analysis using pre-event historic records of meteorological station. The concept of return period is used in designing the life cycle of engineering structure, and the return period for most protected objects is usually 50-100 years and seldom over 200 years in Taiwan. Therefore, the occurrence probability of rainfall intensity less than 1/50 annually might be defined as extreme event which could cause engineering structure failure (Soil and Water Conservation Bureau 2014). As for landslide, there is no concrete data to judge the life cycle of a slope land due to uncertainty. One can adopt an adequate return period for the simulation of an extreme event. A 200-year maximum daily rainfall was selected in this study as the index of return period for consistently simulating the event of Typhoon Morakot (Figure 6(A)). The index of return period was then under the process of linear reversal transformation and normalization to display the variations of environmental vulnerability ( Figure 6(B)).

a. Road density
Road could intercept slopeland run-off, and the intercepted volume of run-off is usually much larger than the drainage capacity of the road during extreme rainfall. As a result, the drained water, especially at the turning section of the road, is susceptible to being spilled and down-slope failure . Hence, road density of subdivision could be one of the vulnerability indices to assess landslide risk. Road density was defined as the ratio of road length in a subdivision to the subdivision area. Figure 7(A) is the spatial distribution of normalized road density.

b. Spatial dispersion of a road
In addition to density, most roads were constructed in a trans-boundary way. The spatial distribution of a road should be considered as the exposure index to modify the vulnerability of road density. The SDI developed by Weng and Tsai (2006) was used to describe the exposure of a road. The formula of SDI is expressed as in the following equation: Equation (4) is based on the statistical index of deviation from the mean; d i;j ð Þ is dispersion degree of two measured objects, a i % + a j % is weight of d i;j ð Þ , SS and i 6 ¼ j is calculated degree of dispersion for one object to the other objects, 1 2 nÀ1 ð Þ is number of repeated calculation and 1 ffiffi ffi A p is a correction for size of statistic unit. The more the objects concentrated, the smaller the SDI derived (Figure 7(B)). Finally, the index of road development (Figure 7(C)) could be calculated from the multiplication of normalized road density and SDI. 2.4.2.3. Index of green deterioration. Satellite imageries of dry and wet seasons were selected to extract the variations of NDVI (Normalized Difference Vegetation Index) (Figure 8(A,B)). Since NDVI value is between -1 and +1, the NDVI has been linear transformed to the value range of 0-1 in order to overlay NDVI with other indices in this study. The NDVI calculation (Rouse et al. 1974) can be obtained as follows: where NIR is the reflectance of the near-infrared waveband and R is the reflectance of the visible red waveband of the satellite radiometer. The higher the NDVI value, the better the photosynthesis activity and the greater the vegetation covers. The average linear transformation of NDVI calculated from dry and wet seasons should be reversely linear transformed to the index of green deterioration to indicate environmental vulnerability (Figure 8(C)).

Vulnerability calculation.
A composite indicator integrated the three indices using overlay calculation and then normalized them to be the unmodified index of vulnerability (Figure 9(A)), which should be further modified under the process of multiplying by the normalized gradient index (Figure 9(B)) to derive the index of vulnerability (Figure 9(C)). The necessity of modification with gradient is to remove the possibility of flat-failure and to reflect that the area with steeper slope will become more unstable and increase additional effect in vulnerability.

Risk calculation
Landslide risk evaluated just before the event of Typhoon Morakot could be derived from multiplying each other using hazard (Figure 4(B)) and vulnerability (Figure 9(C)). Results show that the higher risks of landslides were located at the southern and central mountain areas (Figure 10).

Index of potential landslide scale
Several devastating large-scale landslides such as dip-slope failure at national highway no. 3, the collapses of Hsiao-lin village, Cao-ling and Chiu-fen-er-shan have occurred over the past few decades in Taiwan. Sites of river concave (Figure 11(A)), headward erosion (Figure 11(B)) and dip slope (Figure 11(C)) are most likely to have large-scale landslides during extreme rainfall. The driving force could be attributed to footing failure caused by human or natural action. Two main mechanisms (lateral/vertical erosion and headward erosion) of footing failure under exogenetic process  (Lee and Lin 2015) should be further explored to better understand the sites of potential large-scale landslide.

Index of river concave
The concave river bank is usually susceptible to scouring and may cause footing failure due to centrifugal force (L evy et al. 2012). Generally, compared to sites without concave, places of river concave are usually protected by the natural bedrock and/or artificial strengthened bank which could result in less failure frequency and increasing collapse scale once failure. The index of river concave could be established using variations of tangent slope along a river (Figure 12(A)). The larger the  gradient of curvature, the greater scouring it experienced. Variation of tangent slope (u) could be expressed as Equation (6), and the spatial distribution of the index can then be calculated (Figure 12(B)):

Index of headward erosion
Headward erosion is the erosion at the origin of a stream channel, which causes the origin to move back away from the direction of the stream flow, and thus causes the stream channel to lengthen (Judson and Kauffman 1990). The principles of diminishing returns (Case and Fair 1999) used in the field of economics should be considered for understanding the functionality of plants in slope stability. Generally, vegetation is beneficial at the initial succession; while accompanying the plant growth, the huge biomass increasing loads could result in slope failure at certain sites especially those which occurred at the dip slope due to effects of strong wind and high infiltration during typhoon events (Hung et al. 2012). Besides, the well growth natural forest located at the headwater areas with low soil erosion will ordinarily trap sediments at the slopeland and increase the potential of headward erosion. Since the sites with such characteristics are vulnerable to slope instability and will yield large-scale landslide under the occurrence of high return period event, a catastrophic disaster of thirty thousand hectare landslide and one billion cubic meters of sediment occurred at forestry compartment in southern Taiwan during Typhoon Morakot can be an unforgettable experience. Since headward erosion can easily occur at the origin of river course, stream order developed by Strahler (1952) was applied to delineate the probability of headward erosion. A subdivision of watershed using rules of stream order was illustrated (Figure 13(A)), and the normalized potential headward erosion (Figure 13(B)) could be calculated as the index of headward erosion.

Index of dip slope
A dip slope is a topographic surface which slopes in the same direction, and often by the same amount, as the true dip or apparent dip of the underlying strata (Jackson et al. 2005;Allaby 2008). Dip slope consisting of sandstone-shale alternation is very common in Taiwan. Sandstone with characteristics of thicker, harder and more permeable properties compared to shale could resist weathering, and slope of sand interbedded shale is susceptible to friction drop due to the difference in erosion and/or percolation. As a result, a dip-slope failure could be triggered and cause large amounts of debris to remove its footing support (Wang et al. 2013), and could be an index of largescale landslide. Since failure location affects the amount of landslide scale, topographic wetness index (TWI) shown in Figure 14(A) computed from Equation (7) (Beven and Kirkby 1979) was employed to modify the sites of dip slope (Figure 14(B)) for deriving the index of dip slope (Figure 14(C)): TWI ¼ ln As tan u where As is the upslope area per unit contour length and tan u is the local gradient.

Landslide scale
A comprehensive indicator using overlay computation from indices of river concave, headward watershed and dip slope could be provided to effectively reflect hotspots of potential large-scale landslide (Figure 15(A)). However, the indicator needs to be modified by a computation of multiplying with vegetation index (Figure 15(B)) for adjusting the potential landslide scale (Figure 15(C)) because the amount of soil in a slope could be indirectly detected from vegetation status.

Danger grade
For testing the accuracy of the models, the estimated risk/scale of each subdivision was classified into three categories (high, medium and low) using K-means cluster analysis. A danger grade with five levels (very low, low, medium, high and very high) computed from the array of the risk-scale matrix ( Figure 16) could be established. The collapse ratio and average of landslide ratio for each level were extracted from the map of landslide triggered by Typhoon Morakot (Figure 17) to explore the relationship between danger grade and actual landslide risk/scale of the delineated subdivisions. Collapse ratio is defined as the number of subdivision with collapse to the number of subdivision, and landslide ratio is defined as the area ratio of landslide to subdivision. Figure 18 shows that there exists a parabolic relationship between collapse risk and failure scale which was derived from 3601 subdivisions of Taiwan during Typhoon Morakot. It reveals that the relationship follows the law of diminishing returns; when the collapse risk is greater than 0.7, it will be less likely to have landslide scale occurrence due to frequent disturbance.

Verification of risk analysis
The relationship between risk index and collapse ratio of the subdivisions with certain risk range (Figure 19(A)) shows highly significant correlation (R 2 = 0.9914 ÃÃÃ ). This implies that the model developed in this study could be applied to estimate the collapse ratio of interested subdivisions under the calculation of risk analysis. Moreover, the landslide ratio of the subdivisions (Figure 19 (B)) also depicts that there is a highly significant correlation with risk index (R 2 = 0.917 ÃÃÃ ); and once the risk index is greater than 0.7, the landslide ratio decreases gradually. This might be due to over slope modification at steeper areas in vulnerability analysis.   Figure 20(A,B) indicates that the index of landslide scale shows significant correlation with collapse and/or landslide ratio (R 2 = 0.7801 ÃÃÃ /R 2 = 0.677 Ã ), respectively. This implies that the subdivisions with higher index of landslide scale could cause the risk of greater collapse/landslide ratio, and some subdivisions with high index of landslide scale but showing less collapse ratio need to be concerned of catastrophic disaster for debris accumulation. Figure 21(A,B) depicts that danger grade is showing highly significant correlation with collapse ratio (R 2 = 0.9934 ÃÃÃ ) and significant correlation with landslide ratio (R 2 = 0.9544 Ã ). The more dangerous    a subdivision is, the higher collapse/landslide ratio it will have. The classification of danger grade could be effectively applied in screening and management of hotspots of potential large-scale landslide.

Conclusion
The risk-based models have been established for potential large-scale landslide monitoring and management. The sites of road development, headward erosion, river concave and dip slope adopted in this study can be classified as the places without update necessity, otherwise the event rainfall, vegetation status and return period should be re-evaluated and dynamically updated after an event. Collapse risk could be substituted by vulnerability, and could be used with failure scale to judge the danger degree of the subdivision due to unpredictable hazard at ordinary time. Danger grade can be coupled with real-time rainfall derived from meteorological stations for dynamic landslide risk assessment and potential large-scale landslide monitoring during extreme rainfall. The models of danger grade could also be used to understand the orientation relationship of potential landslide hotspots and protected objects for further disaster management.