(1)
H. Y. Lee
Prof.,
Dept. of Civil Engr., National Taiwan Univ., Taipei, China.1 Section 4
Rooservelt Road Taipei 106, Taiwan
,China，Daytime
phone number: 886-223630231-2404-17，Fax
number: 886-223631558，E-mail
address: hylee@ce.ntu.edu.tw
H. M. Hsieh
Associate Researcher , Water Resources Research Center , Taichuang, Taiwan, China 5f-3 122-19 Section 2 Taichunggun Road Taichung 407, Taiwan ,，Daytime phone number: 886-42555988，Fax number: 886-425551586，E-mail address: hmhsieh@ms3.hinet.net
Abstract: A numerical model, which is capable of simulating scouring and deposition behaviors in a channel network, is developed in this study. The model treats suspended load and bed load separately ,and hence is able to simulate the depositional behavior of the suspended sediment under a nonequilibrium process. The model solves the de Saint Venant equation, and thus can be applied to unsteady flow conditions. An internal boundary condition based on the sediment transport capacity was proposed to distribute the incoming sediment load into the downstream links. The application of this model to the Tanhsui River system in Taiwan gave very convincing results.
Keywords: channel network , scour and deposition , alluvial channel simulation
During the past decade, several numerical models have been developed. A comprehensive literature survey has been conducted by the authors in their previous paper [6], and will not be repeated herein. A new model, named NETSTARS, is developed in this article to simulate the channel bed evolutions of a channel network. It treats the suspended load and bed load separately ,and hence is able to simulate the depositional behaviors of the suspended sediment under a nonequalibrium process. An internal boundary condition based on the sediment transport capacity was proposed to distribute the incoming sediment load into the downstream links. The model is applied to simulate the bed evolution of the Tanhsui River system, which is an esturine channel network to demonstrate its engineering applicabilities.
The NETSTARS model is an uncoupled sediment routing model. It consists of two major parts, namely hydraulic routing and sediment routing. Suspended load and bed load are treated separately in sediment routing. A network algorithm was proposed to resolve the nodal point problem. A node is assumed to be a virtual section that could not accumulate water and sediment, i.e. no scour and no deposition. A nodal continuity equation was proposed in this paper. Double sweep method is applied to solve the discharges and the water stages at those nodes. The allocation of the suspended load at the node is assumed to be proportional to the allocations of the flow discharges, and the net flux of the suspended sediment due to longitudinal dispersion is assumed to be zero at a nodal point.
The NETSTARS model avails some good ideas of CHARIMA[4], BRALLUVIAL[5] and GSTARS[7] models to develope a powerful tool for resolving the problems of unsteady sediment process in a channel network. Since it is basically a 1-D model, secondary current and local scour can not be simulated by this model.
The de Saint Venant equations are used in the unsteady flow computation. These include a continuity equation and a one-dimensional momentum equation.
(1)
(2)
where A = channel
cross-sectional area ; Q = flow discharge ; t = time ; x = coordinate in the
flow direction ; q = lateral inflow/outflow discharge per unit length ; b=momentum
correction coefficient; g = gravitational acceleration ; y= water surface
elevation ;
= friction slope;
= channel conveyance;
= roughness coefficient of Manning's formula ,and
= hydraulic radius.
The governing equations include a sediment continuity equation, a sediment concentration convection-dispersion equation and a bed load equation. The Rouse number W/κu*, where W = sediment fall velocity , κ= Von Karman's constant and u* =shear velocity, is used to distinguish between bed load and suspended load. Particles with W/κu*＞5 are treated as bed load and particles with W/κu*≦5 are treated as suspended load. The sediment continuity equation is given as:
(3)
where
= bed load transport rate and
= depth-averaged concentration of the suspended sediment of size fraction
.
= amount of sediment scouring/deposition per unit length and
=channel bed porosity.
The concentration
is calculated using the
convection-dispersion equation shown as:
(4)
where
= longitudinal dispersion
coefficients., and
= source term of the suspended
sediment of size fraction k. The second term of equation (4) in the left hand
side is the convection term. The first, second, and third terms of equation (4)
in the right hand side are longitudinal dispersion, and reaction terms,
respectively. According to Van Rijn [8] and Holly and Rahuel [3], the source
term
is the combination of deposition
and resuspension.
The bed composition computation procedure is proposed by Bennett and Nordin [1]. According to this procedure the thickness of the active layer is set equal to a preselected parameter, Nt, times the geometric mean of the largest size class used in the simulation. The active layer is defined as the bed material layer that can be worked or sorted through by the action of the flowing water in the time step, t, to supply the volume of material necessary for erosion. The bed armoring is achieved by the erosion of a particular size of material being limited to the amount of that size available in the active layer during a time step.
The simulation processes consist of two parts in every time step, i.e., flow computations, and sediment routing. Flow computation is performed first. Then, sediment routing is performed to calculate the amount of channel bed variations. Eqs.1 and 2 are transformed into difference equations using a Preissmann four point finite difference scheme, and solving by double sweep method. The difference equation for the sediment continuity equation, i.e., equation (3), is shown as:
(5)
where ΔZi
= variation of the bed elevation for every size fraction and Pi =
wetted parameter. The concentration Ck is predetermined by solving
the convection-dispersion equation, i.e., equation (4).
The split operator approach is used in solving this equation. The
governing equation is separated into three portions, i.e., convection,
longitudinal dispersion, and reaction. They are solved subsequently in one time
step. The Ck and CXk
values, where
, obtained in the previous computation are served as the known values for the
next computation. Holly-Preissmann two-point four-order scheme, Tee scheme, and
analytical method are used to solve the Ck and CXk values.
The flow discharge at the nodal point has to
satisfy the continuity equation, i.e., the summation of the inflow and outflow
discharges from all the tributaries at a nodal point must be zero, or there is
no storage at the nodal point. The nodal solution uses de St. Venant equations
along the links between nodes to express the
values in terms of corrections to
water surface elevations at the nodes. The solution algorithm comprises three
phases for each iteration: namely, link forward sweep, node matrix loading and
link backward sweep.
The bed elevation at a nodal point is assumed to be in an equilibrium condition, i.e., there is no net scouring and sediment deposition at the nodal point and total sediment flux at the nodal point equals to zero. Assuming the incoming sediment is fully mixed at the nodal point and then is distributed, according to the sediment transport capacity, to the downstream links. The bed load and suspended load are redistributed and treated separately at the nodal point.
A bed load
equation is used to calculate the bed load transport capacity at the first
section downstream of the nodal point and a rating curve to describe the
relationship between the flow discharge and the bed load transport capacity,
was established.
Where
=the bed load transport capacity, and
and
=constants to be calibrated. The outflowing bed load transport rate is
distributed according to the bed load transport capacity obtained from the
rating curve.
The suspended load is calculated by solving the sediment convection-dispersion equation using the split operator method. There are three portions, namely , Convection , longitudinal dispersion and reaction, involved in the split operator method and each portion has to be treated differently at the nodal point. The split operator method had been applied in numerical solution of the one- and two- dimensional convection-dispersion equation by Cunge, Holly, and Verwey, 1980.[2] The algorithm used to split the sediment load at each node was originally proposed by Holly, Yang, et al. , 1990[5].
(1)
Convection: The sediment concentration
and corresponding concentration
gradient are assumed to be fully mixed at the nodal point and then redistributed
according to the weight of the flow discharge.
(2) Longitudinal dispersion: Assuming net flux of the suspended sediment due to the longitudinal dispersion effect equals to zero at a nodal point.
(3) Reaction: Since the suspended sediment is assumed fully-mixed at the nodal point and hence transverse dispersion and reaction treatments are not needed at the nodal point.
The Tanhusi River system passes through Taipei area, and is one of the most important river in Taiwan. It consists of of three branches, namely, Keelung River, Hsindan Creek , Tahan Creek, and one flood by pass channel, is a typical channel network. The location map is shown in Fig.1. This network consists of 5 links, 6 nodes and 88 crosssections and is within the estuarine area. There are one gage station, namely, Taipei bridge (link 3, Sec.13) in the study area.These data, including water stages and channel bed variations, can be used to calibrate and verify the model. Field data from 1989, including geometric cross-sectional data and bed material data, are used as the inital condition. Data from 1989 to 1990 were used to calibrate the model and data from 1990 to 1991 were used for verification. The upstream inflow suspended sediment concentrations versus the inflow water discharge rating curves obtained by the Taiwan Provincial Water Resources Department are used for the upstream boundary conditions. At the downstream boundary, which is located at Tudigonbi (Link 1, Sec.1), the measured stage hydrographs are used as the downstream boundary conditions. The model has to be "spinned" before real computation proceeds. A flow condition was assigned and let the model run until the water stages and discharges are smoothly-connected in the channel network, and this condition serves as the initial condition. The initial conditions for the sediment concentrations were then determined by the corresponding rating curves. The data of 1990 were used to calibrate the model. A time interval of 1 hour was used in the simulation and the total simulation period is 1489 hours. The Nt value of the thickness of the active layer is set to 1.0. The ranges of size fraction are from 0.016 mm to 16.0 mm.
The measured water stage data from five different gage stations were used to calibrate the Manning's n value of the model. The agreements are very good. The simulation results of Taipei bridge are shown in Fig.2. The simulated bed elevations for links 4 are shown in Fig.3. The agreements are satisfactory.
Using the parameters determined in the previous analysis, the model is applied to simulate the channel bed evolutions of Tanhsui River System from 1990 to 1991. The total simulation period is 672 hours. The comparisons of the longitudinal water surface elevations, and longitudinal bed profiles are shown in Figs.4, and 5 respectively. The overall accuracy is good.
A numerical model, which is capable of simulating scouring and deposition behaviors in a channel network under an unsteady flow condition, is developed in this study. A state-of-the-art sediment routing algorithm which is capable of simulating suspended and bed loads seperately was developed in this model. An internal boundary condition based on the sediment transport capacity was proposed to distribute the incoming sediment load into the downstream links. The proposed measure is proved to be feasible. The model's performance and applicability have been demonstrated through an application to the Tanhsui River system.
Acknowledgement
This study is financially supported by the Sinotech Research Foundation. Helpful comments from Mr. Hsu, R.L., Mr. Chen, L.M. and Mr. Hsieh, K.C., of the Sinotech Consultant Company and Dr. Yang, C.T. from U.S. Bureau of Reclamation are appreciated.
References
[1] Bennet, J.P. and Nordin, C.F. (1977), Simulation of Sediment Transport and Armouring , Hydrological Sciences Bulletin, XXII.
[2] Cunge,J.A., Holly,F.M.,Jr, and Verwey, A.(1980), Practical Aspects of Computational River Hydraulics, Pitman Publishing Ltd., London,, 00.312-349.
[3] Holly, F.M. ,and Rahuel, J.L. (1990), New Numerical/Physical Framework for Mobile-Bed Modelling, Part I, Journal of Hydraulic Research, Vol.28, NO.4.
[4] Holly, F.M., Yang, J.C., Schovarz, P., Scheefer J., Hsu, S.H. and Einhellig, R. (1990), CHARIMA Numerical Simulation of Unsteady Water and Sediment Movements in Multiply Connected Networks of Mobile-Bed Channels, IIHR Report No.0343, The University of Iowa, Iowa City, Iowa, U.S.A..
[5] Holly, F.M., Yang, J.C. and Spasojevic, M.(1985), Numerical Simulation of Water and Sediment Movement in Multi-Connected Networks of Mobile Bed, Iowa Institute of Hydraulic Research, Limited Distribution Report No.131, The University of Iowa, Iowa City, Iowa, U.S.A..
[6] Lee, H.Y., Hsieh, H.M., Yang, J.C. and Yang, C.T. (1997), 'Quasi-Two-Dimensional Simulation of Scour and Deposition in Alluvial Channels, Journal of Hydraulic Engineering, ASCE, Vol.123, No.7, July, pp.600-609.
[7] Molinas, A.M. and Yang, C.T. (1986), Computer Program User's Manual for GSTARS, U.S. Department of Interior Bureau of Reclamation, Engineering and Research Center, Denver, Colorado.
[8] Van Rijn, L.C. (1984), Sediment Transport, Part II: Suspended Load Transport, Journal of Hydraulic Engineering, ASCE, Vol.110, No.11, 1613-1641.
Fig. 1 The location map of Tanhsui river
Fig. 2 The calibrated water surface elevations at the Taipei bridge
Fig. 3 The calibrated longitudinal bed elevations for link 4
Fig. 4 The verified water surface elevations at the Taipei Bridge
Fig. 5 The verified longitudinal bed elevations for link 4