Application of simulation technique on debris flow hazard zone delineation : a case study in the Daniao tribe , Eastern Taiwan

Typhoon Morakot struck Taiwan in August 2009 and induced considerable disasters, including large-scale landslides and debris flows. One of these debris flows was experienced by the Daniao tribe in Taitung, Eastern Taiwan. The volume was in excess of 500 000 m 3, which was substantially larger than the original design mitigation capacity. This study considered large-scale debris flow simulations in various volumes at the same area by using the DEBRIS-2D numerical program. The program uses the generalized Julien and Lan (1991) rheological model to simulate debris flows. In this paper, the sensitivity factor considered on the debris flow spreading is the amount of the debris flow initial volume. These simulated results in various amounts of debris flow initial volume demonstrated that maximal depths of debris flows were almost deposited in the same area, and also revealed that a 20 % variation in estimating the amount of total volume at this particular site results in a 2.75 % variation on the final front position. Because of the limited watershed terrain, the hazard zones of debris flows were not expanded. Therefore, the amount of the debris flow initial volume was not sensitive.


Introduction
Debris flows are frequent phenomena in Taiwan.To minimize the possible hazards caused by debris flows, the common countermeasures include the construction of dams, the limitation of land use, and habitant evacuation.Two of the Correspondence to: Y. C. Hsu (d92521021@ntu.edu.tw)common uncertainties during the planning of countermeasures are the hazard zone area and the path of debris flows.Several empirical formulas may be used to obtain part of the information required in the designing processes.However, empirical formulas may be inaccurate in complicated geographic regions, even for the order of magnitude.A superior method to obtain the required information is to use numerical simulation.
Non-Newtonian fluid models were often used in previous studies of debris flow.Early analytical or numerical studies of Bingham-like fluids were limited mainly to one-or two-dimensional spreading on an inclined plane.Liu and Mei (1989) presented a two-dimensional theory for the unidirectional slow flow of Bingham fluid on a slope.Huang and Garcia (1998) worked on the same problem for Herschel-Bulkley fluid.For three-dimensional flows, the slow and steady spreading of mud released from a point source on a plane was investigated by Hulme (1974) with a Bingham model, and by Coussot and Proust (1996) and Wilson and Burgess (1998) with a Herschel-Bulkley model.The static problem of the final deposit on an inclined plane was studied experimentally by Coussot et al. (1996) and by Osmond and Griffiths (2001).For a horizontal plane bottom, Coussot (1997), Griffiths (2000) and Mei, Liu and Yuhi (2001) derived analytical and numerical solutions for the radically symmetric evolution.For high-speed flows, Liu and Mei (1994) and Ng and Mei (1994) examined the nonlinear formation of roll waves for a Bingham fluid and a power-law fluid, respectively.Similar problems regarding the avalanche of dry granules flowing down an inclined plane were reported by Wieland et al. (1999) and Pouliquen and Forterre (2002).Most debris-flow models focused on laboratory scale experiments or slow debris flow motion in a Published by Copernicus Publications on behalf of the European Geosciences Union.regular channel.However, debris flows occurring in the field markedly differ from those in a controlled environment.It is difficult to simulate debris flow both numerically and experimentally.Major and Iverson (1999) used a two-phase flow model to simulate debris flow moving from a large flume (2 m wide by 95 m long) to a wide deposition basin.O'Brien and Julien (1997) used a quadratic constitute to simulate high concentration flows.The DEBRIS-2D program was developed by Liu and Huang (2006) for field debris flow simulation.This model was verified by a 1-D analysis solution, laboratory testing and a field case.In addition, this model was used in numerous practical applications in Taiwan.
Although numerical simulation is considered a superior approach, the challenges for real engineering projects is in the uncertainties of numerous input data.The geographical data is available, but is rarely in high precision.Two main problems are the total amount of available soil that can be eroded or mobilized during heavy rainfall and the rheological properties that can correctly represent the field material.If these parameters are not precisely resolved, any modelling would be incorrect.However, for engineering purposes, a 20 % error in these data may be common; therefore, we must determine the errors that will be induced in the final result.If the parameters are sensitive, efforts to precisely locate the parameter must be emphasized.Conversely, for insensitive parameters, an approximate estimate may satisfy engineering purposes.Therefore, this paper focused on only a few main factors.
The DEBRIS-2D model (Liu and Huang, 2006;Liu and Hsu, 2008;Liu et al., 2009) was used to identify the sensitivity of various initial volumes relating to the final spreading of real large-scale debris flows.These results are useful for engineering designs and estimates to mitigate effectiveness.A field case of debris flow that occurred during Typhoon Morakot in August 2009 at the Daniao tribe in Eastern Taiwan was used in this study.The Typhoon Morakot produced 740.5 mm of rainfall in 62 h, and induced considerable landslides and debris flow with volumes exceeding 500 000 m 3 .The authors performed a simulation for the Daniao tribe by using DEBRIS-2D in 2007 under various design capacities (Liu et al., 2009).However, the area of influence for this event, with mitigation measures constructed, was almost the same as that predicted with no countermeasures.This verified the simulation ability of DEBRIS-2D; however, it also incurred questions regarding the consistent result.The challenge in determining the answers is the uncertainty of the input data.The geographical data was available; however, it was not highly precise.The two main problems were the total amount of available soil that may be eroded or mobilized during heavy rainfall and the properties that may correctly represent the field material.

Description of Daniao tribe debris flow
The Daniao tribe debris flow occurred in the debris flow potential stream DF097 in Taitung (Eastern Taiwan).Stream DF097 has a high debris flow disaster potential according to the information from the Soil and Water Conservation Bureau in Taiwan.The watershed area of DF097 is approximately 0.86 km 2 .A total of 71.7 % of the area has a slope greater than 15 • , 18.6 % of the area has a slope between 15 • and 6 • , and 28.3 % has a slope less than 6 • (Fig. 1).
Field investigations in September 2006 discovered several locations with large amounts of deposit.The total solid volume exceeded 19 943 m 3 distributed on the hillside and streambed of the watershed (Fig. 2).A total of 63.1 % of this material was located in regions with a slope greater than 15 • , and 7.8 % of the material was located in regions with slope less than 6 • .The formation of the mixture was composed of slate, mudstone, sandstone and weathered gravel, which are all easily movable under external forces.
The aerial photography after the disaster is shown in Fig. 5. Field investigations after the disaster revealed that almost 17.2 % (0.1485 km 2 ) of the watershed was buried, the total volume of debris flow exceeded 500 000 m 3 , and almost 200 000 m 3 flowed out of the valley.The aerial photography before and after Typhoon Morakot are shown in Figs. 6 and  7, respectively, for comparison.

DEBRIS-2D model
The original O'Brien and Julien (1985) rheological relationship from a one-dimensional version was extended to a threedimensional version, as follows: where τ ij is the shear stress tensor and γij is the strain rate tensor.τ 0 is the yield stress, µ d is the dynamic viscosity and µ c is the turbulent-dispersive coefficient.τ ij and γij represent the second invariant of the shear stress and strain rate tensor, respectively.Liu and Lai (2000) defined the portion of debris flow with stress greater than the yield stress as the boundary layer.The depth ratio between the boundary layer and the main debris flow was verified as small.This implied    that most of the flow region was in a weak stress condition, that is, the plug region.The corresponding constitutive law is shown in Eq. ( 2), which can be expressed as follows: where the x-axis coincides with the averaged bottom of the channel and is inclined at angle θ with regard to the horizon.The y-axis is in the transverse direction and the z-axis is perpendicular to both the x-and y-axes.u, v, w are the velocity components in the x, y, z directions, respectively.Because debris flow in a lab or in the field is usually considered long waves, that is, the depth scale is substantially smaller than the horizontal length scales; it can be obtained from Eq. ( 4) by neglecting the small terms, as follows: This implies that the portion of debris flow near the free surface at which the stress free condition applies is a two-dimensional plug flow [that is,u = u(z) and v = v(z)].17 where ρ is density of debris flow, and p is pressure.The stress free condition is applied at the free surface z = h(x,y,t).The upper boundary of the thin boundary layer near the bottom is defined as z = B(x,y,t)+δ(x,y,t) , where the natural bottom of the debris flow is z = B(x,y,t).Because the thickness of the boundary layer is small compared to the flow depth, the natural bottom can be used as the boundary for the plug flow.Equation ( 7) leads to static pressure in z.
Integrating Eqs. ( 5) and ( 6) in z from the bottom to the free surface obtains the results in conservative form as follows: The depth integration of continuity equation gives where H=h(x,y,t)-B(x,y,t) is the flow depth.Equations ( 9), (10), and (11) may be used to solve the three unknowns, H,u and v.This study used the Adams-Bathforth 3rd order scheme for time and central differences and the upwind scheme in space.The upwind method was used for convective terms.The central difference method was used for all other terms.Mathematically, one condition each for H,u and v is required in the physical boundary.For debris flow simulations in the field, it is usually necessary to obtain a computational domain, which contains the whole reach of the debris flows.In real applications, a large computation domain may be selected to prevent debris flows from reaching the domain boundary.Thus, the boundary conditions are as follows: H (x,y,t) = 0,u = 0,v = 0 on all boundaries (12) If debris flows are restrained in a fixed domain, such as a flume, a normal flux condition would not be used on all physical boundaries.However, Eq. ( 12) applies to the front and tail of the debris flow.The tracking of the points with a velocity near zero is crucial.Corrections of overshooting the physical quantities were performed during every time step.The initial condition is the depth contour in the computation domain with all possible debris flow sources.The value of the rheological properties is also required, which must be obtained from field sample measurements.

Daniao tribe debris flow simulation
From the 2006 survey, we discovered that at least 20 000 m 3 solid sources of debris flow were deposited on the triggering areas of the Daniao tribe watershed.However, for a heavy rainfall event, more material can be created.Therefore, we used a different approach by using accumulated rainfall to estimate the volume of debris flow in this site.An equilibrium concentration conceptual of Takahashi (1980) was applied to estimate the debris flow volume.Takahashi (1980) derived the equilibrium concentration, as follows: where C ∞ is the equilibrium concentration, ρ s is the solid density of debris, ρ w is the liquid density of debris, φ is the rest angle of solids, and θ is the slope.In general, ρ w is the water density, ρ s and φ may be measured from field samples, and θ may be calculated from the Digital Topographic Model (DTM).In the Daniao tribe debris flow watershed (Taitung DF097 watershed), the slope was calculated from the average creek bottom slope as 23.2 % (≈13 • ), as shown in Fig. 8. From the field samples, the solid density was ρ s = 2.6 g cm −3 , and φ ≈ 30 • .With the liquid density of ρ w = 1.0 g cm −3 , the equilibrium concentration was calculated as C ∞ = 41.6 % from Eq. ( 13).The amount of water required  to induce debris flow in this watershed was estimated with this equilibrium concentration.If the amount of water was insufficient to mobilize all the source material to form a debris flow, the volume of the debris flows would be smaller.This study assumed that the water that can induce debris flow originates from rainfall in the area with loose deposition and in the correct flow direction.After slope and direction analysis, we discovered that 71.7 % of the Daniao tribe watershed (equal to 0.63 km 2 , as shown in Fig. 1) satisfies this requirement.For the particular event, the rainfall accumulated to 740.5 mm in 62 h before the debris flow occurred.The records from the past 12 yr revealed that this value was      and flowed out of the valley region (in the down stream), as shown in Fig. 9a, b.The maximal velocity was in excess of 20 m s −1 during the start of the debris flow, however, it began to slow rapidly when the debris flow passed the watershed gap (maximal velocity less than 3 m s −1 ), as shown in Fig. 10b, c, d and e.The results indicated that at approximately 10 min, the debris flow reached upstream of the Daniao tribe, as shown in Figs.9d and 10d.As the velocity of the debris flow slowed to less than 0.5 m s −1 , as shown in Fig. 10f, the debris flow front peak continued to maintain the same depth (≈15 m), as shown in Fig. 9d, e and f.The final deposition fronts for all 4 cases are almost identical.The final deposition area for both the simulation and the real event are shown in Fig. 11.The simulation results were nearly consistent to the field measurements.However, drainage ditches were constructed on both sides of the village after the simulation in 2006, therefore, a part of the debris flow spread along the ditches (Fig. 11).Consequently, the front travelled a shorter distance than that in the simulation because the only difference for various total volumes is the front location.The location for the maximum as predicted by numerical simulation is 3 m away from the real location, as shown in Fig. 11.A depth near the front of the field is available by the estimation from rescuers.The estimated depth is between 12 m and 13 m at the location marked by a red star in Fig. 11.The simulation result is 13.06 m for 508 417 m 3 .
The debris flow volume is one of the crucial parameters in torrent hazard, however, it is difficult to forecast in reality because numerous uncertainties occur in a watershed.In this study, the DEBRIS-2D model was applied to assess a hazard zone with total amounts of 200 000 m 3 , 300 000 m 3 , 400 000 m 3 and 500 000 m 3 in four cases.We used 200 000 m 3 volumes as the design basis in 2006, and these volumes are also cases that we simulated.A yield stress of the debris flow of 250 dyne cm −2 was measured in the field.A time step of 0.01 s was set up, and the computational grid size was 5 m × 5 m DTM of the Daniao tribe watershed.The initial debris sources were distributed at the head of the Taitung DF097 stream.Therefore, the sensitivity analysis in various volumes is described as follows: Fig. 12 shows the final simulated deposit contour maps used in various volumes of debris flow sources.Therefore, the spread of the simulation from all four total volumes are equally acceptable.However, the maximal depth of debris flow for the final deposition was 15 m as measured in the field, and is 15.14 m for the 500 000 m 3 simulation and 15.02 m for the 200 000 m 3 simulation.The location of the maximum predicted by numerical simulation is 3 m away from the real location, as shown in Fig. 11.A depth near front of the field is available by the estimation from rescuers.The estimated depth is between 12 m and 13 m at the location marked by a red star in Fig. 11.The simulation result is 13.00 m for 500 000 m 3 and 11.74 m for 200 000 m 3 .
This study compared hazard zones in various volumes, as shown in Fig. 13, and the front positions change for various volumes, as shown in Table 1.This study revealed that a 20 % variation in estimating the volume results in a 2.76 % variation on the final deposition front.This case study of the Daniao tribe debris flow provides support for the usability of numerical simulations in real engineering detailed designs.

Figure 9 .
Figure 9.The debris flow depth contour maps at various times.3

Fig. 9 .
Fig. 9.The debris flow depth contour maps at various times.

Figure 10 .Fig. 10 .
Figure 10.The debris flow velocity vector maps at various times.3 Fig.10.The debris flow velocity vector maps at various times.

Fig. 11 .Figure 12 .Fig. 12 .
Fig. 11.Region marked in red is the area affected by Typhoon Morakot, and the region marked in blue is the simulation result.The red star indicates where field depth estimation is available.The depth estimated by rescuers is between 12 m and 13 m.The simulated result for 508 417 m 3 is 13.06 m.1

Table 1 .
Changes in front positions for various volumes.