- Email: [email protected]
Comparative analysis of some explicit-implicit streamflow models
Comparative analysis of some explicit-implicit streamflow models Ify. L. Nwaogazie
Department of Civil Engineering, University of Port Harcourt, P.M.B. 5323, Port Harcourt, Rivers State, Nigeria The numerical performance and cost analysis of four streamflow models are presented. These models are constructed from the well-known equations of Saint-Venant. The Galerkin finite element method and the Newton-Raphson iterative method were utilized for the models' solution of depth and velocity of flow. Because of the large computer storage associated with complete solutions of Saint Venant equations, particularly for floods of long durations, the approximate models (kinematic- and diffusion-waves) were introduced to investigate the savings which could be made by using each of these models. Model predictions were contrasted with previously found solutions of a flood-wave in an idealized geometry. The effects of large time step (At) and the timeweighting factor, (0) for implicit models on cost and the numerical distortion of the hydrographs were examined. Results indicate that the computer cost is a direct function of the level of model approximation. For field application of the implicit models, regression equations relating: (1) Manning's roughness coefficient and discharge; (2) area and depth of flow; and (3) top-width and depth of flow were utilized. Relative to idealized channel the degree of nonlinearity for the natural channel not only affected model performance but increased the amount of computations and thus the cost. Convergence and unconditional stability were observed for implicit models for 0 in the range of 0.55 and 1.00. For the explicit model, instability restricted the choice of At. Key Words: streamflow explicit-implicit models, finite element, Newton-Raphson, timeweighting factor INTRODUCTION The process of determining the water depths, velocities and discharge in channels, rivers or reservoirs under unsteady conditions arising from flood motions is commonly referred to as flood routing ~. Interest in flood routing and in part in unsteady flows in water resources stems from the need to plan, design, regulate and manage our flood prone areas and other water resources systems. Presently, two basic techniques for unsteady flow simulation are in use, namely: (1) methods which approximate the basic equations of unsteady flow; and (2) methods which solve the basic equations directly. The basic equation as used herein refers to the complete equations of Saint Venant 2. The first methods are sometimes referred to as 'hydrologic' routing methods and the second methods 'hydraulic' routing methods 3. Basically, the hydrologic routing method neglects the momentum equation and solves the continuity equation which relates: (1) inflow into a given reach; (2) outflow from that reach; (3) time period for the flow to travel through the reach; and (4) the change in storage during the time period in the reach. A graphical relationship established between storage and outflow 4 yields a feasible solution to the continuity equation. This solution approach, though analytical is greatly simplified by the averaging technique. Accepted September 1986. Discussion closes August 1987. 0309-1708/87/020069-0952.00
© 1987 Computational Mechanics Publications
A number of hydrologic routing methods have emerged over the years including various coefficient routing techniques such as the Muskingum method 5. These methods are handy when data for hydraulic routing are not available. For detailed discussion on hydrologic routing methods, interested readers are referred to a basic text on hydrology 14. The importance of the hydraulic routing method has come to the limelight with the development of modern high-speed digital computers, allowing numerical solution of unsteady partial differential equations 2 that have no known analytical solutions to be obtained. Numerous unsteady flow phenomena such as surges, effects of tidal fluctuations, back-water resulting from channel junctions or reservoirs, and normal flood waves from excessive rainfall can be analysed using the hydraulic routing and numerical methods. In this paper, four numerical models namely: (1) explicit kinematic finite element model (EKFEM)~; (2) weighted-implicit kinematic finite element model (WlKFEM)6; (3) weighted-implicit diffusion finite element model (WIDFEM)7; and (4) weighted- implicit complete finite element model (WlCFEM) s'9 are presented. The purpose of comparing these is to enable the model-user to assess the individual model performance relative to their level of approximations and solution methodology. Given that the models are solved using the same procedure - weighted-implicit finite
Adv. Water Resources, 1987, Volume 10, June
69
Streamflow
models: I. L. Nwaogazie
method and iterative Newton-Raphson element algorithm - any observable difference in the predicted results for a typical flood routing problem may be attributed to the level of model approximation and solution technique (explicit or implicit). Often, implicit models that are unconditionally stable may conveniently employ large time- and distance-steps with the advantage of reduced computer cost*-“. On the other hand, the trade-off in the use of large time steps becomes too obvious in loss of numerical accuracy particularly the degree of attenuation of peak flow and the precise time of peaking. For instance, consider a flood wave of 3 hours duration in a channel of constant geometry, where the time to peak is to be predicted. Assuming the actual time to peak is 30 minutes, then a selected time step (At) of 600 seconds as compared to 60 seconds will not only dampen the peak flow but will yield fewer plotting points than those of At =60 seconds. In these circumstances, the precise time of peaking may be difficult to isolate. These are points of special interest in this paper. To establish confidence in the capability of any proposed numerical model, it must be validated using three levels of testing”. Level-I testing represents the model validation against analytic solutions. Level-II testing compares the model against at least an existing numerical model. And Level-III testing consists of validating the model using field data. In this paper only Levels-II and III are examined. Effects of time steps and time-weighting factors are considered along with a hypothetical flood problem in a rectangular channel. Cost and model performance are contrasted for individual models in field application and idealized channel geometry, respectively. STREAMFLOW Governing
MODELS
equations
The distribution of discharge, depth of flow, and velocity of flow in a stream are represented in Fig. 1. The mathematical models that predict the depth and velocity on a space and time basic are the continuity and momentum equations’ respectively:
~+?.;+c~-q(x,t)=o
(1)
~+~~+g~+~q(_x,t)-g(so-s,)=O (2) ’
?
Q=t’A
(3)
dx,t)
I
c
Y//T\
I /I&\
//A\
Streamjlow
AX
Fig.
1.
70
Adv. Water Resources,
4
Y
A
”
Y/A\
element
1987,
Volume 10, June
w//K
in which _Vrepresents the depth of flow, r the mean velocity of flow across the section, x the distance along the channel, and t time. q(x, t) represents the lateral inflow into the channel reach (negative for outflow), g the gravity acceleration, sI the friction slope, so the longitudinal bottom slope of the channel, Q the discharge and A the area of flow. The two dependent variables in equations (1) and (2) are the depth of flow, ~~(x,r) and the velocity of flow, t’(x, t). Equations (l), (2) and (3) can represent both an idealized channel geometry and natural rivers. For a natural channel, the geometry is specified by the area of flow, A(x); the top-width, B(x) and the longitudinal slope, so = so(x). The variables A(x) and B(x) are modelled using a fourth-order regression equation -nodal areas and topwidths as functions of depth of flow. Examples of the regression equations are illustrated in this paper in a case study of Illinois River in Oklahoma, USA. The transient values of ,4(x) and B(x) are used to update the nodal discharge values as well as the individual Jacobian terms during iterative simulation. Simplified and complete
flow models
From equations (1) and (2) three distinct streamflow models are possible by completely retaining equation (1) while simplifying assumptions are made in equation (2). By making the friction slope, s/, the subject of the formula and setting q(x, t) equal to zero, equation (2) becomes:
av c au 1 dv s~=so--L___-__ ax gsx g at Two simplified models and a complete flow model are generated using equation (4). These models are presented in the order of increasing complexity as follows: 1. Kinematic flow model. If the last three terms on the right-hand side of equation (4) are small compared to so, the following uniform flow equation6 results: n2c2 ‘/= ‘0 = 2.22 R4’3 in which n denotes Manning’s roughness coefficient and R the hydraulic radius. Equation (5) implies that gravity is balanced by resistance and when coupled to equation (1) is known as the kinematic flow model. The kinematic flow model is best suited for streamflow routing in steep gradient of the order of 10 feet per mile’ (1.9 m per km) or more and overland flows”. 2. Diffusion flow model. If the last two terms on the right-hand side of equation (4) are small compared to so, the steady nonuniform flow equation’ results:
Simultaneous solution of equations (1) and (6) give rise to what is known as the diffusion flow model, which has been successfully applied to very flat gradients’. 3. Complete frow model. If all the terms on the righthand side of equation (4) are retained, the unsteady nonuniform flow equationssl Ois obtained. The Kinematic and diffusion flow models are two extreme cases of slopes (steep and flat) that are frequently
Streamflow models: I. L. Nwaogazie
encountered in the overland flow on watersheds and natural routing of a flood wave in streams. Conceivably, there are possible intermediate values of slope for which all the four slope terms in equation (4) would be appreciable. For instance, the magnitude and timing of backwater effects require all slope terms. This is a case where the complete flow model is employed.
eqlaation, equation (1) is employed to demonstrate the following: (1)derivation of explicit streamflow model, and (2) derivation of the implicit-weighted streamflow models. Inserting equation (1) into equation (9) yields:
Initial and boundary condition
Evaluating the individual terms of equation (10) using linear shape functions over the length L yields the following element equation:
The intial condition describes the depth and velocity of flow at all locations in the stream at time, t =0. This applies to all the models (simplified and complete). However, the boundary conditions refer to the depth and velocity of flow at the upstream and downstream nodes or at other specified nodes in the stream at all times, t > 0. For all models, the upstream boundary adopted was that of a prescribed discharge hydrograph: (7)
V1(t) = QI(t)/B(X)yl(t)
in which Vl(t), Q1 (t) and Yl (t) represent upstream velocity, discharge and depth of flow, respectively, at all times, t > 0. Theoretically, the kinematic model does not require any downstream boundary condition. However, for the diffusion and complete flow models, respectively, a loop rating curve was adopted. The loop rating curve as used herein implies that the relation between depth of flow and velocity is not a unique function, dependent only on the values of depth and velocity of flow but is also a function of the friction slope, s :. If s/is approximated as So - 8y/ Ox (see equation (6)), then the downstream boundary condition for a loop rating curve is given by the following expression:
l--' 1 F( vk+,= 1.486{R213"~J+
,,..
~y'~2]j + 1
LV°- x) A,<-,
DEVELOPMENT OF NUMERICAL MODELS
vz--4v 1
21"!1
2t,Z + V l ] t y , l 4/"2--u1J (Y2)
in which dot denotes differentiation with respect to time. Before the solution algorithm is initiated, element properties expressed in local coordinates, equation (11) need to be transformed into global coordinates. With the node-to-node relationship (linear discretization), it is possible to generate an overall element property matrix for the entire domain. The process is called the assembling of element equations. It is important to emphasize that the concept of discretization is based on the idea that a domain with varying geometric and hydraulic properties can be treated as a collection of sub-domains interacting. Assuming the elements are of variable lengths and that there are N nodes, the assembled global matrix equation for the continuity equation, equation (11) becomes equation (12). The expanded form of equation (12) is given in Appendix 1. [A](y) + [ B ] ( y ) - (C)= 0
(12)
in which [A] and [B] denote global metrices and (C) the global vector. The explicit streamflow model
The explicit finite element streamflow model is generated by replacing the time derivative y in equation (12) by the forward difference sheme: [a](y) j+' = [ B ] ( y ) I + A t [ B ] ( y ) J - A t ( C ) j
(13)
The right-hand side of equation (13) is known at a time level, j. Therefore, the solution of the depth of flow at various nodes proceeds forward in time with the left-hand side evaluated at the current time level, j + 1. Equation (13) can be expressed in a more compact form [A](y) j+ l = (X) j
The Galerkin weighted residual method (GWRM) is used to derive finite element equations. The weighted' residual method utilizes the concept of orthogonal projections of a residual of a differential equation onto a subspace spanned by a certain weighting function. The principle may be expressed mathematically as: 0
zJ i ,g
1V
as:
The Galerkin weighted residual method
Nke dR
1-](~;,1)
(11)
(=)
in whichj + 1 is the current time step and N, and N - 1 are the boundary node and its preceeding node, respectively. According to Freadl 3, the water surface slope is greater than the channel bottom slope, So, when the flow is increasing and less than s o when the flow is decreasing, the loop rating is thus below the single value rating on the rising limb and above it on the recession limb. If s: is determined from the momentum equation (2), terms other than So, 8y/Sx are known to affect s:," however, these terms are usually small in comparison to So, especially for typical flood routing problems and can be safely neglected.
fg
L[2
(9)
Equation (9) implies that the shape (basis) function N k must be orthogonal to the error e over the region R. For implementation of equation (9), the continuity
(14)
in which (X) denotes the sum of the right-hand side of equation (12). Because matrix [A] is of linear and tridiagonal type, a direct solution algorithm is easily applied. For equation (5) and (14) (Kinematic-wave model), the resulting computer programme uses the compact tridiagonal algorithm proposed by Varga 1s. The coded explicit finite element scheme exhibits dynamic instability unless restrictions are made on the size of the time step. This drawback is inherent in an explicit model regardless of the finite element approach. The stability problem is overcome, however, by use of the weighted-implicit streamflow model.
Adv. Water Resources, 1987, Volume 10, June
71
Streamfiow models: I. L. Nwaogazie channel having a depth of flow of 6 ft (1.8 m). The system is subjected to an upstream increase in flow to 2000 cfs (56.7 m3/s) in a period of 20 minutes, and then it decreases uniformly to the initial depth of flow in an additional period of 40 minutes. The channel has a bottom slope of 0.0015ft/ft (0.0015m/m) and an estimated Manning's Coefficient, n of 0.02. To allow comparison with the example of Viessman et al.l 4, distance steps Ax of 528 ft (162m) were used throughout the present simulations, although the proposed models have the option to accept a variable distance steps. A total of four models proposed in this paper were tested individually for this flood routing problem in a rectangular channel in order to assess their respective numerical performance in line with the following: (a) their level of approximations, (b) their explicit and implicit nature, and (c) their degree of numerical distortion with large time steps.
ill
E
-F •
-r
I I I I I L
l-l,j+i
-T
T+,--, j-;r-
l,j+i
l,J
t-l,j
l+l,j
i I I I -I
I I
<~
_L
-]
t-i,J
I
[~1-I ~--- AX
~-I,j- I
L distance (x)
q
O =~-,,t'I-
Fig. 2. Discretization of space and time domains and definition of 0
The implicit-weighted streamflow model The weighted-implicit finite element streamflow model is constructed by combining: (a) the nonlinear continuity matrix equation, equation (12); (b) the dimensionless time weighting factor, 0; and (c) the forward difference scheme for the time derivative. The resulting model yields the following relationship: [A + AtOB]{y} j+ l - OAt{C} ~+' =[A +At(1-0)B]{y}S+At(1-0){C} i 0 = At'/At
(15) (16)
in which At denotes the specified time-step and At' the adjustable time-step for the selection of 0 (see Fig. 2; 0~< At'~< At). The simultaneous solutions of equations (15) and (5) or the finite element equivalents of (6) or (2) are obtained by using the generalized functional iterative method known as the Newton-Raphson method l'6-9'a2'aa. Finally, a iinearized form of the nonlinear weighted-implicit kinematic (or diffusion or complete flow) finite element model is generated and solved by the tridiagonal 6 or bitridiagonal algorithms ~'79, respectively. For the iterative method, the guess values of the dependent variables y and v, required to initiate the solution, are the initial uniform flow depths, Y0 and velocity, %. The initial nodal depths and velocities of flow prior to the floodwave into the channel are good guess values to use.
Hydrographs for explicit and weighted-implicit models Simulated discharge hydrographs using the EKFEM for time step of 2 seconds; WIKFEM, WlDFEM and WlCFEM for time step of 60 seconds and 0= 1.00 are contrasted with the Viessman et al. 14 solution, Fig. 3. As illustrated in the figure, the hydrographs of the two kinematic-wave models, EKFEM and WIKFEM are not only identical but their peaks are much higher than those of the WIDFEM, WlCFEM and Viessman's solution. Also, evident in the results is one unique characteristics of explicit numerical modelling demonstrated by the EKFEM, that of instability with large time steps. Similar to the explicit finite difference model of Viessman et al.14, the EKFEM becomes unstable for time steps, At, larger than 2 seconds. This restriction in the use of the explicit scheme applies to both finite lement and finite difference techniques. From Fig. 3, the degree of approximations in the models is depicted from their respective predictions. For instance, the kinematic-wave by its very nature is expected to merely translate a flood wave while the diffusion-wave and the complete flow models translate as well as attenuate peak flows. Predicted hydrographs
/\
7~00
UPSTREAM4
~
'
~
. DOWNSTREAM
a,
lfiO0
1400
i
o
lO00
/
./~
L'
.
,'/'~'
\
\
V~E~.~ ~o~.~io.
TESTS IN A RECTANGULAR CHANNEL Problem description The proposed models were verified individually using a hypothetical flood in a prismatic rectangular channel as presented by Viessman et al. 14, who used an explicit finite difference scheme. The example problem considers a 2mile (3.2 kin) long and 20 feet (6.1 m) wide rectangular
72
Adv. Water Resources, 1987, Volume 10, June
800
10
20
30 /.,0 TIME , MN I UTES
i
50
Fig. 3. Discharge hydrographs using W I K F E M , W I D F E M , WICFEM and solution
i ~
70
EKFEM, V iessman's
Streamflow models: I. L. Nwaogazie confirmed these model characteristics. The WIDFEM and WlCFEM predicted peak flows about the same as those of the Viessman's solution in which the continuity and momentum equations are solved without any approximation.
- ' , ~ - - WIKFEM FOR' ,",t -"Z2-- WIKFEM FOR At
= 30 SEC. = 120 SEC.
--O-- WIKFEM FOR ,'.t = 300 SEC. -VIESSMAN SOLUTION FOR At = 2 SEE. 9
2O0O
=
1.00
UPSTREAM~ ,
/~
~ MID REACH ~,~f
Hydrographs and numerical distortion Experience has shown that hydrographs predicted using numerical models often times suffer from numerical distortion 1. Herein, numerical distortion is taken as the dispersion of hydrograph and attenuation of the peak flow. Factors controlling numerical distortion and the extent to which it occurs are the subject of this section. For effective assessment of these factors, the weightedimplicit kinematic flow model (WlKFEM) and an idealized channel geometry were selected for the following reasons: (a) Because it is a well-established fact that a kinematicwave only translates, any attenuation in the peak flow is accounted for by the factors controlling numerical distortion, and (b) The use of a rectangular channel greatly simplifies the degree of nonlinearity arising from natural geometries of river channel. Thus, the result of (a) above is easily appreciated. The WIKFEM was run for time steps, At=30, 60, 120 and 300 seconds with various values of weighting factor 8 of 0.55, 0.75 and 1.00. The resulting discharge hydrographs are compared with those of Viessman et al. 14. 1. Effects of lar#e time-steps Fig. 4 compares the predicted flows for lower limiting value of 8=0.55 and a range of time steps from small to large (30 to 300 seconds). The figure indicates that the model converges to the same solution and about the same peak flow. Evidently, the lower 0-value of 0.55 minimizes the numerical attenuation arising from large At. Thus, the kinematic-wave has merely translated the channel reach
- ' O ' - WiKFEM FOR .".! = 30"SEC. - " ~ ' - WIKFEM FOR At =120 SEC-'-C)'-- WIKFEM FOR ,".t =300 SEC. - VIPCKMAN SOLUTION FOR At = 2 SEC
8 : 0.55
/MIDREACH
2000
UpS TREAM ~ } 2 ~ / ' ~ V ~#l ~~f " L j '~f . . -
DOWNSTREA M
1800
~, f
DOWNSTREAM
1,900 P
L)
' ~
~
,
1400 u~
F= 1200
/
£,,/ /¢
tOO0
800 0
i
i
10
20
i
i
30 40 TIME , MINUTES
i
i
50
50
.Fig. 5. Discharge hydrographs using W I K F E M t = 30, 120 and 388 seconds and 8 = 1.00
70
for
of 3.2 km from upstream node to downstream node. In contrast to Fig. 4, Fig. 5 shows the discharge hydrograph using the upper limiting value of 0 -- 1.00 and At = 30, 120 and 300 seconds. Two important points are established by these illustrations. Firstly, for a large time step of 300 seconds and 8 of 1.00, numerical attenuation of 10.4 percent compared to that associated with a time step of 30 seconds is observed at the downstream node. Secondly, although the peak flows from Viessman's solution can compare favourably with those predicted by the WIKFEM, the peak flow predictions can be in error because of the combined effects of large time steps and use of an upper limiting value of 8. Thus, correct peak flow predictions are possible by combined use of lower limiting value of 0 and large time steps. Figs 4 and 5 illustrate the influence of the lower and upper limits of 8=0.55 and 1.00, respectively on numerical distortion for large time steps. Stated concisely, numerical distortion occurs because of large time steps, but is minimized at lower limiting value of 8 = 0.55 as demonstrated in Fig. 4. Although, large (Ax) distance steps accounts for numerical distortion, a constant Ax of 528ft (162m) was used throughout the runs as in Viessman's solution for ease of comparison.
o 1400
12110
8OO
i 10
/ 20
i 30
i 40
i 50
i G0
7O
TIME , MINUTES
Fig. 4. Discharge hydrographs using W I K F E M t = 30, 120 and 300 seconds and 0 = 0.55
for
2. Effects of time-weighting factor, 8 Fig. 6 compares the predicted flows for 0 = 0.55, 0.75, and 1.00 with a time step of 300 seconds. The influence of 8 on the numerical distortion is clearly evident from the figure. In fact, the predicted peak flow is 8.8 percent for 8 = 0.55 and 4.4 percent for 8=0.75 more than the Viessman solution at the downstream node. For 8=I.00, the predicted peak flow is about the same as that of the Viessman's solution. The most accurate solution in Fig. 6 is the hydrograph for 8=0.55. This is true because, a kinematic-wave translates with no peak flow attenuation. The numerical attenuation for 8 = 0.75 and 1.00 is due to a large time step being used in the computation.
Adv. Water Resources, 1987, Volume 10, June
73
Streamflow models." I. L. Nwaogazie an iterative solution algorithm with a prolonged convergence will translate into appreciable computer cpu time and cost. For a Newton-Raphson iterative method, good initial guess values for nodal depths and velocities (flow conditions before the flood wave reaches the upstream node) will accelerate convergence and subsequently reduce cost.
2OOO
t800
c~
1600
TESTS IN A N A T U R A L C H A N N E L
d < 1400 rE}
1200
1000
0
10
20
30 ~ TIME j MINUTES
5O
60
70
Fig. 6. Discharge hydrooraphs using WIKFEM for t=300 seconds and 0=0.55, 0.75 and 1.00
Table 1. Models,CPU timeand cost variationsfor rectangularchannel Time-weightingfactor, 0=0.55 Time step, At 60 seconds Models
CPU +
WIKFEM
3.59
WIDFEM
9.05
WICFEM
10.68
EKFEM
Cost -+
300 seconds CPU ÷
Cost -+
0.91 1.44 0.40 ( - 65.0%) ( - 64.9%) 2.03 2.67 0.69 ( - 21.2%) ( - 39.5%) 2.60 4.54 1.14 (o.o%) (0.0%) At = 2 seconds CPU=20.87 Cost = 5.01 (+92.7%)
Problem description A typical example of a field application of the proposed implicit models ( W I K F E M , W I D F E M and WICFEM) is given here, using the natural channel of the Illinois River in Oklahoma, USA. The Illinois river between Watts a n d Tahlequah gauging stations (sta. 1955 and 1965, respectively) shown in Fig. 7 was selected. The flood data of April 10, 1979, recorded at the upstream gauging stations, Watts, and downstream station, Tahlequah; which are 50.4 miles (80.6 km) apart, and the Flint Creek which is a tributary approximately 13.2 miles (21 km) downstream of the Watts station, were used in the flood analysis. Prior to the arrival of the flood wave at Watts station, the initial discharge values of 482 cfs (13.6 m3/s) at Watts and 596cfs (16.9 m3/s) at Tahlequah stations were recorded. In a 28 hour interval, the discharge hydrograph at the Watts station increased from 482 cfs (13.6 m3/s) to 22 980 cfs (651.0 m3/s) and then decreased to 1722cfs (49.0m3/s) in the next 68 hours. The initial discharge values at the two gauging stations and their corresponding depths of flow were used with linear interpolation to set the intermediate nodal values. These initial estimates of nodal depths of flow and discharges were used as a first approximation for initiating the Newton-Raphson iterative solution procedure. A total of 26 nodes with Ax of 2 miles (3.2 kin) for the first 24 reaches, and 2.4 miles (3.9 kin) for the last reach were employed. Other pertinent data for the simulation are the geometric cross-sections, longitudinal slope, Manning's roughness coefficient, time-weighting factor, 0, the time
+ unit of cpu time is seconds + unit of cost is dollars
4. Computer time and cost analysis for idealized channel Comparisons of the proposed models and their corresponding computer central processing unit (CPU) time and cost are presented in Table I for a given timeweighting factor, 0; time step, At; and hydraulic data. Table 1, shows that the simplified models ( W l K F E M and W I D F E M ) have less cpu time and cost as compared to the complete flow model, W l C F E M . The W l C F E M is used as the base model for cost comparison. All costs are expressed in percentages and shown in parenthesis. On the percentage cost, the negative values signify reduced cost with respect to the base cost while the positive values indicate increased cost over and above the base values. Numerical computations associated with each of the flow models are a direct function of computer cpu time and cost. Thus, Table 1 clearly illustrates the level of computations associated with each model. In fact,
74
Adv. Water Resources, 1987, Volume 10, June
TI
¢
TE RE
14
Fig. 7.
Map of Illinois River Basin
Streamflow models." I. L. Nwaogazie steps, At, and the lateral inflow. Averaged cross-sectional data (composite channel section) developed from topographic maps were used to generate the following fourth-order polynomials using a least square fitting programme 17 (a fourth-order polynomial yielded the best fit from analysis of variance): Ai= - 32.0132 + 84.6053yi + 5.4734y 2 + 0.91215y~ - 0.00453y~
16 14-
O ,.J
(17)
u_
12
14. O
Pi = 10.0824 + 57.6701y~- 4.9013y 2 + 0.34127y 3 - 0.00601y~
Watts {observedl
20
(18)
in which A~and Pi denote the average area and toP-width at node i, respectively. An increase of the order of the polynomial beyond the fourth-order was necessary to see if the increase in the degree of the polynomial significantly improved the fit of the regression. Such statistics as the sum of squares due to deviation defined as the difference between total sum of squares (SST) and sum of squares due to regression (SSR) and the goodness-of-fit defined as SSR / SSr were used for assessment. For natural rivers, equations (17) and (18), and their partial derivatives with respect to the depth of flow, y, were used for updating nodal discharge values and the nodal rate of change of the hydraulic radius. Typically, the volumetric flow rate at any node i is calculated by substituting the values of the area A in equation (3) with A i of equation (17). The initially estimated values of the Manning roughness coefficient were fitted against the discharge values for the Illinois river using the following third-order polynomial regression: n(Qi)=O.03713 +0.14097 × l O- SQi +0.41739 x 10-X°Q~-0.23004 x 10-14Q 3 (19) in which n(Qi) denotes Manning's roughness coefficient as a function of discharge, Q~, evaluated at node i. An average longitudinal slope of s0=4.5ft/mile (0.85 m/km) obtained from topographic map was used. The range of the time-weighting factor for which the implicit models are unconditionally stable is 0.55 to 1.00. The lateral inflow hydrograph at Flint Creek, a tributary of the Illinois river, was imposed as a function of time. The main channel reach corresponding to 13.2 miles (21 km) from Watts stations was allowed to receive the lateral inflow from Flint Creek. The unit of the inflow is represented as cubic feet per second per unit length of the reach within the vicinity of the confluence of Flint Creek. Predicted versus observed hydroyraphs Three weighted-implicit finite element models (WIKFEM, WIDFEM and WICFEM) were used to simulate the flood-wave on Illinois River. Each model employed a time step, At of 900 seconds and a timeweighting factor, 0 of 0.55. The predicted stage (depth) hydrographs are contrasted with the observed flow hydrograph (Fig. 8). A number of observations are evident from the figure. First, the simulated hydrographs deviate from the observed at the early portion of the rising limb and the recession limb of the hydrographs. Second, the predicted time to peak flow for the WIDFEM and WICFEM models are the same as the observed while that of the WIKFEM is about 2 hours earlier. Third, the peak
¢
~
~
(
observed)
212 I'--
o_
8
w a
0 T~lequah ( Prcdictcd, WIVFEM) F
oi
I
i
2O
,
i
4O
i
I
6O
i
~
i
!
I)0
TIME ~ HOURS
Fig. 8. Observed and simulated stage hydrographs at Tahlquah Station for t = 900 seconds and 0 = 0.55 flow predicted using WIKFEM is 6.8 ~o higher than the observed. Essentially, kinematic-wave model only translates and any observable attenuation in peak flow is accounted for by factors as presented in section 4.3. For WIDFEM and WICFEM, the predicted peak values are slightly lower than the observed. Significant in any flood predictions are the peak flow and time of peaking. In part, the time of arrival of the flood-wave at Tahlequah station if predicted correctly affords an adequate warning period for evacuation. Furthermore, the peak flow gives a measure of the level of inudation at the station and corresponding level of damage. Thus, the proposed weighted-implicit finite element models have predicted with good accuracy the two essential flood parameters (Fig. 8). For the model users, the choice of anyone of the proposed models will depend on the associated computer costs for a given flood duration, channel hydraulic and geometric properties; and model performance. The cost analysis for the field problem is presented in the section following. For model performance refer to earlier section on Hydrographs and Numerical distortion. Computer time and cost analysis for natural channel Tabulated results in Table 2 are very much in line with those of Table 1 for the rectangular channel. Contrasting both tables results in two distinct observations: Firstly, the associated cost reductions in percentages for the simplified models for the case of natural channel are slightly less than those of the rectangular channel. This is due to the additional numerical computations resulting from the increase in the nonlinearity of the channel geometric and hydraulic data from rectangular to natural type. Secondly, in either Table I or 2, increasing the time step, At by a constant factor does not cause computer cpu time and cost to decrease by the same factor. In other words, cost response to At is nonlinear and is accounted for by the following factors: (a) nonlinearity in the governing equation (b) nonlinearity in the geometric data (area, top-width, etc.) and
Adv. Water Resources, 1987, Volume 10, June
75
Streamflow models: I. L. Nwaogazie Table 2.
Models, CPU time and cost variations for natural channel Time-weighting factor. 0 = 0.550 Time step, At 900 seconds
Models WIKFEM
CPU +
Cost ±
14.76
5.80
1800 seconds CPU ÷
Cost -+
-
-
(- 70.8%) WIDFEM
WICFEM
66.7
86.50
15.86
(-23.0%) 20.57 (0.0%)
52.58
12.53
16.94
(-21.4%) 15.94 (0.0%)
+ unit of cpu time is seconds +__ unit of cost is dollars (USA)
(c) nonlinearity in hydraulic data (Manning's roughness coefficient, etc.). The advantage of the simplified models (WIKFEM and WIDFEM) in cost economy has to be weighed against the accuracy of the predicted results relative to those of the complete flow model (WlCFEM). However, some restrictions on the choice of the individual models must be highlighted. The kinematic-wave approximation has been used successfully in routing flows in rivers with steep gradient 6 and/or overland flow 18. The diffusion-wave approximation has been tried with good results for rivers with flat slopes 7. And the complete flow model encompasses the steep and fiat slopes 8. CONCLUDING REMARKS An evaluation of four deterministic streamflow models using explicit and weighted-implicit finite element methods has been carried out. The assessment was based on: numerical performance and computer storage and cost analysis, and the effects of different time steps and the time-weighting factor were investigated. Comparisons were made with previous model solutions and observed flows to enable the model user to appreciate the level of approximation associated with the individual models and their unique characteristics. Based on the results of ihis study, the following conclusions are drawn: 1.
2.
76
The predicted discharge hydrographs from the proposed models (EKFEM, WlKFEM, WlDFEM and WlCFEM) are in good agreement with those presented by Viessman et al) 4 for a flood problem in rectangular channel geometry. The two longitudinal slopes, So, of 0.0015ft/ft (0.0015m/m) and 0.00085 ft/ft (0.00085 re~m) adopted for test problem - 1 and 2 examplified approximations of steep and flat gradients for which the application of kinematicwave and diffusion-wave models are generally accepted. Also, the WICFEM performs equally well in either steep or flat slopes. The explicit kinematic finite element model (EKFEM) like the explicit finite difference scheme of Viessman becomes unstable for At larger than 2 seconds for the rectangular flood routing problem tested. This restriction (on the choice of time steps) in the use of the explicit scheme applies regardless of the numerical solution technique (finite element or difference) chosen.
Adv. Water Resources, 1987, Volume 10, June
3.
Numerical distortion (dispersion and attenuation of the computed peak flows) occurs because of large At; and is minimized at a lower limiting value of 0 = 0.55. This 0-value is considered 'conservative' and is recommended for a wide range of time and distance steps for which the weighted-implicit models remain stable. 4. The weighted-implicit models have been applied to a field problem of the Illinois River in Oklahoma for a flood observed on April 10, 1979. Predicted results using A t = 9 0 0 seconds and 0=0,55 are in close agreement with the observed stage hydrographs at Tahlequah station 50.4 miles (80.6 km) downstream. 5. The proposed models are sensitive to the variations of Mannings's roughness coefficient, lateral inflow, large time and distance steps and time-weighting factor. Proper calibration of the model is necessary for reliable predictions. For instance, raw data on nodal cross-sectional areas and top-widths were fitted with regression equations with depths of flow as well as those of Manning's roughness coefficient versus discharge. This is a necessary part of data preparation and provides essential smoothing of nonlinear data for computational ease during the modelling. 6. The level of a model's approximations are a direct function of the amount of savings in total cost per run and computer central process unit (cpu) time. The most simplified model using the implicit technique (WlKFEM) afforded the least cost with some degree of loss in accuracy. However, the performance of WIDFEM and WICFEM are identical, but costs differ. This implies that approximations of the level of diffusion or higher, better represented the complete solutions of Saint Venant equations. ACKNOWLEDGEMENTS This research was supported by the Allotment Grant No. A-103-0K at the Oklahoma State University. The author wishes to acknowledge the assistance of Mr Phil Weigant of the National Weather Service, Tulsa, Oklahoma in providing streamflow data on Illinois River. Also, sincere thanks are due to Dr A. K. Tyagi of the Department of Civil Engineering, Oklahoma State University for some valuable advice. NOTATION The following Symbols are used in this paper A B q Y x
L t
sf v So n
R
Q Nr
channel cross-sectional area width of channel water surface lateral inflow in the channel reach depth of flow distance along the channel length of an element time friction slope mean velocity across the section longitudinal bottom channel slope Manning's roughness coefficient hydraulic radius volumetric flow rate across the crosssection transpose of a basis function
Streamflow models: I. L. Nwaogazie 0 g {c} k [A] k
EKFEM WIKFEM WlDFEM WlCFEM
dimensionless time-weighting factor acceleration due to gravity column vector for nodal residual values at previous time level k linearized banded tridiagonal or bitridiagonal Jacobian matrix evaluated at time level k specified error criterion for depth explicit kinematic finite element model weighted-implicit kinematic finite element model weighted-implicit diffusion finite element model weighted-implicit complete finite element model
unpublished manuscript presented at Conference o[ the North Atlantic Division, Corps of Engineers, US Army, 1938 Nwaogazie, I. L. An approximate river flood forecasting model based on Galerkian finite element: Part 1, Int. J. Dev. Tech. 1985a, 3(1), 34-45 Nwaogazie, I. L. An approximate river flood forecasting model based on Galerkian finite element: Part 2, Int. J. Dev. Tech. 1985b, 3(2), 107-122 Nwaogazie, I. L. WlCFEM - a Fortran program for solutions of Saint-Venant Streamflow equations, Adv. Eng. Sc~(tware 1985c, 7(4), 182-198 Nwaogazie, I. L. and Tyagi, A. K. Unified streamflow routing by finite elements, Proc. ASCE, Journal of the Hydraulics Div. 1984, II0(HY7), 1595-1611 Cooley, R. L. and Moin, S. A. Finite element solution of SaintVenant equations, Proc. ASCE, Journal oJ" the Hydralics Div., 1976, 102(HY6), 759-776 Tyagi, A. K. and Nwaogazie, I. L. Testing and validation of solute transport models for groundwater systems, Tech. Report No• 15, Hydrologic Systems, Inc., Stillwater, Oklahoma, 1983 Amein, M. and Fang, C. S• Implicit flood routing in natural channels, Proc. ASCE, Journal of the Hydraulics Div. 1970 96(HY12), 2481-2500 Fread, D. L. Theoretical development of implicit dynamic routing model, presented at the Dynamic Routing Seminar, Lower Mississippi River Forecast Center, Slidell, Louisiana, 1976 Viessman, W., Jr et al. Introduction to Hydrology, 2nd edition, Harper and Row Publishers, New York, 1972 Varga, R. S. Matrix Iterative Analysis. Ist edition, Prentice-Hall, New York, 1962 Fread, D. L. Effects of time step size in implicit dynamic routing, Water Resources Bulletin, AWRA, 1973, 9(2), 339-351 Davis, J. C. Statistics and Data Analysis in Geology, 1st edition, John Wiley and Sons, New York, 1973 Ross, B. B. and Shanholtz, V. O. A one-dimensional finite element for modelling the hydrology of small upland watersheds, Proc. Hydrologic Transport Modelling Symposium, American Society of Agricultural Engineers, 1979, 42-59
6 7 8 9 10 11 12
REFERENCES 1
13
Nwaogazie, I. L. Finite element modelling ofstreamflow routing, PhD Thesis Oklahoma State University, Stillwater, Oklahoma, 1982 Barr~ de Saint-Venant. Theory of unsteady water flow with application to floods and to propagations of tides in river channels, Proc. the French Academy of Science 1871, 73, 148154, 237-240. Translated in the English by G-edding, W. W., Jr. Waterways Experiment Station, Publication No. 49-9, Vicksburg, Mississippi, 1949 Thomas, W. A. Water Surface Profiles, Hydrologic Engineering Methods for Water Resources Development, HEC, Corps of Engineers, US Army, Davis, California, 1975, 6 Puls, L. G. Flood regulation of the Tennessee river, House Document No. 185, 70th Congress, Ist session, 1928 McCarthy, G. T. The unit hydrograph and flood routing,
2
3 4 5
14 15 16 17 18
A P P E N D I X 1: E X P A N D E D F O R M OF E Q U A T I O N (12) 2L t
La
0
0
0
0
0
Yl
L~
2(Lt + L2)
L2
0
0
0
0
Y2
0
L2
2 ( L 2 + L3)
L3
0
0
0
Y3 • , ,
0
0
0
Li
2(Li+Li+l)
Li+l
0
Yi
0
0
0
0
0
LN-1
2LN_ t
YN
0
0
0
0
Yl
0
0
0
Y2
2Va + V3
0
0
Y3
2VI+2 + Vi+ 1
0
Yi
VZ--4V 1
2VZ+V l
--Vz--2v 1
l?3--v 1
2V3+V 2
0
-- U3 -- 2V2
V4 -- V2
0
0
-vi+l - 2 v i
0
0
0
+ Ui+2--U
0
i
- -
U N
- -
2vu-1
4VN -- VN-1
~ YN
qtL1 qlL1 +q2L2 q2L2 +q3L3
=o
-3
(20)
qiLi + qi+ ILl+ 1 •
.
.
q~ - 1LN - 1
Adv. Water Resources, 1987, Volume 10, June
77
Department of Civil Engineering, University of Port Harcourt, P.M.B. 5323, Port Harcourt, Rivers State, Nigeria The numerical performance and cost analysis of four streamflow models are presented. These models are constructed from the well-known equations of Saint-Venant. The Galerkin finite element method and the Newton-Raphson iterative method were utilized for the models' solution of depth and velocity of flow. Because of the large computer storage associated with complete solutions of Saint Venant equations, particularly for floods of long durations, the approximate models (kinematic- and diffusion-waves) were introduced to investigate the savings which could be made by using each of these models. Model predictions were contrasted with previously found solutions of a flood-wave in an idealized geometry. The effects of large time step (At) and the timeweighting factor, (0) for implicit models on cost and the numerical distortion of the hydrographs were examined. Results indicate that the computer cost is a direct function of the level of model approximation. For field application of the implicit models, regression equations relating: (1) Manning's roughness coefficient and discharge; (2) area and depth of flow; and (3) top-width and depth of flow were utilized. Relative to idealized channel the degree of nonlinearity for the natural channel not only affected model performance but increased the amount of computations and thus the cost. Convergence and unconditional stability were observed for implicit models for 0 in the range of 0.55 and 1.00. For the explicit model, instability restricted the choice of At. Key Words: streamflow explicit-implicit models, finite element, Newton-Raphson, timeweighting factor INTRODUCTION The process of determining the water depths, velocities and discharge in channels, rivers or reservoirs under unsteady conditions arising from flood motions is commonly referred to as flood routing ~. Interest in flood routing and in part in unsteady flows in water resources stems from the need to plan, design, regulate and manage our flood prone areas and other water resources systems. Presently, two basic techniques for unsteady flow simulation are in use, namely: (1) methods which approximate the basic equations of unsteady flow; and (2) methods which solve the basic equations directly. The basic equation as used herein refers to the complete equations of Saint Venant 2. The first methods are sometimes referred to as 'hydrologic' routing methods and the second methods 'hydraulic' routing methods 3. Basically, the hydrologic routing method neglects the momentum equation and solves the continuity equation which relates: (1) inflow into a given reach; (2) outflow from that reach; (3) time period for the flow to travel through the reach; and (4) the change in storage during the time period in the reach. A graphical relationship established between storage and outflow 4 yields a feasible solution to the continuity equation. This solution approach, though analytical is greatly simplified by the averaging technique. Accepted September 1986. Discussion closes August 1987. 0309-1708/87/020069-0952.00
© 1987 Computational Mechanics Publications
A number of hydrologic routing methods have emerged over the years including various coefficient routing techniques such as the Muskingum method 5. These methods are handy when data for hydraulic routing are not available. For detailed discussion on hydrologic routing methods, interested readers are referred to a basic text on hydrology 14. The importance of the hydraulic routing method has come to the limelight with the development of modern high-speed digital computers, allowing numerical solution of unsteady partial differential equations 2 that have no known analytical solutions to be obtained. Numerous unsteady flow phenomena such as surges, effects of tidal fluctuations, back-water resulting from channel junctions or reservoirs, and normal flood waves from excessive rainfall can be analysed using the hydraulic routing and numerical methods. In this paper, four numerical models namely: (1) explicit kinematic finite element model (EKFEM)~; (2) weighted-implicit kinematic finite element model (WlKFEM)6; (3) weighted-implicit diffusion finite element model (WIDFEM)7; and (4) weighted- implicit complete finite element model (WlCFEM) s'9 are presented. The purpose of comparing these is to enable the model-user to assess the individual model performance relative to their level of approximations and solution methodology. Given that the models are solved using the same procedure - weighted-implicit finite
Adv. Water Resources, 1987, Volume 10, June
69
Streamflow
models: I. L. Nwaogazie
method and iterative Newton-Raphson element algorithm - any observable difference in the predicted results for a typical flood routing problem may be attributed to the level of model approximation and solution technique (explicit or implicit). Often, implicit models that are unconditionally stable may conveniently employ large time- and distance-steps with the advantage of reduced computer cost*-“. On the other hand, the trade-off in the use of large time steps becomes too obvious in loss of numerical accuracy particularly the degree of attenuation of peak flow and the precise time of peaking. For instance, consider a flood wave of 3 hours duration in a channel of constant geometry, where the time to peak is to be predicted. Assuming the actual time to peak is 30 minutes, then a selected time step (At) of 600 seconds as compared to 60 seconds will not only dampen the peak flow but will yield fewer plotting points than those of At =60 seconds. In these circumstances, the precise time of peaking may be difficult to isolate. These are points of special interest in this paper. To establish confidence in the capability of any proposed numerical model, it must be validated using three levels of testing”. Level-I testing represents the model validation against analytic solutions. Level-II testing compares the model against at least an existing numerical model. And Level-III testing consists of validating the model using field data. In this paper only Levels-II and III are examined. Effects of time steps and time-weighting factors are considered along with a hypothetical flood problem in a rectangular channel. Cost and model performance are contrasted for individual models in field application and idealized channel geometry, respectively. STREAMFLOW Governing
MODELS
equations
The distribution of discharge, depth of flow, and velocity of flow in a stream are represented in Fig. 1. The mathematical models that predict the depth and velocity on a space and time basic are the continuity and momentum equations’ respectively:
~+?.;+c~-q(x,t)=o
(1)
~+~~+g~+~q(_x,t)-g(so-s,)=O (2) ’
?
Q=t’A
(3)
dx,t)
I
c
Y//T\
I /I&\
//A\
Streamjlow
AX
Fig.
1.
70
Adv. Water Resources,
4
Y
A
”
Y/A\
element
1987,
Volume 10, June
w//K
in which _Vrepresents the depth of flow, r the mean velocity of flow across the section, x the distance along the channel, and t time. q(x, t) represents the lateral inflow into the channel reach (negative for outflow), g the gravity acceleration, sI the friction slope, so the longitudinal bottom slope of the channel, Q the discharge and A the area of flow. The two dependent variables in equations (1) and (2) are the depth of flow, ~~(x,r) and the velocity of flow, t’(x, t). Equations (l), (2) and (3) can represent both an idealized channel geometry and natural rivers. For a natural channel, the geometry is specified by the area of flow, A(x); the top-width, B(x) and the longitudinal slope, so = so(x). The variables A(x) and B(x) are modelled using a fourth-order regression equation -nodal areas and topwidths as functions of depth of flow. Examples of the regression equations are illustrated in this paper in a case study of Illinois River in Oklahoma, USA. The transient values of ,4(x) and B(x) are used to update the nodal discharge values as well as the individual Jacobian terms during iterative simulation. Simplified and complete
flow models
From equations (1) and (2) three distinct streamflow models are possible by completely retaining equation (1) while simplifying assumptions are made in equation (2). By making the friction slope, s/, the subject of the formula and setting q(x, t) equal to zero, equation (2) becomes:
av c au 1 dv s~=so--L___-__ ax gsx g at Two simplified models and a complete flow model are generated using equation (4). These models are presented in the order of increasing complexity as follows: 1. Kinematic flow model. If the last three terms on the right-hand side of equation (4) are small compared to so, the following uniform flow equation6 results: n2c2 ‘/= ‘0 = 2.22 R4’3 in which n denotes Manning’s roughness coefficient and R the hydraulic radius. Equation (5) implies that gravity is balanced by resistance and when coupled to equation (1) is known as the kinematic flow model. The kinematic flow model is best suited for streamflow routing in steep gradient of the order of 10 feet per mile’ (1.9 m per km) or more and overland flows”. 2. Diffusion flow model. If the last two terms on the right-hand side of equation (4) are small compared to so, the steady nonuniform flow equation’ results:
Simultaneous solution of equations (1) and (6) give rise to what is known as the diffusion flow model, which has been successfully applied to very flat gradients’. 3. Complete frow model. If all the terms on the righthand side of equation (4) are retained, the unsteady nonuniform flow equationssl Ois obtained. The Kinematic and diffusion flow models are two extreme cases of slopes (steep and flat) that are frequently
Streamflow models: I. L. Nwaogazie
encountered in the overland flow on watersheds and natural routing of a flood wave in streams. Conceivably, there are possible intermediate values of slope for which all the four slope terms in equation (4) would be appreciable. For instance, the magnitude and timing of backwater effects require all slope terms. This is a case where the complete flow model is employed.
eqlaation, equation (1) is employed to demonstrate the following: (1)derivation of explicit streamflow model, and (2) derivation of the implicit-weighted streamflow models. Inserting equation (1) into equation (9) yields:
Initial and boundary condition
Evaluating the individual terms of equation (10) using linear shape functions over the length L yields the following element equation:
The intial condition describes the depth and velocity of flow at all locations in the stream at time, t =0. This applies to all the models (simplified and complete). However, the boundary conditions refer to the depth and velocity of flow at the upstream and downstream nodes or at other specified nodes in the stream at all times, t > 0. For all models, the upstream boundary adopted was that of a prescribed discharge hydrograph: (7)
V1(t) = QI(t)/B(X)yl(t)
in which Vl(t), Q1 (t) and Yl (t) represent upstream velocity, discharge and depth of flow, respectively, at all times, t > 0. Theoretically, the kinematic model does not require any downstream boundary condition. However, for the diffusion and complete flow models, respectively, a loop rating curve was adopted. The loop rating curve as used herein implies that the relation between depth of flow and velocity is not a unique function, dependent only on the values of depth and velocity of flow but is also a function of the friction slope, s :. If s/is approximated as So - 8y/ Ox (see equation (6)), then the downstream boundary condition for a loop rating curve is given by the following expression:
l--' 1 F( vk+,= 1.486{R213"~J+
,,..
~y'~2]j + 1
LV°- x) A,<-,
DEVELOPMENT OF NUMERICAL MODELS
vz--4v 1
21"!1
2t,Z + V l ] t y , l 4/"2--u1J (Y2)
in which dot denotes differentiation with respect to time. Before the solution algorithm is initiated, element properties expressed in local coordinates, equation (11) need to be transformed into global coordinates. With the node-to-node relationship (linear discretization), it is possible to generate an overall element property matrix for the entire domain. The process is called the assembling of element equations. It is important to emphasize that the concept of discretization is based on the idea that a domain with varying geometric and hydraulic properties can be treated as a collection of sub-domains interacting. Assuming the elements are of variable lengths and that there are N nodes, the assembled global matrix equation for the continuity equation, equation (11) becomes equation (12). The expanded form of equation (12) is given in Appendix 1. [A](y) + [ B ] ( y ) - (C)= 0
(12)
in which [A] and [B] denote global metrices and (C) the global vector. The explicit streamflow model
The explicit finite element streamflow model is generated by replacing the time derivative y in equation (12) by the forward difference sheme: [a](y) j+' = [ B ] ( y ) I + A t [ B ] ( y ) J - A t ( C ) j
(13)
The right-hand side of equation (13) is known at a time level, j. Therefore, the solution of the depth of flow at various nodes proceeds forward in time with the left-hand side evaluated at the current time level, j + 1. Equation (13) can be expressed in a more compact form [A](y) j+ l = (X) j
The Galerkin weighted residual method (GWRM) is used to derive finite element equations. The weighted' residual method utilizes the concept of orthogonal projections of a residual of a differential equation onto a subspace spanned by a certain weighting function. The principle may be expressed mathematically as: 0
zJ i ,g
1V
as:
The Galerkin weighted residual method
Nke dR
1-](~;,1)
(11)
(=)
in whichj + 1 is the current time step and N, and N - 1 are the boundary node and its preceeding node, respectively. According to Freadl 3, the water surface slope is greater than the channel bottom slope, So, when the flow is increasing and less than s o when the flow is decreasing, the loop rating is thus below the single value rating on the rising limb and above it on the recession limb. If s: is determined from the momentum equation (2), terms other than So, 8y/Sx are known to affect s:," however, these terms are usually small in comparison to So, especially for typical flood routing problems and can be safely neglected.
fg
L[2
(9)
Equation (9) implies that the shape (basis) function N k must be orthogonal to the error e over the region R. For implementation of equation (9), the continuity
(14)
in which (X) denotes the sum of the right-hand side of equation (12). Because matrix [A] is of linear and tridiagonal type, a direct solution algorithm is easily applied. For equation (5) and (14) (Kinematic-wave model), the resulting computer programme uses the compact tridiagonal algorithm proposed by Varga 1s. The coded explicit finite element scheme exhibits dynamic instability unless restrictions are made on the size of the time step. This drawback is inherent in an explicit model regardless of the finite element approach. The stability problem is overcome, however, by use of the weighted-implicit streamflow model.
Adv. Water Resources, 1987, Volume 10, June
71
Streamfiow models: I. L. Nwaogazie channel having a depth of flow of 6 ft (1.8 m). The system is subjected to an upstream increase in flow to 2000 cfs (56.7 m3/s) in a period of 20 minutes, and then it decreases uniformly to the initial depth of flow in an additional period of 40 minutes. The channel has a bottom slope of 0.0015ft/ft (0.0015m/m) and an estimated Manning's Coefficient, n of 0.02. To allow comparison with the example of Viessman et al.l 4, distance steps Ax of 528 ft (162m) were used throughout the present simulations, although the proposed models have the option to accept a variable distance steps. A total of four models proposed in this paper were tested individually for this flood routing problem in a rectangular channel in order to assess their respective numerical performance in line with the following: (a) their level of approximations, (b) their explicit and implicit nature, and (c) their degree of numerical distortion with large time steps.
ill
E
-F •
-r
I I I I I L
l-l,j+i
-T
T+,--, j-;r-
l,j+i
l,J
t-l,j
l+l,j
i I I I -I
I I
<~
_L
-]
t-i,J
I
[~1-I ~--- AX
~-I,j- I
L distance (x)
q
O =~-,,t'I-
Fig. 2. Discretization of space and time domains and definition of 0
The implicit-weighted streamflow model The weighted-implicit finite element streamflow model is constructed by combining: (a) the nonlinear continuity matrix equation, equation (12); (b) the dimensionless time weighting factor, 0; and (c) the forward difference scheme for the time derivative. The resulting model yields the following relationship: [A + AtOB]{y} j+ l - OAt{C} ~+' =[A +At(1-0)B]{y}S+At(1-0){C} i 0 = At'/At
(15) (16)
in which At denotes the specified time-step and At' the adjustable time-step for the selection of 0 (see Fig. 2; 0~< At'~< At). The simultaneous solutions of equations (15) and (5) or the finite element equivalents of (6) or (2) are obtained by using the generalized functional iterative method known as the Newton-Raphson method l'6-9'a2'aa. Finally, a iinearized form of the nonlinear weighted-implicit kinematic (or diffusion or complete flow) finite element model is generated and solved by the tridiagonal 6 or bitridiagonal algorithms ~'79, respectively. For the iterative method, the guess values of the dependent variables y and v, required to initiate the solution, are the initial uniform flow depths, Y0 and velocity, %. The initial nodal depths and velocities of flow prior to the floodwave into the channel are good guess values to use.
Hydrographs for explicit and weighted-implicit models Simulated discharge hydrographs using the EKFEM for time step of 2 seconds; WIKFEM, WlDFEM and WlCFEM for time step of 60 seconds and 0= 1.00 are contrasted with the Viessman et al. 14 solution, Fig. 3. As illustrated in the figure, the hydrographs of the two kinematic-wave models, EKFEM and WIKFEM are not only identical but their peaks are much higher than those of the WIDFEM, WlCFEM and Viessman's solution. Also, evident in the results is one unique characteristics of explicit numerical modelling demonstrated by the EKFEM, that of instability with large time steps. Similar to the explicit finite difference model of Viessman et al.14, the EKFEM becomes unstable for time steps, At, larger than 2 seconds. This restriction in the use of the explicit scheme applies to both finite lement and finite difference techniques. From Fig. 3, the degree of approximations in the models is depicted from their respective predictions. For instance, the kinematic-wave by its very nature is expected to merely translate a flood wave while the diffusion-wave and the complete flow models translate as well as attenuate peak flows. Predicted hydrographs
/\
7~00
UPSTREAM4
~
'
~
. DOWNSTREAM
a,
lfiO0
1400
i
o
lO00
/
./~
L'
.
,'/'~'
\
\
V~E~.~ ~o~.~io.
TESTS IN A RECTANGULAR CHANNEL Problem description The proposed models were verified individually using a hypothetical flood in a prismatic rectangular channel as presented by Viessman et al. 14, who used an explicit finite difference scheme. The example problem considers a 2mile (3.2 kin) long and 20 feet (6.1 m) wide rectangular
72
Adv. Water Resources, 1987, Volume 10, June
800
10
20
30 /.,0 TIME , MN I UTES
i
50
Fig. 3. Discharge hydrographs using W I K F E M , W I D F E M , WICFEM and solution
i ~
70
EKFEM, V iessman's
Streamflow models: I. L. Nwaogazie confirmed these model characteristics. The WIDFEM and WlCFEM predicted peak flows about the same as those of the Viessman's solution in which the continuity and momentum equations are solved without any approximation.
- ' , ~ - - WIKFEM FOR' ,",t -"Z2-- WIKFEM FOR At
= 30 SEC. = 120 SEC.
--O-- WIKFEM FOR ,'.t = 300 SEC. -VIESSMAN SOLUTION FOR At = 2 SEE. 9
2O0O
=
1.00
UPSTREAM~ ,
/~
~ MID REACH ~,~f
Hydrographs and numerical distortion Experience has shown that hydrographs predicted using numerical models often times suffer from numerical distortion 1. Herein, numerical distortion is taken as the dispersion of hydrograph and attenuation of the peak flow. Factors controlling numerical distortion and the extent to which it occurs are the subject of this section. For effective assessment of these factors, the weightedimplicit kinematic flow model (WlKFEM) and an idealized channel geometry were selected for the following reasons: (a) Because it is a well-established fact that a kinematicwave only translates, any attenuation in the peak flow is accounted for by the factors controlling numerical distortion, and (b) The use of a rectangular channel greatly simplifies the degree of nonlinearity arising from natural geometries of river channel. Thus, the result of (a) above is easily appreciated. The WIKFEM was run for time steps, At=30, 60, 120 and 300 seconds with various values of weighting factor 8 of 0.55, 0.75 and 1.00. The resulting discharge hydrographs are compared with those of Viessman et al. 14. 1. Effects of lar#e time-steps Fig. 4 compares the predicted flows for lower limiting value of 8=0.55 and a range of time steps from small to large (30 to 300 seconds). The figure indicates that the model converges to the same solution and about the same peak flow. Evidently, the lower 0-value of 0.55 minimizes the numerical attenuation arising from large At. Thus, the kinematic-wave has merely translated the channel reach
- ' O ' - WiKFEM FOR .".! = 30"SEC. - " ~ ' - WIKFEM FOR At =120 SEC-'-C)'-- WIKFEM FOR ,".t =300 SEC. - VIPCKMAN SOLUTION FOR At = 2 SEC
8 : 0.55
/MIDREACH
2000
UpS TREAM ~ } 2 ~ / ' ~ V ~#l ~~f " L j '~f . . -
DOWNSTREA M
1800
~, f
DOWNSTREAM
1,900 P
L)
' ~
~
,
1400 u~
F= 1200
/
£,,/ /¢
tOO0
800 0
i
i
10
20
i
i
30 40 TIME , MINUTES
i
i
50
50
.Fig. 5. Discharge hydrographs using W I K F E M t = 30, 120 and 388 seconds and 8 = 1.00
70
for
of 3.2 km from upstream node to downstream node. In contrast to Fig. 4, Fig. 5 shows the discharge hydrograph using the upper limiting value of 0 -- 1.00 and At = 30, 120 and 300 seconds. Two important points are established by these illustrations. Firstly, for a large time step of 300 seconds and 8 of 1.00, numerical attenuation of 10.4 percent compared to that associated with a time step of 30 seconds is observed at the downstream node. Secondly, although the peak flows from Viessman's solution can compare favourably with those predicted by the WIKFEM, the peak flow predictions can be in error because of the combined effects of large time steps and use of an upper limiting value of 8. Thus, correct peak flow predictions are possible by combined use of lower limiting value of 0 and large time steps. Figs 4 and 5 illustrate the influence of the lower and upper limits of 8=0.55 and 1.00, respectively on numerical distortion for large time steps. Stated concisely, numerical distortion occurs because of large time steps, but is minimized at lower limiting value of 8 = 0.55 as demonstrated in Fig. 4. Although, large (Ax) distance steps accounts for numerical distortion, a constant Ax of 528ft (162m) was used throughout the runs as in Viessman's solution for ease of comparison.
o 1400
12110
8OO
i 10
/ 20
i 30
i 40
i 50
i G0
7O
TIME , MINUTES
Fig. 4. Discharge hydrographs using W I K F E M t = 30, 120 and 300 seconds and 0 = 0.55
for
2. Effects of time-weighting factor, 8 Fig. 6 compares the predicted flows for 0 = 0.55, 0.75, and 1.00 with a time step of 300 seconds. The influence of 8 on the numerical distortion is clearly evident from the figure. In fact, the predicted peak flow is 8.8 percent for 8 = 0.55 and 4.4 percent for 8=0.75 more than the Viessman solution at the downstream node. For 8=I.00, the predicted peak flow is about the same as that of the Viessman's solution. The most accurate solution in Fig. 6 is the hydrograph for 8=0.55. This is true because, a kinematic-wave translates with no peak flow attenuation. The numerical attenuation for 8 = 0.75 and 1.00 is due to a large time step being used in the computation.
Adv. Water Resources, 1987, Volume 10, June
73
Streamflow models." I. L. Nwaogazie an iterative solution algorithm with a prolonged convergence will translate into appreciable computer cpu time and cost. For a Newton-Raphson iterative method, good initial guess values for nodal depths and velocities (flow conditions before the flood wave reaches the upstream node) will accelerate convergence and subsequently reduce cost.
2OOO
t800
c~
1600
TESTS IN A N A T U R A L C H A N N E L
d < 1400 rE}
1200
1000
0
10
20
30 ~ TIME j MINUTES
5O
60
70
Fig. 6. Discharge hydrooraphs using WIKFEM for t=300 seconds and 0=0.55, 0.75 and 1.00
Table 1. Models,CPU timeand cost variationsfor rectangularchannel Time-weightingfactor, 0=0.55 Time step, At 60 seconds Models
CPU +
WIKFEM
3.59
WIDFEM
9.05
WICFEM
10.68
EKFEM
Cost -+
300 seconds CPU ÷
Cost -+
0.91 1.44 0.40 ( - 65.0%) ( - 64.9%) 2.03 2.67 0.69 ( - 21.2%) ( - 39.5%) 2.60 4.54 1.14 (o.o%) (0.0%) At = 2 seconds CPU=20.87 Cost = 5.01 (+92.7%)
Problem description A typical example of a field application of the proposed implicit models ( W I K F E M , W I D F E M and WICFEM) is given here, using the natural channel of the Illinois River in Oklahoma, USA. The Illinois river between Watts a n d Tahlequah gauging stations (sta. 1955 and 1965, respectively) shown in Fig. 7 was selected. The flood data of April 10, 1979, recorded at the upstream gauging stations, Watts, and downstream station, Tahlequah; which are 50.4 miles (80.6 km) apart, and the Flint Creek which is a tributary approximately 13.2 miles (21 km) downstream of the Watts station, were used in the flood analysis. Prior to the arrival of the flood wave at Watts station, the initial discharge values of 482 cfs (13.6 m3/s) at Watts and 596cfs (16.9 m3/s) at Tahlequah stations were recorded. In a 28 hour interval, the discharge hydrograph at the Watts station increased from 482 cfs (13.6 m3/s) to 22 980 cfs (651.0 m3/s) and then decreased to 1722cfs (49.0m3/s) in the next 68 hours. The initial discharge values at the two gauging stations and their corresponding depths of flow were used with linear interpolation to set the intermediate nodal values. These initial estimates of nodal depths of flow and discharges were used as a first approximation for initiating the Newton-Raphson iterative solution procedure. A total of 26 nodes with Ax of 2 miles (3.2 kin) for the first 24 reaches, and 2.4 miles (3.9 kin) for the last reach were employed. Other pertinent data for the simulation are the geometric cross-sections, longitudinal slope, Manning's roughness coefficient, time-weighting factor, 0, the time
+ unit of cpu time is seconds + unit of cost is dollars
4. Computer time and cost analysis for idealized channel Comparisons of the proposed models and their corresponding computer central processing unit (CPU) time and cost are presented in Table I for a given timeweighting factor, 0; time step, At; and hydraulic data. Table 1, shows that the simplified models ( W l K F E M and W I D F E M ) have less cpu time and cost as compared to the complete flow model, W l C F E M . The W l C F E M is used as the base model for cost comparison. All costs are expressed in percentages and shown in parenthesis. On the percentage cost, the negative values signify reduced cost with respect to the base cost while the positive values indicate increased cost over and above the base values. Numerical computations associated with each of the flow models are a direct function of computer cpu time and cost. Thus, Table 1 clearly illustrates the level of computations associated with each model. In fact,
74
Adv. Water Resources, 1987, Volume 10, June
TI
¢
TE RE
14
Fig. 7.
Map of Illinois River Basin
Streamflow models." I. L. Nwaogazie steps, At, and the lateral inflow. Averaged cross-sectional data (composite channel section) developed from topographic maps were used to generate the following fourth-order polynomials using a least square fitting programme 17 (a fourth-order polynomial yielded the best fit from analysis of variance): Ai= - 32.0132 + 84.6053yi + 5.4734y 2 + 0.91215y~ - 0.00453y~
16 14-
O ,.J
(17)
u_
12
14. O
Pi = 10.0824 + 57.6701y~- 4.9013y 2 + 0.34127y 3 - 0.00601y~
Watts {observedl
20
(18)
in which A~and Pi denote the average area and toP-width at node i, respectively. An increase of the order of the polynomial beyond the fourth-order was necessary to see if the increase in the degree of the polynomial significantly improved the fit of the regression. Such statistics as the sum of squares due to deviation defined as the difference between total sum of squares (SST) and sum of squares due to regression (SSR) and the goodness-of-fit defined as SSR / SSr were used for assessment. For natural rivers, equations (17) and (18), and their partial derivatives with respect to the depth of flow, y, were used for updating nodal discharge values and the nodal rate of change of the hydraulic radius. Typically, the volumetric flow rate at any node i is calculated by substituting the values of the area A in equation (3) with A i of equation (17). The initially estimated values of the Manning roughness coefficient were fitted against the discharge values for the Illinois river using the following third-order polynomial regression: n(Qi)=O.03713 +0.14097 × l O- SQi +0.41739 x 10-X°Q~-0.23004 x 10-14Q 3 (19) in which n(Qi) denotes Manning's roughness coefficient as a function of discharge, Q~, evaluated at node i. An average longitudinal slope of s0=4.5ft/mile (0.85 m/km) obtained from topographic map was used. The range of the time-weighting factor for which the implicit models are unconditionally stable is 0.55 to 1.00. The lateral inflow hydrograph at Flint Creek, a tributary of the Illinois river, was imposed as a function of time. The main channel reach corresponding to 13.2 miles (21 km) from Watts stations was allowed to receive the lateral inflow from Flint Creek. The unit of the inflow is represented as cubic feet per second per unit length of the reach within the vicinity of the confluence of Flint Creek. Predicted versus observed hydroyraphs Three weighted-implicit finite element models (WIKFEM, WIDFEM and WICFEM) were used to simulate the flood-wave on Illinois River. Each model employed a time step, At of 900 seconds and a timeweighting factor, 0 of 0.55. The predicted stage (depth) hydrographs are contrasted with the observed flow hydrograph (Fig. 8). A number of observations are evident from the figure. First, the simulated hydrographs deviate from the observed at the early portion of the rising limb and the recession limb of the hydrographs. Second, the predicted time to peak flow for the WIDFEM and WICFEM models are the same as the observed while that of the WIKFEM is about 2 hours earlier. Third, the peak
¢
~
~
(
observed)
212 I'--
o_
8
w a
0 T~lequah ( Prcdictcd, WIVFEM) F
oi
I
i
2O
,
i
4O
i
I
6O
i
~
i
!
I)0
TIME ~ HOURS
Fig. 8. Observed and simulated stage hydrographs at Tahlquah Station for t = 900 seconds and 0 = 0.55 flow predicted using WIKFEM is 6.8 ~o higher than the observed. Essentially, kinematic-wave model only translates and any observable attenuation in peak flow is accounted for by factors as presented in section 4.3. For WIDFEM and WICFEM, the predicted peak values are slightly lower than the observed. Significant in any flood predictions are the peak flow and time of peaking. In part, the time of arrival of the flood-wave at Tahlequah station if predicted correctly affords an adequate warning period for evacuation. Furthermore, the peak flow gives a measure of the level of inudation at the station and corresponding level of damage. Thus, the proposed weighted-implicit finite element models have predicted with good accuracy the two essential flood parameters (Fig. 8). For the model users, the choice of anyone of the proposed models will depend on the associated computer costs for a given flood duration, channel hydraulic and geometric properties; and model performance. The cost analysis for the field problem is presented in the section following. For model performance refer to earlier section on Hydrographs and Numerical distortion. Computer time and cost analysis for natural channel Tabulated results in Table 2 are very much in line with those of Table 1 for the rectangular channel. Contrasting both tables results in two distinct observations: Firstly, the associated cost reductions in percentages for the simplified models for the case of natural channel are slightly less than those of the rectangular channel. This is due to the additional numerical computations resulting from the increase in the nonlinearity of the channel geometric and hydraulic data from rectangular to natural type. Secondly, in either Table I or 2, increasing the time step, At by a constant factor does not cause computer cpu time and cost to decrease by the same factor. In other words, cost response to At is nonlinear and is accounted for by the following factors: (a) nonlinearity in the governing equation (b) nonlinearity in the geometric data (area, top-width, etc.) and
Adv. Water Resources, 1987, Volume 10, June
75
Streamflow models: I. L. Nwaogazie Table 2.
Models, CPU time and cost variations for natural channel Time-weighting factor. 0 = 0.550 Time step, At 900 seconds
Models WIKFEM
CPU +
Cost ±
14.76
5.80
1800 seconds CPU ÷
Cost -+
-
-
(- 70.8%) WIDFEM
WICFEM
66.7
86.50
15.86
(-23.0%) 20.57 (0.0%)
52.58
12.53
16.94
(-21.4%) 15.94 (0.0%)
+ unit of cpu time is seconds +__ unit of cost is dollars (USA)
(c) nonlinearity in hydraulic data (Manning's roughness coefficient, etc.). The advantage of the simplified models (WIKFEM and WIDFEM) in cost economy has to be weighed against the accuracy of the predicted results relative to those of the complete flow model (WlCFEM). However, some restrictions on the choice of the individual models must be highlighted. The kinematic-wave approximation has been used successfully in routing flows in rivers with steep gradient 6 and/or overland flow 18. The diffusion-wave approximation has been tried with good results for rivers with flat slopes 7. And the complete flow model encompasses the steep and fiat slopes 8. CONCLUDING REMARKS An evaluation of four deterministic streamflow models using explicit and weighted-implicit finite element methods has been carried out. The assessment was based on: numerical performance and computer storage and cost analysis, and the effects of different time steps and the time-weighting factor were investigated. Comparisons were made with previous model solutions and observed flows to enable the model user to appreciate the level of approximation associated with the individual models and their unique characteristics. Based on the results of ihis study, the following conclusions are drawn: 1.
2.
76
The predicted discharge hydrographs from the proposed models (EKFEM, WlKFEM, WlDFEM and WlCFEM) are in good agreement with those presented by Viessman et al) 4 for a flood problem in rectangular channel geometry. The two longitudinal slopes, So, of 0.0015ft/ft (0.0015m/m) and 0.00085 ft/ft (0.00085 re~m) adopted for test problem - 1 and 2 examplified approximations of steep and flat gradients for which the application of kinematicwave and diffusion-wave models are generally accepted. Also, the WICFEM performs equally well in either steep or flat slopes. The explicit kinematic finite element model (EKFEM) like the explicit finite difference scheme of Viessman becomes unstable for At larger than 2 seconds for the rectangular flood routing problem tested. This restriction (on the choice of time steps) in the use of the explicit scheme applies regardless of the numerical solution technique (finite element or difference) chosen.
Adv. Water Resources, 1987, Volume 10, June
3.
Numerical distortion (dispersion and attenuation of the computed peak flows) occurs because of large At; and is minimized at a lower limiting value of 0 = 0.55. This 0-value is considered 'conservative' and is recommended for a wide range of time and distance steps for which the weighted-implicit models remain stable. 4. The weighted-implicit models have been applied to a field problem of the Illinois River in Oklahoma for a flood observed on April 10, 1979. Predicted results using A t = 9 0 0 seconds and 0=0,55 are in close agreement with the observed stage hydrographs at Tahlequah station 50.4 miles (80.6 km) downstream. 5. The proposed models are sensitive to the variations of Mannings's roughness coefficient, lateral inflow, large time and distance steps and time-weighting factor. Proper calibration of the model is necessary for reliable predictions. For instance, raw data on nodal cross-sectional areas and top-widths were fitted with regression equations with depths of flow as well as those of Manning's roughness coefficient versus discharge. This is a necessary part of data preparation and provides essential smoothing of nonlinear data for computational ease during the modelling. 6. The level of a model's approximations are a direct function of the amount of savings in total cost per run and computer central process unit (cpu) time. The most simplified model using the implicit technique (WlKFEM) afforded the least cost with some degree of loss in accuracy. However, the performance of WIDFEM and WICFEM are identical, but costs differ. This implies that approximations of the level of diffusion or higher, better represented the complete solutions of Saint Venant equations. ACKNOWLEDGEMENTS This research was supported by the Allotment Grant No. A-103-0K at the Oklahoma State University. The author wishes to acknowledge the assistance of Mr Phil Weigant of the National Weather Service, Tulsa, Oklahoma in providing streamflow data on Illinois River. Also, sincere thanks are due to Dr A. K. Tyagi of the Department of Civil Engineering, Oklahoma State University for some valuable advice. NOTATION The following Symbols are used in this paper A B q Y x
L t
sf v So n
R
Q Nr
channel cross-sectional area width of channel water surface lateral inflow in the channel reach depth of flow distance along the channel length of an element time friction slope mean velocity across the section longitudinal bottom channel slope Manning's roughness coefficient hydraulic radius volumetric flow rate across the crosssection transpose of a basis function
Streamflow models: I. L. Nwaogazie 0 g {c} k [A] k
EKFEM WIKFEM WlDFEM WlCFEM
dimensionless time-weighting factor acceleration due to gravity column vector for nodal residual values at previous time level k linearized banded tridiagonal or bitridiagonal Jacobian matrix evaluated at time level k specified error criterion for depth explicit kinematic finite element model weighted-implicit kinematic finite element model weighted-implicit diffusion finite element model weighted-implicit complete finite element model
unpublished manuscript presented at Conference o[ the North Atlantic Division, Corps of Engineers, US Army, 1938 Nwaogazie, I. L. An approximate river flood forecasting model based on Galerkian finite element: Part 1, Int. J. Dev. Tech. 1985a, 3(1), 34-45 Nwaogazie, I. L. An approximate river flood forecasting model based on Galerkian finite element: Part 2, Int. J. Dev. Tech. 1985b, 3(2), 107-122 Nwaogazie, I. L. WlCFEM - a Fortran program for solutions of Saint-Venant Streamflow equations, Adv. Eng. Sc~(tware 1985c, 7(4), 182-198 Nwaogazie, I. L. and Tyagi, A. K. Unified streamflow routing by finite elements, Proc. ASCE, Journal of the Hydraulics Div. 1984, II0(HY7), 1595-1611 Cooley, R. L. and Moin, S. A. Finite element solution of SaintVenant equations, Proc. ASCE, Journal oJ" the Hydralics Div., 1976, 102(HY6), 759-776 Tyagi, A. K. and Nwaogazie, I. L. Testing and validation of solute transport models for groundwater systems, Tech. Report No• 15, Hydrologic Systems, Inc., Stillwater, Oklahoma, 1983 Amein, M. and Fang, C. S• Implicit flood routing in natural channels, Proc. ASCE, Journal of the Hydraulics Div. 1970 96(HY12), 2481-2500 Fread, D. L. Theoretical development of implicit dynamic routing model, presented at the Dynamic Routing Seminar, Lower Mississippi River Forecast Center, Slidell, Louisiana, 1976 Viessman, W., Jr et al. Introduction to Hydrology, 2nd edition, Harper and Row Publishers, New York, 1972 Varga, R. S. Matrix Iterative Analysis. Ist edition, Prentice-Hall, New York, 1962 Fread, D. L. Effects of time step size in implicit dynamic routing, Water Resources Bulletin, AWRA, 1973, 9(2), 339-351 Davis, J. C. Statistics and Data Analysis in Geology, 1st edition, John Wiley and Sons, New York, 1973 Ross, B. B. and Shanholtz, V. O. A one-dimensional finite element for modelling the hydrology of small upland watersheds, Proc. Hydrologic Transport Modelling Symposium, American Society of Agricultural Engineers, 1979, 42-59
6 7 8 9 10 11 12
REFERENCES 1
13
Nwaogazie, I. L. Finite element modelling ofstreamflow routing, PhD Thesis Oklahoma State University, Stillwater, Oklahoma, 1982 Barr~ de Saint-Venant. Theory of unsteady water flow with application to floods and to propagations of tides in river channels, Proc. the French Academy of Science 1871, 73, 148154, 237-240. Translated in the English by G-edding, W. W., Jr. Waterways Experiment Station, Publication No. 49-9, Vicksburg, Mississippi, 1949 Thomas, W. A. Water Surface Profiles, Hydrologic Engineering Methods for Water Resources Development, HEC, Corps of Engineers, US Army, Davis, California, 1975, 6 Puls, L. G. Flood regulation of the Tennessee river, House Document No. 185, 70th Congress, Ist session, 1928 McCarthy, G. T. The unit hydrograph and flood routing,
2
3 4 5
14 15 16 17 18
A P P E N D I X 1: E X P A N D E D F O R M OF E Q U A T I O N (12) 2L t
La
0
0
0
0
0
Yl
L~
2(Lt + L2)
L2
0
0
0
0
Y2
0
L2
2 ( L 2 + L3)
L3
0
0
0
Y3 • , ,
0
0
0
Li
2(Li+Li+l)
Li+l
0
Yi
0
0
0
0
0
LN-1
2LN_ t
YN
0
0
0
0
Yl
0
0
0
Y2
2Va + V3
0
0
Y3
2VI+2 + Vi+ 1
0
Yi
VZ--4V 1
2VZ+V l
--Vz--2v 1
l?3--v 1
2V3+V 2
0
-- U3 -- 2V2
V4 -- V2
0
0
-vi+l - 2 v i
0
0
0
+ Ui+2--U
0
i
- -
U N
- -
2vu-1
4VN -- VN-1
~ YN
qtL1 qlL1 +q2L2 q2L2 +q3L3
=o
-3
(20)
qiLi + qi+ ILl+ 1 •
.
.
q~ - 1LN - 1
Adv. Water Resources, 1987, Volume 10, June
77