OHD/HL - papers: 1-dim mud

ONE-DIMENSIONAL ROUTING OF MUD/DEBRIS FLOWS USING NWS FLDWAV MODEL

            
            ABSTRACT
        
               A one-dimensional 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 non-Newtonian
        fluids.  This technique involves determining the friction slope of mud/debris flows
        based on a semi-empirical rheological power-law equation and a wave-front
        tracking technique.  Three similar techniques are compared for model performance
        on three real-case mud/debris flow simulations and with some model sensitivity
        studies.
        
        INTRODUCTION
        
          Mud/debris floods, such as those caused by a landslide-induced mud/debris flow
        or those emanating from the dam-break-failure 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
        one-dimensional dynamic unsteady flow equations by adding an additional friction
        slope term in the momentum equation according to the rheological properties of
        flowing mud/debris-water 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 non-Newtonian fluid is used.
        
          The NWS FLDWAV model is a generalized dynamic flood routing model based
        on an implicit weighted four-point, nonlinear, finite-difference solution of the 
        one-dimensional unsteady flow (Saint-Venant) 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 power-law equation of a non-
        Newtonian fluid.  Also, a new wave-front tracking scheme is developed for
        modeling the mud/debris flow situations where a steep-fronted 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 wave-front 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 one-dimensional Saint-Venant unsteady flow equations used in FLDWAV
        as modified to include the mud/debris flow friction slope term, Si, 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 cross-sectional area of flow, A0 is the
        inactive (off-channel storage) cross-sectional 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, Sf is the friction slope due to turbulent boundary
        shear stress and determined by Manning's equation, Se is the slope due to local
        expansion-contraction (large eddy loss), Si is the friction slope associated with
        internal viscous dissipation of non-Newtonian mud/debris fluids, L is the
        momentum effect of lateral flow, Wf is the wind term, and B is the channel flow
        width.
        
          The additional friction slope term, Si, in Eq.(2) is obtained by applying the
        rheological power-law equation of non-Newtonian fluids to a two-dimensional
        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 power-law component of the shear stress, y is the
        yield shear strength and  is the apparent viscosity,   is the bulk density of the fluid
        mixture, Si is the friction slope, and Si=S0 in which S0 is the channel bottom
        slope.  Equation (3) can be solved for the depth mean velocity V=f(y,,y,,,Si)
        by integrating over flow depth y and assuming a parabolic velocity distribution in
        combination with a uniform velocity for y>z>y-y/(Si) (Chen, 1983); however, the
        resulting equation for V is so complicated that the friction slope, Si, cannot be
        derived explicitly and therefore this approach does not lend itself for unsteady flow
        routing purposes.  Instead, an alternative semi-empirical 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 D0=y/(Si) 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 Si can 
        be derived from Eq.(4). The derived equation for Si can be written as:
        
                                              (5)

          In this study, the following two equations for Si 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).
        
        WAVE-FRONT TRACKING TECHNIQUE
        
          Equations (1) and (2), together with one of the equations for Si (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 steep-fronted leading edge associated with the mud/debris flood
        wave.  FLDWAV contains a new wave-front 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 front-edged wave between xj and xj+1 to move.  Extensive tests show
        that this technique is excellent in simulating the moving steep-fronted 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 dam-break 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, depth-averaged 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 water-debris mixture in the reservoir is
        about 8.45 105 m3; top width of the reservoir at dam is 245 m and height of the
        dam is 21.7 m; and the dam-break induced flow lasted less than 5 minutes. A
        reservoir with a final rectangular-shape dam breach of width of 240m and a 1
        minute time for breach failure is modeled in the FLDWAV model.  It is assumed
        that cross-sections 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/m2 (0.044 lb·s/ft2), y =38 N/m2 (0.80 lb/ft2), and
         =15700 N/m3 (100 lb/ft3).  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 Si (Eqs.5-7) 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 coal-mining operation liquefied and flowed down a
        slope of about 12° (S0=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 dam-break 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 dam-failure
        hydrograph.  The mud/debris mixture properties used are: =958 N·s/m2 (20
        lb·s/ft2), y=4794 N/m2 (102 lb/ft2),  =17640 N/m3 (112 lb/ft3).  Initial flow is
        zero and the new wave-front tracking option is used.
        
          Figure 4 shows the computed wave-front travel times of the mud/debris flood. 
        These results are obtained by using the three different equations for Si.  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 Landslide-Induced Mudflow
        
          The FLDWAV model is also applied to a 1983 landslide-induced Rudd Creek
        mud/debris flow case which occurred in Davis County, Utah.  This case has been
        successfully simulated by using a two-dimensional unsteady mud/debris flow model
        (O'Brien, Julien, and Fullerton, 1993).  This landslide-induced 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 two-dimensional, it is found that the one-dimensional 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
        cross-sections is set up.  The width of the cross-sections 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 cross-sections.

        
        
          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/m2 (20 lb·s/ft2), y=956 N/m2 (20 lb/ft2), and
         =15750 N/m3 (100 1b/ft3).  The flood hydrograph used as the upstream boundary
        condition is shown in Fig.10 (x=0). Initial flow is zero and the wave-front tracking
        option is used.  The computed peak discharge and mud/debris depth profiles from
        the three expressions for Si 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 1415m3/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: Q5,6(%)
        or  Q5,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), y5,6(%) or
         y5,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 11-14 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 one-dimensional Saint-Venant unsteady flow equations can be applied to
        simulate non-Newtonian mud/debris flows if an additional friction slope (Si)
        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 wave-front tracking technique
        can be a useful tool in unsteady mud/debris flow analysis associated with  
        landslide-induced or dam-break 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., McGraw-Hill, New York, 10.1-10.36.
        
        Fread, D.L. (1988). The NWS DAMBRK model: Theoretical background and user
        documentation, HRL-258, 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), 262-265.
        
        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), 150-189.
        
        O'Brien, J.S., Julien, P.Y., and Fullerton, W.T. (1993). "Two-dimensional water
        flood and mudflow simulation," J. Hydraulic Eng., ASCE, 119(2), 244-261.
        
        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/G-85/03, Utah Water Research Laboratory,
        Utah State University, 260-278.
        
        Schamber, D.R., and MacArthur, R.C. (1985). "One-dimensional model for mud
        flows," in Hydraulics and Hydrology in the Small Computer Age, Vol.1, ASCE,
        1334-1339.
        

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.