Estimating the topographic predictability of debris flows

Geomorphology 207 (2014) 114–125

Contents lists available at ScienceDirect

Geomorphology journal homepage: www.elsevier.com/locate/geomorph

Estimating the topographic predictability of debris flows Nele Kristin Meyer a,b,⁎, Wolfgang Schwanghart c, Oliver Korup c, Bård Romstad d, Bernd Etzelmüller b a

International Centre for Geohazards (ICG) c/o NGI, Sognsveien 70, 0855 Oslo, Norway University of Oslo, Department for Geosciences, Sem Sælands vei 1, 0371 Oslo, Norway c University of Potsdam, Institute of Earth and Environmental Science, Karl-Liebknecht-Str. 24-25, 14476 Potsdam-Golm, Germany d Center for International Climate and Environmental Research (CICERO), Gaustadalléen 21, 0349 Oslo, Norway b

a r t i c l e

i n f o

Article history: Received 1 May 2013 Received in revised form 28 October 2013 Accepted 29 October 2013 Available online 6 November 2013 Keywords: Weights-of-Evidence Debris flows Susceptibility Slope–area plot Process domains Norway

a b s t r a c t The Norwegian traffic network is impacted by about 2000 landslides, avalanches, and debris flows each year that incur high economic losses. Despite the urgent need to mitigate future losses, efforts to locate potential debris flow source areas have been rare at the regional scale. We tackle this research gap by exploring a minimal set of possible topographic predictors of debris flow initiation that we input to a Weights-of-Evidence (WofE) model for mapping the regional susceptibility to debris flows in western Norway. We use an inventory of 429 debris flows that were recorded between 1979 and 2008, and use the terrain variables of slope, total curvature, and contributing area (flow accumulation) to compute the posterior probabilities of local debris flow occurrence. The novelty of our approach is that we quantify the uncertainties in the WofE approach arising from different predictor classification schemes and data input, while estimating model accuracy and predictive performance from independent test data. Our results show that a percentile-based classification scheme excels over a manual classification of the predictor variables because differing abundances in manually defined bins reduce the reliability of the conditional independence tests, a key, and often neglected, prerequisite for the WofE method. The conditional dependence between total curvature and flow accumulation precludes their joint use in the model. Slope gradient has the highest true positive rate (88%), although the fraction of area classified as susceptible is very large (37%). The predictive performance, i.e. the reduction of false positives, is improved when combined with either total curvature or flow accumulation. Bootstrapping shows that the combination of slope and flow accumulation provides more reliable predictions than the combination of slope and total curvature, and helps refining the use of slope–area plots for identifying morphometric fingerprints of debris flow source areas, an approach used outside the field of landslide susceptibility assessments. © 2013 Elsevier B.V. All rights reserved.

1. Introduction Hillslope mass wasting processes account for up to 14% of the loss of lives due to natural hazards worldwide (Aleotti and Chowdhury, 1999). National costs related to landslides in the U.S., Japan, and the Alpine countries vary between US$ 1 and 5 billion annually (Aleotti and Chowdhury, 1999; Sassa and Canuti, 2008), and yet much of this damage may remain underestimated (Petley, 2012). A large fraction of the world's terrestrial landslide occurrences are concentrated in tectonically active mountain belts, fault zones, and volcanic island arcs (Korup, 2012). Yet mass wasting in mountainous terrain along passive continental margins may also cause substantial problems. The Norwegian traffic network, for example, is impacted by thousands of landslides, avalanches, and debris flows each year that incur annual economic losses of the order of 12 million € (Bråthen et al., 2008; Bjordal and Helle, 2011). Rapid debris flows, i.e. mixtures of unconsolidated sediment and water may travel at speeds N 10 m s−1 (Hungr et al., 2008), are ⁎ Corresponding author at: International Centre for Geohazards (ICG) c/o NGI, Sognsveien 70, 0855 Oslo, Norway. Tel.: +47 45185079. E-mail address: [email protected] (N.K. Meyer). 0169-555X/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.geomorph.2013.10.030

the most devastating among these mass wasting processes, and cause damage to transport infrastructure of 0.8 million € per year on average (Romstad, 2013). Yet routing traffic arteries through the deeply dissected fjords of western Norway leaves little options such that knowledge about potential debris flow source areas and runout paths is highly desirable. A large toolbox of methods for identifying such areas susceptible to landslides and debris flows has become available (Guzzetti et al., 2006). Susceptibility maps show the spatial propensity or proneness to landslides without quantifying any probability of occurrence, as opposed to hazard maps (Van Westen et al., 2003; Fell et al., 2008). Direct susceptibility mapping of debris flow source areas usually relies on detailed field mapping and is thus restricted to small study areas. For regional-scale studies, indirect and less field-work intensive susceptibility mapping methods use statistical modelling (Van Westen et al., 2003; Guzzetti et al., 2006; Thiery et al., 2007). This approach assumes that the spatial distribution of slope instability is not stochastic but determined by environmental conditions that also control future landslide occurrences (Fabbri et al., 2003; Fabbri and Chung, 2008). Indirect susceptibility mapping uses a set of environmental indicators (= predictor variables) and the spatial distribution of past landslide

N.K. Meyer et al. / Geomorphology 207 (2014) 114–125

events (= response variables) to identify and characterize potentially prone locations via a degree of susceptibility for a given unit area (Guzzetti et al., 2006; Blahut et al., 2010). A special type of bivariate models concerns the concept of process domains delineated in plots of local slope gradient and contributing catchment area, also known as slope–area plots. Such process domains are largely empirical in nature and entail hillslopes, unchannelled valleys, debris flow dominated (or hillslope-controlled) channels, bedrock and alluvial (self-formed) channels (Montgomery and Foufoula-Georgiou, 1993; Montgomery, 1999, 2001; Stock and Dietrich, 2003; Brardinoni et al., 2012; Williams, 2012). The transition from hillslopes to unchannelled or debris flow dominated channels is often characterized by a kink in slope–area plots: Slope increases with drainage area on hillslopes, though decreases at the point of channel or debris flow initiation (Montgomery and Foufoula-Georgiou, 1993; Brardinoni and Hassan, 2006). Debris flow initiation requires steep slopes, accumulation of unconsolidated sediment, and sufficient porewater pressures (Innes, 1983), and these factors can be approximated by the terrain metrics of slope, total curvature, and flow accumulation. Flow accumulation quantifies the upstream contributing area that serves as proxy of accumulated water runoff in humid areas (Blahut et al., 2010; Kappes et al., 2011). Total (or mean) curvature is the average curvature of any two orthogonal normal sections of the terrain surface, describes the overall curvature of the surface, regardless of slope direction (Wilson and Gallant, 2000), and has been used as a proxy for local soil depth or sediment thickness (Montgomery, 1999; Minasny and McBratney, 2001). Weights-of-Evidence (WofE) is a data-driven statistical bivariate method using a log-linear form of Bayes' theorem to determine the weight (or importance) of evidence (Bonham-Carter, 1994). In several recent studies (Pradhan et al., 2010; Oh and Lee, 2011; Prasannakumar and Vijith, 2012; Quinn et al., 2012; Lee, in press; Regmi et al., in press) WofE is used for regional-scale landslide susceptibility assessments for areas 101 to 102 km2 in size, based on landslide inventories with up to several hundreds of entries. The choice of predictor variables is usually based on prior (or a priori) knowledge on the causes of landslides, though often restricted by data availability. The number of predictor variables included in published WofE models differs between studies and ranges from seven to 20 variables (Mathew et al., 2007; Regmi et al., 2010; Xu et al., 2012). Using multiple predictor variables, however, may violate the assumption of conditional independence between these predictors, which is required for an unbiased susceptibility estimate (Thiery et al., 2007). The derivation of unique condition units (UCU) has often allayed this problem (Piacentini et al., 2012). However, such UCUs do not allow quantifying the relative importance of single factors, which is a distinctive advantage of the WofE model (Süzen and Doyuran, 2004). While the majority of the many hundreds of published landslide susceptibility studies achieve a prediction accuracy of N 80% (e.g. Dahal et al., 2007; Regmi et al., 2010; Sterlacchini et al., 2011; Pourghasemi et al., 2013; Prasannakumar and Vijith, 2012), this seeming success is unrelated to the type of model or the number of predictor variables (Korup and Stolle, in press). More complex models with more predictor variables often fail to significantly improve the predictive accuracy obtained by simpler models with a limited number of predictor variables. Sterlacchini et al. (2011) test different WofE models for the same study area and obtain similar predictive performance for all of them despite significant spatial disagreement between the different resulting susceptibility maps. Therefore, susceptibility maps should always be probed for uncertainties related to each susceptibility class. Acknowledging these constraints, our strategy of predicting topographic susceptibility to debris flow initiation in western Norway follows a minimalist approach with regard to the number of predictors, which should allow maximum transparency when interpreting our predictions. Numerous susceptibility studies show that slope, total curvature, and flow accumulation are among the most important predictor

115

variables of mass wasting (Lee and Choi, 2004; Dahal et al., 2007; Regmi et al., 2010; Fischer et al., 2012), and we further test this proposal here. We acknowledge that, among others, lithology, precipitation, vegetation cover or human activity, may also influence debris flow initiation (Van Westen et al., 2003). Yet these remain difficult to quantify at high spatial resolution (Fischer et al., 2012). Moreover, Fabbri et al. (2003) and Guinau et al. (2007) found that landslide susceptibility maps based on topographic metrics outperformed maps based on other thematic layers or a combination of both. To quantify the uncertainties arising from previous knowledge about debris flows in our area, we aim at deriving from training data a model with a minimum number of predictor variables. These should be both readily available and sufficiently representative of the physical processes and environmental controls of debris flow initiation. We then quantify model uncertainties arising from different classification schemes and random subsets of observed debris flows, before computing the errors of our prediction on separate test data. 2. Study area and data Many parts of Norway are exposed to rapid mass movements (Jaedicke et al., 2009), which in parts have been documented in a national database (Skrednett, 2013; Fig. 1). Within this database 710 entries are debris flows reported by the national road and railway authorities, and consulting agencies from 1900 to 2008. Each entry features the timing and approximate coordinates of where debris flows had deposited or blocked a traffic line. We focus on a study area in western Norway, extending from Ålesund in Møre og Romsdal to Notodden in Telemark (close-up in Fig. 1). The Precambrian basement dominates the largest part of the study area. In the southern and central parts, nappes from the Caledonian orogenesis form the bedrock. The lithology is dominated by granites and gneisses (Bøe et al., 2010; Fischer et al., 2012). Topographically the area is characterized by an extensive plateau at elevations between 1200 and 1500 m, which is deeply incised by fjords and U-shaped valleys. The treeline ranges between 600 and 1050 m with a strong west–east gradient, rising with increasing distance from the coast (Rössler et al., 2008). Within our study area 429 debris flows are recorded along the transport routes that closely follow valley bottoms or lower slopes of the steep fjord sidewalls (Fig. 1). Since the source areas of these events are unknown, we checked each reported coordinate, and relocated each event to its point of initiation identified on aerial photographs using 3D visualization (Norkart Virtual Globe, http://www.virtual-globe.info). Initiation points were set by tracking debris flow paths upstream to the deposition site where channels or bare hillslopes were clearly discernible. This manual mapping procedure enabled us to validate each debris flow event, and classify the events into open-slope (38%) and channelized (46%) debris flows. For 16% of the events it was difficult to identify a distinct flow path from the optical imagery since significant land cover changes had occurred after the debris flow occurrence and before the acquisition of the aerial photograph. We selected 263 debris flows (~60%) documented prior to 2005 as training data, and the remaining 166 debris flows (~40%) as testing data. The higher data density of the inventory in recent years results from a reporting bias rather than increasing debris flow activity (Jaedicke et al., 2009). We used a digital elevation model (DEM) with 10 m grid resolution for extracting three topographic candidate predictors of susceptibility to debris flow initiation. The DEM is derived from 5 m contour lines (FKB-H5) in settled areas, and 20 m contour maps (N50) in remote areas. In settled areas the vertical error is 2 m while in remote areas it is 6 m (Statens Kartverk, 2011). We resampled the DEM to 20 m resolution to reduce a possible resolution bias in remote areas (Fischer et al., 2012), and clipped the DEM to account for documentation gaps of debris flows in remote areas. Lastly, we extracted only those first-order catchment areas draining towards roads or railways, and excluded fjords and lakes (Fig. 1).

116

N.K. Meyer et al. / Geomorphology 207 (2014) 114–125

Fig. 1. Distribution of documented debris flows in Norway from 1979 to 2008 (Skrednett, 2013), and close-up of study area along transport route corridors. Point and star symbols show the spatial distribution of training and test data, respectively.

N.K. Meyer et al. / Geomorphology 207 (2014) 114–125

3. Methods 3.1. Weights-of-Evidence We apply WofE to debris flows and calculate the probability of debris flow initiation per unit area by updating its prior probability P(D) with additional evidence from one or more predictor variables, where D represents the presence or absence of a debris flow. We compute the prior probability from the spatial density of reported debris flows: NfDg P ðDÞ ¼ N fT g

ð1Þ

where N{D} is the number of raster cells containing a debris flows, where each event is considered to be limited to one cell, and N{T} is the total number of cells in the study area. By integrating additional evidence derived from a classified predictor variable we compute an updated conditional probability, or posterior probability. To this end, we combine the prior probability with positive and negative weights calculated from the conditional probabilities given the presence or absence of a binary predictor variable, B (e.g. a certain slope gradient interval), and D. To obtain the posterior probability, we first transform the prior probability to logit (D) using the odds formulation O(D). The logit is the natural logarithm of the odds, which is in turn the ratio of debris flow occurrence probability to the probability that it will not occur. OðDÞ ¼

P ðDÞ 1−P ðDÞ

ð2Þ

Logit ðDÞ ¼ loge ðOðDÞÞ

ð3Þ

The WofE method requires discrete predictor variables. The simplest case is a binary pattern representing the presence or absence of an environmental factor, B, say a topographic metric such as slope gradient, total curvature, or contributing catchment area. For the method to work, the value range of B needs to be classified. For each class of B we compute the conditional probability of being present or absent given that debris flow D occurred (Table 1). From these conditional probabilities we obtain positive and negative weights indicating the likelihood of each class B to be associated with debris flows. Positive weights W+ are calculated as P ðBjDÞ W ¼ log
 P BjD

!

þ

ð4Þ

where W+ N 0, debris flow occurrences are positively correlated with the predictor variable, and negatively correlated otherwise. No correlation exists if W+ = 0; in this case the posterior probability equals the prior probability, and no further predictive information is gained from the predictor. We set W+ = 0 if there is no debris flow data available within B. Negative weights W− are calculated as W

−


! P BjD  : ¼ log
P BjD

ð5Þ

Table 1 Conditional probabilities of a classified predictor variable B being present or absent for a documented debris flow occurrence D. Topographic predictor variable

Debris flow

Present

Absent

Present

g P ðBjDÞ ¼ NNfB∩D f Dg


 NfB∩Dg P BjD ¼ NfDg

Absent


 NfB∩Dg P BjD ¼ NfDg


 NfB∩Dg P BjD ¼ NfDg

117

For W− N 0 debris flow occurrences are negatively correlated with predictor variable B, and positively correlated otherwise. No correlation exists if W− = 0. The difference between W+ and W− is called the contrast C, and indicates the spatial association of the weight with each class of a predictor variable. þ

C ¼ W −W

−

ð6Þ

Since we estimate each weight assigned to a class from a finite set of observations, the value of each weight is subject to uncertainty. The reliability of the weight as a predictor can be calculated as standard deviation S(C) of the contrast. SðC Þ ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi s2 ðW þ Þ þ s2 ðW − Þ

ð7Þ

where the variances of W+ and W− are calculated as follows   þ ¼ W

1 1 þ
 P ðDjBÞ P DjB

ð8Þ

1 1 2 − s ðW Þ ¼
þ
: P DjB P DjB

ð9Þ

2

s

One of the advantages of the formulation of weights over calculating the posterior probability directly is the possibility to quantify and visualize the uncertainty associated with the final susceptibility map. The use of logits is also computationally more advantageous, since very low probabilities may introduce round-off errors in floating point arithmetics. To obtain the posterior logit the local weight is added to the prior logit given the presence or absence of B. The posterior probabilities can be derived by retransformation Logit ðDjBÞ ¼ Logit ðDÞ þ W

þ=−

ð10Þ

OðDjBÞ ¼ expðLogit ðDjBÞÞ P ðDjBÞ ¼

ð11Þ

OðDjBÞ : 1 þ OðDjBÞ

ð12Þ

The logit formulation also allows combining the conditional probabilities of several predictor variables by adding the weights of each predictor variable. Hence, the logit formulation for n variables can be extended to Logit ðDjBÞ ¼ Logit ðDÞ þ

Xn i¼1

þ=−

Wi

:

ð13Þ

Since each binary pattern of a predictor variable has a unique spatial distribution (if conditionally independent of all other predictors), the spatial variation in posterior probabilities has a maximum of 2n different values. Eq. (9) can be extended to use WofE with one to several multiclass predictor variables by encoding each multiclass variable with k classes into k independent binary variables. The posterior logits are calculated by summing the weights of the binary patterns of all multiclass variables:   Xn Xk j ð j¼1 to kÞ þ=− Logit DjBi ði¼1 to nÞ ¼ Logit ðDÞ þ W ij i¼1 j¼1

ð14Þ

where Bji is the jth class of the ith predictor variable. The posterior logit is transformed back to the posterior probability using Eqs. (11) and (12).

118

N.K. Meyer et al. / Geomorphology 207 (2014) 114–125

conditionally dependent, their joint use has to be rejected. Two classes are conditionally independent if:

3.2. Predictor classification schemes, and bootstrapping Since WofE requires categorical data, continuous variables must be classified. This classification is critical since subjective choice of the number of classes entails a trade-off between retaining high information content and obtaining statistically significant results. While binary classifications lead to high information loss, too many class breaks tend to create sparsely populated or even empty classes, for which it is difficult to detect conditional dependencies (Bonham-Carter, 1994). Over-fitting can be avoided by multi-class predictor variables with an adequate number of classes instead of binary maps, although the computational effort is larger (Mathew et al., 2007). We tested the sensitivity of WofE to univariate classification of the predictor variables with two different classification schemes. The first scheme uses arbitrary class breaks guided by visual interpretation of the estimated probability density functions of the predictor variables, whereas the second scheme uses deciles as class breaks. We further tested the robustness of the decile classification by repeatedly using randomly selected subsets of the input data. Robust weights are calculated with bootstrapping by re-evaluating the model 1000 times, and using a random selection of 2/3 of the training data set for each run.

3.3. Test of conditional independence WofE is based on the assumption that conditional independence exists between the different predictor variables (Bonham-Carter, 1994; Regmi et al., 2010). Conditional independence is required because the final weights are calculated as the sum of the single weights obtained for each predictor variable. Weights attributed to conditionally dependent factor maps would overestimate the posterior probability where both predictor variables are present. Agterberg et al. (1990) and Bonham-Carter (1994) pointed out cases where too many predictor variables compromised the assumption of conditional independence at the expense of minute changes in model prediction. All multiclass maps have to be tested pair-wise and the contingency table must include a χ2-value for all possible combinations of the classes of two predictor variables (Dunlop, 2010). If two predictors are

NfB1 ∩Dg  N fB2 ∩Dg : NfDg

ð15Þ

The left-hand side of this equation is the observed number of debris flows in the overlap region where both patterns B1 (e.g. a given slope class) and B2 (e.g. a given curvature class) are present; and the righthand side is the calculated expected number of deposits in this overlap zone. The individual χ2-value using Yates correction for small frequencies (Bonham-Carter, 1994) with differences N0.5 is 2

x ¼

X4 ðjObserved −Expected j−0:5Þ2 i i : i¼1 Expectedi

ð16Þ

The null hypothesis states that the weights in B1 and B2 are independent from each other. If the calculated χ2-value exceeds the tabulated value, then the hypothesis has to be rejected, and both classes have to be regarded as conditionally dependent (Bonham-Carter, 1994). To reduce the computational effort we adopted the method of Dunlop (2010) calculating the sum of χ2-values for all class combinations of two predictor maps, and compared it to the tabulated value with the appropriate degrees of freedom. 3.4. Model evaluation The success rate relates the model to the debris flow data that was used to calibrate the model, while the prediction rate assesses the relation between model and validation or test data set, thus indicating the predictive performance of a model. We use receiver operating characteristic (ROC) curves to analyse the relation between correctly classified debris flow source cells (true positives or sensitivity) and falsely classified cells (false positives or 1 − specificity). The area under curve (AUC) is a measure of the likelihood that the obtained model classifies a phenomenon correctly. An AUC value of 0.5 (diagonal straight line) implies that the predictive performance is equal to random sampling, and an AUC value of 1 indicates a perfect prediction (Fawcett, 2006).

C

B

Weights

Weights

A

NfB1 ∩B2 ∩Dg ¼

10

20

30

40

50

Slope (°)

60

70

80

-10 -5 -1 -0.5 -0.1 0 0.1 0.5 1

5 10

1

10 5

Curvature (m-1)

100 50

1000 10000 100000 1000000 500 5000 50000 500000

Flow accumulation (no. of cells)

Fig. 2. Weights for the arbitrarily set classes of the predictor variables (A) slope, (B) total curvature, and (C) flow accumulation, for the training data set. Weights N0 and b0 indicate positive and negative correlations with debris flow initiation points, respectively. Error bars are standard deviations of the contrast (Eq. (7)) that serves as a measure of uncertainty related to the calculated weights.

A

B

C

Curvature (m-1)

Slope (°)

119

% of debris flow data population

% of debris flow data population

N.K. Meyer et al. / Geomorphology 207 (2014) 114–125

Flow accumulation (no. of cells)

Fig. 3. Calculated class breaks based deciles of debris flow data for (A) slope, (B) total curvature, and (C) flow accumulation, based on bootstrap sampling using 2/3 of the training data set. Boxes contain 50% of the data showing the median as vertical bar. Whiskers show the lower and upper quantiles. Outliers are indicated as points (see Fig. 4).

4. Results The prior probability of debris flow initiation anywhere in the study area is 2.8 × 10−6. In the following we refer to cells as susceptible if their posterior probability exceeds their prior probability, and are thus characterized by positive weights. Results on the manual classification are based on nine equal-size classes of 10° width each for slope angle (Fig. 2A), and 12 classes of variable width for total curvature. We chose bin widths of N0.1 m−1 for planar surface segments, and larger ones (b5 m−1) for highly convex and concave slopes (Fig. 2B). We classified the log-transformed flow accumulation into 14 equal-size classes (Fig. 2C). Weights are positive for slope angles between 30° and 60°, convex total curvatures between −5 and −0.05 m−1, and flow accumulations of 50 to 5000 cells (equal to a contributing area of 0.2 to 2 km2, respectively; Fig. 2). The standard deviation of the contrast is not equally distributed among the different classes. Classes containing no debris flows have zero contrast. The largest errors are tied to classes covering only small portions of the study area, and containing few debris flows. In contrast, the data-driven classification uses equal portions of the documented debris flows distributed over the variable space

B

C

Weights

Weights

A

(Fig. 3), where we defined class breaks from the bootstrapped median of each decile. We calculated an individual weight for each decile class, and find that slopes from 22° to 61°, convex total curvatures from −0.73 to −0.03 m−1, and flow accumulations from 44 to 5065 cells (i.e. a contributing area of ~0.02 to ~2 km2) have positive weights. All three predictors were conditionally independent with an error probability b0.05 for the manual classification scheme. However, pairwise tests reveal that some of the class combinations feature only few debris flows, resulting in a weak conditional test of independence (Agterberg and Cheng, 2002). A more reliable test of the class combinations is given by the decile classification as each class features nearly the same debris flow density and uncertainties related to the weights are approximately similar (Fig. 4). Here, the test of conditional independence shows that flow accumulation and total curvature are conditionally dependent (Table 2). Analysis of the pair-wise tests reveals that six class combinations exceed the tabulated value at 0.95 confidence level (Fig. 5). Convex surfaces correlate with low flow accumulations while concave surfaces correlate with higher flow accumulations. This relation has been previously attributed to the notion that total curvature is a measure of the local accumulation of the topography while flow

17 22 28 32 34 38 41 45 50 61

-0.73 -0.19 -0.08 -0.03 0 0.21 -0.26 -0.14 -0.05 -0.01 0.04

Slope (°)

Curvature (m-1)

1

2

6

12 18 44 86 253 732 5065

Flow accumulation (no. of cells)

Fig. 4. Weights for the percentile classes of predictors (A) slope, (B) total curvature, and (C) flow accumulation, calculated for the training data set. Weights N0 and b0 indicate positive and negative correlations with debris flow occurrence, respectively. Error bars are standard deviations of contrast (Eq. (7)), and serve as a measure of uncertainty associated with the weights.

120

N.K. Meyer et al. / Geomorphology 207 (2014) 114–125

Table 2 Conditional test of independence for multiclass data (Dunlop, 2010) for decile classes. Sum of calculated chi-square values (bold) for all combinations of the predictor variables compared to the tabulated value at 0.95 confidence level. Degrees of freedom are given in brackets.

Flow accumulation Curvature

Slope

Curvature

90.43 b103.01 (81) 39.38 b103.01 (81)

129.02 N103.01 (81)

accumulation refers to its global accumulation (Schorghofer and Rothman, 2002). Based on these results we refrained from using the manual classification scheme, and excluded total curvature and flow accumulation from the same model. The distributions of weights calculated for slope, total curvature and flow accumulation depict the individual performances of the three predictor variables (Fig. 6A–C). The success rate of the estimated weights for the slope model is 88%, and the predictive rate is 86%. Randomly sampling from the study area ~2 × 106 cells within N 30 m of each other, we find that 37% of the entire area has weights N0. In contrast, the success (predictive) rates of total curvature and flow accumulation are only 61% (66%) and 58% (58%), respectively. The area classified susceptible to debris flows is commensurately lower, i.e. 21% and 24%, respectively. Since we cannot exclude conditional dependence between some curvature and flow accumulation classes, we checked combinations of slope with either total curvature or flow accumulation (Fig. 6D–E). Both combinations have success rates (76% and 80%) and predictive rates (78% and 80%) that are lower than those

obtained by using slope alone, but higher than those solely from total curvature or flow accumulation. Also, the area attributed with positive weights (21% and 24%) is much less than the classification by slope alone. The ROC curves and the corresponding AUC values underline the better overall performance of the WofE model using the combination of two predictor variables only (Fig. 7). The AUC values indicate that the predictor combination of slope and total curvature predicts slightly better than the combination of slope and flow accumulation. To measure the uncertainty related to model performance we bootstrapped training and test data sets (n = 100), and analysed the scatter of ROC curves related to the three variables slope, total curvature, and flow accumulation (Fig. 8). The relative ranking of the variables with regard to their predictive capabilities remains the same, but the AUC values of the total curvature have a higher standard deviation, and thus higher prediction uncertainty than the other two variables. When locating the recorded debris flow data together with a randomly sampled subset of the study area in a slope–area plot (Fig. 9), we find that locations with the highest weights are congruent with the debris flow-dominated process domain. Vice versa, locations with low weights occur within the domain of unchannelled valleys, and the channel network; 44% of the debris flow source locations and 12% of the study area are within the domain of debris flow-dominated channels (Table 3) if adopting the domain mapping proposed by Brardinoni and Hassan (2006). In addition, we also observe a densely populated hillslope domain at slopes N0.5 m−1 (~11.3°), which we interpret as a previously undocumented process domain of debris flow dominated hillslopes. This new domain covers 38% of the debris flow initiation points and 15% of the study area, and agrees especially well with the

860 440 180 120 60 10

20

Flow accumulation (no. of cells)

2530 7320

3.84 (tabulated critical value at 0.95 confidence interval with 1 degree of freedom

-0.73

-0.26

-0.19

-0.14

-0.08

-0.05

Curvature

-0.03

-0.01

0

0.04

0.21

(m-1)

Fig. 5. Results of pair-wise conditional independence test for the predictor variables flow accumulation and total curvature. Class combinations with test scores exceeding the tabulated value at 0.95 confidence level have to be judged as conditional dependent. The grey ellipsis indicates tendency for conditional dependence.

N.K. Meyer et al. / Geomorphology 207 (2014) 114–125

Sucess rate: 88% Predictive rate: 86% Susceptible area: 37%

B Density

Sucess rate: 61% Predictive rate: 66% Susceptible area: 21%

C

Sucess rate: 58% Predictive rate: 57% Susceptible area: 24%

D

Sucess rate: 76% Predictive rate: 78% Susceptible area: 21%

E

Sucess rate: 80% Predictive rate: 80% Susceptible area: 24%

Density

A

121

Weights Random area sample: n = ~2 x 10

6

Training data: n = 263 Test data: n = 166

Weights Fig. 6. Probability density estimates of single weights for (A) slope, (B) total curvature and (C) flow accumulation, and combined weights for (D) slope and total curvature, and (E) slope and flow accumulation for ~2 × 106 randomly sampled locations in the study area (black); the training data set (blue, n = 263); and the test data set (red, n = 166). Negative and positive weights indicate below- and above-average susceptibility levels of debris flow initiation, respectively. Success and prediction rates, and susceptible area refer to a cut-off value at weight 0.

initiation locations that we previously classified as open-slope debris flows (Fig. 9). We defined those areas that were assigned positive weights as susceptible, and thus had a higher posterior probability of debris flow initiation than the prior probability. Yet no distinction was made within this

“susceptible” class, and class boundaries have to be found for visualization. We arbitrarily defined four susceptibility classes “high”, “medium”, “low” and “very low”, containing 25%, 50%, 75% and 100% of the debris flow data, respectively. Posterior probabilities assigned to area cells range over three orders of magnitude from 1.4 × 107 to 2.46 × 104.

bootstrapping (n = 100)

Predictive performance (AUC)

Predictive performance (AUC) Slope (0.80) Curvature (0.76) Flow accumulation (0.70) Slope/ Curvature (0.86) Slope/ Flow accumulation (0.85)

Slope mean: 0.82, sd: 0.016 Curvature mean: 0.75, sd: 0.021 Flow accumulation mean: 0.71, sd: 0.016

1 - Specificity 1 - Specificity Fig. 7. ROC curves for different predictors and predictor combinations. Sensitivity corresponds to the true positive rate (= debris flow affected area correctly classified), while specificity is the true negative rate (= area not affected by debris flows correctly classified). Numbers in brackets refer to calculated areas under the curves (AUC). The closer the curve is to the upper left corner (and the higher AUC), the better the predictive performance.

Fig. 8. ROC curves for assessing the predictive performance of slope, total curvature, and flow accumulation, based on n = 100 bootstrap iterations of randomly subsets of debris flow observations, highlighting the bandwidth of uncertainty tied to the initial data set to which the WofE method is applied. Means and standard deviations (sd) refer to the areas under the curves (AUC values).

122

N.K. Meyer et al. / Geomorphology 207 (2014) 114–125

2.00 0.50 1.00

Table 4 Area, class boundaries and posterior probabilities calculated for the four susceptibility classes. The number in brackets in the field “% of area” shows the area distribution of the generalized map (see Fig. 10).

0 -3

Susceptibility class

0.10

0.20

High Medium Low Very low

hillslopes

alluvial channels

unchanneled valleys

0.02

0.05

Local slope (m/m)

Weights 6

debris flow dominated channels

0.001

0.01

0.1

1

10

100

Contributing area (km²) Fig. 9. Slope–area plot showing study-area sample points (grey) and documented debris flow initiation points (colour ramp indicates modelled weights assigned to the debris flow locations). Superimposed process domains are separated by dashed lines proposed by Brardinoni and Hassan (2006). The steeper portion of the hillslope domain can be regarded as debris flow prone hillslopes (N0.5 m−1). Symbols specify source locations: squares are open-slope, and triangles are channelized debris flows.

However, the boundaries separating the four different classes are within one order of magnitude (Table 4). More than 90% of the area is classified as very low or low susceptible, and b 1% is in the high susceptibility class. The ratio of the percentage of debris flows captured by one class (25%), and the area occupied by the same class indicates the medium and high susceptibility classes as strong predictors. Especially the high susceptibility class seems to give a good indication of potential debris flow sources (Fig. 10A,C). However, the division into susceptibility classes is artificial and lacks physical justification. Thus, we additionally applied a terrain classification according to the process domains outlined above (Fig. 10B,D). Almost all hillslopes are occupied by the two debris flow dominated domains. Debris flow dominated channels are well indicated but are too abundant to be of value for further practical application in subsequent risk analysis.

5. Discussion We used the WofE method for estimating the susceptibility to debris flow initiation in western Norway using an inventory of historic events, and three topographic predictor variables that we classified both manually and by deciles. Our results are in good agreement with threshold values reported by Fischer et al. (2012) for entire Norway. In their study critical slopes for debris flow initiation were found between 25° and 45° and contributing areas exceeding a minimum threshold varying between 0.03 and 0.1 km2. However, our prime interest was not so much the final result in terms of a regional susceptibility map based

% of area 0.4 6.1 14 79.5

Boundary weights 2.7–4.5 1.5–2.7 0.3–1.4 −3–0.3

Posterior probability N4.04 N1.27 N3.82 b3.82

× × × ×

5

10 105 106 106

Ratio 62.5 4.1 1.8 0.3

on posterior probabilities, but rather ways of gauging the predictability of potential debris flow source locations from the three fundamental terrain attributes of slope gradient, total curvature, and contributing catchment area (flow accumulation). Moreover, our results highlight a number of potential pitfalls that may compromise the interpretation of WofE in the context of regional susceptibility to debris flow initiation. These pitfalls are tied to (a) the choice and classification of predictors; (b) the combination of potentially interacting (= correlated) predictors; and (c) the choice and size of documented debris flow data for building the WofE model. We find that the decile classification scheme is more robust than a manual, and inadvertently arbitrary, classification of a given candidate predictor. This is because the uncertainty tied to the weights is evenly distributed among all classes, and the equal distribution of events over the classes allows for a more reliable test of conditional independence. In our case, the equal distribution of debris flows between classes ensured large enough sample sizes to increase the statistical power of the conditional independence test. The majority of studies using WofE or related statistical methods rarely justify their discretization of parameter space, and the same applies to aptly testing for variable interactions. In this context, the pairings of the decile classes showed conditional dependencies for the predictor variables total curvature and flow accumulation. The ROC curves indicate that the model combining slope and flow accumulation as predictors has a high predictive capacity and is preferable to any other tested model. Although we obtain several models with a prediction rate of ~80%, the true predictive performance and accuracy may be much less. This shows that ROC curves and prediction rates may point at model performances, but may contain little information about the precision, sharpness, or accuracy of the obtained models. Fig. 8 underscores how both the composition of the debris flow data set and choice of predictor variables influence the ROC curves and associated AUC values for the seemingly simple case of only three terrain metrics and several hundreds of documented events. To our knowledge, so far very few studies on debris flow or landslide susceptibility have explicitly considered the quality of their training and testing data sets. We find that models with similar predictive power can have large differences regarding their predictive accuracies. Nevertheless, the consistency of our results with those from other, more empirical based bivariate approaches is highly encouraging. Tests for conditional independence revealed that a combination of slope and contributing catchment area is preferable over a combination of slope and total curvature. This finding paves the way for comparing our findings with previous work on slope–area relationships in the context of catchment geomorphometry. Plotting the predictions of our model into slope–area space (Fig. 9) helps refine earlier predictions of debris flows from geomorphometric domains thought to reflect

Table 3 Spatial extent of process domains over the entire study area and captured debris flow types. Numbers in brackets indicate the ratio between the portion of debris flows and the portion of the entire study area in the corresponding process domain. Ratios N 1 indicate a stronger spatial relationship than normal distribution, ratios b 1 indicate less correlation. Process domain

% of study area

% of all debris flows

% of open-slope debris flows

% of channelized debris flows

Unchannelled valleys Hillslopes b 0.5 m−1 Hillslopes N 0.5 m−1 Alluvial channels Debris flow dominated channels

14.5 56.7 14.8 2.4 11.9

3.7 (0.26) 12.8 (0.23) 37.8 (2.55) 1.6 (0.69) 44.1 (3.71)

3.7 (0.25) 18.3 (0.32) 40.9 (2.76) 0.6 (0.26) 36.6 (3.08)

1.5 (0.1) 4.0 (0.07) 37.9 (2.56) 2.0 (0.85) 54.5 (4.59)

N.K. Meyer et al. / Geomorphology 207 (2014) 114–125

123

A

B

C

D

Fig. 10. Susceptibility map (A) with coloured classes “low”, “medium” and “high” and map showing the process domains (B) introduced in Fig. 9. (C) and (D): Close-ups of a small valley where debris flows were documented.

124

N.K. Meyer et al. / Geomorphology 207 (2014) 114–125

characteristic processes (Montgomery and Foufoula-Georgiou, 1993; Brardinoni and Hassan, 2006). We find that the WofE model using slope and flow accumulation as predictor variables does not only have a high predictive performance, but also is consistent with the geomorphic concepts of landscape partitioning, process domains and channel initiation. Fig. 9 demonstrates that the highest weights obtained by the WofE method correspond to debris flow source locations that cluster well within the proposed morphometric domain of debris flow dominance. The majority of locations clustering in this domain are channelized debris flows as opposed to those that originated on hillslopes. These open-slope debris flows are more difficult to predict since their occurrence is only loosely determined by flow accumulation, e.g., they may occur at all slope positions. In this context, we introduce a new process domain of debris flow dominated hillslopes, and underscore the importance of open-slope events that have not been captured in slope–area plots before (Table 3). Slope–area plots are bivariate models not originally aimed at quantifying debris flow susceptibility, but they offer a simple yet previously untested way of cross-validating purely terrainbased predictions of debris flow sources. The key question remains of whether high predictive performance may be gained by including additional predictors. We have intentionally focused on few and readily available terrain metrics. The variance of ROC curves in Fig. 8 indicates that the selection of predictors is likely to influence AUC values more than using a random subset of recorded events. Predictive performance thus seems more amenable to correctly choosing and classifying predictors rather than simply to increase their number and thus the potential for variable interaction, which may lead to a biased representation of debris flow susceptibility. Also, predictive performance cannot solely be measured by true positives, i.e. the fraction of debris flow occurrences that were correctly predicted. Given the rarity of such events, prediction of non-events becomes unduly straightforward, and may distort performance statistics. The ultimate answer to the optimal choice and classification of topographic predictors of debris flow source locations will be dictated by data availability and processing constraints eventually. Our results indicate that using fewer predictors can achieve predictive performances that rival those reported from the literature, though based on many more variables that we used here. 6. Conclusions We summarise the key findings of our study as follows: • Classification of predictor variables for correctly using the WofE method is essential. For our particular study we found that a simple quantile-based classification provided a more robust approach to defining class boundaries than did a subjective classification, simply because it assured sufficiently balanced sample sizes for each class and, thus, features a similar uncertainty related to each cell weight. • We found that for our study area in western Norway the two terrain variables of total curvature and flow accumulation are conditionally dependent, and thus excluded one of them in joint models using slope gradient. Repeated random subsetting of the debris flow inventory data set for training the WofE model affected the predictive performance less than did the choice of predictor variables. • The model using slope and flow accumulation as predictor variables shows a high predictive performance (AUC ~ 0.85) and is most robust concerning its predictive accuracy. This minimalist predictor pair is also consistent with previously postulated geomorphic process domains that are delineated on slope–area plots. We find that the highest weights derived from the WofE method cluster within the debris flow-dominated domain, and coincide mostly with channelized debris flows. We introduce a new morphometric process domain of open-slope debris flows that consistently complements previously proposed domains in slope–area plots.

Acknowledgements The study was conducted within the Norwegian research project InfraRisk (http://infrarisk.ngi.no). The authors were funded by the International Centre of Geohazards (ICG), the Norwegian Research Council, the Norwegian Public Roads Administration, the Norwegian National Rail Administration, and the Potsdam Research Cluster for Georisk Analysis, Environmental Change and Sustainability (PROGRESS). We computed statistics using the R software environment (www.r-project.org). References Agterberg, F.P., Cheng, Q., 2002. Conditional independence test for Weights-of-Evidence modeling. Nat. Resour. Res. 11, 249–255. Agterberg, F.P., Bonham-Carter, G.F., Wright, D.F., 1990. Statistical pattern integration for mineral exploration. In: Gaal, G., Merriam, D.F. (Eds.), Computer Applications in Resource Estimation Prediction and Assessment for Metals and Petroleum. Pergamon Press, Oxford, pp. 1–21. Aleotti, P., Chowdhury, R., 1999. Landslide hazard assessment: summary review and new perspectives. Bull. Eng. Geol. Environ. 58, 21–44. Bjordal, H., Helle, T.E., 2011. Skred og flom på veg. Report of the Norwegian Public Road Authority. SVV, Oslo. Blahut, J., Van Westen, C.J., Sterlacchini, S., 2010. Analysis of landslide inventories for accurate prediction of debris flow source areas. Geomorphology 119, 36–51. Bøe, R., Fossen, H., Smelror, M., 2010. Mesozoic sediments and structures onshore Norway and in the coastal zone. Nor. Geol. Unders. Bull. 450, 15–32. Bonham-Carter, G.F., 1994. Geographic Information Systems for Geoscientists, Modelling with GIS. Pergamon Press, Oxford. Brardinoni, F., Hassan, M.A., 2006. Glacial erosion, evolution of river long profiles, and the organization of process domains in mountain drainage basins of Coastal British Columbia. J. Geophys. Res. 111 http://dx.doi.org/10.1029/2005JF000358 (F01013). Brardinoni, F., Church, M., Simoni, A., Macconi, P., 2012. Lithologic and glacially conditioned controls on regional debris flow sediment dynamics. Geology 40, 455–458. Bråthen, S., Husdal, J., Rekdal, J., 2008. Samfunnsøkonomisk verdi av rassikring. Rapport 0801. Møreforsking, Molde. Dahal, R.K., Hasegawa, S., Nonomura, A., Yamanaka, M., Masuda, T., Nishino, K., 2007. GISbased weights-of-evidence modelling of rainfall-induced landslides in small catchments for landslide susceptibility mapping. Environ. Geol. 54, 311–324. Dunlop, S., 2010. Changing Climate: Establishing Relationships Between Meteorological Conditions and Rockslides in Southwestern Norway for the Purposes of Developing a Hazard Analysis. (Master thesis) Queen's University, Kingston, Canada. Fabbri, A.G., Chung, C.-J., 2008. On blind tests and spatial prediction models. Nat. Resour. Res. 17, 107–118. Fabbri, A.G., Chung, C.-J.F., Cendrero, A., Remondo, J., 2003. Is prediction of future landslides possible with a GIS? Nat. Hazards 30, 487–503. Fawcett, T., 2006. An introduction to ROC analysis. Pattern Recogn. Lett. 27, 861–874. Fell, R., Corominas, J., Bonnard, C., Cascini, L., Leroi, E., Savage, W.Z., 2008. Guidelines for landslide susceptibility, hazard and risk zoning for land-use planning. Eng. Geol. 102, 99–111. Fischer, L., Rubensdotter, L., Sletten, K., Stalsberg, K., Melchiorre, C., Horton, P., Jaboyedoff, M., 2012. Debris flow modeling for susceptibility mapping at regional to national scale in Norway. In: Eberhardt, et al. (Eds.), Landslides and Engineered Slopes: Protecting Society through Improved Understanding. Taylor & Francis Group, London, pp. 723–729. Guinau, M., Vilajosana, I., Vilaplana, J.M., 2007. GIS-based debris flow source and runout susceptibility assessment from DEM data — a case study in NW Nicaragua. Nat. Hazards Earth Syst. Sci. 7, 703–716. Guzzetti, F., Reichenbach, P., Ardizzone, F., Cardinali, M., Galli, M., 2006. Estimating the quality of landslide susceptibility models. Geomorphology 81, 166–184. Hungr, O., McDougall, S., Wise, M., Cullen, M., 2008. Magnitude–frequency relationships of debris flows and debris avalanches in relation to slope relief. Geomorphology 96, 355–365. Innes, J.L., 1983. Debris flows. Prog. Phys. Geogr. 7, 469–501. Jaedicke, C., Lied, K., Kronholm, K., 2009. Integrated database for rapid mass movements in Norway. Nat. Hazards Earth Syst. Sci. 9, 469–479. Kappes, M.S., Malet, J.P., Remaitre, A., Horton, P., Jaboyedoff, M., Bell, R., 2011. Assessment of debris flow susceptibility at medium-scale in the Barcelonette Basin, France. Nat. Hazards Earth Syst. Sci. 11, 627–641. Korup, O., 2012. Landslides in the Earth system. In: Clague, J.J., Stead, D. (Eds.), Landslides: Types, Mechanisms, and Modelling. Cambridge University Press, pp. 10–23. Korup, O., Stolle, A., 2014. Landslides and machine learning. Geol. Today (in press). Lee, S., 2013. Landslide detection and susceptibility mapping in the Sagimakri area, Korea using KOMPSAT-1 and weight of evidence technique. Environ. Earth Sci. http:// dx.doi.org/10.1007/s12665-013-2385-0 (in press). Lee, S., Choi, J., 2004. Landslide susceptibility mapping using GIS and the weight-ofevidence model. Int. J. Geogr. Inf. Sci. 18, 789–814. Mathew, J., Jha, V.K., Rawat, G.S., 2007. Weights of evidence modelling for landslide hazard zonation mapping in part of Bhagirathi. Curr. Sci. 92, 628–638. Minasny, B., McBratney, A.B., 2001. A rudimentary mechanistic model for soil formation and landscape development: II. A two-dimensional model incorporating chemical weathering. Geoderma 103, 161–179. Montgomery, D., 1999. Process domains and the river continuum. J. Am. Water Resour. Assoc. 35, 397–410.

N.K. Meyer et al. / Geomorphology 207 (2014) 114–125 Montgomery, D., 2001. Slope distributions, threshold hillslopes, and steady-state topography. Am. J. Sci. 301, 432–454. Montgomery, D., Foufoula-Georgiou, E., 1993. Channel network source representation using digital elevation models. Water Resour. Res. 29, 3925–3934. Norkart Virtual Globe, e. available on http://www.virtual-globe.info (accessed 29.04.2013). Oh, H.J., Lee, S., 2011. Landslide susceptibility mapping on Panaon Island, Philippines using a geographic information system. Environ. Earth Sci. 62, 935–951. Petley, D., 2012. Global patterns of loss of life from landslides. Geology 40, 927–930. Piacentini, D., Troiani, F., Soldati, M., Notarnicola, C., Savelli, D., Schneiderbauer, S., Strada, C., 2012. Statistical analysis for assessing shallow-landslide susceptibility in South Tyrol (south-eastern Alps, Italy). Geomorphology 151 (152), 196–206. Pourghasemi, H.R., Pradhan, B., Gokceoglu, C., Mohammadi, M., Moradi, H.R., 2013. Application of weights-of-evidence and certainty factor models and their comparison in landslide susceptibility mapping at Haraz watershed, Iran. Arab. J. Geosci. 6 (7), 2351–2365 http://dx.doi.org/10.1007/s12517-012-0532-7. Pradhan, B., Oh, H.J., Buchroithner, M., 2010. Weights-of-evidence model applied to landslide susceptibility mapping in a tropical hilly area. Geomatics Nat. Hazards Risk 1, 199–223. Prasannakumar, V., Vijith, H., 2012. Evaluation and validation of landslide spatial susceptibility in the Western Ghats of Kerala, through GIS-based weights of evidence model and area under curve technique. J. Geol. Soc. India 80, 515–523. Quinn, P.E., Hutchinson, D.J., Diederichs, M.S., Rowe, R.K., 2012. Regional-scale landslide susceptibility mapping using the weights of evidence method: an example applied to linear infrastructure. Can. J. Civ. Eng. 37, 905–927. Regmi, N.R., Giardino, J.R., Vitek, J.D., 2010. Modeling susceptibility to landslides using the weight of evidence approach: Western Colorado, USA. Geomorphology 115, 172–187. Regmi, A.D., Devkota, K.C., Yoshida, K., Pradhan, B., Pourghasemi, H.R., Kumamoto, T., Akgun, A., 2013. Application of frequency ratio, statistical index, and weights-ofevidence models and their comparison in landslide susceptibility mapping in Central Nepal Himalaya. Arab. J. Geosci. http://dx.doi.org/10.1007/s12517-012-0807-z (in press).

125

Romstad, B., 2013. A Model for Multi-hazard Risk Analysis to Transport Infrastructure (personal communication). Rössler, O., Bräuning, A., Löffler, J., 2008. Dynamics and driving forces of treeline fluctuation and regeneration in central Norway during the past decades. Erdkunde 62 (2), 117–128. Sassa, K., Canuti, P., 2008. Landslides: Disaster Risk Reduction. Springer, BerlinHeidelberg. Schorghofer, N., Rothman, D.H., 2002. Acausal relations between topographic slope and drainage area. Geophys. Res. Lett. 29, 11.1–11.4. Skrednett, 2013. http://www.skrednett.no (accessed 29.04.2013). Statens Kartverk, 2011. Terrengmodell – Landsdekkende digital terrengmodell i 10 eller 20 meters rutenett. Product sheet published by the Norwegian Mapping Authority (in Norwegian). Sterlacchini, S., Ballabio, C., Blahut, J., Masetti, M., Sorichetta, A., 2011. Spatial agreement of predicted patterns in landslide susceptibility maps. Geomorphology 125, 51–61. Stock, J., Dietrich, W.E., 2003. Valley incision by debris flows: evidence of a topographic signature. Water Resour. Res. 39, 1089 http://dx.doi.org/10.1029/ 2001WR001057. Süzen, M.L., Doyuran, V., 2004. A comparison of the GIS based landslide susceptibility assessment methods: multivariate versus bivariate. Environ. Geol. 45, 665–679. Thiery, Y., Malet, J.-P., Sterlacchini, S., Puissant, A., Maquaire, O., 2007. Landslide susceptibility assessment by bivariate methods at large scales: application to a complex mountainous environment. Geomorphology 92, 38–59. Van Westen, C.J., Rengers, N., Soeters, R., 2003. Use of geomorphological information in indirect landslide susceptibility assessment. Nat. Hazards 3, 399–419. Williams, K.M., 2012. Hillslope to Fluvial Process Domain Transitions in Headwater Catchments. (Dissertation) Montana State University, Bozeman, Montana. Wilson, J.P., Gallant, J.C., 2000. Primary topographic attributes. In: Wilson, J.P., Gallant, J.C. (Eds.), Terrain Analysis. Principles and Applications. Wiley and Sons, New York, pp. 51–86. Xu, C., Xu, X., Dai, F., Xiao, J., Tan, X., Yuan, R., 2012. Landslide hazard mapping using GIS and weight of evidence model in Qingshui River watershed of 2008 Wenchuan earthquake struck region. J. Earth Sci. 23, 97–120.