Application of Geomorphologic Factors for Identifying Soil Loss in Vulnerable Regions of the Cameron Highlands

Abstract

The main purpose of this study is to propose a methodology for identifying vulnerable regions in the Cameron Highlands that are susceptible to soil loss, based on runoff aggregation structure and the energy expenditure pattern of the natural river basin, within the framework of power law distribution. To this end, three geomorphologic factors, namely shear stress and stream power, as well as the drainage area of every point in the basin of interest, have been extracted using GIS, and then their complementary cumulative distributions are graphically analyzed by fitting them to power law distribution, with the purpose of identifying the sensitive points within the basin that are susceptible to soil loss with respect to scaling regimes of shear stress and stream power. It is observed that the range of vulnerable regions by the scaling regime of shear stress is much narrower than by the scaling regime of stream power. This result seems to suggest that shear stress is a scale-dependent factor, which does not follow power law distribution and does not adequately reflect the energy expenditure pattern of a river basin. Therefore, stream power is preferred as a more reasonable factor for the evaluation of soil loss. The methodology proposed in this study can be validated by visualizing the path of soil loss, which is generated from the hillslope process (characterized by the local slope) to the valley through a fluvial process (characterized by the drainage area as well as the local slope).

1. Introduction

Frequent extreme storms, such as heavy rainfall and typhoons, are causing landslides and slope failures, which result in massive soil loss now prevailing all over the world [1]. In the past, most landslides were concentrated in mountainous areas, and their frequency was so low that they did not receive much attention, compared to the other disasters like floods and earthquakes [1,2]. However, the changing climate has taken its toll on local rainfall patterns and characteristics, contributing to increased local rainfall amount and intensity [2]. Thus, occurrence of landslides is more frequent, and their severity has been threatening human life. In addition, rapid industrialization and urbanization have expanded our urban periphery to hillslopes, and there have been several tragedies for developments in this critical zone, without a proper engineering evaluation of slope stability [3].

Landslides occur simultaneously in a wide spatial range, and they are challenging to deal with once initiated, so it is recommended to minimize the effects of landslides by installing appropriate preventive facilities in the landslide-prone areas. For this purpose, the slope assessment for the prediction of landslides has to be conducted, in order to generate landslide risk maps. In fact, landslide risk maps have been prepared in many countries, to prevent and manage landslides at a national level. In United States, the standardized landslide risk maps have been provided and operated under the governance of the United States Geological Survey (USGS). In Malaysia, the Department of Mineral and Geosciences (DMG) and the Centre of Remote Sensing (MACRES) are responsible for carrying out research and informing the government of landslide-prone areas, while the Public Works Department (PWD) plays a major role in slope remedial works and the development of landslide assessment and management. However, there is still ambiguity in the landslide assessment system provided by the respective authorities, because there is room for the subjective judgment by the evaluator in measuring risk, by estimating the probability weighting factors for landslide occurrence. Therefore, in order to overcome these limitations, it is necessary to develop more scientific and clear landslide assessment techniques.

It is very important to know in advance the locations that are susceptible to erosion due to landslides and surface flows, for efficient management of the watershed in dealing with problems related to soil loss. A watershed is a natural system that provides a spatial basis for the hydrologic cycle. It is defined as the range of area that contributes to runoff at a single outlet located on any section of the river reach, due to precipitation [4]. Generally, a watershed is described as a complex transformation system, consisting of a channel network and hill slope, which is constantly evolving to efficiently drain inflows. The gradual evolution occurring in a watershed involves the complex interactions of numerous phenomena, ranging from macroscopic tectonic uplifts to fluvial erosion. It has been pointed out that explaining those complicated interactions using basic physical governing equations, such as Newton’s law of motion, is beyond the reach of present knowledge and technology [5]. Fractal geometry [6] and the theory of self-organized criticality (SOC) [7] have emerged as a scientific approach to analyze the complex phenomena of watersheds, in which huge systems with large spatial degrees of freedom are seen. Empirical relations, such as Horton’s Law and Hack’s law, have been widely adopted to describe the complexity of the geomorphological characteristics of a watershed, and Kim et al. [8] regarded this as explicit evidence to prove that the watershed is one of the most complex systems.

Power law distribution has been known as the most representative characteristic of a complex system, notably for scale invariance systems. Gutenberg-Richter’s law, Pareto’s law, and Zip’s law, which are frequently mentioned as typical examples of the power law distribution, have been widely adopted in various academic fields, such as geophysics, economics, and linguistics. Although these laws are depicted in different notations and applied in different fields, in fact, all of them are based on power law distribution [7]. In the field of hydrology, Rodriguez-Iturbe et al. [9] combined the random mass aggregation system and the concept of fractional Brown motion to describe a basic relationship for mass (or runoff) aggregation and energy expenditure structure in a watershed [10,11]. Since then, the power law distribution over the watershed has been constructed as a theoretical system, and has been actively studied as one of the representative geomorphological features [12,13,14,15].

The main objective of this study is to present a methodology for identifying soil loss in vulnerable regions, based on the analysis of runoff aggregation structure and the energy expenditure pattern of the natural watershed within the framework of power law distribution. For this purpose, the drainage area for each point in the watershed is selected as the major topographic factor, and energy dissipative factors of the watershed, such as shear stress and stream power, are defined as geomorphologic factors; these factors were analyzed using power law distribution. The soil loss in vulnerable regions in the Susu River Basin, which is located in the Cameron Highlands of Malaysia, was identified based on the relationship between local slope and drainage areas for assessing its slope stability condition.

2. Theoretical Background

2.1. Runoff Aggregation Structure of a River Basin

2.1.1. Random Walk Aggregation System

Takayasu et al. [10] demonstrated that an aggregation structure with constant mass injection tends to follow a power law distribution. Let us assume a one-dimensional lattice receiving particles in unit mass on every site at each time step. Random walk aggregation of those particles results in a dendritic trajectory over time steps (FIG.2 in Takayasu et al. [10], p. 3111, where a vertical axis represents time domain). This one-dimensional aggregation system corresponds to a two-dimensional stochastic model of river networks [16], by converting the vertical axis into a spatial one ($L_e$). Thereby, the magnitude of each aggregated mass should be proportional to its own drainage area (A) having an outlet at which two random walk trails collide, originating from the same source. A plan form of A can be defined as a closed curve area surrounded by both random walk trails (FIG.4 in Takayasu et al. [10], p. 3112).

According to the statistic feature of random walk, the boundary of A (FIG.4 in Takayasu et al. [10], p. 3112) has the same probabilistic properties as Brownian motion [17]:

$$p(x,L_e|x_0)=\frac{1}{\sqrt{4\pi D{L}_e}}{e}^{-\frac{{(x-x_0)}^2}{4D{L}_e}}$$

(1)

$$f(L_e)=\frac{\left|x_c-x_0\right|}{\sqrt{4\pi D{L}_e^3}}{e}^{-\frac{{(x_c-x_0)}^2}{4D{L}_e}}\approx \frac{\Delta x}{\sqrt{4\pi D{L}_e^3}}\propto L_e^{-\frac{3}{2}}$$

(2)

where $x$ is a location of an aggregated mass particle and $D$ denotes Brownian diffusivity. Equation (1) addresses the random behaviour of Brownian motion, starting at $x_0$ with $L_e=0$, expressed in the same form of the probability density function as normal distribution, with the mean $\mu =x_0$ and the variance ${\sigma }^2=2D{L}_e$ ($\sigma \propto L_e{}^{1/2}$). Equation (2) is a probability density function of the first passage of time ($L_e$), to determine when an aggregated mass particle will arrive at a specific location $x_c$ for the first time. In the case of $L_e\propto {(\Delta x)}^2/4D$ ($\Delta x=\left|x_c-x_0\right|$), $f(L_e)$ can be replaced with either the second (approximation) or the third term (proportional expression) on the right-hand side of Equation (2) [18]. To rewrite Equation (2) in terms of $A$, Takayasu et al. [10] assumed a geometric relation for the shape of the drainage area, on the basis of Equation (1):

$$A\propto W{L}_e\propto L_e{}^{\frac{3}{2}}$$

(3)

where $W$ is a width of $A$, being proportional to $L_e{}^{1/2}$. Therefore, a substitution of Equation (3) into Equation (2) gives rise to a probability density function of aggregated mass distribution with respect to the drainage area $f(A)$ and its corresponding exceedance probability, $F(A)(=1-{\displaystyle \int f(A)dA})$.

$$f(A)=f(L_e)\frac{d{L}_e}{dA}\propto A^{-\frac{4}{3}}$$

(4)

$$F(A)\propto A^{-\frac{1}{3}}$$

(5)

By conducting numerical experiments of the random walk aggregation system, Takayasu et al. [10] obtained identical exponent values as those in the right-hand sides of Equations (4) and (5).

2.1.2. Runoff Aggregation System

Rodriguez-Iturbe et al. [9] combined the random walk aggregation system by Takayasu et al. [10] with a concept of fractional Brownian motion [11], in order to provide a power law distribution for the drainage area of river basins:

$$p(x,L_e|x_0)=\frac{1}{\sqrt{4\pi D_H{L}_e}}{e}^{-\frac{{(x-x_0)}^2}{4D_H{L}_e}}=\frac{1}{\sqrt{4\pi D{L}_e{}^{2H}}}{e}^{-\frac{{(x-x_0)}^2}{4D{L}_e{}^{2H}}}$$

(6)

$$f(L_e)\propto L_e{}^{-(H+1)}$$

(7)

where $H$ is the Hurst exponent, and $D_H(=D{L}_e^{2H-1})$ represents fractional Brownian diffusivity. In fact, Equations (6) and (7) address the stochastic behaviour of fractional Brownian motion, each having a similar meaning to Equations (1) and (2), respectively. However, being a geometric relation to $A$, Equation (3) should be transformed into Equation (8), due to an effect of $D_H$ for random walk trajectory:

$$A\propto W{L}_e=L_e{}^{H+1}$$

(8)

where $W$ is a width of $A$, similarly to Equation (3) but proportional to $L_e{}^H$. By the same rationale as Equations (4) and (5), a probability density function for drainage area of river basin $f(A)$ and its exceedance probability distribution $F(A)$ can be derived as:

$$f(A)\propto A^{-{\alpha }_A};{\alpha }_A=\frac{2H+1}{H+1}$$

(9)

$$F(A)\propto A^{-\left({\alpha }_A-1\right)}$$

(10)

where α_A is the exponent of the power law distribution for the drainage area, and $H$ is the Hurst exponent, which defines the behaviour of the fractional Brownian motion of the divide of the catchment plan form, or the trajectory of the stream line. Equations (4) and (5), or Equations (9) and (10), represent the fundamental relationships for such power law distributions as Gutenberg-Richter law and Zip’s law.

2.2. Energy Expenditure within a River Basin

2.2.1. Energy Expenditure Process of an Open Channel

Let us consider water flow in an open channel, with length being $L$ and bottom slope being $S$(S $=\tan \theta$, where $\theta$ denotes the bottom slope angle), having a rectangular cross-section with width $w$. In the case of uniform flow, with water depth being $d$, a dynamic balance can be written between gravity and friction force in terms of equilibrium condition:

$$\rho g\cdot wdL\cdot S=\tau \left(w+2d\right)L$$

(11)

where $\rho$ is a density of water, and $g$ denotes gravitational acceleration. If $\theta$ is small enough ($\tan \theta \approx \sin \theta$), the left-hand side of Equation (11) represents the flow directional component of gravity. And $\tau$ in the right hand side of Equation (11), meaning shear stress between water body and bottom surface, which can be expressed as follows under the assumption of incompressible turbulent flow:

$$\tau =\rho gRS=C_f\rho v^2$$

(12)

where $C_f$ is a non-dimensional resistance coefficient, and $v$ represents average flow velocity. After some manipulations with Equations (11) and (12), $S$ can be expressed as:

$$S=\frac{C_f{v}^2}{Rg}$$

(13)

where $R$($=wd/(w+2d)$) denotes the hydraulic radius. By comparing Equation (13) with the Darcy-Weisbach law, relation $S$ can be interpreted as a kind of energy loss per weight per unit length within an open channel. Based on Equations (12) and (13), general relations for shear stress and energy expenditure in open channel water flow can be formulated in the form of proportional expression:

$$\tau \propto RS$$

(14)

$$P\propto QS$$

(15)

where P denotes energy expenditure by water flow with discharge $Q$, implying stream power and the time rate of energy expenditure per unit length.

2.2.2. Inferential of Geomorphological Factors from DEM

A geomorphologic factor I, which quantifies the properties of a watershed based on DEM, can be expressed as:

$$I=k{Q}^m{S}^n$$

(16)

where k is the constant, and m and n are the exponents. Equation (16) can be transformed into various forms depending on the properties of the watershed in which one is interested. Equation (17) is the relationship which is often applied to the delineation of river networks:

$$Q_i\propto A_i$$

(17)

where A_i can be interpreted as the surrogate variable for Q at the ith pixel. One can substitute the surrogate variable obtained in Equation (17) into Equation (16), to infer the geomorphologic factor at ith pixel. Thus, the shear stress and stream power at the ith pixel can be inferred using Equations (18) and (19):

$${\tau }_i\propto A_i^b{S}_i$$

(18)

$$P_i\propto A_i{S}_i$$

(19)

where the value of exponent b can be inferred from the hydraulic geometry relationships [19] in the form of power functions:

$$w=c_1{Q}^a;d=c_2{Q}^b;v=c_3{Q}^e$$

(20)

where C₁, C₂, and C₃ are the constants and a, b, and e are the exponent. With the consideration of Equation (17) and the relationship between R, α, and d, the b in Equation (18) can be thought to be the same as the b in Equation (20). Leopold et al. [19] also presented the exponent values as a = 0.5, b = 0.4, and e = 0.1 (a + b + e = 1). Assumptions have been made in previous studies [9,20,21,22] that a b of 0.5 should be adopted to assess the shear stress.

It is widely known that there exists an empirical relation between the channel slope and drainage area within river basins [23,24].

$$S\propto A^{-\xi }$$

(21)

where $\xi$ is an exponent. Equation (21) explains the scaling property of the channel slope with the magnitude of the drainage area. By using Equations (18), (19), and (21), the probability density function for shear stress and stream power can be defined as following:

$$f(\tau )=f(A)\frac{dA}{d\tau }\propto {\tau }^{-{\alpha }_{\tau }};{\alpha }_{\tau }=\frac{1}{b-\xi }\left({\alpha }_A+b-\xi -1\right)$$

(22)

$$f(P)=f(A)\frac{dA}{dP}\propto P^{-{\alpha }_P};{\alpha }_P=\frac{1}{1-\xi }\left({\alpha }_A-\xi \right)$$

(23)

Equations (9), (22), and (23) are the basic relationships of the power law distribution, and Rodriguez-Iturbe et al. [9] found that the relationships between A and P are universally established for several watersheds with different characteristics of geological condition, climate, vegetation, etc.

2.3. Power Law Distribution

The power law distributions in proportional expressions like Equations (9), (22), and (23) can be converted into a more complete form of power function [25,26]:

$$f(x)=C{x}^{-\alpha }$$

(24)

where $C$ is a normalization constant, and $\alpha$ being an exponent. Equation (24) goes to infinity as x approaches zero, implying the existence of the lower limit x_min. Within the range of $x_{min}<x<\infty$, the complete form of $f(x)$ can be derived as:

$$f(x)=\frac{\alpha -1}{x_{min}}{\left(\frac{x}{x_{min}}\right)}^{-\alpha }$$

(25)

Equation (25) is a full mathematical expression of the probability density function of power law distribution, which gives rise to a closed form of statistical moments for x:

$$E\left[x^z\right]={\displaystyle {\int }_{x_{min}}^{\infty }{z}^{z}f(x)dx}=\{\begin{array}{cc}\infty & ;\alpha <z+1\\ \left(\frac{\alpha -1}{\alpha -1-z}\right)x_{min}^z& ;\alpha >z+1\end{array}$$

(26)

where z represents the order of statistical moment. It can be seen that Equation (26) depends on the magnitude of $\alpha$. All of the statistical moments for x go to infinity when $\alpha$ is less than 2. This is a typical feature of the complex system without a representative scale. In addition, an analytical expression for exceedance probability $F(x)$ also can be derived from Equation (25):

$$F(x)={\left(\frac{x}{x_{min}}\right)}^{-(\alpha -1)}$$

(27)

It is interesting that $F(x)$ is also in the form of power function, as is $f(x)$, and both of them are in a linear form on a double logarithmic paper.

3. Materials and Methodology

3.1. Overview of Target Area

The Susu River Basin is located at Cameron Highlands, Malaysia. The Susu reservoir is located downstream of the river basin, for hydropower generation and flood control purposes. There are three major rivers within the river basin: the Bertam River, the Telom River, and the Lemoi River. A gravity dam 88 m high and 500 m long has been constructed to impound the Bertam River, leading to the creation of the Susu reservoir. The contributing area of the Susu reservoir is about 230 km², as shown in Figure 1. The DEM data used in this study is from the Shuttle Radar Topography Mission (SRTM) DEM produced by the National Aeronautics and Space Administration (NASA) and the National Geospatial-Intelligence Agency (NGA), with grid resolution of 30.7 m × 30.7 m. Arc GIS was used to conduct the terrain analysis of the basin. First and foremost, the DEM data was processed by executing the fill function for removing any undesired sink, to ensure proper delineation of basins and streams. Secondly, the flow direction function was executed to determine the direction of flow from every cell in the raster, and an eight-direction (D8) approach was adopted. Consequently, the local slope of each cell within the basin was computed by dividing the differences of elevation by flow path length between two cells. Finally, a flow accumulation function was carried out to calculate the accumulated flow in each cell, whereby cells with high flow accumulation could be recognized as stream channels, while cells with zero accumulation could be recognized as ridges. Figure 1 shows the geographical location, elevation map, and three-dimensional plot for the Susu River Basin.

3.2. Parameter Estimation for Power Law Distribution Using Maximum Likelihood

The power functions can be generally fitted to the observations by the least square. However, in the case of power law distribution, the least-square-based fitting methods are likely to result in biased estimates [25,26]. As an alternative, Clauset et al. [26] proposed fitting maximum likelihood instead of least-square. Therefore, we adopted the approach of maximum likelihood in this study, in order to estimate the parameter of A_min (drainage area), τ_min (shear stress), P_min (stream power), and α, as presented in Equation (25). The procedures of parameter estimation using maximum likelihood are elaborated as following. Firstly, let us consider a set containing n observations for x (x_i $\ge$ x_min, i = 1, 2, …, n), where the observations follow power law distribution. Based on Equation (25), the likelihood of A, τ, and P can be estimated as:

$$L(x|\alpha )={\displaystyle \prod _{i=1}^n\frac{\alpha -1}{x_{min}}{\left(\frac{x_i}{x_{min}}\right)}^{-\alpha }}$$

(28)

where ${\displaystyle \prod }$ is product operator. Taking the logarithm of both sides of Equation (28), setting $\partial \ln L(x|\alpha )/\partial \alpha =0$, and solving for $\alpha$, the maximum likelihood estimate of $\widehat{\alpha }$ can be derived as follows:

$$\widehat{\alpha }=1+n{\left[{\displaystyle \sum _{i=1}^n\ln \left(\frac{x_i}{x_{min}}\right)}\right]}^{-1}$$

(29)

The application of Equation (29) to real-world observations requires an a priori determination of x_min. To this end, Clauset et al. [26] suggested Kolmogorov-Smirnov (KS) statistic D for the determination of x_min:

$$\min D=\underset{A\ge A_{min}}{\max }\left|S(x)-F(x)\right|$$

(30)

where S(x) and F(x) denote exceedance probabilities, the former from the observation and the latter from Equation (27). We are assuming a certain x_min can produce a pair of $\widehat{\alpha }$ and D by Equations (29) and (30). Therefore, if this estimation process is repeatedly performed on all of the assumptions for x_min or all of x_i, a final pair of $\widehat{\alpha }$ and x_min can be obtained by minimizing the D. Clauset et al. [26] demonstrated the applicability of Equations (29) and (30) to power law distribution fitting by numerical experiments. In this study, we adopted the methodology as discussed above to fit the drainage area, shear stress, and stream power of the study area to power law distribution.

3.3. Classification of Critical Regions Based on Mechanisms of Channel Initiation

The channel head (source of a stream) has been known to occur near the geomorphologic threshold, which is the boundary between diffusive erosion and incisive erosion [22]. This implies that the source of a stream is located at the transition point of those erosion processes. The watershed would automatically distribute the mass (precipitation) and energy exerted from the exterior in a separate way, eventually resulting in different forms of the landscapes.

McNamara et al. [22] proposed a methodology for classifying the watershed into multiple regions, based on the above-mentioned principle of channel initiation by relating local slope s (the steepest downward slope on eight-directional pixel windows) to A and subsequently extracting points susceptible to soil losses. The FIG. 7 in McNamara et al. [22] (p. 153) shows the relationship between s and A, in which each respective region is classified by superimposing the transition points (thresholds) of s (s₁), A (A₁, A₂, A₃, A₄), and energy index EID (EID₁, EID₂), using the high-resolution DEM data of the Pang Khum watershed in Thailand. McNamara et al. further suggested that the thresholds of the watershed correspond to changes in the complementary cumulative distribution curves.

Table 1. Stability condition of slope–area behavioural regimes based on channel initiation.

Behavioural Regimes	Energy	Slope	Area	Stability Condition
1a	Insufficient	-	-	Divergent hillslope
1b	Insufficient	-	-	Convergent valley
2a	Sufficient	Flat	Small	Stable zone
2b	Sufficient	Steep	Small	Landslide active
2c	Sufficient	Steep	Large	Overland flow active
2d	Sufficient	Flat	Large	Groundwater active
3	Excessive	-	-	River channels

summarizes the results of FIG. 7 in McNamara et al. [22] (p. 153), and we applied this mechanism to our study area, to identify the regions that are susceptible to soil losses. However, there is a notable point in McNamara et al.’s study that requires further consideration, prior to adopting their methodology in our study area. The authors applied Equation (18) to estimate the energy index, and it was discussed in Section 2.2 that Equation (18) actually represents shear stress rather than energy. Therefore, there are some contradictions in describing the state of energy, as depicted in Table 1, if Equation (18) is used to define the energy index (EID) of a watershed. To this end, besides using Equation (18), we are also using Equation (19), which represents stream power as the energy index, and both of these geomorphological factors are assessed in order to determine which energy index is more suitable to be used as the energy threshold for the classification of the behavioural regimes.

4. Case Study of the Susu River Basin in Cameron Highlands

4.1. Relationship between Local Slope and Drainage Area

Based on the topographic analysis results obtained from Figure 1 and using the procedures as explained in Section 3.1, the relationship between the local slope and drainage area is illustrated in Figure 2. As mentioned frequently in previous studies [21,27], there is a tremendous scatter in this plot. Figure 3 portrays the slope–area plot after the data was binned and grouped, in order to reduce the scattering characteristic observed in the original data. Four distinguishable regions can be observed from the figure as the data has been filtered: Region 1 indicates that the local slope increases with drainage area until the turnover is achieved at threshold A₁. Divergent topography is seen in this region, which is typically dominated by diffusive erosion process, such as rain splash, soil creep, or sheet wash [22]. Local slope decreases with drainage area in Region 2, representing convergent topography until threshold A₂ is achieved. The slope–area relationship in Region 3 resembles the transition section between landscape surfaces and river channels, where soil loss is active. The local slope decreases with the drainage area in Region 4, following a power law with an exponent of 0.28, and Flint [24] has recognized this region as one with fluvial channels.

Special attention has to be paid to the transition region between Region 1 and Region 2. There should be a turnover in threshold A₁ indicating a transition from a divergent landscape to a convergent one. However, turnover cannot be observed in Figure 2. Even though turnover is observed after the data has been binned, as illustrated in Figure 3, there is only one data point located in Region 1, leading to the ambiguous slope–area relationship in this region. This phenomenon is due to the resolution of the DEM (about 30 × 30 m) used in this study. In other words, it is impossible for us to accurately identify the position of the turnover if the scale of A₁ is smaller than the area of a single pixel (about 946 m²). Kim et al. [28] plotted the relation of local slope and drainage area for the Choyang Creek Basin, using DEM data with resolution of 90 × 90 m; in addition, no turnover has been observed in the transition region between Region 1 and Region 2, even though the data has been binned. Therefore, it is suggested to adopt DEM data with a high resolution, in order to obtain the prominent trend for the relationship between local slope and drainage area. As there were insufficient data points observed from our data to compare with the slope–area plot derived in previous studies [22,29], we decided to set the turnover at drainage area A₁ to 1893 m², as a positive slope had been observed from 946 m² to 1893 m², and a negative slope was observed when the drainage area was greater than 1893 m². It is also illustrated in Figure 3 that the threshold area A₂ for Region 2 was 31,464.8 m²; A³ for Region 3 was 202,510.3 m²; and the mean slope s₁ calculated in Region 3 was 0.2174.

4.2. Analysis of Behaviour of Geomorphological Factors

Figure 4 is the result of fitting the complementary cumulative distribution of the geomorphologic factors of the Susu River Basin: A, τ, and P, to the power law distribution. τ and P were computed using Equations (18) and (19), where 0.5 has been applied for the exponent value of b in Equation (18). It is noted that all three distribution curves show linear characteristics at the central part of the curves in a log–log space, but they behave differently (non-linearly) in the vicinity of the origin and at the tail of the curves. In this study, we selected the points which those factors deviate from the linear characteristic at the boundary of each region, as illustrated in Figure 4.

Table 2. Parameter estimation of power law distribution.

I	α	xmin
A	αA = 1.50	Amin = 23184.6 m2
τ	ατ = 2.85	τmin = 156.6 m
P	αP = 1.66	Pmin = 4444.3 m2

shows the parameter estimation of the power law distribution for A, τ, and P using the maximum likelihood approach. The maximum likelihood approach has been adopted for parameter estimation of power law, as it avoids the subjective assumption seen in conventional least-square approximation, used for determining the initiation point of linear relationship (x_min). We can observe from the estimated α value that A and P have similar characteristics, while τ behaves differently. According to the basic properties of power law distribution, the statistical moment cannot be defined if the estimated α value is less than 2 [25,26]. The results revealed that both α_A and α_P were less than 2. Thus, the statistical moment cannot be defined for A and P, for all orders, while α_τ is greater than 2, implying that the first order of the statistical moment can be defined for τ. Therefore, it can be concluded that geomorphologic factors for A and P are scale-invariant, but τ is a scale-dependent (finite scale) geomorphologic factor.

Kim et al. [28] plotted a complementary cumulative distribution curves for the Choyang Creek Basin in China, using DEM data with a resolution of 90 × 90 m. They similarly concluded that shear stress is a scale-dependent geomorphologic factor, as the α value obtained for shear stress in their study is 4.42 (greater than 2). The α values estimated for the drainage area and stream power in Kim et al. [28] were 1.50 and 1.78, respectively, which is quite similar to the values estimated in the current study. It can be seen that the resolution of the DEM data influences the shape of the complementary cumulative distribution curve, most notably at the tail of the curve. It can be easily assessed from the shape of cumulative distribution curve for the shear stress of the Choyang Creek Basin, in that it shows a scale-dependent trend, as the tail of the curve deviates upward independently from the characteristic of scale. However, we adopted DEM data with a resolution of 30.7 × 30.7 m in this study, and found that the shape of the cumulative distribution curve for shear stress of the Susu River Basin shows a finite scale effect on scale-independent variables. Thus, it can be concluded that the shape of the cumulative distribution curve for shear stress might depend on the resolution of the DEM adopted for topographic analysis. In other words, it is difficult to determine that whether shear stress is a scale-dependent factor just by assessing the shape of the curve. The α values listed in Table 2 provide us with the quantitative approach for determining the scale dependency, whereby a factor is considered as scale-dependent if an α value of greater than 2 is estimated and, as elaborated in previous paragraph, if the α value for shear stress obtained in this study is greater than 2. Therefore, shear stress is recognized as a scale-dependent factor, even though its shape of cumulative distribution curve is obeying the characteristic of scale independence.

4.3. Validation of Exponent Estimates

The exponent values in Table 2 were validated based on the outcomes drawn from previous research. In this study, α_A was estimated to be about 1.50, which is somewhat larger than the value estimated in the previous study. Indeed, Rodriguez-Iturbe et al. [9] suggested an α_A of 1.45. Kim et al. [8] inferred that this discrepancy might stem from the difference between the methodologies used to estimate α_A, in which the least-squares method had been applied to fit the power law distribution, until the emergence of maximum likelihood method proposed by Clauset et al. [26]. In order to examine the variations induced by different fitting approaches, the least-squares method was applied to the data points greater than A_min in Figure 4a, resulting in Figure 5, where the black dotted line is the fitting curve (upper regression), excluding the points at the tail part, while the solid orange line (lower regression) includes the points at the tail part. The effect of the tail part can be confirmed through the value of the exponent. In particular, we can estimate α_A of 1.47 from the upper regression. This is very similar to the value from previous studies by least squares, and thus it validates the estimations presented in this study.

Rodriguez-Iturbe et al. [9] proposed an α_P of 1.90, which can be easily deduced from Equation (23). This value of α_P can be obtained if α_A and ξ are assumed to be 1.45 and 0.5, respectively. Tarboton et al. [30] confirmed this result by suggesting an ξ value of 0.5 from link-based terrain analysis, which is consistent with the information mentioned above. In this study, α_P was estimated as 1.66, as shown in Table 2, which is significantly smaller than the value obtained from previous studies. By using the α_A value of 1.50 obtained in this study, an ξ of 0.24 can be obtained according to Equation (23), which is very similar to the exponent of regression for Region 4 in Figure 3. More interestingly, ξ of 0.23 can be obtained from Equation (22), by applying an α_τ value of 2.85, as shown in Table 2, further confirming the validity of the exponents estimated in this study.

4.4. Classifications of Behavioural Regions

Figure 6 illustrates a slope–area plot that is superimposed with the slope threshold (s), area thresholds (A), and energy threshold (τ and P). Each behavioural region has been classified according to the stability conditions stated in Table 1, and the geomorphologic thresholds are listed in

Table 3. Geomorphologic thresholds in the space of local slope and drainage area.

EID	τ1	156.6 m
	τ2	1604.1 m
	P1	4444.3 m2
	P2	2,250,964.0 m2
s	s1	0.2174
A	A1	1893 m2
	A2	27,324.7 m2
	A3	202,510.3 m2
	A4	100,608,888 m2

. It has been mentioned in Section 3.2 that shear stress (τ) and stream power (P) were used as the energy index (EID) in our study, and two thresholds have been determined for each factor, as presented in Figure 4b (τ₁ and τ₂) and Figure 4c (P₁ and P₂), according to the parameter obtained in Table 2. As mentioned previously, the slope threshold (s₁) has been computed by averaging the local slope in Region 3 for the data, as shown in Figure 3. Four area thresholds have been identified: A₁ and A₃ were identified from the boundaries determined for Region 1 and Region 3 in Figure 3; A₄ was identified through the cumulative area distribution curves, as shown in Figure 4a, where A₄ represents the point of the curve starting to deviate from the linear relationship; and A₂ was determined according to the approach proposed by McNamara et al. [22] and Ijjasz-Vasquez el al. [27], by averaging the A_min in Table 2 and the boundary determined for Region 2 in Figure 3. The statements below summarize the classification process, as well as the elaborations for the respective behavioural regions:

(1)

The energy thresholds of shear stress (τ₁ and τ₂) and stream power (P₁ and P₂) have divided the slope–area plot (Figure 6a,b) into three distinct regions. The energy state in each respective region has been described in Table 1, such that energy was insufficient for Region 1 but sufficient in Region 2. Energy in Region 3 was excessive, which typically corresponds to river channels. From here, we could determine that soil loss of the ground surface was active in Region 2, where there was sufficient energy for soil motion.

(2)

Region 2, as elaborated in Equation (1), was further divided into two sub-regions, which were separated by the slope threshold (s₁). It can be interpreted from Table 1 that the lower sub-Region 2 represents river channels, as its slopes are lower than the slope threshold. Therefore, soil loss of ground surface will only be able to occur at the upper sub-Region 2.

(3)

The slope–area plot was also divided into five regions, by four area thresholds (A₁, A₂, A₃, and A₄). Figure 6 summarises the behavioural regions, which have been classified using the geomorphologic thresholds as discussed above. It can be concluded that the parts which are susceptible to soil loss are located in Regions 2b and 2c, where 2b is steep-slope dominant and 2c is overland-flow dominant.

We have reiterated in the previous section that we have been adopting two types of EID, instead of a single EID in this study. The behavioural regions classified using τ were different from the ones classified using P as the energy threshold. It can be visually observed from the comparison between Figure 6a and Figure 6b that the Region 2 classified using τ was much narrower than the one classified using P. Furthermore, sub-Region 2a, as proposed in Table 1 and FIG. 7 in McNamara et al. [22] (p. 153), was not even defined in Figure 6a, and no data was found in sub-Region 2b. A similarity has been noticed when comparing Figure 6b (with P as the energy threshold) to FIG. 7 in McNamara et al. [22] (p. 153).

The inconsistency of results obtained from using τ as the energy threshold might be due to the findings discussed in Section 4.2, that τ is a scale-dependent, factor which does not follow a power law distribution. In fact, the tiny Pang Khum watershed area considered in McNamara et al. [22] was just about 0.973 km²_, with an energy threshold of τ₁ = 14 m and τ₂ = 87 m. However, the watershed area considered in this study (the Susu River Basin) was about 230 km², which is about 200 times greater than the area of the Pang Khum watershed; in addition, the energy threshold of τ computed from this study is about 10 times greater than that of the Pang Khum watershed. This implies that τ, which is presented by McNamara et al. [22] as the EID, does not adequately reflect the magnitude of the distribution of energy in the watershed under the framework of the power law distribution. Therefore, in this study we suggest that the geomorphologic factor of stream power P as the suitable EID, rather than shear stress, as well as τ for the evaluation of soil loss under the framework of the power law distribution.

4.5. Illustration of Regions Susceptible to Soil Losses

Figure 7 and Figure 8 show the spatial distribution of regions susceptible to soil loss for the Susu River Basin, overlaid with DEM based on the geomorphological factors of shear stress and stream power, as presented in Figure 6. In Figures Figure 7a and Figure 8a, each lattice represents the behavioural region to which it belongs, while Figure 7b and Figure 8b illustrate the regions susceptible to landslides (Region 2b: cyan) and active overland flow (Region 2c: red). It can be seen that the extent of critical regions (Regions 2b and 2c) observed in Figure 7 is smaller than the one observed in Figure 8, and no data point has been observed in sub-Region 2b for characteristic regime classification using shear stress. This phenomenal is due to the distribution characteristic of τ (Figure 4b), as described in previous section. Assuming that τ was distributed according to the power law distribution and overestimation of the lower limit of energy threshold, τ₁ would eventually lead to underestimation of the extent of soil loss-vulnerable regions. Therefore, the energy threshold of stream power (P) is recognized as the appropriate factor, rather than shear stress (τ) as proposed in McNamara et al. [22], for evaluation and identification of soil losses vulnerable regions.

Figure 9 shows the three-dimensional illustration of particularly vulnerable soil loss regions in the Susu River Basin, as presented in Figure 8b. There is an interesting point that can be observed from this three-dimensional illustration: the movement path of the soil loss generated from Region 2b (cyan pixels), which is susceptible to landslides to the valley through Region 2c (red pixels), is seen clearly, validating the methodology proposed in this study for identification of critical regions susceptible to soil loss. In Figure 9, it can been seen that the enlarged images indicate the movement paths of debris flow (red pixels), which carry soil particles generated from landslide-prone regions through the accumulated overland flow.

5. Conclusions

In this study, the runoff aggregation structure and energy distribution pattern (according to power law) within the watershed of the Susu River Basin has been studied, and the outcomes have been used to estimate and extract the critical regions that are susceptible to landslides and slope disasters. In other words, this study presented a new methodology for assessing the vulnerable soil loss regions within a watershed, using the geomorphologic factors of drainage area, shear stress, and stream power as the threshold parameters. GIS has been extensively adopted in this study to extract the drainage area, as well as the geomorphologic factors of shear stress and stream power, which defines the hydrodynamic force of each pixel within the watershed of the Susu River Basin. The magnitude distribution of these geomorphologic factors was fitted to power law distribution, and the results obtained were used to define the boundary thresholds of respective characteristic regimes. The outcomes of this study can be summarised as follows:

(1)

It can be clearly visualized that the critical regions classified by using shear stress are narrower than those classified by using stream power as the indicator. This is due to the fact that shear stress (τ) is a scale-dependent factor, which does not follow the power law distribution. Thus, shear stress is recognized as the geomorphologic factor, which does not adequately reflect the scale of the natural watershed energy in the framework of the power law distribution. Therefore, this study proposed the identification of vulnerable soil loss regions based on the geomorphologic factor of stream power, as it was found to be scale-independent and reflects the scale of the natural watershed energy more appropriately than shear stress.

(2)

The selected geomorphologic factors—namely drainage area, shear stress, and stream power—of each pixel within the Susu River Basin were fitted to the power law distribution, and the approach of maximum likelihood has been adopted to estimate the parameters of α and x_min, respectively. The use of maximum likelihood provides a useful alternative for parameter estimation of power law, as it avoids the subjective assumption used in conventional least-square approximation for determining the initiation point of the linear relationship (x_min).

(3)

Regions that are landslide-active (2b) and overland flow-active (2c) have been identified by adopting stream power as the threshold parameter. In contrast, regions of landslide activity (2b) could not be identified when shear stress was adopted as the threshold parameter. This might be due to the explanations mentioned in the outcome of Equation (1). The critical regions of landslide and overland flow activity within the watershed have been identified and illustrated in three dimensions, as shown in Figure 9, and the validity of the methodology proposed in this study has also been visually verified through the movement paths of the soil particles generated from the vulnerable soil loss regions to the valley, through the accumulation of overland flow.

Lastly, we realize that accuracy assessment with data collected on site is required to verify the critical regions of landslide activity and overland flow actives as identified and shown in Figure 9. Therefore, the direction of future study will focus on the accuracy assessment of the results concluded from this study, with the field work data and the results computed from other methods.