Statistical approach to storm event-induced landslides susceptibility

For the interpretation of the storm event-induced landslide distribution for an area, deterministic methods are frequently used, while a region’s landslide susceptibility is commonly predicted via a statistical approach based upon multi-temporal landslide inventories and environmental factors. In this study we try to use an event-based landslide inventory, a set of environmental variables and a triggering factor to build a susceptibility model for a region which is solved using a multivariate statistical method. Data for shallow landslides triggered by the 2002 typhoon, Toraji, in central western Taiwan, are selected for training the susceptibility model. The maximum rainfall intensity of the storm event is found to be an effective triggering factor affecting the landslide distribution and this is used in the model. The model is built for the Kuohsing region and validated using data from the neighboring Tungshih area and a subsequent storm event – the 2004 typhoon, Mindulle, which affected both the Kuohsing and the Tungshih areas. The results show that we can accurately interpret the landslide distribution in the study area and predict the occurrence of landslides in the neighboring region in a subsequent typhoon event. The advantage of this statistical method is that neither hydrological data, strength data, failure depth, nor a long-period landslide inventory is needed as input.

In landslide susceptibility analysis (LSA), it has been most common to use a statistical approach, where landslide inventories and causative factors are utilized to build a susceptibility model for mapping or delineating areas prone to landslides.Many different methods and techniques for assessing landslide hazards have been proposed and tested.These have already been systematically compared and their advantages and limitations outlined in the literature (Carrara, 1983;Varnes, 1984;Carrara et al., 1995;Hutchinson, 1995;Mantovani et al., 1996;Aleotti and Chowdury, 1999;Chung and Fabbri, 1999;Guzzetti et al., 1999;Wang et al., 2005;Chung, 2006;van Westen et al., 2006).Most of these approaches require multi-temporal landslide inventories, so that the susceptibility model can predict landslide occurrence for a given time period (Guzzetti et al., 1999).In previous statistical models the triggering factors have seldom been emphasized.
In recent years, Dai and Lee (2003) used the rolling 24h rainfall as an independent variable for the building of a storm-induced shallow landslide probabilistic model.Chang et al. (2007) used the maximum 3-h rainfall and rainfall duration in their logit model to model the rainfall conditions critical for triggering landslides.Dahal et al. (2008) used extreme 1-day rainfall records in their weights-of-evidence model to predict the rainfall-induced landslide hazard.It has become a trend to incorporate rainfall as an independent variable into storm event-induced landslide modeling, but this method still needs validated and further defined for different environments.
Multi-temporal landslide inventories are not yet available for Taiwan.However, several typhoons (tropical cyclones) strike Taiwan every year, which has given us the opportunity to map an event-based landslide inventory for each major typhoon.In this study we use an event-based landslide inventory, considering rainfall as an independent variable in our multivariate statistical susceptibility model.We carefully compare pre-event and post-event landslide inventories to produce an event-based landslide inventory.We select data for the shallow landslides triggered by the 2002 typhoon Toraji in the Kuohsing region in central western Taiwan to train the susceptibility model.The susceptibility model is then validated using data from the neighboring Tungshih region, as well as data from a subsequent storm event -the 2004 typhoon, Mindulle.The advantages and shortcomings, as well as practical use of this statistical method are discussed.

Regional settings and storm events
The study area comprises the Kuohsing and Tungshih quadrangles (coincide with the 1 to 50 000 scale map of Taiwan; 705.3 km 2 and 703.9 km 2 in area), which are located in central western Taiwan (Fig. 1).Kuohsing is east of the Taichung Basin where metropolitan Taichung is located, and west of the Central Range.Tungshih is located immediately north of Kuohsing.Geologically, the area falls partly in the Western Foothills geologic province (western part) and partly in the Hsueshan Range sub-province of the Central Range geologic province (eastern part).The Western Foothills are characterized by fold-and-thrust Neogene sedimentary strata; the Hsueshan Range is typically characterized by a Paleogene slate belt of argillite and quartzitic sandstone (Ho, 1975).The Shuilikeng Fault forms a boundary between these two provinces (see Fig. 1).
Geomorphologically, the study area can be divided into two types of terrains by the Hsuangtung Fault: western hilly terrain and eastern mountainous terrain.The elevation in the hilly terrain is generally less than 500 m with rocks consisting of weakly cemented Pliocene and Pleistocene mudstone, sandstone and conglomerate.The elevation in the mountainous terrain is generally greater than 500 m with many peaks exceeding 1000 m with rocks consisting of better-indurate Eocene, Oligocene and Miocene, sandstone, shale, argillite and quartzite.The geomorphic appearance in the hilly terrain is more fragmented with short slopes, while the mountainous terrain is solider and has longer slopes.Typical slope gradients in the hilly terrain are between 10 • and 30 • , with a mode of about 16 • , while in the mountainous terrain, typical slope gradients are between 10 • and 45 • , with a mode of about 26 • .
The study area has a subtropical climate with average annual precipitation of about 2400 mm, and on average, 159 rainy days per year.Approximately three typhoons strike this area each year, mostly between July and October.The heaviest rainfall typically occurs in June, July and August; with a monthly average during that period of about 400 mm.In the dry season, October through March, the monthly average precipitation is about 80 mm with a monthly average of about 8 rainy days.The sloped lands are usually green and covered with vegetation.Gentle or moderate slope with gradients less than 30% are usually cultivated while steep slopes tend to be covered by bushes and shrubbery.Underneath the hill surface, the slopes in the study region are for the most part mantled by permeable colluvium soils.Although the slopes in the study area are green, shallow landslides are common during a major earthquake or a typhoon event.
A disastrous earthquake -the M w 7.6 Chi-Chi earthquake -occurred on 21 September 1999 in central western Taiwan (Ma et al., 1999;Kao and Chen, 2000).A huge thrust ruptured along the Chelungpu fault about 3 km west of the Kuohsing quadrangle, extending into the southwest corner of the Tungshih quadrangle.The Kuohsing quadrangle and most of the Tungshih quadrangle are located on the hanging wall of the thrust fault.These areas suffered from severe shaking during the main shock and aftershocks triggering 9272 large landslides (with areas greater than 625 m 2 ), with a total area of 127.8 km 2 (Liao and Lee, 2000).
After the 1999 Chi-Chi earthquake disaster, two major typhoons struck Taiwan, reactivating many landslides, especially in central western Taiwan, where slopes had already been damaged or destabilized by the earthquake (Liao et al., 2002).Typhoon Toraji passed through Taiwan from 28 July to 31 July 2001.Typhoon Mindulle passed over Taiwan from 28 June to 3 July 2004.The properties of the rainfall and losses due to these two typhoons are listed in Table 1.Storminduced landslides became a big problem after the 1999 Chi-Chi earthquake.Debris flows and sediment transport became a major problem for governmental authorities to cope with.Although typhoon Mindulle brought a maximum total precipitation of 2125 mm to Taiwan, it caused less loss than Toraji did.This may have been due to a more developed hazard mitigation program and evacuation of people by the government before the typhoon hazard occurred.

Methodology
This study comprises one element of an investigation of landslides in Taiwan carried out by the Central Geological Survey (CGS) of Taiwan (Lee et al., 2005).The overall aim of the CGS investigation is to produce a set of landslide susceptibility maps (scale: 1/25 000) for all of Taiwan.Since this is a target-oriented project, the methodology for producing the landslide susceptibility maps must be objective and simple enough to meet the needs of future work.The specific  This third problem was solved by selecting a mature analytical method with the whole of Taiwan being divided into several terrain zones.An LSA model was trained for each specific terrain (Lee et al., 2005).To solve the first and second problems, however, we were forced to consider using an event-based landslide inventory and event-based landslide susceptibility analysis (EB-LSA).We tested the EB-LSA for different regions and different events in Taiwan before publishing an EB-LSA for earthquake-induced landslides in the Kuohsing area (Lee et al., 2008).

Methods and working procedure
The methods and working procedure utilized in the present study of storm-induced landslides generally follow those of Lee et al. (2008).The first step includes image and data collection, after which an event-based landslide inventory is established (in this case for a storm event).In parallel with this, the causative factors of the landslides are processed and the triggering factors determined.These factors are then statistically tested, and the effective factors selected for susceptibility analysis.Each selected factor is rated, and their weighting analyzed.
Discriminant analysis allows us to determine the maximum difference for each factor between the landslide group and the non-landslide group, as well as the apparent weights of the factors.A linear weighted summation of all factors is used to calculate the landslide susceptibility index (LSI) for each grid point.The LSIs are used to establish a landslide ratio to LSI curve and determine the spatial probability of landslide at each grid point.The landslide ratio used here is the ratio of landslide pixels to total pixels in an LSI interval (Lee et al., 2005(Lee et al., , 2008)), called the proportion of landslide cells in Jibson et al. (2000).The spatial probability of landslides is then used for landslide susceptibility mapping.For a more detailed description of EB-LSA, please refer to Lee et al. (2008).

Selection of causative factors
There are more than fifty different landslide-related factors commonly used (both in Taiwan and worldwide) for LSA (Lin, 2003).In the present storm event-induced landslide study, we selected fourteen of the most frequently used, based on data abundance and availability.These causative factors were lithology, slope gradient, slope aspect, terrain roughness, slope roughness, total curvature (Wilson and Gallant, 2000), local slope height, total slope height, topographic index (Kirkby, 1975), distance from a road, distance from a fault, distance from a river head, distance from a river bend, and the normalized differential vegetation index (NDVI, Paruelo et al., 2004).
These factors were further tested, including the normality of each factor, correlation coefficient between any two factors, and calculation of standardized differences (Davis, 2002) between the landslide group and non-landslide group for each factor.A final selection of effective factors was decided based upon the evaluation and test results.

Selection of a triggering factor
It is well known that rainfall plays an important role in triggering landslides during a storm event.Most shallow slope failures during or after a rainstorm are triggered by an increase in pore pressure and a corresponding reduction in effective stress in the soil (Dai and Lee, 2003).Shallow soils are often underlain by a relatively low permeability layer or an impeding layer.The raising of the soil-water level above the impeding layer during a rain storm may be influenced by antecedent rainfall, the condition of the soil surface and its vegetation, the soil properties, such as its porosity and hydraulic conductivity, and the upslope drainage area (Kirkby, 1975;Chow et al., 1988).Because typhoon rainfall in the Taiwan region is usually heavy, the antecedent rainfall may be neglected as just having a minor effect.Vegetation cover is reflected by the NDVI.Instead of slope length and upslope area, we used a total-slope-height factor and the topographic index.The soils are loose and permeable colluvium.We may assume that the colluvium soils within a lithologic unit will have similar properties but that there will be variations between different lithologic units.Therefore, the lithology factor is also related to hillslope hydrology.
To include the above-mentioned causative factors in the statistical susceptibility model may connect rainfall to the triggering of slope instability, and help in our interpretation of the physical meaning.Rainfall intensity and duration or cumulative rainfall are the most commonly used factors to delineate landslide occurrence (Zezere and Rodrigues, 2002;Guzzetti et al., 2007).In this study, the maximum hourly rainfall (maximum rainfall intensity), the rolling 24-h rainfall, and the total rainfall of the storm event were considered as candidate triggering factors.These factors were statistically tested, and evaluated to make a final decision.

Assessment of model performance and validation
The error matrix (Stehman, 1997), the receiver operating characteristic (ROC) (Swets, 1988), and the prediction rate curve (Chung and Fabbri, 2003) are the most commonly used methods to assess a model's performance in landslide and other types of studies.For both the error matrix and the ROC, the classification of true positives and false positives are needed; a probability value of 0.5 is used in the logit model to determine whether the model has made a correct prediction (>0.5) or not (<0.5); a discriminant index λ 0 is used in the discriminant analysis to determine whether the model has made a correct prediction (>λ 0 ) or not (<λ 0 ).This clear cut boundary is deemed to be not as friendly as the prediction rate curve, which needs only successive classes to reflect the landslide potential.Furthermore, we used a discriminant function to develop the LSIs; λ 0 was not used in the discriminant model.Therefore, the prediction rate curve was selected in the present study.
Landslide data used to establish the model were first grouped into several classes according to their LSIs.The number of landslide pixels in each class was then divided by the total number of pixels in that class, and a cumulative curve was plotted.The area under the curve (AUC) is between 0 and 1; a higher value indicates a higher prediction rate, whereas a value near 0.5 means the prediction is no better than a random guess (Chung and Fabbri, 2003).In the assessment, we classify AUC 0.9 as excellent, 0.9>AUC 0.8 as good, 0.8>AUC 0.7 as fair, 0.7>AUC 0.6 as poor, AUC<0.6 as very poor.
The prediction rate curve may be used to assess the prediction performance of the model, or used to validate a model using different data sets.When the prediction rate curve is used to assess model performance, we use the same data set as for building the model to compute a success rate curve (Chung and Fabbri, 1999).The success rate curve is used to indicate how well the model fits the data.

Data acquisition and processing
The basic data utilized in this study included a 40 m×40 m grid digital elevation model (DEM), SPOT images, 1/5000 photo-based contour maps, 1/50 000 geologic maps, and hourly rainfall data.The DEMs were collected by the Aerial Survey Office of Taiwan's Forestry Bureau.They were transferred to a color-shaded image and were visually checked.Defects were replaced by re-digitizing from a 1/5000 scale photo-based contour map.Other abnormal points were corrected using a median filter.Finally the DEMs were interpolated to a 20 m×20 m grid cells using cubic spline interpolation.
Selected SPOT images taken before and after the typhoon events are shown in Table 2.All SPOT images were received, processed and rectified by the Center for Space and Remote Sensing Research, National Central University, Taiwan.Both multi-spectral (XS) and panchromatic (PAN) images were used.A fusing technique (Liu, 2000) was utilized to produce a higher resolution false-color composite image to facilitate landslide recognition.The pixel resolution after fusing was 6.25 m.
1/50 000 geological maps were collected from the CGS.Each map was overlaid with a shaded DEM and visually inspected in a Geographic Information System (GIS).Some abnormal boundaries, mostly associated with alluvial and terrace deposits, were corrected.The ERDAS IMAGINE system (ERDAS, 1997) was used to transform the geologic vector map to a raster image of 20 m×20 m grid cells.
The hourly rainfall data in and around the study region were collected from the Central Weather Bureau and the Water Resources Agency, Taiwan.These data were first plotted and visually inspected to compare individual records for consistency with neighboring gauge stations; abnormal data was deleted.Rainfall data were finally interpolated into 20 m×20 m grid cells data as will be explained in Sect.4.3.All later processing and analysis for each susceptibility factor and for the EB-LSA are based on the 20 m×20 m grid-cells unit.

Event-based landslide inventories
An event-based landslide inventory is difficult to extract using aerial photographs alone.Complete sets of photographs taken just before and soon after an event are rarely available.An event-based landslide inventory, extracted from satellite imagery, is often more practical and valid.A major disadvantage of using satellite imagery for extracting landslides is that the spatial resolution may be less than that of an aerial photograph, so that some small landslides may be missed.However, not all of the landslide data are needed to establish   a susceptibility model (Weirich and Blesius, 2007); missing some small landslides is not an important issue.SPOT images taken before and after the 2 typhoon events were selected for each quadrangle (Tables 2 and 3) from which two event-based landslide inventories for each of the two quadrangles were prepared.A pre-event landslide inventory was prepared for each event from a SPOT image taken before the event; a post-event landslide inventory was prepared from an image taken after that event, and finally an event-based landslide inventory was derived by comparing the pre-event and post-event inventories.An event-induced landslide could be absent from the pre-event landslide inventory, or present in both inventories (re-activated landslide).Landslides found in both inventories were examined very carefully for changes in tone and/or enlargement of extent (Pan et al., 2004).
We found that landslide deposit areas should be included in the above-mentioned landslide mapping procedure.Since potential landslide sources are of interest in susceptibility analysis, only source areas may be used to train the susceptibility model.Therefore, we had to differentiate between source areas and actual deposit areas.Landslide deposits were identified by comparing the GIS landslide layer with the 1/5000 scale photo-based contour map.The slope angle or concentration of contour lines was used to differentiate deposits from sources.
Each landslide inventory was completed using a geographic information system (GIS).The attributes include detailed descriptions of the date/event, and the size and type of each landslide object.Four event-based landslide inventories were developed: Toraji event -Kuohsing, Toraji event -Tungshih, Mindulle event -Kuohsing, and Mindulle event -Tungshih.The spatial distribution of landslides triggered    by the two events is summarized in Fig. 2, and listed in Table 4.Only source areas of shallow landslides (including rock falls), were used in the susceptibility analysis.Deep seated slides, rock avalanches (located outside the present study area) and debris flows were excluded in the present study.

Landslide causative factors
The data set for each terrain was divided into a landslide group and a non-landslide group.Tests of normality for each factor, the calculation of the correlation coefficient between any two factors, and the calculation of standardized differences between the landslide group and non-landslide group for each factor were performed in order to select effective factors for the discriminant analysis.Eight causative factors -lithology, slope gradient, slope aspect, terrain roughness, slope roughness, total curvature, total slope height, and NDVI -were judged to be effective factors to be used in the discriminant analysis.Six of them have also been used for earthquake-induced landslide susceptibility analysis by Lee et al. (2008), whereas the NDVI factor and the total-slopeheight factor are additional factors.
The NDVI is an environmental factor indicating the abundance of vegetation cover at a specific pixel point, and may be used as an indicator at land-use.A higher NDVI value indicates denser vegetation.Denser vegetation increases precipitation interception, and decreases the amount of infiltration, whereas barer land areas will have higher infiltration rates, which are associated with higher soil water levels, which in turn contributes to slope failure.To obtain accurate land-use information for establishing the susceptibility model, the NDVI should be calculated from an image taken prior to a storm event.
The total-slope-height for a given point on a slope is defined as the height of the upper slope above the point (height from point to crest) plus the slope height below the point (height from toe to point) (F+G in Fig. 3).The total-slopeheight factor may be physically related to the magnitude of the stress and the pore-water pressure in the lower slope; for long slopes the surface and subsurface water is more likely to be concentrated in the lower slope, causing instability.It has been observed that storm-induced landslides tend to occur relatively low on the slopes (Chang and Hsu, 2004).
The topographic index, which is defined as the area draining through a point from upslope divided by the local slope gradient at that point (Kirkby, 1975), is an important hydrological factor related to landslides.However, it was found to be ineffective for discriminating between the landslide and non-landslide group data.This unusual result may be due to the weakening of ridges in this study area by the coseismic shaking of the Chi-Chi earthquake, so that many landslides occurred on the upper-slopes close to the ridge tops.The occurrence of these landslides was due to weak soil strength and open cracks on the slope which allowed overland flows to seep into the soil.
All eight causative factors, except for total-slope-height, were processed by the ERDAS IMAGINE system (ERDAS, 1997).The total-slope-height values were calculated from a series of along-slope sections by a FORTRAN program   developed for this study.The spatial distribution of the totalslope-height and NDVI factors are shown in Fig. 4; those for the slope gradient, slope aspect, topographic roughness, slope roughness, total curvature and lithology are similar to those in Fig. 4 in Lee et al. (2008) so it is not necessary to repeat them here.The descriptive statistics for the causative factors are listed in Tables 5 and 6.Plots of the frequency distributions of the landslide and non-landslide groups in the different terrains (hilly terrain and mountainous terrain) are shown in the left and middle columns of Fig. 5.The righthand column shows the landslide ratios with respect to the factor values.There is a positive correlation between the factor and corresponding landslide ratio for both hilly terrain and mountainous terrain for slope gradient, topographic roughness, slope roughness, total curvature and total-slopeheight (Fig. 5c, i, l, o, r).A straight line may be fitted between the lower and upper bounds.There is a negative correlation between NDVI and landslide ratio for both hilly terrain and mountainous terrain (Fig. 5u).A straight line may also be fitted between the lower and upper bounds.All plots were visually inspected and evaluated.The landslide ratio of the slope aspect shows a sinusoidal curve (Fig. 5f), and a sinusoid is fitted.
A factor was assigned a rating according to its landslide ratio.The rated values were then normalized to be between 0 and 1, with 0 being a factor with a value less than or equal to its lower threshold and 1 being a factor value greater than or equal to its upper bound.Although the lithology factor data are categorical in nature, a landslide ratio can also be found for each lithologic unit, and a normalized score between 0 and 1 assigned.

Landslide triggering factors
In this study we used data from 68 rain gauge stations located within the 12 map quadrangles in the study area: 5 stations in the Kuohsing quadrangle and 11 in the Tungshih quadrangle (Fig. 6).Each rainfall record plotted was checked visually, and abnormal data were not used.Maximum rainfall intensity and total rainfall values were calculated station by station and then interpolated for each grid point using the ordinary Kriging method (Goovaerts, 1997).Because the duration of rainfall at the rain gauge stations in the Kuohsing quadrangle ranged from 19 to 23 h, neither maximum daily rainfall nor rolling 24-h rainfall were considered.
The frequency distribution of the landslide and nonlandslide groups as well as landslide ratio for these factors were tested, plotted, and the results examined (Fig. 7a-f).Both total rainfall and the maximum rainfall intensity factors were tested statistically for effectiveness at discriminating between the landslide group and the non-landslide group, based on which the maximum rainfall intensity was selected as the triggering factor for the storm-induced LSA.
The maximum rainfall intensity was also assigned a rating according to the landslide ratio and considering a lower threshold, similar to that for a causative factor.We also  needed to assign a slightly higher upper bound than the actual value encountered in the event.An EB-LSA model with a credible upper bound for the triggering factor would be of benefit when using this model for prediction, as higher rainfall intensity may be encountered in some scenarios.

Results and evaluation
Data in the Kuohsing quadrangle for shallow landslides triggered by typhoon Toraji were used for the EB-LSA.All data sets from the landslide group (5575 pixels in hilly terrain; 11 540 pixels in mountainous terrain) and a randomly selected non-landslide data set of similar size were used in the discriminant analysis.Gentle slopes, i.e., a slope gradient of less than 10%, were regarded as stable, and were not included in the discriminant analysis.These areas were also not considered in the evaluation of the results or the calculation of the success rate or prediction rate curves.
The results of the analysis included a coefficient (apparent weight) for each factor for both the hilly and mountainous terrain (Table 7).Among the 9 factors used, the slope gradient factor has the highest coefficient and a large percentage of the weighting.The environmental factor -NDVI, also carries a large percentage of the weighting.As compared with the above-mentioned two factors, the rainfall intensity factor is not weighted as highly and needs further discussion.
These various weights were used to calculate the LSI for each grid point, as with Eq. ( 5) in Lee et al. (2008).LSIs were used to calculate the landslide ratio for each LSI class.The spatial probability of landslide occurrence is indicated Nat.Hazards Earth Syst.Sci., 8, 941-960, 2008 www.nat-hazards-earth-syst-sci.net/8/941/2008/  by the relation between the probability of failure (landslide ratio) and the LSI (Fig. 8).The probability of failure curves does show a strong trend of increasing probability of landslide occurrence with increasing susceptibility values.The relation generally follows a Weibull distribution (Lee, 2006;Lee et al., 2008).There are some fluctuations in the landslide ratio when the susceptibility value is higher, simply because there is less data available.The spatial probability of a landslide can then be used to map the susceptibility classes, as shown in Fig. 9.

Success rate for the Kuohsing quadrangle
We now compare the distribution of actual landslides triggered by typhoon Toraji with the event-based susceptibility map (Fig. 9).This comparison shows that the landslide pattern generally agrees with the pattern of the high susceptibility classes.We further examine how the results fit the data, using the prediction rate curve method (Chung and Fabbri, 2003).Landslide data used to establish the model were grouped into several classes according to their LSIs, the number of landslide pixels in each class was divided by the total number of pixels in that class, and a cumulative curve was plotted.Since it is the full set of data from the landslide group used to train the model and calculate the rates, we essentially calculate a success rate.The success rate curves for the two terrains are plotted in Fig. 10.The AUC values were also calculated and are shown in Fig. 10.The results for the hilly terrain are excellent (AUC=0.9343)while those for the mountainous terrain (AUC=0.8859)(in the Kuohsing quadrangle for the Toraji typhoon event).The exclusion of slopes less than 10% from these calculations, as noted above, may reduce the bias in the success rate calculation, leading to a more conservative result than if these slopes were included.

Validation in the Tungshih quadrangle
The same data sources and the same procedure were used to validate the susceptibility model.We processed factors from data from the same sources for the Tungshih quadrangle, and used the factor weights from the Kuohsing quadrangle to calculate the LSI for each grid point in the Tungshih quadrangle.The LSIs were then transferred to indicate probability of failure using the equations shown in Fig. 8.The map of the susceptibility classes for the grid points in this quadrangle is shown in Fig. 11.
The results of the prediction rate calculation are shown in Fig. 10.Again, flat areas, with slopes of less than 10%, were not considered in the calculation of the prediction rate curve; the results are again conservative.The prediction rate for hilly terrain (AUC=0.8171)was less than that for the Kuohsing quadrangle, whereas for mountainous terrain, the prediction rate for both quadrangles was similar (AUC=0.8859vs. 0.8903).This is a fairly good result.The lower prediction rate for the hilly terrain may be explained as due to the most widespread rock type -the Houyenshan conglomerate, where numerous landslides were triggered by the Chi-Chi earthquake and also reactivated during the typhoon Toraji, not being present in the Tungshih quadrangle.

Validation using data from a subsequent event
We looked at landslides induced by Typhoon Mindulle for additional validation.The same data sources and procedure were used.Most factors are time-invariant, but this invariance cannot be extended to environmental factors (Guzzetti et al., 1999), such as rainfall and NDVI.We used rainfall data gathered during typhoon Mindulle, from 28 June to 3 July 2004.NDVIs were calculated from a SPOT image taken just before typhoon Mindulle.The landslides used for validation were those in the event-based inventory that had been derived from SPOT images taken before and after the Mindulle event.Factor weights from the model of the Toraji event in the Kuohsing quadrangle were adopted.The LSI was then calculated for each grid point in both quadrangles.The LSIs were then transferred to indicate probability of failure using the equations given in Fig. 8.The maps of the susceptibility classes for the grid points in both quadrangles are shown in Fig. 12.The prediction rate curves are calculated and shown in Fig. 13.The prediction rate for the subsequent event in the Kuohsing quadrangle was less than the success rate by about 8 and 5 percent (in terms of the AUC) for the hilly and mountainous terrains, respectively.For the Tungshih quadrangle, the prediction rate for the subsequent event in the hilly terrain was slightly lower than the success rate in the Kuohsing quadrangle by about 2 percent; however, in the mountainous terrain, it was higher than the success rate by about 5 percent (Table 8).
In summary, the application of the susceptibility model to a subsequent event, whether in the training site or a neighboring site, was good, as shown by the good prediction rate.The performance in the Tungshih quadrangle was especially good, with a prediction rate better than the success rate.

Discussion
The present study provides a statistical approach for the interpretation of the storm event-induced landslide distribution and for the mapping of regional landslide susceptibility.There are some interesting aspects which need to be discussed.

Methodology and landslide triggering factors
The present EB-LSA differs from traditional statistical LSAs in two ways.First, instead of a multi-temporal landslide inventory, an event-based landslide inventory is used, and second, the triggering factor is emphasized.It is a very important to use both an event-based landslide inventory and a triggering factor.Without the event-based landslide inventory, the triggering factor is not significant; without the triggering factor, the single-period landslide inventory provides no insight into temporal changes in the landslide distribution.
Of the two key issues in EB-LSA, the selection of a triggering factor is more important.To select the best triggering factor, we compared the effectiveness of the total rainfall and the maximum rainfall intensity of a storm event to discriminate between the landslide group and the non-landslide group.Based on these results, the maximum rainfall inten-sity was selected as a discriminator.Hill slopes in the study region are for the most part mantled by permeable colluvium soils.Shallow landslides on such slopes are often correlated to short-period rainfall intensity (Aleotti and Chowdury, 1999).On the other hand, if a slope is mantled by residual soils or clayey soils, it may be the total rainfall that controls the soil water level, making this a better discriminator.However, as mentioned in Sect.5, the percentage of weighting for rainfall intensity factor is not as high as expected, when compared to the percentage of weighting for Arias intensity factor in the earthquake-induced landslide study (Lee et al., 2008): 0.419 and 0.242 for hilly terrain and for mountainous terrain, respectively.
The reason for the relatively low weights for rainfall intensity in the present case may be due to the inclusion of many reactivated landslides in the landslide inventory used in training the susceptibility model.Reactivated landslides tend to occur under smaller rainfall than new ones.This reduces the effectiveness of the rainfall intensity factor and lowers the weighting percentage.Recent studies by our team in the catchments area of the Shihmen Reservoir shows a high percentage of weighting for rainfall intensity factor, apparently due to the fact that the landslides used to build the model are mostly new (i.e., few reactivated landslides).If this is true, it is recommended that reactivated landslides not be considered in the establishment of a landslide susceptibility model in future studies.
When selecting a triggering factor, the dependency between factors must be examined.Generally, inter-factor dependencies are common, as can be seen from the correlation coefficients in Table 9.The correlation coefficients between the factors used in this study ranged from -0.119 to 0.542.The triggering factor -maximum rainfall intensity -also had some dependency on other factors, however the correlation coefficients were relatively small (0.024∼0.392) so it could be used as an independent factor.
Is the maximum hourly rainfall the right explanatory variable?Would the maximum 3-h, 6-h, or 12-h rainfall do?This question needs to be tested with real data in further study.If one has investigated the land cover, soil depth and permeability, this problem may be solvable using hillslope hydrology theory (Freeze, 1978;Dietrich et al., 1995;Borga et al., 2002;Wilkinson et al., 2002;Malet et al., 2005)  Note: Name of each factor is same as that in Table 7.
hydraulic conductivity of colluvium soils in Taiwan is about 10 −3 -10 −4 m/s, or 10 1 -10 2 m/day (first author personal experience from many engineering projects).It may take about 1 to several days for infiltration from the upper slope to be transmitted to the lower slope to reach a steady-state hydrologic condition.However, colluvium soils may be loosened by subsurface storms (Whipkey and Kirkby, 1978) so the flow may be much faster.Field measurements of subsurface flow or test data are absolutely needed for verification.Prior to this verification, it may be more practical to include both the rainfall intensity and duration, such as done by Chang et al. (2007), or the maximum 1-hour and maximum rolling 24-h rainfall, in a statistical model.
To use more accurate rainfall distribution in the modeling is also influential.Although the density of rain gauge stations is high, interpolation of rainfall data is still necessary.In this study we used the ordinary Kriging method to interpolate rainfall data for each grid point.In future studies, multivariate geostatistical methods and the incorporation of more auxiliary variables, such as surface elevation, slope gradient, slope aspect, and radar data, into the analysis may be considered to improve the quality of the interpolation and minimize the estimation error minimized (Goovaerts, 2000;Diodato and Ceccarelli, 2005;Chang, 2007;Haberlandt, 2007).
Proper data selection is important to the building of a susceptibility model.Dai and Lee (2003) used two event-based landslide inventories (for two different years) to build a susceptibility model.Only one specific event was used in training the susceptibility model in the present study.If a specific storm event gives a wide range of rainfall intensity in the study area, then the model would be good for prediction.Otherwise, it is better to select several storm events in the study area to include a wide range of rainfall intensities.To use more event data and to plan a good data selection scheme may be considered in a further study.

Use of the susceptibility model for mapping
A landslide susceptibility map can be used for regional planning, site selection, and policy making for hazard mitigation.It should faithfully reflect the relative hill-slope stability and danger zones.The direct product of a landslide susceptibility map from the present EB-LSA allows for interpretation of the landslide distribution after the specific event, but this may be event-dependent.For example, if there is high rainfall intensity present at some portion of the map, then that portion would show high susceptibility or high landslide probability.A landslide susceptibility map of this kind does not properly represent the relative slope stability under general rainfall conditions in the region.
To prepare better landslide susceptibility maps using the EB-LSA model we should replace the event rainfall data by the rainfall data for a particular return period, such as 10, 20, 50 or 100 years so as to produce a more uniform temporal probability map of the region, such as that done by Dai and Lee (2003).The temporal probability of a storm event may be obtained by frequency analysis of the rainfall (Chow et al., 1988).By combining the spatial probability and the temporal probability, we could establish a probabilistic landslide hazard model.Lee et al. (2005) proposed some examples of landslide spatial probability for certain return-periods of maximum rainfall intensity.

Use of the susceptibility model for prediction
Similar to the deterministic models, an EB-LSA storm model is capable of landslide prediction, provided that the rainfall distribution is known.In EB-LSA, an LSI can be calculated from a pre-trained landslide susceptibility model and the value of the triggering factor of a scenario event.The spatial probability of a landslide at any given point can then be derived from a diagram of the probability of failure vs. the LSI (Fig. 8).
When using the EB-LSA model for prediction, the value of the triggering factor is better within the range that the model was trained.This is the case here.The range of maximum rainfall intensity used in training the model is 43-131 mm/hr, the range in the Tungshih quadrangle is 53-102 mm/hr; well within this range.During the typhoon Mindulle event, the range in the Kuohsing quadrangle is 67-138 mm/hr, the range in Tungshih quadrangle is 36-125 mm/hr; either within this range or not exceeding it by much.This is why the prediction rates in the validations are good.A range of predicted rainfall intensity similar to that of the trained model is preferable.Extrapolation or over range prediction is not so controlled.
When used for prediction, the inherent past condition (like land-use and soil moisture) must be preserved so that the causative factors used in establishing the model are representative of the background conditions for prediction.Of the causative factors, NDVI is most sensitive, so it should be derived from satellite images made close to the time of prediction during a season which has similar weather conditions as when the model was trained.
In Taiwan, hill slopes are commonly covered by loose and permeable colluvium.Antecedent rainfall and land-use do not affect the soil moisture very much.Land-use does affect the amount of interception of rainfall, but the influence is small during a heavy rainstorm.Factors derived from the DEM may also remain unchanged, for shallow landslides do not change the elevation too much (within the accuracy of a DEM).The causative factors should preserve past conditions with the NDVI being renewed, for the triggering factor -maximum rainfall intensity to be effective.Therefore, the susceptibility model may be used for prediction.

Perspectives
Up to now, the proposed EB-LSA has been capable of finding the shallow landslide probability for certain return-periods of storm rainfall, but we still regard this as a susceptibility analysis, because we cannot predict the size or run out distance of a potential landslide (Burton and Bathurst, 1998;Guzzetti et al., 2005;Claessens et al., 2007b).This means that another important step is needed to upgrade the EB-LSA for landslide hazard analysis (LHA).
The grid-cell units used in the present study are judged to be valid, and their use is efficient and correct.However, further study should take the interaction of the nearest neighbors into consideration, to better predict the size of the landslide.Slope units (Carrara et al., 1991;Guzzetti et al., 1999;Xie et al., 2004) are one frequently used susceptibility mapping unit.Parameters derived from a slope unit may reflect the overall characteristics of the slope where a landslide occurs.In the present study, we have used a total-slope-height factor to reflect the longitudinal dimension of a slope unit.If a slope unit is used, then the transverse dimensions of a slope would also be known.A method for slope-unit based grid-cell analysis for EB-LSA is under active study by our group.In the second-phase of the CGS landslide project, slope-unit based landslide danger zone mapping will be employed.
The frequency of landslide occurrence is often dominated by the triggering agent (Aleotti and Chowdury, 1999).Therefore, EB-LSA models need to be closely linked to hydrological frequency analysis for storm-induced landslide prediction.Because the density of rain gauge stations is seldom ideal, interpolation of rainfall data is always necessary.Improvement by incorporating radar data (Chiang and Chang, 2008) and/or using multivariate geostatistical methods (Chang, 2007) is needed.
The present discussion has been limited to shallow landslides.Other types of landslides, such as rock-falls, deepseated slides, and debris flows, require different methods (mainly the incorporating of different mapping units and different causative factors).The movement of some types of landslide can range from very slow to very rapid, meaning that the danger level will be very different.This feature should be carefully considered in any future LHA.

Conclusions and recommendations
A method of event-based landslide susceptibility analysis (EB-LSA) was applied to data from a storm event-induced shallow landslide study in the Kuohsing quadrangle.Results show that our EB-LSA for storm event-induced landslides is effective as confirmed by careful validation on a neighboring region and on a subsequent event.The present methodology and working procedure are feasible for both storm-induced and earthquake-induced landslide studies.
EB-LSA uses an event-based landslide inventory derived from a pre-event and a post-event remotely sensed image, a set of environmental variables and a triggering factor to train the susceptibility model.The combination of an event-based landslide inventory and a triggering factor in the model allows for effective interpretation of the event-induced landslide distribution.Landslide susceptibility mapping of this kind is event-dependent; landslide probability at a given point is known only when the triggering factor is given.For actual mapping, it is recommended that the landslide susceptibility of a region be represented by a susceptibility map using certain return-period rainfall values.
The maximum rainfall intensity (maximum hourly rainfall) is found to be an effective factor in the interpretation of the event-induced landslide distribution in the present study.However, if hill slopes are mantled by clayey soils, then the story could be different; the rainfall duration or the total rainfall may control the result.An optimum selection of a rainfall factor or other factors needs to be determined for a given region.
The present approach is feasible for regional susceptibility mapping and is capable of prediction for a scenario event without hydrological data, strength data, failure depth, or a www.nat-hazards-earth-syst-sci.net/8/941/2008/Nat.Hazards Earth Syst.Sci., 8, 941-960, 2008 long-period landslide inventory.It may be further developed so that landslide size and run-out distance could be included. Fig.1

Fig. 1 .
Fig. 1.Location and geology of the study sites: (a) geologic map of Tungshih quadrangle; (b) geologic map of Kuohsing quadrangle; (c) geologic provinces of Taiwan and index map; (d) legend for geologic maps.

Fig. 2 .
Fig. 2. Spatial distribution of landslides in the study region.(a) Typhoon Toraji, (b) typhoon Mindulle.Upper quadrangle is the Tungshih site, lower quadrangle is the Kuohsing site; black polygons indicate landslides, white polygons indicate no information is available from the SPOT images.

Fig. 3 .
Fig. 3. Schematic map showing the definition of factors.(A) elevation of crest, (B) horizontal distance to drainage, (C) height relative to riverbed, (D) elevation of toe, (E) total slope height, (F) height relative to crest, (G) height relative to toe, (H) horizontal distance to crest, (I) horizontal distance to toe, (J) horizontal distance between crest and toe, (K) slope length.

Fig. 9 Fig. 9 .
Fig. 9Fig.9. Landslide susceptibility map of the Kuohsing quadrangle developed using susceptibility model trained with Toraji inventory at the Kuohsing quadrangle.Red indicates high susceptibility, yellow moderately high, green moderate, cyan low, and gray stable.Landslides triggered by the typhoon Toraji are shown.

Fig. 11 .
Fig. 11.Landslide susceptibility map of the Tungshih quadrangle developed using susceptibility model trained with Toraji inventory in the Kuohsing quadrangle.Landslides triggered by the typhoon Toraji are shown.

Table 1 .
Rainfall characteristics and losses during the typhoon Toraji and the typhoon Mindulle.
* 1: At the Nantienchi gauge station in Kaohsiung County, southern Taiwan.* 2: At the Shihnan gauge station in Kaohsiung County, southern Taiwan.* 3: Name in Parentheses is name of gauge station in the quadrangle.

Table 2 .
SPOT images used in the Kuohsing quadrangle.

Table 2
SPOT images used in the Kuohsing quadrangle

Table 3 .
SPOT images used in the Tungshih quadrangle.

Table 2
SPOT images used in the Kuohsing quadrangle

Table 4 .
Landslide numbers and areas in the Kuohsing and Tungshih quadrangles.

Table 5 .
Descriptive statistics of causative factors other than lithology.

Table 6 .
Descriptive statistics of lithology factor.
The formation names in 8-10 are local names for the Shuilikeng Formation as shown in Fig.1.

Table 7 .
Results of coefficients of the discriminant function for each terrain in the Kuohsing quadrangle.TerrainLitho Slope Slope Asp Topo Rou Slope Rou Total Curv Total Slope NDVI Rain Note: Litho: Lithology.Slope: Slope gradient.Slope Asp: Slope aspect.Topo Rou: Topographic roughness.Slope Rou: Slope roughness.Total Curv: Total curvature.Total Slope: Total slope height.NDVI: Normalized differential vegetation index.Rain: Maximum rainfall intensity.

Table 8 .
AUC of success and prediction rates for the Kuohsing and Tungshih quadrangles.

Table 9 .
Correlation coefficients between susceptibility factors.