Numerical Simulation of a Debris Flow on the Basis of a Two-Dimensional Continuum Body Model

Abstract

A two-dimensional debris and mud flow model that considers both laminar and turbulence flow was developed. Subsequently, the model was applied to a debris flow that occurred in Asaminami, Hiroshima, Japan in August 2014. The applicability of the model and the debris flow characteristics are discussed. The calculated horizontal distribution of sediment deposited in the Asaminami residential area was in good agreement with the horizontal distribution of the deposited large rocks and driftwood. This result indicates that the fine material in the downstream area was transported by water flow resulting from heavy rain that occurred after the debris flow. The scale of the initial debris flow was small; however, it increased with time, because eroded bed material and water were entrained to it. Therefore, it is important to reproduce the development process of debris flows to predict the amount of sediment produced, the deepest flow depth, the maximum flow velocity, and the inundation area. The averaged velocity of the simulated debris flow was about 9 m/s, and the velocity at the entrance to the residential area was about 8 m/s. This kind of information can be used to design sediment deposition dams. The travel time of the simulated debris flow from the upstream end of the main channel to the entrance of the residential area was 96 s. This kind of information can be used for making evacuation plans. Valley bed steps can suppress the deepest flow depth which is very important for the design of check dams; therefore, the high-resolution elevation data and fine numerical grids that reproduce step shapes are required to accurately calculate the deepest flow depth and maximum flow velocity.

1. Introduction

Sediment disasters caused by debris and mud flows that originate from surface landslides generated by high intensity rainfall events happen frequently, for example, in Japan, Miyagi in 2019, Hiroshima in 2014 and 2018, Fukuoka in 2017, Hokkaido in 2016, Kagoshima in 2015, and Izu Oshima in 2013. According to an IPCC (Intergovernmental Panel on Climate Change) report [1], because of global warming that will result in extreme weather characteristics, it is expected that the frequency of high intensity rainfall events will increase even if average precipitation does not change. When rainfall intensity is high, even if the rainfall duration is short and the total precipitation is small, the debris and mud flows that originate from surface landslides are easily generated. Therefore, the frequency of debris and mud flows caused by weather characteristics resulting from global warming are expected to increase, and effective countermeasures based on the results of numerical analysis of the debris and mud flows are required. The information required from these numerical analyses includes velocity, depth, travel time to residential regions, inundation regions, and sediment production. Travel time to residential regions and inundation regions can be used for the resident evacuation action. Velocity, depth, and sediment production can be used to design hard countermeasures (e.g., a check dam, channel work).

Debris flow and snow avalanches are multiphase mass movement on mountain slope and similar governing equations are used to reproduce them by numerical analysis. Therefore, the same ideas can be applied to both phenomena under some conditions and phenomena.

Many previous studies (e.g., [2,3,4]) have suggested that the entrained material volume from valley beds is much larger than landslide volume and many field surveys support this idea (e.g., [3,5,6]). These findings suggest that entrained material volume should be accounted for when simulating debris flow and snow avalanche behavior. Despite the importance of the material entrainment from the valley bed estimating debris flow and snow avalanche volume and the temporal change of their transport rate, it tends to be neglected. In recent years, several models for simulating debris flow and snow avalanche have been developed (e.g., [7,8,9,10,11,12,13,14]). However, some of these models assume a constant volume of debris flow and snow avalanche, and only simulate their deformation process along valleys. When the unstable sediment and snow cover on the base rock is thin and the land slope is mild, the entrainment can be neglected (e.g., [15,16,17]). However, when the unstable sediment and snow cover on the base rock is thick and the land slope is steep, the entrainment should be considered.

Some studies (e.g., [4,18,19]) have used the idea that the thickness of the bed surface layer where the bed shear stress is larger than the bed resistance is erosion depth, and many bed erosion rate equations are proposed (e.g., [10,13,20,21,22,23]). These models can calculate the suitable bed erosion rate using suitable parameters. However, the calculated results are well affected from the many parameters used in the analysis. Additionally, some of them cannot reproduce equilibrium bed conditions.

When the cohesive characteristics of the bed layer on bed rock can be neglected and the equilibrium bed condition that bed erosion rate is the same as sediment and snow deposition rate exists, erosion rate can be estimated using only flow velocity, solid material concentration in the bed layer, and relationship between bed slope and equilibrium bed slope [24]. Incidentally, when the solid material size is large, the debris flow and the snow avalanche must be laminar flow, because turbulence is not well developed in the flow depth scale. In contrast, when the solid material size is small (e.g., fine granite, volcanic ash, dry snow, etc.), a turbulence flow layer is formed on the laminar flow layer near the bed, which becomes mud flow and mixed-motion avalanche. The formation of both laminar and turbulence flows reduces the equilibrium bed slope and increases the bed erosion rate. Therefore, considering the formation of both laminar and turbulence flow layers is key to estimating the bed erosion rate.

In this study, a two-dimensional debris and mud flow model that considers both laminar and turbulence flows was developed. The model was subsequently applied to debris flows that occurred in Asaminami, Hiroshima, Japan in 2014 and the applicability of the model and the flow characteristics of the debris flows are discussed.

2. The Basic Equations of the Numerical Debris Flow Simulation Model on the Basis of a Two-Dimensional Continuum Body

Debris and mud flows comprise a mixture of flowing water and sediment, which is treated as a one phase continuum fluid body. The mass conservation equation of the water and sediment mixture is as follows [24]:

$$\frac{\partial h}{\partial t}+\frac{\partial hu}{\partial x}+\frac{\partial hv}{\partial y}=\frac{E}{c_{*}}$$

(1)

where, t is the time, h is the flow depth, while u and v are depth averaged velocities in the x and y directions, respectively. The term on the right hand side indicates the sink and source of the mass and expresses the development and decrescence of a debris flow by exchange of the mixture of water and sediment between bed surface and debris flow, where c_* is the concentration of sediment in the static deposition layer (bed layer), and E is the erosion rate of the bed. When the cohesive characteristics of the bed layer on the bed rock can be neglected and the equilibrium bed condition, wherein the bed erosion rate is the same as sediment deposition rate, exists, the bed erosion rate can be estimated using the following equations [24]:

$$\frac{E}{\sqrt{u^2+v^2}}=c_{*}\tan \left(\mathsf{\theta }-{\mathsf{\theta }}_e\right)$$

(2)

where θ, the bed slope in the flow direction, is calculated using the following equation:

$$\sin \mathsf{\theta }=\frac{u\sin {\mathsf{\theta }}_x+v\sin {\mathsf{\theta }}_y}{\sqrt{u^2+v^2}}$$

(3)

where θ_x is the bed slope in the x direction and θ_y is the bed slope in the y direction, and θ_e is the equilibrium bed slope. Figure 1 shows the relationship between bed deformation and bed slope. When the bed slope is steeper than the equilibrium bed slope, the bed surface is eroded, and the eroded materials are introduced into the debris flow. When the bed slope is equal to the equilibrium bed slope, the bed surface is not deformed, and the mass of the debris flow is unchanged. When the bed slope is milder than the equilibrium bed slope, the materials in the debris flow are deposited on the bed and the mass of the debris flow is reduced.

As shown in Figure 2, when the sediment size is large, the laminar flow layer is dominant and the flow becomes a debris flow. When the sediment size is small, the turbulence flow layer is formed on the laminar flow layer near the bed, which becomes the mud flow. Considering the formation of both the laminar and turbulence flows, and referring to the averaged depth and the sediment concentration $\overline{c}$, the equilibrium bed slope θ_e in the flow direction is obtained as follows:

$$\tan {\mathsf{\theta }}_e=\frac{\left(\mathsf{\sigma }/\mathsf{\rho }-1\right)\overline{c}}{\left(\mathsf{\sigma }/\mathsf{\rho }-1\right)\overline{c}+1}\frac{h_s}{h}\tan {\mathsf{\varphi }}_s$$

(4)

where φ_s is the angle of repose. The laminar flow layer thickness h_s is estimated by the following equation [25]:

$$\frac{h_s}{d_m}=\frac{1}{\overline{c}\cos \mathsf{\theta }\left(\tan {\mathsf{\varphi }}_s-\tan \mathsf{\theta }\right)}{\mathsf{\tau }}_{*m}$$

(5)

where d_m is the mean diameter of the debris flow, and τ_*m is the non-dimensional shear stress of the mean diameter. When the laminar flow layer is thin, the equilibrium bed slope is small. Consequently, debris flow can be transported to the milder slope area.

The mass conservation of sediments in the sediment and water mixture is

$$\frac{\partial \overline{c}h}{\partial t}+\frac{\partial \overline{c}hu}{\partial x}+\frac{\partial \overline{c}hv}{\partial y}=E$$

(6)

The momentum conservation equations are as follows:

$$\frac{\partial hu}{\partial t}+\frac{\partial huu}{\partial x}+\frac{\partial huv}{\partial y}=-gh\frac{\partial z_b}{\partial x}-\frac{1}{{\mathsf{\rho }}_m}\frac{\partial P}{\partial x}-\frac{{\mathsf{\tau }}_{bx}^{}}{{\mathsf{\rho }}_m}$$

(7)

$$\frac{\partial hv}{\partial t}+\frac{\partial huv}{\partial x}+\frac{\partial hvv}{\partial y}=-gh\frac{\partial z_b}{\partial y}-\frac{1}{{\mathsf{\rho }}_m}\frac{\partial P}{\partial y}-\frac{{\mathsf{\tau }}_{by}^{}}{{\mathsf{\rho }}_m}$$

(8)

where g is the gravity acceleration and z_b is the bed elevation. The pressure P is assumed to be the static pressure. The density of debris flow ρ_m is as follows:

$${\mathsf{\rho }}_m=\left(\mathsf{\sigma }-\mathsf{\rho }\right)\overline{c}+\mathsf{\rho }$$

(9)

where ρ is the water density and σ is the sediment density, while τ_bx and τ_by are the shear stress in the x and y directions, respectively. The energy dissipation of both solid and fluid motion bases are considered to estimate the shear stress. When the turbulence flow region is dominant, the shear stress is as follows:

$${\mathsf{\tau }}_{bx}^{}=\left\{{\mathsf{\tau }}_y+\mathsf{\rho }{f}_b\left(u^2+v^2\right)\right\}\frac{u_{}}{\sqrt{u_{}^2+v_{}^2}}$$

(10)

$${\mathsf{\tau }}_{by}^{}=\left\{{\mathsf{\tau }}_y+\mathsf{\rho }{f}_b\left(u^2+v^2\right)\right\}\frac{v_{}}{\sqrt{u_{}^2+v_{}^2}}$$

(11)

where τ_y is the yield stress, which uses the following relation for non-cohesive materials:

$${\mathsf{\tau }}_y={\left(\frac{\overline{c}}{c_{*}}\right)}^{\frac{1}{5}}\left(\mathsf{\sigma }-\mathsf{\rho }\right)\overline{c}g{h}_s\cos \mathsf{\theta }\tan {\mathsf{\varphi }}_s$$

(12)

where f_b is the resistance coefficient. The following relation is used:

$$f_b=\frac{4}{25}\left\{k_f\frac{{\left(1-\overline{c}\right)}^{\frac{5}{3}}}{{\overline{c}}^{\frac{2}{3}}}+k_d\frac{\mathsf{\sigma }}{\mathsf{\rho }}\left(1-e^2\right){\overline{c}}^{\frac{1}{3}}\right\}{\left(\frac{h_s}{d_m}\right)}^{-2}+9{\alpha }^2\left(1-\frac{h_s}{h}\right)$$

(13)

The first and the second terms on the right-hand side are from the laminar flow layer, and the third term is from the turbulence flow layer. The first term is from the turbulence of the liquid among the solid particles. Small scale turbulence can be formed in a laminar flow layer and must be considered when estimating energy dissipation. The second term is from the inelastic collisions of sediment particles, where k_f=0.16, k_d=0.0828, and e are the inelastic coefficients of the particles herein, and d_m is the mean particle size of the sediment in the debris flow. α is as follows:

$$\alpha =\frac{\mathsf{\kappa }}{6}$$

(14)

The bed elevation equation is as follows [24]:

$$\frac{\partial z_b}{\partial t}=-\frac{E}{c_{*}}$$

(15)

The developed numerical analysis model is installed into the free public domain river analysis interface iRIC [26], which is developed by a research group including the first author. The name of the numerical analysis mode is Morpho2DH and can be used for the two-dimensional water flow and bed deformation analysis of bed material load [27], debris flow, and mud flow.

3. The Sediment Disaster at Asaminami, Hiroshima, Japan in 2014

As shown in Figure 3, Asaminami is located in the northern Hiroshima City. The average land slope is 20 degrees. Figure 4 presents the geological map around Asaminami. The red circle in Figure 3 and Figure 4 indicates the research area. The geology around Asaminami is mapped as granite and accretionary prism. However, most parts of the land surface are covered by granite, based on our field observations. The residential area is located at the foot of a slope, and 41 casualties resulted from debris flows.

Figure 5 shows the horizontal distribution of rainfall intensity at 03:00 on August 20, 2014 obtained by Japanese Ministry of Land, Infrastructure, Transport and Tourism. As shown in Figure 5, a high rainfall intensity area is formed around the black circle in the direction northeast from southwest, which is called a line-shaped precipitation system, and Asaminami received heavy rain from August 19 to 20, 2014. Figure 6 shows the temporal change of the rainfall intensity at Asaminami. The rainfall intensity was measured by X-band radars operated by Japanese Ministry of Land, Infrastructure, Transport and Tourism. Rainfall intensity is measured every minute and the measured rainfall intensity is converted to one-hour rainfall intensity. The rainfall intensity intermittently exceeded 100 mm for about two hours between 01:35 and 03:30, and the three debris flows were transported to the residential area at about 03:00 on August 20. Therefore, the heavy rain continued for 30 min after the debris flows occurred.

Figure 7 shows the slope eroded by the surface landslides in the (a) main valley and (b) tributary, and (c) and (d) the beds eroded by the transportation of debris and water flows. At least three debris flows occurred. They started from (a) the upstream end of the main valley, (b) the upstream end of the tributary, and (d) the midstream area of the main channel. Photo (c) in Figure 7 shows a high step in the valley bed, and many other steps were observed in the valley. The slope surface is composed of well weathered granite, which can be easily broken by hand, as shown in Figure 8. The permeability in the granite is very high because many joints are formed as shown in Figure 9. Hence, the slope failures do not occur as a result of a high level of accumulated precipitation with low rainfall intensity; the failures of the surface slope occurred because of high rainfall intensity.

Figure 10 shows the confluence of the main valley and the tributary, and photo (b) shows the right bank (on the western side of the valley) at the confluence having been eroded by the debris flow from the tributary, and exposing the base rock. Contrastingly, the left bank (the eastern side of the valley) at the confluence was eroded by a debris flow from upstream of the main valley, but the land surface is covered by fine grained granite after the disaster (Figure 10c). This indicates that the last debris flow was from the tributary, and fine-grained granite from the tributary was deposited on the left, inner bank (the eastern side).

Figure 11 shows the damaged residential area. Three-storied reinforced concrete buildings, shown in (b), and wooden houses are located at the downstream end of the valley. While many of the wooded houses located along the road collapsed, all the reinforced concrete buildings remained intact despite being in the path of the debris flows. Even though the glass windows of the reinforced concrete buildings were broken, and sediment flowed into the first-floor rooms, damage to the reinforced concrete buildings was minor. Figure 12 shows the no-inundation areas of the debris flow. Because of the existence of the reinforced concrete buildings, the debris flows were concentrated along north–south roads and did not spread in the east-west direction.

4. Calculation Conditions

Figure 13 shows the calculation domain and polygons of houses and buildings. The land elevation data (DEM), indicated by the colored contours, was measured by the Japanese Ministry of Land, Infrastructure, Transport, and Tourism in 2008. Houses and buildings left after the disaster were considered. No heavy rain events associated with landslides or debris flows occurred between 2008 and 2014. The horizontal resolution of the elevation data is 1 m. To express the shapes of houses, roads, and steps in the valley bed, each cell size of the numerical analysis is 2 m × 2 m.

The surface landslides, which are the origins of the three debris flows, considered in these numerical analyses are indicated by the three brown circles in Figure 12, showing the location of the landslides. They are located at (a) the upstream end of the main valley, (b) the upstream end of Z tributary and (d) the midstream area of the main channel (see Figure 7). As explained in Section 3, the debris flow from the tributary was the last one. Therefore, the debris flows from the upstream end of the main valley and the midstream area of the main channel, respectively, were started first in the analysis. After 100 s, the debris flow from the upstream end of the tributary was started. The sizes of the landslides were based on field survey data. The sediment volume of each landslide are 70 m³ at the upstream end of the main valley, 90 m³ at the upstream end of tributary, and 60 m³ at the midstream area of the main channel, respectively. Friction slope is 34°. Sediment specific weight is 2.65. Mean sediment diameter is 1 cm, which is obtained by field observation. 0.85 is used for the inelastic coefficients of the particles e.

Two analysis cases were performed. For Case 1, the houses and buildings left after the disaster as shown in Figure 13b were considered. For Case 2, the houses and buildings were neglected. Houses and buildings were treated as non-permeable non-destructive structures and behave as obstacles.

5. Results and Discussion

Figure 14 shows the horizontal distribution of the land deformation thickness in the residential area obtained by the numerical analysis of Case 1. The brown line indicates the region where large rocks and driftwoods were deposited. As shown in the photos in Figure 14, sediment sorting was observed in the residential area. Large rocks and driftwoods were observed near the exit of the valley and many wooden houses were seriously damaged. However, only fine-grained granite was deposited downstream of the residential area. As shown in Figure 14, the land deformation area is in good agreement with the deposition area of the large rocks and driftwoods. This model was developed to reproduce the development, transportation, and deposition processes of the debris flow; the erosion, transportation, and deposition of fine-grained sediment after the deposition of sediment transported as a debris flow was not considered. As shown in Figure 6, the rainfall intensity was high after the sediment transported by the debris flows was deposited. Hence, the deposited sediment in the residential region was eroded by the heavy rain and the only fine-grained material was transported downstream. Therefore, the two-layer model based on Equations (4) and (5) did reproduce the debris flow land deformation region.

Figure 15 shows the spatiotemporal change in debris flow depth in Case 1. Initially, the masses of the debris flows were small, but they increased with time, because eroded bed material and water were added. The debris flow development process was considered by the terms in the right-hand side of Equation (1). Therefore, reproducing this development process must be important for predicting the amount of sediment production and the inundation area.

The average velocity of the calculated debris flow in the valley was 9 m/s, and the calculated velocity upstream of the residential area was 8 m/s. The calculated travel time of the debris flow from the upstream end of the main channel to the residential area was 96 s. The debris flows happened around 03:00 on August 20, and the rainfall intensity exceeded 100 mm/h. The residents were unable to move to an evacuation shelter after the occurrence of the landslide because the debris flow travel times were very short, they occurred very early in the morning, and the rainfall intensity was high.

Figure 16 shows the spatiotemporal change in debris flow depth in the residential area including houses and buildings (Case 1). Figure 17 shows the spatiotemporal change in debris flow depth in the residential area without the inclusion of houses and buildings (Case 2). In Case 1, debris flows concentrate along the north-south road and cannot spread in the east-west direction because of houses and buildings. In contrast, in Case 2, the debris flow spreads widely and forms a fan. However, the debris flow in Case 1 does not flow further to the south, even though the flow is concentrated on the north-south road because sediment was deposited upstream of the houses and buildings, and the debris flow discharge decreased along the north-south road. Comparing the inundation areas in Case 1 with those in Case 2, there are safe areas located downstream of the houses and buildings that were not destroyed. These areas, especially those downstream of the reinforced concrete structures, can be used for effective evacuation areas because they are located in the residential area and the residents could easily evacuate very quickly. When the travel time of debris flow to the residential area is too short and inhabitants cannot move to non-inundation area, evacuation to second or third floor can increase survival rate. As shown in Figure 16, the flow depth in the downstream of inundation area is shallower than 2.5 m and the debris flow cannot reach second floor. As shown in Figure 11e, the first floors of houses are destroyed, however, the second floors have no damage.

As mentioned above, numerical simulation of debris flows can be a useful tool for the planning of evacuation sites.

Numerical simulation of debris flows can be a also useful tool in taking hard countermeasures against sediment disasters. Figure 18 shows the longitudinal distribution of the temporally accumulated sediment transport rate, which increased up to 880 m away from the upstream end of the calculation region. This reduced after 880 m due to deposition. The dense residential area was located 1008 m from the upstream end of the calculation region, and the temporally accumulated sediment transport rate was 33,504 m³, which is in good agreement with the measured value of 33,000 m³ from the DEM data obtained by the Japanese Ministry of Land, Infrastructure, Transport, and Tourism. The longitudinal distribution of the temporally accumulated sediment transport rate is very useful for the design and locating of check dams. If accessible, a check dam could be constructed in the upstream region at 200 m; a small check dam or a steel net as shown in Figure 19a able to deposit 6000 m³ would be sufficient. Locating a check dam at 600 m that would deposit sediment from both the main valley and the tributary able to deposit 25,000 m³ would also be sufficient.

To design each check dam, the deepest flow depth and the maximum flow velocity at each location would be the most useful information. These values are very affected by the geometry of the land upstream of the check dams. Figure 20 shows the longitudinal distribution of the deepest flow depth of each cross-section, and as Figure 18 shows, the temporally accumulated sediment transport rate increased up to 880 m from the upstream end of the calculation area. However, the flow depth did not substantially increase from 100 m onwards, only changing from 2 to 7 m (Figure 20), and started to decrease from 900 m because of deposition. The steps shown in Figure 7c suppress the development of flow depth. Figure 21 shows the temporal change in the flow depth upstream of the main valley. There are several steps up to 10 m high in this reach, and, as shown in Figure 2, the flow depth is shallower at the steps and deeper upstream of the steps and the flow velocity is faster at the steps and slower upstream of steps. As a result, the debris flow is extended along the flow direction around the step, and the flow depth cannot be much deeper because of the extensions. In the analysis, fine numerical grids were used to reproduce the shape of the steps in valley rock bed. Therefore, high resolution elevation data and fine numerical grids are needed to reproduce the shape of steps when calculating deepest flow depth and maximum flow velocity.

As mentioned previously, the numerical simulation of a debris flow can be an essential tool for estimating the deepest flow depth and maximum flow velocity when designing a check dam. The numerical simulation of a debris flow can also be an essential tool when designing a channel work in residential area (Figure 19b). If there is a large river downstream of the residential area, channel work through the residential area can drain the debris flow to the large river. Width, depth, slope, and horizontal shape of a channel work, which can suppress debris flow inundation, can be designed with the use of the numerical analysis.

The numerical simulation of a debris flow can also be used for future large debris flow disasters. It is very difficult to predict the timing and location of landslide. However, the largest debris flow can be predicted, when the entrainment of the material from the valley bed is much larger than the landslide volume. When landslides occur at the upstream end of all valleys in the basin and the all debris flows from each valley merge, the largest debris flow is formed. Numerical simulation of the largest debris flow can be useful information to make soft and hard countermeasures against sediment disasters.

The horizontal distribution of house destruction is important information to reproduce the spatial distribution of debris flow inundation. When the strength of each house and building is surveyed, it can be judged if those houses and buildings are destroyed or non-destructed by the use of the information on flow depth and flow velocity in the residential area.

6. Conclusions

A two-dimensional debris and mud flow model considering both laminar and turbulence flows was developed. The model was then applied to the debris flows that occurred in Asaminami, Hiroshima, Japan in 2014. The results are summarized as follows:

(1)

In the two-dimensional debris and mud flow models under a continuum body treatment, the formation of both laminar and turbulent flows was considered to estimate the equilibrium bed slope and energy dissipation.

(2)

The calculated horizontal distribution of the sediment deposit area in the residential region was in good agreement with the horizontal distribution of the deposition of large rocks and driftwood. This indicates that the fine-grained material in the downstream area was transported by water flow resulting from the heavy rain that continued after the debris flows.

(3)

The debris flow was initially small. However, it increased in size over time, because eroded bed material and water were incorporated. Hence, reproducing the development process of a debris flow is important for predicting the amount of sediment production and the inundation area.

(4)

The averaged velocity of the simulated debris flow was about 9 m/s, and the simulated velocity at the entrance to the residential area was about 8 m/s. This kind of information can be used to design sediment deposition dams.

(5)

The travel time of the simulated debris flow from the upstream end of the main channel to the entrance of the residential area was 96 s. This kind of information can be used for evacuation planning.

(6)

Steps in the valley bed can suppress the deepest flow depth, which is important to know when designing a check dam. Therefore, high resolution elevation data and fine numerical grids for the reproduction of step shapes are required to accurately calculate the deepest flow depth and the maximum flow velocity.