Numerical Analysis of the Effects of Saturation Overland Flow and Spatial Variability of Soil Strength on Sediment Movement Processes in Headwaters

Numerical Analysis of the Effects of Saturation Overland Flow and Spatial Variability of Soil Strength

on Sediment Movement Processes in Headwaters

Masayuki Hata ^1* , Taro Uchida ² , Junichi Kanbara ² , and Masahiro Kaibori ³

1Chuden Engineering Consultants Co., Ltd (Deshio, Minami-ku, Hiroshima, Hiroshima 7348510, Japan) 2National Institute for Land and Infrastructure Management (Asahi, Tsukuba, Ibaraki 3050804, Japan)

3Graduate School of Integrated Arts and Sciences, Hiroshima University (Kagamiyama, Higashihiroshima, Hiroshima 7398521, Japan)

*Corresponding author. E-mail: hata@cecnet.co.jp

Sediment movement in headwaters includes a variety of processes, such as landslides, debris flow, and surface erosion. For predictions of the occurrence of shallow landslides, several physically based models have been proposed. However, in most of these models, the effects of saturation overland flow and the spatial distribution of soil strength have not been taken into consideration. Thus, the accuracy of sediment yield prediction is not adequate in relatively gentle slopes (i.e., 15°–30°). Therefore, we confirmed the effects of saturation overland flow and spatial variability of soil strength, respectively, and found considering them simultaneously is efficient to improve the accuracy of sediment yield prediction in headwaters. We suggest using the advanced H-SLIDER method with a simple classification of cohesion according to the scaling relation of local slope and upslope areas. This method is efficient, especially in the areas where the upslope drainage area is relatively large or the local slope gradient is not very gentle.

Key words: shallow landslide, sediment movement, saturation overland flow, destabilization, advanced H-SLIDER method

1. INTRODUCTION

Prediction of volume of sediment yield is one of the key issues for planning countermeasures against sediment disasters, such as debris flow and landslides. Many empirical relations between sediment yield volume and topography and geology have been proposed. However, to date, no widely used physically based model for prediction of sediment yield volume is available.

Hillslopes and channels have different characteristics in terms of runoff processes and soil layer formation processes. In channels, even if no perennial overland flow occurs, overland flow is generated at times of heavy rain. On the other hand, in hillslopes, subsurface flow is a dominant contributor to runoff generation (Fig. 1) [e.g., Uchida et al., 2003].

In hillslopes, the soil is generated by weathering at the original position and is modulated by soil creep [e.g., Heimsath et al., 1997]. In channels,

except for exposed rock areas, the soil consists mainly of deposition from upper hillslopes, and this soil is movable [e.g., Tsukamoto, 1973].

Thus, different concepts for processes of sediment yield in hillslopes and channels have been proposed. In hillslopes, the shallow landslide is a dominant process of sediment yield and is modeled by combining slope stability analysis and a subsurface flow model [e.g., Okimura et al., 1985;

Hiramatsu et al., 1990; Montgomery and Dietrich, 1994; Wu and Sidle, 1995; Kosugi et al., 2002;

Uchida et al., 2009]. These models ignored the effects of saturation overland flow, and the accuracy of sediment yield prediction is not adequate in relatively gentle slopes (i.e., 15°–30°). Takahashi [1991] proposed an equation for evaluation of the effects of saturation overland flow depth on initiation of debris flow in channels. However, this equation regarded the effects of soil cohesion as negligible.

In headwater catchments, there are areas known

as zero-order basins and unchanneled hollows [e.g., Tsukamoto et al., 1973; Montgomery et al., 2009].

Runoff processes in soil formation in zero-order basins might show intermediate characteristics between hillslopes and channels (Fig. 1). For example, the generation of overland flow is rare in both zero-order basins and hillslopes [e.g., Uchida et al., 2003], and soil consists not only of weathered materials at the original point but also of deposited materials delivered from upper hillslopes [e.g., Reneau and Dietrich, 1991].

Therefore, we connected two concepts of sediment yields, shallow landslide on hillslopes and initiation of debris flow in channels, to assess the occurrence of shallow landslide and debris flow

simultaneously [Hata et al., 2014]. We tested the applicability of the method proposed by [Hata et al., 2014] (the advanced H-SLIDER method) to predict sediment yield in headwaters.

2. METHODS AND NUMERICAL FEATURES Based on infinite slope stability analysis, sediment movement processes can be classified into three types, depending on the magnitude of soil thickness and slope gradient, as shown in Fig. 2.

One type is shallow landslide occurrence caused by subsurface flow (Type 1 of Fig. 2). The other two types are sediment yield caused by not only subsurface flow but also saturation overland flow (Types 2 and 3, respectively, of Fig. 2). In a relatively steep hillslope, once the flow depth reaches critical value, the soil layer moves simultaneously, as in Type 1 and Type 2 of Fig. 2, but in a relatively gentle slope. Even though flow depth reaches a critical value, the soil layer moves gradually from the surface to deeper zones (Type 3 of Fig. 2).

In the study, to assess the effects of subsurface flow on shallow landslide occurrence (Type 1 of Fig. 2), we applied a method proposed by [Uchida et al., 2011]. The method combines infinite slope stability analysis and a subsurface flow model, assuming rainfall intensity is steady. Thus, we calculated a safety factor for Type 1 (Fs ₁ ) using the following equation:

Fig. 2 Classification of transition to unstable state of soil layer Fig. 1 Characteristics and sediment yield prediction method

of hillslopes, channels, and hollows

(1)

where c is the cohesion, γ is the unit weight, the subscripts h and w indicate hillslope and water, respectively, h is the soil depth, θ is the slope angle, r is the rainfall intensity, A is the drainage area for unit counter length, K s is the saturated hydraulic conductivity, and φ is the friction angle of soil.

In a case where the rainfall intensity (r) is zero and the Fs 1 is less than 1.0, we interpreted the grid cells as “unconditionally unstable cells,” similar to previous studies [e.g., Montgomery and Dietrich, 1994]. Hereafter, destabilization of this type is referred to as Type 0.

To assess the effects of saturation overland flow on hillslope stability (Types 2 and 3 in Fig. 2), we used the equation proposed by [Takahashi, 1991].

However, while [Takahashi, 1991] regarded the effects of cohesion as negligible, we added the effects of cohesion to his equation. Then, if a grid cell satisfies the equation

(2) even if the soil mantle is fully saturated, Fs ₁ becomes more than 1.0. In this case, we calculated the safety factor for Type 2 (Fs ₂ ) and Type 3 (Fs ₃ ) as follows:

(3)

(4)

where C * is the sediment concentration of soil, γ p is the unit weight of a soil particle, h w is the depth of saturation overland flow, and d is the grain diameter of the surface of a hillslope.

Moreover, we classified Types 2 and 3 using the following equation:

(5) where γ _s is the unit weight of saturated hillslope. If the slope angle satisfied Equation (5), we used Equation (3) described as Type 2. If not, we used Equation (4).

3. APPLICATION TO A REAL HEADWATER 3.1 Study area and rainstorm

Our study area is a headwater in the Tsurugi River, in Hofu City, Japan. The region is humid and temperate. The mean annual precipitation in this region is approximately 1,600 mm, and the mean temperature is approximately 16 °C. The headwater is mainly underlain by granite, and the drainage area is 6.4 ha (Fig. 3). The site is deeply incised and dominated by hillslopes. Slope gradients range from 0° to 56°, and the average is 32°. Slope lengths range from 10 to 270 m.

On July 21, 2009, many shallow landslides occurred in this region, and these landslides triggered many debris flows. Eight of these shallow landslides occurred in the study area. The landslide scarred, the range from 3 to 30 m in width. The total rainfall amount and maximum rainfall intensity of the triggering event was 275 mm and 63 mm/h, respectively, measured at the Hofu rainfall observing station where the distance to the study site was 5 km. As shown in Fig. 3, most of the sediment movement occurred in zero-order basins.

3.2 Parameters setting

The parameters set in the calculation model are shown in Table 1.

3.2.1 Topographical parameters

We made the 5-m digital elevation model (DEM) of the ground surface using the LiDAR data.

      Basin divide        Stream 

      Observed erosion area         Zero‑order basin in above area        Penetration test point        Soil test point 

Fig. 3 Study area

θ θ γ

θ φ θ γ

γ

sin cos

tan tan cos ²

1 h

K h rA

c Fs

h

s

h w  

 



 −

+

=

( ) ^{( )}

[ ]

( )

( ^{γ γ γ} ) ^{γ γ ( ) θ φ}

θ θ γ

γ γ

tan tan / 1

sin cos /

2 1

h h C

C

h h h C

Fs c

w w w p

w p

w w w p

+ +

− + −

+ +

= −

∗

∗

∗

      

( ) ( )

[ ]

( )

( ^{γ γ} ) ^{γ γ γ} ( ) ^{θ φ}

θ θ γ

γ γ

tan tan / 1

sin cos /

3 1

d h C

C

d d h C

Fs c

w w w p

w p w w w p

+ +

− + −

+ +

= −

∗

∗

∗

       c h

s w

− >

+ cos (sin cos tan ) tan

cos ² θ φ γ θ θ θ φ

γ

γ φ γ

θ γ ^tan

tan

s w s −

≥

※Contour interval : 1m 

Then, we used a kriging interpolation scheme to make a 5-m DEM of the bedrock surface from a 5-m ground surface DEM and the soil depth data. Thus, we calculated local slope angle and upslope drainage area data sets for the bedrock surface topography. Also, we set the soil depth where sediment movement occurred. We added the difference of the before and after of the disaster at each grid point to the value.

We calculated the local slope angle and upslope area using the D-infinity Flow Direction method (D-infinity) [Tarboton, 1997]. D-infinity is an algorithm based on proportioning flow between two downslope pixels and represents flow direction as a single angle taken as the steepest downward slope on the eight triangular facets centered at each grid point. Upslope area is then calculated by proportioning flow between two downslope pixels according to how close this flow direction is to the direct angle of the downslope pixel.

We used the soil depth data measured by Public Works Research Institute, using a cone penetrometer (Portable dynamic cone penetration test) with a cone diameter of 25 mm, a weight of 5 kg, and a fall distance of 50 cm. N _d represents the number of blows required for 10 cm of penetration. The penetration tests were conducted at intervals of about 15 m in 2010-2011. In all, the penetration tests were conducted at 151 points, 140 of which were in the study area, 11 of which were outside of it (Fig. 3). In this study, we assumed the depth of the layer with N _d < 10 as the soil depth, according to the research of [Osanai et al., 2005]. In this study area, the soil depth wasn’t observed in the northwest area, so we excluded the northwest area from the examination.

3.2.2 Soil parameters

We set the friction angle and the unit weights of saturated and unsaturated (wet) soil by measuring undisturbed soil samples. The soil tests were also conducted by Public Works Research Institute in 2010. The sampling was conducted at 7 points (shown as ‘ ’ in Fig. 3) of 2 depths (50 and 100 cm), using a sampler with a diameter of 75 mm, and a height of 150 mm. The cohesion was back calculated using soil depth data and the DEM using the method proposed by [Uchida et al., 2009 ]. We assumed that if the soil mantle was unsaturated, the factor of safety of most grid cells would be larger than 1.0. Therefore, we used Equation (1) to calculate the minimum cohesion, including root strength, required for all grid cells to remain stable under unsaturated condition. Effective saturated hydraulic conductivity in a hillslope is often different from that measured in a small undisturbed

soil sample [e.g., Uchida et al., 2003]. Therefore, we used the value from hydrologic observation in the Hiroshima western mountains, where geologic features are the same as the study area.

3.3 Case setting

We summarized the calculation cases in Table 2. In cases A-1, A-2, and B-1 through B-3, we ignored the effects of saturation overland flow, similar to most shallow landslide prediction models [e.g., Okimura et al., 1985; Montgomery and Dietrich, 1994]. That is, we calculated only Fs 1

using Equation (1). In these cases, if the condition described as Equation (2) were satisfied, we interpreted such grid cells as “unconditionally stable” with regard to shallow landsliding.

Table 1 Parameters setting

Items Setting method Setting value Rainfall

condition

Rainfall intensity

Equivalent to a 6 h average rainfall intensity before expected collapse time

40 mm/h Topographic

conditions

Soil depth Value of before the disaster, complemented by kriging interpolation scheme using soil depth of after the disaster, and added the difference of the elevation between the before and after disaster LIDAR data

−

Local slope

angle，

Upslope drainage area

Calculated by D-infinity Flow Direction method using bedrock elevation estimated by LIDAR data and soil depth

−

Soil conditions

Cohesion Value back calculated by soil depth and terrain, or regarded 0 as same as the debris flow prediction method

6.5 kN/m ² 0.0 kN/m ²

Friction angle

Soil test

35°

Saturated unit weight

of soil

Soil test

18.1 kN/m ³

Wet unit

weight of soil

Soil test

15.2 kN/m ³ Grain

diameter

Grain diameter to pass the mass of 90% from the smallest diameter

10 mm Saturated

hydraulic conductivit y

Value of hydrologic observation in Hiroshima western mountain where geologic features are the same as the study area

0.05 cm/s

Sediment concentrati

on of soil

Value back calculated by the saturated hydraulic conductivity

0.53 The others Roughness

coefficient

Value used as a standard in

many hillslopes 0.1 m ^-1/3 ・s

Unit weight

of water

− 9.8 kN/m ³

Mesh size − 5 m

In cases C-1 through C-3, we considered the effects of saturation overland flow to be similar to the equation proposed by [Takahashi, 1991]. In these cases, we judged stability of grid cells based not only on Fs 1 but also on Fs 2 and Fs 3 .

In cases A-1 and A-2, we set the single values as the cohesion of soil layers. In other cases, we considered the spatial variability of cohesion of soil layers. Here, we simply assumed that the cohesion of in-situ weathered materials was different from the cohesion of deposited materials. We assumed that the cohesion of deposited materials might be negligible, similar to the concept of the [Takahashi, 1991] equation.

Previous studies indicated that by use of a characteristic break in the scaling relation of local slope and upslope areas, the catchment can be classified into some parts such as hillslopes, valley heads, colluvial and alluvial [e.g., Montgomery, 2001]. Montgomery [2001] indicated the transition from the hillslope domain, where gravitational processes dominate, to the channel, where fluvial processes dominate, can be determined from a DEM. According to these studies, we simply classified the catchment into two parts, hillslopes, where in-situ weathered materials dominate, and channels, where deposited materials dominate, by use of the characteristic break in the scaling relation of local slope and upslope areas of the study area (Fig. 4). In Fig.4, we plotted the average value of 20 grid cells. We selected 20 grid cells according to the rules below.

1. The values of upslope areas are similar.

2. The values of local slopes are similar.

Because we found three breaks in the scaling relation of local slope and upslope areas, we used three upslope areas as boundaries between hillslope and channel (Table 2).

Table 2 Case of calculation Case

*

Effects of saturation

overland flow

Cohesion (kN/m ² )

Boundary upslope

area of channel

(m ² ) Hillslopes Channels

A-1 Ignored 6.5 6.5 - A-2 Ignored 0.0 0.0 - B-1 Ignored 6.5 0.0 100 B-2 Ignored 6.5 0.0 900 B-3 Ignored 6.5 0.0 30,000 C-1 Included 6.5 0.0 100 C-2 Included 6.5 0.0 900 C-3 Included 6.5 0.0 30,000

*Case A, B: H-SLIDER method Case C: Advanced H-SLIDER method

3.4 Evaluation of calculated result

To evaluate results of the calculation, we determined three values, hit ratio, cover ratio, and threat score, as follows (Table 3).

Table 3 Calculation method of hit ratio, cover ratio, and threat score

Observed sediment movement

Move Not move

Predicted sediment movement

Move a b

Not move c d

Hit ratio = a/(a+b) Cover ratio = a/(a+c) Threat score = a/(a+b+c)

4. RESULTS

We show the calculation results in Fig. 5.

In case A-1, setting the cohesion (c = 6.5 kN/m ² ) uniformly all over the basin, sediment movements where slope gradients were relatively steep were predicted as Type 1, but the grid cells where slope gradient are relatively gentle were predicted to be stable. However, in some of these grid cells, sediments were yielded in the storm of July 2009. In this case, although the hit ratio was relatively high (0.65), the cover ratio was small (0.25). In case A-2, the cohesion was set as zero for all of the grid cells in the watershed. Destabilization Types 0 and 1 were calculated for most of the grid cells, although they were the areas where sediments movement did not occur in the storm of July 2009.

In this case, although the cover ratio was relatively high (0.62), the hit ratio was very small (0.12).

The B cases, setting the classified cohesion, show intermediate results between case A-1 and A-2. The hit ratios of the B cases were smaller than

0

100 m ²   900 m ²   30,000 m ²  

Loca l s lop e (° ) 

0       100    1,000     10,000   100,000  10

20 30 40 50 60

Fig. 4 Scaling relation of local slope and upslope area

   Upslope area (m ² )