Presented at First International Conference
on DebrisFlow Hazards Mitigation: Mechanics,
Prediction, and Assessment
USGS
San Francisco, California
August 79, 1997
ONEDIMENSIONAL ROUTING OF MUD/DEBRIS FLOWS USING NWS FLDWAV MODEL
Ming Jin
D. L. Fread
Office of Hydrology
NOAA/National Weather Service
1325 EastWest Highway
Silver Spring, Maryland 20910
ABSTRACT
A onedimensional unsteady mud/debris flow modeling technique is being
incorporated into the National Weather Service (NWS) FLDWAV dynamic flood
routing model enhancing its capability to model unsteady flows of nonNewtonian
fluids. This technique involves determining the friction slope of mud/debris flows
based on a semiempirical rheological powerlaw equation and a wavefront
tracking technique. Three similar techniques are compared for model performance
on three realcase mud/debris flow simulations and with some model sensitivity
studies.
INTRODUCTION
Mud/debris floods, such as those caused by a landslideinduced mud/debris flow
or those emanating from the dambreakfailure of a tailings or a debris dam, are a
unique unsteady flow phenomenon in which the flow changes rapidly and the
properties of moving fluid from the mixture of mud/debris and water are very
different from pure water. One method of modeling this special flow is to use the
onedimensional dynamic unsteady flow equations by adding an additional friction
slope term in the momentum equation according to the rheological properties of
flowing mud/debriswater mixtures. The derivation of the friction slope term of the
mud/debris flow depends on which rheological model (constitutive equation) for
shear stress of a nonNewtonian fluid is used.
The NWS FLDWAV model is a generalized dynamic flood routing model based
on an implicit weighted fourpoint, nonlinear, finitedifference solution of the
onedimensional unsteady flow (SaintVenant) equations. FLDWAV combines the
capabilities of the popular NWS DAMBRK and DWOPER models (Fread,1993)
and provides some additional features. A recent enhancement of the FLDWAV
model is a new mud/debris flow routing technique in which the mud/debris flow
friction slope is derived from the shear stress powerlaw equation of a non
Newtonian fluid. Also, a new wavefront tracking scheme is developed for
modeling the mud/debris flow situations where a steepfronted leading edge of the
mud/debris wave propagates along an initially dry channel (zero initial flow), and
the downstream boundary of the unsteady mud/debris flow is the wave front. In
this paper, the new mud/debris flow and wavefront tracking technique is presented
and its performance is tested in modeling three real mud/debris flow case studies.
It is also compared with two existing expressions for the mud/debris flow friction
slope term for three case studies and with some sensitivity studies.
EQUATIONS AND MODEL FORMULATION
The onedimensional SaintVenant unsteady flow equations used in FLDWAV
as modified to include the mud/debris flow friction slope term, S_{i}, are (Fread,
1988, Fread 1993):
(1)
(2)
in which t is time, x is distance along the longitudinal axis of the waterway, h is the
water surface elevation, A is the active crosssectional area of flow, A_{0} is the
inactive (offchannel storage) crosssectional area of flow, q is the lateral inflow or
outflow, is the coefficient for nonuniform velocity distribution within the cross
section, g is the gravity constant, S_{f} is the friction slope due to turbulent boundary
shear stress and determined by Manning's equation, S_{e} is the slope due to local
expansioncontraction (large eddy loss), S_{i} is the friction slope associated with
internal viscous dissipation of nonNewtonian mud/debris fluids, L is the
momentum effect of lateral flow, W_{f} is the wind term, and B is the channel flow
width.
The additional friction slope term, S_{i}, in Eq.(2) is obtained by applying the
rheological powerlaw equation of nonNewtonian fluids to a twodimensional
steady uniform open channel flow of depth, y, as follows:
(3)
in which is the internal shear stress, u=u(z) is the longitudinal velocity in the x
direction, is an exponent of the powerlaw component of the shear stress, _{y} is the
yield shear strength and is the apparent viscosity, is the bulk density of the fluid
mixture, S_{i} is the friction slope, and S_{i}=S_{0} in which S_{0} is the channel bottom
slope. Equation (3) can be solved for the depth mean velocity V=f(y,,_{y},,,S_{i})
by integrating over flow depth y and assuming a parabolic velocity distribution in
combination with a uniform velocity for y>z>y_{y}/(S_{i}) (Chen, 1983); however, the
resulting equation for V is so complicated that the friction slope, S_{i}, cannot be
derived explicitly and therefore this approach does not lend itself for unsteady flow
routing purposes. Instead, an alternative semiempirical equation which produces
an approximate solution to Eq.(3) is proposed:
(4)
in which m=1/ and m=1 represents a Bingham fluid, D is the hydraulic depth,
and D_{0}=_{y}/(S_{i}) can be regarded as the minimum depth for the mud/debris mixture
to move because of the yield shear strength. The difference of the velocity profiles
from Eq.(4) and that from Chen's equation is less than 5%, but an equation for S_{i} can
be derived from Eq.(4). The derived equation for S_{i} can be written as:
(5)
In this study, the following two equations for S_{i} are also tested and compared
with Eq. (5): (1) the equation used in NWS DAMBRK model which comes from a
similar derivation from Eq.(3) in which a parabolic velocity distribution is assumed
(Fread 1988); and (2) the equation based on a linear velocity distribution of laminar
Bingham fluids (Jeyapalan, Duncan and Seed, 1983; Schamber and McArthur,
1985). These equations are expressed, respectively, as:
(6)
(7)
Equation (6) is equivalent to an equation used by O'Brien and Julien (1985) for a
Bingham fluid (m=1). It was slightly modified in a later application (O'Brien,
Julien, and Fullerton, 1993).
WAVEFRONT TRACKING TECHNIQUE
Equations (1) and (2), together with one of the equations for S_{i} (Eq.(5), (6), or
(7)), are solved numerically with appropriate external (upstream/downstream) and
internal (dam/bridge) boundary conditions. One method in routing unsteady mud/
debris flows is to simulate them from an assumed initial mud/debris flow condition
throughout the entire routing reach. There are many cases, however, where the
mud/debris mixture moves over an initially very small water flow or a dry channel
and often has a steepfronted leading edge associated with the mud/debris flood
wave. FLDWAV contains a new wavefront tracking technique in which the model
tracks the moving wave front as its computational downstream boundary and uses
an automatically generated Q=f(y) loop rating as the boundary condition. Moving
of the downstream boundary is controlled by checking, at every time step, the
mud/debris flow volume passed from the current boundary x=xj with the minimum
volume for a frontedged wave between x_{j} and x_{j+1} to move. Extensive tests show
that this technique is excellent in simulating the moving steepfronted waves of
mud/debris flow from a zero (dry bed) or a very small initial flow condition.
APPLICATION CASE STUDIES
Case 1. Anhui Debris Dam Failure Flood
A tailings dam of the Jinshan debris reservoir in Anhui, China, breached in the
early morning of April 30, 1986 (Han and Wang, 1996). The dambreak induced
mud/debris flooding engulfed a village about 0.75 km downstream of the reservoir,
and all of the village residents were killed in the disaster. Measurements of the
inundation area were made after the flooding event.
Han and Wang simulated the unsteady mud/debris flow using a two
dimensional, depthaveraged model and assumed an inflow hydrograph as the
upstream boundary condition. Data provided by these authors was used in the one
dimensional FLDWAV model to simulate the outflow from the breached dam. The
following data were used: total volume of waterdebris mixture in the reservoir is
about 8.45 10^{5} m^{3}; top width of the reservoir at dam is 245 m and height of the
dam is 21.7 m; and the dambreak induced flow lasted less than 5 minutes. A
reservoir with a final rectangularshape dam breach of width of 240m and a 1
minute time for breach failure is modeled in the FLDWAV model. It is assumed
that crosssections are irregular trapezoids with an average width of 210m to 580m
and a channel bottom slope from about 0.012 upstream to 0.00076 downstream.
Values of 0.035 and 0.04 are used for Manning's n. The following Bingham fluid
properties are used: =2.1N·s/m^{2} (0.044 lb·s/ft^{2}), _{y} =38 N/m^{2} (0.80 lb/ft^{2}), and
=15700 N/m^{3} (100 lb/ft^{3}). Since the initial flow is almost zero, the new wave
front tracking option is selected for the routing, and Eq.(5) is used to determine the
friction slope associated with the internal viscous dissipation of the mud/debris
flow.
Figure 1 shows computed mud/debris surface profiles at t=0.005, 0.01, 0.02,
0.05 and 0.09 hours. The dam breaching started at t=0.0 due to an assumed
overtopping failure, and the mud/debris mixture wave front propagated downstream
to a final inundation limit within a total time of about 5 minutes. This agreed with
the site report that the flooding lasted less than 5 minutes. The computed flooding
distance of 1200m compares well with the observed inundation distance of about
1210m. Figure 2 shows the computed discharge hydrographs at three locations
along the reach (x=0 at dam site, x=400m, and x=800m). One characteristic
feature simulated by the model is that the mud/debris flood wave moves with a
steep front and both the discharge and stage hydrographs reach their peak very
quickly. This, along with the greater density of mud/debris floods, contributes to
the fact that even a small debris flood can cause devastating damage in life and
property.
Figure 3 compares the computed peak discharge profiles using the different
equations for S_{i} (Eqs.57) to determine the additional friction slope. In this case
using these different equations produces only small computational differences.
Case 2. Aberfan Coal Waste Dump Failure
The Aberfan mud/debris flow disaster occurred in Wales in 1965 and was
investigated by using an analytical mudflow simulation method based on the laminar
Bingham fluid model (Eq.(7)) (Jeyapahan et al., 1983). In this case, waste material
from a 37 m high dump from a coalmining operation liquefied and flowed down a
slope of about 12° (S_{0}=0.208) over a distance of about 600 m until it inundated a
school and other buildings located in the flow path. One hundred and twenty lives
were lost in the buried portion of the school building. The FLDWAV model is
applied to this event by simulating it as a dambreak mud/debris flow case. The
downstream channel has triangular cross sections with a side slope of about 1
vertical to 2 horizontal, and Manning's n of 0.05 is used in the computation. A
waste reservoir is placed upstream, and the dam is breached at the beginning of
routing. The final breach has the same shape as the triangular channel cross
sections and a time of failure of about 1 minute is used to simulate the damfailure
hydrograph. The mud/debris mixture properties used are: =958 N·s/m^{2} (20
lb·s/ft^{2}), _{y}=4794 N/m^{2} (102 lb/ft^{2}), =17640 N/m^{3} (112 lb/ft^{3}). Initial flow is
zero and the new wavefront tracking option is used.
Figure 4 shows the computed wavefront travel times of the mud/debris flood.
These results are obtained by using the three different equations for S_{i}. An
observed data point is also shown in the figure. The results from using Eqs. (5)
and (6) agree closely with the observed data, and the results from using Eq.(7) are
only slightly different. The computed maximum flow depth profiles are shown in
Fig.5. A small difference is noticed between using Eq.(5) and Eq.(6), while a
larger difference is noted in the use of Eq.(7). Figure 6 shows computed discharge
hydrographs at three locations using Eq.(5).
Case 3. Rudd Creek LandslideInduced Mudflow
The FLDWAV model is also applied to a 1983 landslideinduced Rudd Creek
mud/debris flow case which occurred in Davis County, Utah. This case has been
successfully simulated by using a twodimensional unsteady mud/debris flow model
(O'Brien, Julien, and Fullerton, 1993). This landslideinduced mud/debris flow
submerged a residential block just downstream of the landslide. In this case, the
mud/debris mixture flowed along the hill slope and not in single channel. Although
the phenomena is twodimensional, it is found that the onedimensional model can
also be used successfully. Figure 7 shows a topological map of the site. The
landslide occurred at an elevation of about 4500ft (1372m). The possible flow
paths are drawn in the figure. A nonprismatic channel with specifiled rectangular
crosssections is set up. The width of the crosssections are estimated as the
distance between the possible flow paths and a reasonable extension beyond the
paths according to anticipated flow depth, and the channel bottom slope can be
determined according to an average elevation of the crosssections.
The debris properties are determined according to the an average of about 48%
by volume mud/debris sediment concentration and an empirical relationship from
experimental data (O'Brien and Julien, 1985). The following properties are used in
the computations: =958 N·s/m^{2} (20 lb·s/ft^{2}), _{y}=956 N/m^{2} (20 lb/ft^{2}), and
=15750 N/m^{3} (100 1b/ft^{3}). The flood hydrograph used as the upstream boundary
condition is shown in Fig.10 (x=0). Initial flow is zero and the wavefront tracking
option is used. The computed peak discharge and mud/debris depth profiles from
the three expressions for S_{i} are shown in Figure 8. The results of using Eq.(5)
compares very well with the observed maximum inundation distance, while using
Eq.(6) produced a shorter inundation distance and using Eq.(7) produced a longer
distance. Figure 9 shows the computed times of maximum flow depth from which
the mud/debris flow wave speed can be determined. The average wave speeds
between x=100m and x=300m are 0.93 m/s (using Eq.(5)), 0.67 m/s (using
Eq.(6)) and 1.67 m/s (using Eq.(7)). Eyewitness accounts estimated the wave
speed from 0.6 to 1.2 m/s. Figure 10 shows the upstream hydrograph and
computed hydrographs at two other locations using Eq.(5).
SENSITIVITY STUDIES
The case studies have shown that among the three equations, Eqs.(5), (6), and
(7), Eq. (5) provides the best performance in modeling these mud/debris flows.
Some numerical experiments were also conducted to further investigate the
differences in computational results by using the three equations for a range of
channel slopes and mud/debris flow properties.
A mud/debris flood wave with a peak discharge of 1415m^{3}/s and 0.1 hour time
of rise is routed through a 3218m long prismatic channel with rectangular cross
sections of width of 60m and Maning's n of 0.05. The channel is equally divided
into two reaches by two bottom slopes. The upstream reach has a slope of
S(1)=0.0379 (200 ft/mile) and the slope of the downstream reach, S(2), is changed
to examine the effect of the slope.
In order to evaluate the computational differences in the mud/debris flow peak
profiles, the results from Eq.(5) are compared with those from Eq.(6) and from
Eq.(7), and the following quantities are used to measure the differences: Q_{5,6}(%)
or Q_{5,7}(%) is an average percentage difference in the peak discharge profile
between results of Eq.(5) and Eq.(6) and between Eq.(5) and Eq.(7), y_{5,6}(%) or
y_{5,7}(%) is an average percentage difference in the peak mud/debris flow depth
profile between results of Eq.(5) and Eq.(6) and between Eq.(5) and Eq.(7). The
flow conditions and results are listed in Table 1.
Figures 1114 show some results from one group in which S(2)=0.5 S(1) and
the computed peak discharge and mud/debris flow depth profiles from the three
equations are shown in the figures. It is noticed from these figures that the results
from Eq.(5) are between those from Eq.(6) and those from Eq.(7). Higher values
of the viscosity, causes a larger departure of the results of Eq.(7) from those of
Eq.(5) and Eq.(6), while higher values of the yield shear strength, _{y,} results in a
larger departure of the results of Eq.(6) from those of Eq.(5) and Eq.(7).
These results indicate that the overall difference between Eq.(5) and Eq.(6) or
between Eq.(5) and Eq.(7) are less than 15% in most situations. The largest
difference between Eq.(5) and Eq.(7) occurs under the conditions where both
and _{y,} are large values, while the largest difference between Eq.(5) and Eq.(7)
occurs when the value of _{y,} is large. The results in groups E and F show that the
difference of Eq.(5) and Eq.(6) increases as the fluid departs from a Bingham fluid
( =1) and becomes viscoplastic ( >1). In the FLDWAV model, can be
specified as any value.
CONCLUSIONS
The onedimensional SaintVenant unsteady flow equations can be applied to
simulate nonNewtonian mud/debris flows if an additional friction slope (S_{i})
representing the internal viscous dissipation is appropriately specified. The
excellent computational results in the case studies show that the mud/debris flow
enhanced FLDWAV model which uses Eq.(5) and a wavefront tracking technique
can be a useful tool in unsteady mud/debris flow analysis associated with
landslideinduced or dambreak induced mud/debris flood prediction. The
sensitivity studies suggest that use of Eq.(5) in determining the mud/debris flow
friction slope produces computational flow peak profiles which are between those
from Eq. (6) and those from Eq. (7). The difference between Eq. (5) and Eq. (6)
or between Eq. (5) and Eq. (7) are less than 15% for Bingham fluids but can be 20
30% for viscoplastic mud/debris flows.
APPENDIX I. REFERENCES
Chen, C.L. (1983). "On frontier between rheology and mudflow mechanics," in
Frontiers in hydraulic engineering, Edited by Shen, H.T., ASCE, New York, 113
118.
Fread, D.L. (1993). Chpt. 10: Flow routing, in Handbook of Hydrology, Edited by
Maidment, D.R., McGrawHill, New York, 10.110.36.
Fread, D.L. (1988). The NWS DAMBRK model: Theoretical background and user
documentation, HRL258, Hydrological Research Laboratory, National Weather
Service, Silver Spring, Maryland 20910.
Han, G., and Wang, D. (1996). "Numerical modeling of Anhui debris flow,"
J. Hydraulic Eng., ASCE, 122(5), 262265.
Jeyapalan, J.K., Duncan, J.M., and Seed, H.B. (1983). "Analysis of flow failures
of mine tailing dams," J. Geotechnical Eng., ASCE, 109(2), 150189.
O'Brien, J.S., Julien, P.Y., and Fullerton, W.T. (1993). "Twodimensional water
flood and mudflow simulation," J. Hydraulic Eng., ASCE, 119(2), 244261.
O'Brien, J.S., and Julien, P.Y. (1985). "Physical properties and mechanics of
hyperconcentrated sediment flows," in Delineation of landslide, flash flood and
debris flow hazards in Utah, UWRL/G85/03, Utah Water Research Laboratory,
Utah State University, 260278.
Schamber, D.R., and MacArthur, R.C. (1985). "Onedimensional model for mud
flows," in Hydraulics and Hydrology in the Small Computer Age, Vol.1, ASCE,
13341339.
As required by 17 U.S.C. 403, third parties producing works consisting predominantly of the material appearing in NWS Web pages must provide notice with such subsequently produced work(s) identifying such incorporated material and stating that such material is not subject to copyright protection.
