COMPARISON OF FLOOD ROUTING METHODSbyEDY ANTO SOENTOROSarjana, Bandung Institute of Technology, IndonesiaA THESIS SUBMITTED IN PARTIAL FULFILMENT OFTHE REQUIREMENT FOR THE DEGREE OFMASTER OF APPLIED SCIENCEinTHE FACULTY OF GRADUATE STUDIESDepartment of Civil EngineeringWe accept this thesis as conformingto the required standardTHE UNIVERSITY OF BRITISH COLUMBIANovember, 1991© Edy Anto Soentoro, 1991In presenting this thesis in partial fulfilment of the requirements for an advanceddegree at the University of British Columbia, I agree that the Library shall make itfreely available for reference and study. I further agree that permission for extensivecopying of this thesis for scholarly purposes may be granted by the head of mydepartment or by his or her representatives. It is understood that copying orpublication of this thesis for financial gain shall not be allowed without my writtenpermission.(Signature) Department of C--■^X L- t^EThe University of British ColumbiaVancouver, Canadat^DateDE-6 (2/88)AbstractRecent developments in flood routing have resulted in several numerical modelswith different features. These models produce different results, and the discrepanciesbetween computed and observed flows vary depending on the values of the channel'sfriction coefficient and bed slope. The objective of this study is to find the mostreliable model for a particular range of combinations of channel friction coefficientand bed slope.Results from five flood routing models -- the dynamic wave, the characteristic,the kinematic wave, Muskingum-Cunge and UBC Flow models -- were compared.The first three models are based on a hydraulic approach, but they solve the unsteadyflow equations with different methods. The fourth model is based on the hydrologicalapproach, and the last one is a hybrid model. Whith some modifications, these fivedifferent models were run using the same input parameters, and their results werecompared.The comparison consists of two steps. The first involved applying the models toa set of actual Fraser River flows. The second involved comparisons of the resultswhen the models were applied to an artificial channel with various friction coefficientsand bed slopes. The latter procedure was used because of the difficulties in obtainingsufficient data from natural rivers. The computed Fraser River flows from all theabove models had a good agreement to the observed outflows. However, the mostiiAbstractaccurate result was that from the dynamic wave method, followed by those from thecharacteristic, the kinematic wave, the UBC Flow and Muskingum-Cunge methodsrespectively.There are two main results from this study. First, the most reliable model is thedynamic wave method because of its accuracy and applicability, and second, thediscrepancies between computed and actual flows decrease with smaller frictioncoefficient or steeper bed slope channel.iiiTable of ContentsAbstract ^ iiTable of Contents^ ivList of Figures viList of Symbols^ viiiAcknowledgement xi1. Introduction^ 12. Literature Review2.1. Background^ 42.2. Dynamic Wave method^ 112.3. Characteristic method 122.4. Kinematic Wave method^ 152.5. Muskingum-Cunge method 182.6. UBC Flow model^ 213. The Methods3.1. Dynamic Wave method^ 233.1.1. Boundary and Initial Conditions ^ 263.1.2. The Generalized Newton Raphson Iteration Method^ 283.1.3. Stability Consideration ^ 313.2. Characteristic method 32ivTable of Contents3.2.1. Boundary Conditions ^ 353.3. Kinematic Wave method 373.4. Muskingum-Cunge method^ 393.5. UBC Flow model^ 424. Comparison of the Methods4.1. Comparison Procedure^ 464.2. Results and Discussion 525. Conclusions^ 67Bibliography 71AppendicesA. Flowcharts of the methods^ 77B. Comparison of the computed outflows and their discrepancies^ 88C. Jacobian matrix for the dynamic wave method and example data forthe UBC Flow model (area-discharge and velocity-discharge relationship)^ 103vList of Figures2.1. Energy in gradually varied unsteady flow^ 62.2. Range of influence and domain^ 133.1. Four point weighted implicit finite difference scheme^ 233.2. Explicit finite difference scheme^ 343.3. Finite difference scheme for the kinematic wave method^ 383.4. Flood wave translation diagram^ 434.1. Location of Texas creek and Marguerite hydrometric sta. ^ 484.2. Fraser river hydrographs (May - July 1988) ^ 494.3. Computed outflow hydrographs, Fraser river at Texas creek^ 554.4. Computed outflow hydrographs using the dynamic wave method^ 564.5. Comparison with various Manning's n values ^ 594.6. Comparison of discrepancies with So =0.0010 and n =0.015 ^ 604.7. Comparison with various bed slopes ^ 614.8. Comparison of discrepancies with So =0.0005 and n =0.035 ^ 624.9. Comparison of flood routing discrepancies ^ 656.1. Comparison with various Manning's n values (rt = 0.015 - 0.055) ^ 896.2. Comparison of discrepancies with S o=0.0010 and n =0.015-0.055 ^ 926.3. Comparison with various bed slopes (S0 = 0.0005 - 0.0050) ^ 976.4. Comparison of discrepancies with S o =0.0005-0.0050 and n =0.035 ^ 99viList of Figures6.5. Example data for the UBC Flow model(area -discharge and velocity -discharge relationship) ^ 106viilist of Symbols- a coefficient, for the kinematic wave method- wetted cross-sectional area of the channelcelerity of the flood wavecelerity at points L and R, for the characteristic methodpositive and negative characteristic linesC3^parameters for the Muskingum-Cunge methoddiscrepancy indicatora general function with two independent variablesaACCL, CRC 1 C2C 1 C2 ,DIF ( u,1711(x) F. (Xi)^matrix of F (x) , each row may have different functions, and eachcolumn may also have different independent variables (e.g., x i, x2)g acceleration due to gravityinflow discharge, for the Muskingum-Cunge methodK a general function, for implicit finite difference scheme, expressingthe influence of conditions at the previous time step (or distancestep) to conditions at the next time step (or distance step)K storage parameter, for the UBC Flow and Muskingum-CungemethodKT^K /A tn Manning's friction coefficiento outflow discharge, for the Muskingum-Cunge methodPp^wetted perimeter of the channel at point PviiiList of Symbolsq^- lateral inflowOr^inflow discharge, for the UBC Flow modelQo^outflow discharge, for the UBC Flow modelQOT^translated flow discharge, for the UBC Flow modelQw^flow discharge in the channel per unit width, for the Muskingum-Cunge methodQii+i^discharge at point (i + 1) and time step jRp^hydraulic radius of the channel at point PRTK^correction factor for kinematic wave travel time, for the UBC Flowmethodabsolute storage within the reachSi^friction slope, depends on channel's geometry, friction and dischargef R^friction slope at point R, for the characteristic methodSo^channel's bed slopeSTK^correction factor for storage parameter, for the UBC Flow model7'^top width of water in the channelTT^travel time, for the UBC Flow modelu^velocity of flowUL, UR^velocity at points L and R, for the characteristic methoduij+1^velocity at point (i + 1) and time step jt^timeixList of Symbolsx^weighting factor, for the Muskingum-Cunge methodx^distance, through the longitudinal axis of the channelxN^distance from the upstream boundary to point N, at the downstreamboundary0^parameter, for the kinematic wave methodAX^length of sub reach, for computation purposeA t^- time intervalA^- parameter, for the characteristic methode^- Preissmann's weighting factor, for the dynamic wave methodxAcknowledgementThe author would like to express his sincere appreciation and gratitude to histhesis supervisor, Prof. S.O. Russell, for his advice and guidance throughout theresearch work and in the preparation of this thesis. Thanks are also extended toProf. M.C. Quick for explaining the UBC Flow model, and to Prof. W.F. Caseltonfor checking the manuscript of this thesis.xiChapter 1INTRODUCTIONFlood routing -- a tracing of the dynamic movement of a flood wave movingdown a channel -- is one of the most interesting and challenging of unsteady flowphenomena. It is also one of the most important subjects in hydraulics and hydrology.With an understanding of the physical processes and behaviour of a flood wave in achannel or a natural river, engineers can predict its effects downstream, when andhow high the water level will be. The variation of water levels in a channel is themain input for choosing optimal actions to take in order to protect life and propertiesin the area affected by floods (e.g., any action related to channel improvement, earlywarning system and reservoir operation for flood control). Also, an understandingof flood wave movement could be applied in other areas. For example, one coulddesign a tidal irrigation system by taking advantage of tides for flushing and wateringpadi (rice) field in swamp-land area. If the fluctuation of water levels at a givenpoint can be predicted, that would provide the main input for tidal irrigation as wellas for river management.Mathematical models of flood routing may broadly be classified into two groups-- hydraulic and hydrological models. Hydraulic models are based on the physicalprocesses of water movement in a channel, natural stream or overland. The laws ofconservation of mass and momentum are applied to elements of water. Examples1Chapter 1. Introductionof hydraulic models are the characteristic method, dynamic wave, kinematic wavemethod, finite element, quasi steady and diffusion method. On the other hand,hydrological models (e.g., Muskingum method, Storage Routing method) use asimplified system approach, using only the law of conservation of mass to define therelationship between storage, inflow and outflow within a reach.Since there are many flood routing models, questions have emerged about theapplicability of these models to channels with varying friction factors (Manning's n)and bed slope (S0). Which one of them gives the most accurate result when it isapplied to a particular combination of Manning's n and S0 ? The objective of thisstudy is to find the answers to these questions. For this purpose, five methods arecompared. Ideally, observed data both at the upstream and downstream boundariesare needed as well as So and Manning's n for each channel. However, because of thedifficulties in obtaining enough data from natural rivers, the main comparison is donefor an artificial prismatic channel but a check was also made with the Fraser River.Flows in the Fraser River (near Marguerite) in 1988 were used as inflows. Computedoutflows were compared with one another, and also with recorded flows from thedownstream gauge at Texas Creek.Five models are programmed and the result are compared here. Three of themare based on the hydraulic approach (i.e., dynamic wave, characteristic and kinematicwave method), one model is based on the hydrological approach, and the last one is2Chapter 1. Introductiona hybrid model taken from the UBC Flow model. Each of these models has differentfeatures, and originally needs different parameters for input data. Parameters ofinput data of the models that are based on the hydrological approach usually do notrepresent directly changes in the channel's bed slopes and friction coefficients. Someof the models also require data from both the upstream and the downstreamboundaries (e.g., the dynamic wave method), while others are just from the upstreamboundary. With some modifications, however, it is possible to apply them by usingthe same input parameters and compare their results for the purpose of this study.Background and development of mathematical models of flood routing are givenin Chapter 2. A review of previous studies is also discussed in that chapter.Chapter 3 provides more detailed explanation about the five chosen methods (i.e.,the dynamic wave, the characteristic, the kinematic wave, Muskingum-Cunge, andUBC Flow methods), the equations, and the boundary conditions applied in thisstudy. The comparison procedures, the results and discussion are reported inChapter 4, followed by summary and conclusions in Chapter 5. Flowcharts of themodels are given in appendix A. Appendix B provides the figures of comparisons ina prismatic channel with various friction coefficients and bed slopes.3Chapter 2LITERATURE REVIEW2.1. BACKGROUNDFlood routing has been an object of investigation for centuries, sincescientists wanted to predict the characteristic features of flood wave travel downa river. In earlier time, the main purpose of such investigation was to determinenecessary action to protect life and properties from the effects of flooding (e.g.,early warning, improvement of natural or man-made waterways). Later, with thedevelopment of more accurate predictions, scientists can provide useful data forflood control, tidal irrigation (taking advantage of tidal propagation in a river)and river management.The investigation of flood routing began by observing empirical relationshipsbetween inflow and outflow within a reach (e.g., gauge relationship, lag model).Since then, mathematical techniques have been continually developed to givemore accurate predictions of flood wave movement. The development can bedivided into two group models, based on the hydrological approach and thehydraulic approach.Hydrological models are based on the relationships between inflow, outflow,rate of change of the storage within the reach, and time lag (time to peak).Some example of these models are Muskingum method, Storage Routing method4Chapter 2. Literature Review(Henderson, 1966), Modified Puls method (Chow,1988), Working R & D andLevel Pool Reservoir Routing (Corps of Engr., 1987). In the Muskingum method,there are many variations of the method such as those which were presented byNash (1959), Kulandaiswamy (1966), Diskin (1967), Overton (1966), Cunge(1969, in: Koussis, 1980), Gill (1977, 1978, 1979), Ponce & Yevjevich (1978b),Tung (1985), and Singh (1987).Hydrological models are simpler than hydraulic models, but they have somedisadvantages. Generally, they are limited in application because they needobserved inflow and outflow hydrographs from a reach to determine the routingcoefficients. Back-water effects from tides, significant tributary inflow, dams orbridges are not considered.Hydraulic models (e.g., dynamic wave, kinematic wave method, etc.) arebased on the physical processes of water movement in a channel, natural streamor over land. With contribution from Saint Venant in the 19th century, the basictheory of one-dimensional analysis of flood wave propagation was formulated.However, the application of these models are very limited, because of thecomplexity of the equations (the Saint Venant equations). Significant progressappeared only after if become possible to use a high speed computer todetermine the solution of the Saint Venant equations. Additionally, aconsiderable amount of data is needed to describe the channel geometry.5Ifriction slope Sy IIdePth water slope 8vv+ ax 6.3(;^aybeclah,^Y+^ox''Sojz.29---t -29Chapter 2 Literature ReviewThe derivation of the Saint Venant equations from conservation of mass andmomentum has been given in standard reference books (e.g., Henderson, 1966,Abbot,1979). The result is presented here, including a brief description.Consider a short reach of channel length AX with the flow taking place fromsection 1 to section 2 as shown in Fig.2.1.ElevationWater surface at different timeDistanceFig. 2.1 Energy in gradually varied unsteady flow.Conservation of mass and momentum leads to the continuity equation,a_y + A au^u ay - q 0 ^ (2.1)at T ax^ax Tand the momentum equation,6Chapter 2. Literature Reviewatay + u aau g[ aY^(sf_so)] +^- 0 ^ (2.2)in which :y = depth of flowu = velocity of flowt = timex = distance, through the longitudinal axis of the channelA = wetted cross-sectional areaT = top width of water in the channelst = friction slope, depends on channel's geometry, friction, and dischargeSo = bed slope = Ax/ Ayq = lateral inflowg = acceleration due to gravityTo show how non-uniformity and unsteadiness introduce extra terms into theequations, eq. 2.2 can be rearranged as :_ qu _ ay _ u au _ 1 augA ax g ax g atSteady uniform flow'Steady non uniform flow^Unsteady non uniform flow^Equation 2.1 and eq. 2.2 or eq. 2.3 are partial differential equations of thehyperbolic type. These equations are not easy to solve. Any analytical solutionshould be based on reducing these equations into simpler forms. However, withthis procedure, some important features can be lost and the resulting solutionsmay lack generality. The solution may valid in a particular application, but inother cases it may not valid or it may difficult to find the solution.St So (2.3)7Chapter 2. Literature ReviewOne of the simplest ways to simplify the unsteady flow (Saint Venant)equations is what it is called the monoclinical wave approach. A single wave isassumed to travel down the channel with constant shape and speed. Thisassumption make it possible to study the movement of the wave as a steady flowcondition by referring the motion to a moving coordinate system. An observertravelling at the wave speed in the same direction as the wave would see aconstant profile. (for more detail, see Henderson, 1966, chapter 1, and Abbott1979). With contribution from Lighthill and Whitham (1955) in kinematic wavetheory, and the work of Henderson (1963), the analytical method has hadsignificant development. However, the concept is still limited to the assumptionthat So is much greater than the two other slopes (i.e. , —-Y and^au + 1^)au .' dx^_ug ax g at /Thus it is not true or applicable in all cases.Another method to find solutions of the unsteady flow equations is anumerical method, which does not require any drastic simplification. It uses finitedifferences to substitute derivatives in partial differential equations, so theequations represent simple algebraic equations. However, the procedure is notas simple as it seems. Difficulties involving convergence, stability, accuracy, andefficiency of the calculation emerged. Although the procedure only deals withsimple algebraic equations, the number of equations is extremely large, far largerthan with analytical methods. Due to these difficulties, and by considering8Chapter 2. Literature Reviewphysical processes in natural conditions, some simplification of the Saint Venantequations have been made. These led to the development of methods such askinematic wave, diffusion, quasi steady and dynamic wave methods.In the kinematic wave method, the momentum equation has been simplifiedto become St. - so ---q—u (more detail can be found in section 2.4). While ingA\the diffusion method, it is assumed that the inertia terms (i.e., u au + 1 au—g ax g at 'are insignificant, and can be ignored. The momentum equation (eq. 2.2) can beu aywritten as st - s -- - — . In the Quasi steady method, the factor 1 au0 gA ax^ g atis dropped from the momentum equation. The last is the dynamic wave methodin which the complete unsteady flow equations are used.Solutions of the dynamic wave method can be categorized as direct andindirect or characteristic methods. In the direct method, finite differenceapproximations for derivatives are substituted directly into eq. 2.1 and eq. 2.2.In the characteristic method, the partial differential equations are firsttransformed into characteristic form, and then approximated with finitedifference equations to obtain solutions, either in a characteristic grid mesh orin a rectangular grid mesh. Numerical solution of these models can also beclassified further as either explicit or implicit, depends on the type of finitedifference scheme used in the solution. Implicit models for Saint Venantequations -- either linear or nonlinear -- are generally more complex than explicit9Chapter 2 Literature Reviewmodels. The development of implicit models is to overcome limitations of theexplicit models, e.g. the small time steps for numerical stability. (see Amien et.al.,1970, 1975, Price, 1974).Solution of the complete dynamic method could also be obtained by a finiteelement method, although for the time being, it has not found widespreadapplication in mathematical modelling of river or channel hydrodynamics. As faras one-dimensional models are concerned, the finite element method does notshow any advantage over finite difference or characteristic method. Also, thelegitimacy of its application to time dependent problems (e.g., dynamic flow) isnot always clear (Cunge et al., 1980, chapter:3, and Jayawardena, 1977).Here is a summary of the different terms considered in the hydraulic modelsdiscussed above. The dynamic wave method discussed here refers to thenonlinear implicit finite difference solution of the complete dynamic unsteadyflow equations.Sf = SO_ qu _ ay _ u au _ 1 augA ax g ax g atKinematic Wave method...Diffusion method^ ■•■Quasi steady method^Characteristic methodDynamic Wave method^^Finite Element method^10Chapter 2 Literature Review2.2. DYNAMIC WAVE METHODThe development of the dynamic wave method has been pioneered byPreissmann (1961, in: Samuels, 1990), Abbot & Ionesscu (1967, in Amien &Fang, 1970), Blatzer & Lai (1968), Dronkers (1969), Amien & Fang (1970) andFread (1973). However, the most significant progress occurred after Amien &Fang (1970) introduced the use of Newton Raphson iteration method to findrapid, stable and accurate solutions of the complete Saint Venant equations.Later, Fread (1973) and Amien & Chu (1975) improved the techniques and triedto apply to broader cases. Fread applied the method for forecasting flash floodscaused by a dam failure and floods in rivers with flood-plains and meanders(Fread 1976, and Chow 1988). With some modifications, Ball (1985) applied themethod in urban drainage pipe networks, Schulte & Chaudhry (1987) in channelnetworks. Garbrecht et.al . (1991) reported that DAMBRK (1988 version) hasbeen used as a standard for testing other models. DAMBRK is the NationalWeather Service's dam break flood forecasting model program based on theimplicit complete dynamic wave method, and it has had good success inapplications to real problems.Studies of its numerical stability and accuracy have been reported by Amien(1970), Price (1974), Ponce & Simons (1977)and Schulte & Chaundry (1987).The method has been reported numerically stable with small, medium and large11Chapter 2. Literature Reviewtime step. However, Price (74) found that instability may occur if the time stepwere too large and the weighting factor was not sufficiently weighted to its futuretime step (time: j+1). The time step is also restricted in size for reason ofaccuracy. It is found to depend upon the shape of the wave, Courant condition,distance step, and of course, the type of implicit scheme used. Numericalinstability may also occur due to irregular channel cross section that vary rapidlyin the x - direction.23. CHARACTERISTIC METHODThe method of Characteristics assumes that a flood wave is a disturbance inthe free water surface. Any disturbance occurring at a particular point in openchannel flow propagates along the channel in two directions, downstream andupstream. A pair of paths of a disturbance from point A at time t = 0 movingdown-stream and upstream can be drawn in x-t plane. The paths are calledcharacteristic lines, C1 and C2. The area between a pair of C1 and C2 and pointA at time t > 0 represent the range of influence of the disturbance at point A.Inversely, taking a point D one can define backward in the time domain whichdisturbances can influence the condition at point D, see Fig. 2.2. By using theconcept of a characteristic line, range of influence and domain, solutions of12Chapter 2. Literature Reviewunsteady flow can be found.A^B^C^xFig. 2.2 Range of influence and domainComing back to basic mathematical definition of the total derivative:ay ^at+dy --clx —ax^ dt .dy _ ay dx + ayordt^ax dt^atIn these equations, y is function of x (distance) and t (time). If x and t aresimultaneously varied in some prescribed manner as dxldt, eq. 2.4 will representthe rate of change of y. A similar result would be true for any other parameter,such as velocity, u : du - au dx au+dt^ax dt^at (2.5)With some algebraic manipulation, from eq. 2.1, 2.2, 2.4 and 2.5 one could(2.4)13Chapter 2 Literature Reviewwrite:dx - u ± c^ (2.6)dtdu1 dydt ± —c dt + g(St-So) - 0 ^(2.7)in which : c - , I gA ^ (2.8)\TEquation 2.6, 2.7, and 2.8 are called characteristic forms of the equation ofunsteady flow. Solutions to those equations could be found by using finitedifference approximations, using either a rectangular grid mesh or a characteristicgrid mesh of C1 and C2. More detailed explanations of the characteristic formand characteristic grid mesh including its applications to simple problems areprovided in standard reference books (e.g., Abbott 1979, Cunge et al. 1980, andHenderson 1966).Blatzer (1968), Cunge et al. (1980), and Abbott (1979) reported that flowscomputed by the characteristic method generally indicate good agreement withfield measured flows. Liggett (1967) said that the method of characteristics is themost suitable general method in that it gave a good result over wide range ofparameters. Even Cunge et al. (1980) noted the method could be used as astandard method. It is possible to prove that its solutions may be brought closeto the solution of the basic equations.14Chapter 2 Literature Review2.4. KINEMATIC WAVE METHODThe kinematic wave method has been used to solve the unsteady flowequations in many flood routing models. One of them is the widely used packageprogram HEC-1 that applies the kinematic wave approach for overland flow aswell as for stream flow in a watershed. Some of reasons for using this methodare because it is simple, it does not require downstream boundary conditions forsolving the above equations, and it was believed that its approach approximatesthe natural condition of flood flow. The method assumes that the effects of theinertia and depth slope terms in natural flood flow are small compared with thebed slope term, so that they can be neglected. The friction term in momentumequation mainly depends on bed slope, and eq. 2.2 can be simplified as :Sf - SO (2.9)With this assumption, it is much simpler to find the solution to the unsteady flowequations. Eq. 2.9 leads to a condition of steady uniform flow in whichManning's equation can be applied.Q - 1" ARV 3 S(V2 ^ (2.10)nEq. 2.10 can also be expressed as Q - aA b .or A - f3Q a ^ (2.11)15Chapter 2 Literature Reviewn p2/3Comparing eq. 2.11 with eq. 2.10 : 0 - —1/2^ , and a = 0.6. And usingSovariable Q (discharge) and A (area) to substitute y (depth) and u (velocity),eq. 2.1 will appear as aA^ac.4-^- Qat ax (2.12)Differentiating eq 2.11 to t gives,aAc,_ a p Qa_i aat^at (2.13)and substituting eq. 2.13 to eq. 2.12 gives,a 1 Qa-1a + aQ - 4^ax (2.14)Kinematic wave can also be discussed in term of the method of characteristics.Unlike the dynamic wave method, in the case of the kinematic wave method,there is only one kind of characteristic line C1 in the x-t plane, corresponding toa single direction of wave motion.Using the total derivative for increment in discharge flow :dQ = ao cbc + ac, dt ,ax^ator^ao dt + ao _ dQat cbc ax^dx (2.15)With q = dQIdx = lateral inflow per unit length, eq. 2.14 and eq. 2.15 result :16TP^Sou171 yChapter 2. Literature Reviewdxdt1 a pQa -1 (2.16)From differentiating eq.2.11 and rearranging it, gives :dQ _ ^1 dA^aP Qa -1(2.17)Eq. 2.16 and 2.17 give Ck = dx/dt = dQ/dA, where Ck = kinematic wave velocity.An observer moving along the river bank at velocity Ck = dx/dt would see theflow rate increasing at the rate of dQ/dx = q, and if q = 0, the observer wouldsee a constant discharge.Studying analytically a flood wave in a prismatic channel, Henderson (1963)reported that the flood wave velocity does not differ from that of themonoclinical wave. It does not subside when it travels downstream the channel.If the velocity is materially less than that of the monoclinical wave, it must beascribed to storage effects, either because of channel irregularities or other kindsof off-channel storage.Ponce et al. (1978a) indicated that the wave period must be very long for thekinematic wave method to be applicable with accuracy 95% or greater to a mildslope channel.(2.18)where : 7;P = wave period of a sinusoidal perturbation of steady uniform flow17Chapter 2 Literature Reviewy and u = depth and velocity of steady uniform flowIn general, the steeper the channel bed slope, the shorter the wave periodneeded to satisfy the kinematic wave assumption.Hromadka (1988) said that the selection of ex and e t values may havesignificant impact on the kinematic wave modelling result. Internal checks forselecting AX and A t values are needed in order to achieve an accurate solution.Internal checks should notify program users when inappropriate values of exand A t cause computation error due to channel storage effects. MeanwhileDawdy et al. (1989) reported that kinematic wave models should be usedproperly, especially in choosing the values of AX and A t. With the kinematicwave approach, the leading edge of the flood wave becomes steep with distancedownstream (Hunt,1987). This phenomenon may form a kinematic shock in along open channel, so that a standard kinematic wave is not recommended to beapplied to this kind of problem. The method is most useful in a quick-responding (urban) watershed, where lateral inflow and short channels preventthe formation of the kinematic shock.2.5. MUSKINGUM-CUNGE METHODOne of standard methods of hydrological flood routing models is the well18Chapter 2. Literature Reviewknown Muskingum method. The method is based on conservation of mass whichapplied to storage, inflow and outflow within the reach.Storage equation : S - K( _z'x+ (1 -x) 0) ^ (2.19)Continuity equation : dS/c/t - I - 0 (2.20)in which : s = absolute storage within the reachi = inflow dischargeo = outflow dischargex = a weighting factorK = a gradient of the storage vs. the weighted flow curve, and K isrelated to time lag or travel time of the flood wave throughthe reach.From the equations above, it seems that the Muskingum method is really asimple method. The difficulties are only how to determine proper values of Kand x which lead to an accurate result or prediction. K and x could bedetermined graphically from the values of weighted flows vs. their pertinentstorages (Linsley, 1982), in which x varies from 0.0 to 0.5. K and x can also bedetermined by the least squares method proposed by Gill (1978) and by Aldama(1990). However, all methods mentioned above require a large amount ofhistorical data on storage, inflow and outflow within the reach. It means thateach particular reach (with particular So and Manning's n) should have itshistorical data set. And for other conditions of So or Manning's n, it alwaysneeds another historical data set. Because the author wishes to compare theMuskingum method with other methods for different conditions of bed slope So19Chapter 2 Literature Reviewas well as Manning's n, he need to find another way of determining K and xwhich correspond to the changing of So and Manning's n.Cunge (1969, in: Koussis, 1980), Dooge (1973) and Koussis (1978) proposedother ways to determine K and x based on diffusive wave analogy, and Ponce(1978b), Gill (1979) and Koussis (1980) reviewed and improved their method.From eq. 2.19 and 2.20, it is clear that the way of determining K and x influencesthe solution obtained. The assumption of constant K and x by Dooge (1973) andKousis (1978) makes the solutions depend on the reference values (Qo yo and uo)chosen. Ponce et al. (1978b) proposed their improved Muskingum-Cunge methodwith variable parameters, while Koussis (1980) considered that x varies accordingto time step. In general, however, it is more physically correct to allow both Kand x vary in time and space according to flow variability. For this reason, theauthor chose to use the Muskingum-Cunge method with K and x varying withinthe calculation.Another improvement to the conventional Muskingum method was proposedby Tung (1985) and Singh (1987) by introducing non-linearity in eq. 2.19.However, the methods do still not allow changing So and Manning's n. K and xfor each particular reach must still be determined from a large of historical datasets of inflows, outflows and storages.20Chapter 2 Literature Review2.6. UBC FLOWThe University of British Columbia (UBC) Flow model was originallydeveloped to cope with local problems encountered (i.e., ungauged lateral inflowcoming from snowmelt, steep bed slope, short travel time, etc.) in the FraserRiver, British Columbia, Canada. The method has been tested over many yearson the Fraser river system, and on other rivers (e.g., North Saskachewan River,the upper reaches of Columbia river). The results indicated that the flow modelis flexible, easy to calibrate, and suitable for both large and small rivers (Quickand Pipes, 1975). It is also easy to modify the storage portion of the model forlake and reservoir routing. Lateral inflow is excluded from the routingcalculation within the reach, because it can be added in at downstream boundary.Routing coefficients are determined directly from velocity-discharge and area-discharge relationships.UBC Flow uses a kinematic wave approach, and then adds an extra term toaccount for the influence of channel storage. This extra term is the depth slopeterm, ay/ax . Instead of using a finite difference approximation, the solutionsare determined by decomposing simply into a pure of translation of the floodwave, followed by a decay function. These two separate operations involvecalculating the travel time from kinematic wave velocity, and then by routing theflow through a simple linear reservoir. (for more detail, see Quick and Pipes,21Chapter 2. Literature Review1975, 1976).The UBC Flow model consists of two independent procedures. The firstprocedure is a calculation of the travel time for a particular reach as function ofriver stage (e.g., velocity-discharge and area-discharge relationship). The inflowis translated through the reach using this travel time to give translated outflow.The second procedure is routing. The translated outflow is routed through asimple linear reservoir which represents the storage characteristic of the reach.Lateral inflow coming from snow melt or any additional tributary inflow can beadded at the end of each reach.These two simple and independent procedures mentioned above make theUBC Flow model flexible and easy to fit to a real system. Also, the physicalbehaviour of a flood wave in a channel can be almost completely described bya travel time and the subsidence of the flood wave (modeled by routing througha reservoir). The travel times and reservoir storage are related to the channel'sparameters (i.e., So, Manning's n, bottom width, bank slope) as velocity-dischargeand area-discharge relationships. This condition makes the model accurate andcomparable to other methods that use basic river flow data.22Chapter 3THE METHODS3.1. DYNAMIC WAVE METHODThe dynamic wave method discussed here refers to the non-linear implicitfinite difference approximation to find solutions of the complete dynamic unsteadyflow equations. The method has been modified to fit the condition and purposeof this study. To obtain the results, the method uses an implicit finite differencescheme as seen in figure 3.1,t(time) 0.5(x)t'+'KAt'AtVXI^X14.1^X1+2^X(distance)Fig.3.1. Four point weighted implicit finite difference scheme1^•2^+K = 2-0 [KI:+K1 +1 ] +^( 1-0) [Kii1 +Kij] ^ (3.1)23Chapter 3. The Methods0K^0 rE,,,J+1 „, , , 1-0 r^JiLAxi+1^J+1^-r ^/ Ln.1+1 - 1\7:ax ex AXaKat 1^i+1.2 AX [K1 +1 +Kiwhere 0 is Preissmann's weighting factorThe lateral inflow is removed from eq. 2.1 and 2.2, and added in at thedownstream end of each routing step, in order to simplify the equations as wellas the comparison (e.g., since one of the compared models -- the UBC Flowmodel -- eliminates the lateral inflow from its main equations, and adds it in thecomputation at the end of the reach). Equations 2.1 and 2.2 can be written as,(3.2)(3.3)F(u,y)ay + A au u ay 0at T ax^ax (3.4)^G(u,y)^at^u—ax g axut^au^ay ._ (sf_ so)^0 ^ (3.5)Applying finite difference formula to eq 3.4 and 3.5 leads to,AX •^• -^•F )2e t^(u,y^•^[Y1++11+yr^Y j ]^0 [ A ( j+1 '1+1 )^j+1 j+1i+i --^—u -u •^+u yi+i -y ) ]T^.+ (1-0) [—I^-L. 1:1 )^- 0T 1+1 1AXi ul++11+u1 4. 1 _ uji+i _ uij + 0 E IT( u11:11: _ up' ) ±g yil++3.1 y + 1 ) ]^G(u,y)^A t+ (1-6) [17(u_i +1 -u ji ) +g(yi+ i -yy-1. 1 ) + geoci [5 f-So ] = 0(3.7)(3.6)24Chapter 3. The Methodsin which i and j indicate distance and time step, and ;^= 2 e ^+^(1-0) [uj +u1+1^i 1+12S1 ^1 0 [S j+ 1 +s i.+ 1 ] +^(1-0) [S^+S j.]2^fi+i f.1^2 fAT- 1 0[(-14-)j.+ 1+1+(l)i:'1 ] +^(1-0) [ ( 1 )^+ (4) 1^2^T 1^T 1^2 T 1+1 Tn 2 U2 I U2 ISf (Ri) 4/3Equations 3.6 and 3.7 prevail along the channel, while at both the upstream anddownstream end of the channel, boundary conditions should be applied. Allvariables with superscript (j+1) are unknown, on the other hand, all variables withsuperscript j could be calculated. Within a reach, those equations constitute asystem of two non linear equations with four unknown. If the channel is dividedinto (N-1) reaches (i.e., N node points), there are (N-1) interior rectangular gridswhich consist of 2(N-1) equations with 4(N-1) unknowns. However, two unknownsare common to any adjoining grids, so that the system consists of 2(N-1) equationswith 2N unknowns. Two additional equations needed for determining all theunknown are obtained from the boundary conditions.25Chapter 3. The Methods3.1.1. Boundary and Initial ConditionsBoundary ConditionsConditions at both the upstream and downstream ends of the channel(external boundaries) should be specified to obtain solutions of the completedynamic unsteady flow equations. Either a specified discharge hydrograph or waterelevation time series could be used at the upstream boundary. However, thehydrograph should not be affected by downstream flow conditions at the upstreamboundary.Those conditions for the upstream boundary conditions can also be used forthe downstream boundary conditions. Other conditions which can be applied at thedownstream boundary are the relationship between discharge and water elevation(depth), loop-rating curve based on Manning's normal flow, or critical flow at thedownstream boundary. Location of the downstream boundary should be chosencarefully, so that the flow at that point is not greatly affected by flow conditionfurther downstream. There are always some minor influences on the flow due tochannel irregularities, however, these effects usually can be neglected unless theirregularities are very large such as to cause significant backwater or drawdowneffects.For subcritical and critical flow conditions, the dynamic wave method requirestwo external boundary conditions, one for the upstream-end and the other one for26Chapter 3. The Methodsthe downstream-end boundary. However, some of the compared methods in thisstudy (e.g., the kinematic wave, Muskingum-Cunge method) use only upstreamboundary conditions. To compare all methods here, therefore, boundary conditionsfor the dynamic wave method should be modified, so that the method iscomparable to the other methods. In this case, a hydrograph discharge is used atthe upstream, and a normal flow is assumed at the downstream boundary.The upstream boundary condition is described by an equation below,:^Q3-1. + 1Fl(i,y) — QUP (t) — 0 ^ (3.8)and the downstream boundary condition by an equation below,Sc;12^ (3.9)FN(,, y) : UN-1+1 (Ri+1)213^0 ^Nnin which : QUP ( t) = discharge hydrograph at the upstream boundaryR^= hydraulic radiusIf unsteady flows are routed through the waterways passing through a dam, ora bridge-embankment intruding into the river channel, the flows may vary rapidlyat that point (i.e., location of the dam or the bridge) so that the Saint Venantequations are not applicable. The dynamic wave method copes with theseproblems by providing internal boundary conditions. Internal boundaries arechosen at the upstream and downstream extremities of the section of waterwaywhere rapidly varying flow occurs. The distance between the two internal27Chapter 3. The Methodsboundaries (i.e., upstream to downstream) can be any appropriate value from zeroto the actual measured distance. Because of these internal boundaries, twoadditional equations are required. The first equation should represent an empiricalrapidly varied flow relationship at that point such as flow over spillway, crestcontrol and so forth. The second equation should represent the conservation ofmass at that point with negligible time dependent storage.Initial ConditionEither discharge or water depth for initial conditions at time t = 0 should begiven to obtain solutions of the Saint Venant equations, and these must beaccurate enough to give convergence in the Newton-Raphson iterations.The initial condition may be obtained by any of the following :- discharge or water depth at gauging stations- computed values from a steady flow backwater solution- computed values from a previous unsteady flow solutionDischarge at gauging station is used for initial condition here, while the programis written to have options whether discharge or water depth may be used.3.1.2. The generalized Newton Raphson iteration methodThe application of the generalized Newton Raphson iteration method to28Chapter 3. The Methodsdetermine simultaneously solutions of system of nonlinear equations is describedhere. Let the 2N equations mentioned above (i.e., 2(N-1) eq. from interior pointand 2 eq. from boundary conditions) be represented by a set of equations below,^F,^ = 0^-^X2 ,^, ^ X,,,F2 (xi^xm) = 0F2N ( xi•x2 , x3^xm) - 0Or they can be written as, F IX)^Fi (xj) - ^i - 1,2^ 2Nj - 1,2^Newton Raphson method (after k-th iteration) is defined by :(3.10)Xk+1 - X k- f (2rk) 1 1 17(xk) (3.11)in which .f (z.k) is a Jacobian matrix at k-th iteration step,aFi(x.)f (z) - ax; - 1,2,3, ^ 2Nj - 1,2 3 Eq. 3.11 can be written as, f (xid (xk+1 --xk) -29Chapter 3. The Methodsor, f (ark) Axk + k) - 0 (3.12)Equation 3.12 represents a system of 2N linear equations with 2N unknowns.Assigning initial trial values for z k , the solution can be obtain by the Gaussianelimination method with pivotal condensational strategy to achieve accurateresults. This will gives AlC k , which leads to the determination of xk+1 . Theiteration procedure (i.e., the calculation of channel's variables, Jacobian matrix,Gaussian elimination, etc.) is repeated until the difference of any unknownbetween two consecutive iteration falls below an error tolerance. In this computerprogram, two kinds of error tolerance are set, the first one is YERR for y (flowdepth) and the other is UERR for u (velocity). A stagnation may happen duringthe iteration procedure, in which Axk of any unknown stands at the same valueafter two or more consecutive iterations, and it is greater than the error tolerance,YERR or UERR . To avoid this problem, an algorithm is set to enhance theiteration.An initial condition gives the first values of y (water depth) or Q (discharge)at time t=j, and then u (velocity) is calculated by iteration. Those values at timet=(j+/) are unknown. Assigning initial trial values for the unknown and runningthe iteration mentioned above until convergence will result those values at timestep^ 371++1, u 4+1 and 111+1^17 j41-2 • For the firsttime step, the trial values along the channel can be obtained by using extrapolation30Chapter 3. The Methodsof water depth and velocity at the upstream boundary considering bed slope anddistance from the upstream end. For the next time step, the values from solutionsat the previous time step can be used. The closer the trial values to the exactsolutions, the faster the iteration achieves convergence. If the convergence is stillhard to accomplish, then Preismann's weighting factor should be added with asmall increment until a convergence is obtained for a particular 0.3.1.3. Stability ConsiderationAny numerical method used to find solutions of partial differential equationsmay be affected by truncation, round-off error or any other mistakes. These errorsmay grow as the calculation proceeds and lead to the wrong result. In order tocope with this condition, a stability requirement is needed to prevent thosegrowing errors. A numerical stability requirement for solving hyperbolic partialdifferential equations is given by the Courant-Friedrichs-Lewy (CFL) condition(Anderson, 1984, p.75), which is,I c-L-' I s 1 (3.13)AXThe dynamic wave method uses four point weighted implicit finite differencescheme to find solutions of unsteady flow equations. According to a Fourierstability analysis, the implicit method is usually stable ( Anderson, 1985, p.71-154),31Chapter 3. The Methodshowever, a system of equations must be solved simultaneously at each time step.As long as the time step is not too large, the result will be close to the exactsolutions. In this case, the time step is taken as same as that of the othermethods, but the CFL condition is applied. The spatial mesh spacing (i.e., distancestep) is chosen in a particular way so that the Courant number is close to one. Adouble precision arithmetic is chosen during calculation to limited a growing errorcaused by round-off errors.3.1. CHARACTERISTIC METHODIn the method of Characteristic, equations of unsteady flow are transformedinto characteristic forms before they are solved by finite difference equations. Theequations of unsteady flow from eq. 3.4 and 3.5 are,F(u,y)ay 4_ A au + u ayat T ax^ax 0a t^au^ay + _ , s fm-90) - 0G^u-Ec. + g- I -c.^gk(u,y)^ +In characteristic form, these two equation above can be united to form a linearcombination with an unknown A as follow;H(u,y) : G (u,y) + A F (u,y) - 0H(u,y) : {wau + Tcau (u+1 _A )] + A [ atay +ax^iL7 (u+) + g ( sfm90 ) = 0...32dtdx = u ± c (3.17)Chapter 3. The Methods^(3.14)The first and second term of equation 3.14 corresponds to the total derivative ofdu/dt and dy/dt,du au dx audt^ax dt^atdy - ay dx aydt^ax dt^atand, so that one can rewrite eq.3.14 as,H(u,y) : du +A dy + g(Sf-S0 ) = 0 ^ (3.15)dt^dtdxdt ▪ U+ AT(3.16a)and dx ▪ u +dt^A (3.16b),\1Finding A, _ ±^gT , one can write equations of unsteady flow as below,Adu t 1 dydt^ c dt g (S f-S0 ) - 0 (3.18)in which : c^gAEquations 3.17 and 3.18 are characteristic forms of unsteady flow equationswhich describes the movement of a flood wave through a channel or natural river.The solutions to these equations can be determined numerically by finitedifference approximation.33Chapter 3. The MethodstL^RXFig. 3.2. Explicit finite difference schemeWith explicit finite difference scheme as in figure 3.2 above, the positivecharacteristic line ( C1 ) can be presented as equations below,up - Liz. + 9 (_yp-ydcL+ (St L-SO ) ( tp-td - 0 ^ (3.19)xp-XL - ( uL+cL) ( tp- tL) = 0 ^ (3.20)and negative characteristic line ( C2) can be presented as equations below,UP - UR - cgR (37P--YR ) + (SfR-So ) ( tp-tR) - 0 ^ (3.21)xp-xR - (UR-cd ( tp-tR) = 0 ^ (3.22)Equations 3.20 and 3.22 give,XL - XR + tR ( UR - CR )^ti, ( UL+ CL )tp( UR CR )^( UL+ CL ) (3.23)34Chapter 3. The MethodsXp = XR + ( uR- cR) ( tp- tp) ^ (3.24)and from eq. 3.19 and 3.21, one can determine,YpCL CR 1^YL YR"'^[- (uL — uR) +—+--( tp—tL ) (Sf —S )+(tp—tR ) (Sr —S,,)]+CR g CL CR^L 0^ -^(3.25)UP - UL - g (Yp-YL) - g(StL -So )(tp- tL) ^ (3.26)L3.2.1. Boundary ConditionsSince the method of characteristics uses the same basic equations as that usedin the dynamic wave method (i.e. both of them use complete unsteady flowequations), the boundary and initial conditions are taken to be the same. In thisstudy, the discharge hydrograph is used at the upstream end, and normal flow isassumed at the downstream boundary.Conditions at the upstream boundary can be described in the following equations,xp = 0tp = tRFB : Up- UB- --g- (37p-37B ) +g(S —S )(t —t ) = 0r R 0^P RCRFQ '• up AP — QUP (t) — 0XR UR -CR35Chapter 3. The MethodsaU p A PaF,^aApuayp^P aypConditions at the downstream boundary can be described in the followingequations,XN = XNtp " t i, +FF = Up- UL+-9L (yp-YL ) +g(S f L-S0 ) (t p-tdCL1 2/3 1/2FN "' Up n—Rp SoaF ^aFN 1au p^1^all paFF^g^aFN 2 up ( ailp R app )aYp^CL^aYp - 3 PpRp ayp P aYFin which : xN = distance to end of the reachRp = hydraulic radius at point PPp = wetted perimeter at point PAP = wetted area at point PFrom the equations above, one can find y, and up at the upstream boundary byNewton Raphson iteration, route the flow downstream, and determine y p and upat the downstream boundary. The method uses variable time steps to keep theaFBau paFBaypiCRaFQXN-XLUL+ CL36Chapter 3. The Methodsdomain of dependence of characteristic lines C1 and C2 from point L and R alwaysover point P. Due to this condition, discharge at the upstream boundary is takenfrom hydrograph discharge by linear interpolation. While at the downstreamboundary or any other fixed point P, Lagrang's four point interpolation (Kreyszig,1983, p.'7'79) is used to convert any values of y1, and up at time tp to theappropriate time step the same as the input's time step.3.3. KINEMATIC WAVE METHODThe kinematic wave method assumes that the friction term (Sf) mainlydepends on the bed slope (So). The inertia terms (i.e., —u au + —1 au ) and theg ax g atdepth slope term (i.e.,h'- ) are small and they can be neglected. By theseaxassumptions, it is possible to simplify the unsteady flow equations to the muchsimpler form as in eq. 2.14 (see section 2.4.).The lateral inflow ( q ) can be taken out from the equation and added in atdownstream of the reach as well as the other compared method (e.g. UBC Flow).Equation 2.14 is rewritten as,apQa_i aQ + aQ - 0 ^at^ax (3.27)Solutions of eq. 3.27 is determined by finite difference approximation. A37At•QI DX041Chapter 3. The Methodsbackward difference scheme as in figure 3.3 is used to set up the equation.tt'+'t'X i^X1+1^Xi+2XFig. 3.3 Finite difference scheme for the kinematic wave method.^ao ^-x^AX^ao ^- 01-1+1t^AtQ1 +1^Qii+1 Q 2Applying those equations above, linear finite difference form of equation 3.27 iswritten as,^j ^j+1 a-1^,-, j+1( Q 1 + 1 ^) vi+1 ^) ^ - 011+12 A t^ AX^)^0 ....(3.28)38Chapter 3. The Methodsnij^nj+i a-1[ A t 01+ + a [ 02+1 ( 51 +i + wi )^]n +1 ^AX^ sei+i - ^) + fl.1+3. a-1[ At ÷ a p ( Vi+1^5.".z )AX^2 (3.29)For each ex , the water depth is calculated by iteration to determine p asin equation 2.11. The Courant-Friedrichs-Lewy condition (eq. 3.13.) is applied forselecting AX .AX^ dA1 , in which c = dQ , a kinematic wave velocity.To achieve the most accurate result, ex should be chosen appropriately so thatthe Courant number is close to one.3.4. MUSKINGUM-CUNGE METHODThe Muskingum-Cunge method is based on the conservation of mass withinthe reach (i.e., continuity equation) and the relationship between inflow, outflowand storage (i.e., storage equation).Continuity equation : ds/dt = I - 0Storage equation : S = lc( /x+ (1 -x) 0)in which : s = absolute storage within the reach39Chapter 3. The Methodsi = inflow discharge0 = outflow dischargeX = a weighting factorK = a gradient of the storage vs. the weighted flow curve, and K isrelated to time lag or travel time of the flood wave throughthe reach.Those equations above lead to the determination of outflow discharge as follow, Cl Q ."^c, (21+'^c, Qij+ , (3.30)At/K + 2xC2At/K + 2 (1-x)At/K - 2xAt/K + 2 (1-x)C3 - 2 (1-x) - At/KAt/K + 2 (1-x)Using kinematic wave theory and an assumption of single value depth-dischargerelationship, Cunge developed equation 3.30, and finally he found that K and xcould be determined by the following techniques,K= XC1 OW rix ^-2 cSoexc dQCunge found that parameter K and x can be obtained from the length of the reach(ex), the flood wave celerity (c), and the discharge per unit width ( Qi, ), insteadof determining them from a large amount of historical data. Later, Ponce anddA40Chapter 3. The MethodsYevjevich continued this development by assuming that K and x are variableparameters. K and x are allowed to vary in time and space (distance) followingthe fluctuation of the flow.Flow depth for each point grid is calculated by iteration. Knowing the flowdepth, one can find c and Q,,, for each distance step using a three-point averageof the values at grid points (i,j), (i+/,j) and (i,j+1). From these values, he cancalculate Muskingum's parameters : K, x, C1, C2 and C3. These parameters varyin time and space during the calculation.An accurate result is always the main objective of computation in floodrouting. The result of computed flow should be close to the observed flow. Oneof the efforts to achieve an accurate result comes from the way of choosing Ax(distance step) which influences the values of Muskingum's parameters. Fromhistorical flood data, it was known that parameter x varies between 0.0 and 0.5(Gill, 1978, Hjelmfelt, 1985). Hjelmfelt investigated a negative outflow that maycome out as a result of Muskingum's flood routing. To ensure positive outflow forany positive inflow sequences, he derived a criterion as follow,0. 5 A t X 5^5 (1-x) , for x s 0.50 .From this argument, ex should be limited as follows,(cat - Q„ ) s 46,x s ( cat + ^ (3.31)^cS^ c S0 041Chapter 3. The MethodsIn addition, ex should be chosen appropriately so that the Courant number isclose to one.3.5. UBC FLOWThe University of British Columbia (UBC) Flow model consists of twoindependent procedures. The first procedure is the calculation of travel time fora particular reach, followed by translation of the upstream flow down the reach.It is assumed that the flood wave travels downstream at the kinematic wavecelerity. The travel time, TT, is determined from the velocity-discharge relationshipfor the reach, and then, is corrected for kinematic wave travel time by a constant(i.e., RTK = v/c). RTK is usually assumed equal to 0.67. The length of the reach,L, is chosen appropriately so that the value of TT is just below or equal to onetime interval used for upstream inflow data (e.g., 1 day, or whatever time-intervalis used).Using this travel time, the upstream flow (Q1 ) is translated through the reachto give a translated outflow (Q07. ). To maintain Q0T data at the same datainterval as that of Q., (e.g., 1 day, or any interval of time), an interpolationprocedure is used. Figure 3.4 defines the translation of an upstream hydrographto some downstream point at one interval time later.42TT ( TT < 1 )Yesterday's flow at^ Today's flow atupstream stationAOYesterday's flow Ifat downstream ate.9,7If TT w 1downstream station.^:26,••If TT < 1 , then today's downstream flowIs greater thanyesterday's upstream flowfor a rising hydrographChapter 3. The MethodsKINEMATIC WAVE TRANSLATION OF FLOOD WAVEby ignoring storageFlowTimeFig. 3.4. Flood wave translation diagram(taken from Quick & Pipes, 1975)The upstream hydrograph in figure 3.4. is assumed to be linearly increasingand the interval time of inflow data used is one day. If 77' were exactly equal toone day, the downstream flow today would be equal to yesterday's upstream flow.Using figure 3.4 and defining the upstream flow increment for the previous dayas AQ , the relationship between the translated downstream flow Q07.(i) andyesterday's inflow Qi(j _1) will be as follow,Q0T(J)^( 1 - TT) A Q ^ (3.32)There are two alternative situations provided,If 77' > 1 , qz-u)^(1-TT) [QI(j) - QI(j-1)] ...(3.33)If TT < 1, 0-0T(J.1) a01- (;) + (1-TT) [^(i +1) - QI(j)] ...(3.34)43Chapter 3. The MethodsThe second procedure is to route the translated downstream flow ( _007.(.i) )through a simple reservoir, which represent the storage characteristic of the reach.It is assumed that the storage, S , is linearly related to the outflow, Q0 , witha storage constant K which is continuously updated for every changing ofdischarge.This relationship gives^S - K Q 0 ^ (3.35)Conservation of mass gives^dS/dt - Qi. - Qo ^ (3.36)Quick and Pipes (1975) stated that for a linear process, each day's flows areindependent of their previous flows, and that yesterday's flow is used as referencelevel for both inflow and outflow, e.i.,AO./ - Q'(j+1) - Qou) ^ (3.37)AQ0 - Qo(J+i) - Qou) (3.38)Equation 3.37 stated that the daily increment of inflow is calculated relative toyesterday's outflow.Rearranging eq. 3.37 and 3.38, and applying them to eq. 3.36 gives,dSd t - °Q.r - °Qo (3.39)Combining eq. 3.39 with eq.3.35 gives,K A tQ° ACl/ - AQ0(3.40)44Chapter 3. The MethodsWriting x/A t as KT and solving for Apo gives, 1 A Qo^eQ,1 + K7'^-(3.41)and the routing procedure is described by :Qou,i) - Qo(J) + AQ0or 00(j+1) - QC(J) + 1+ KT1^ ( Q I ( j+i) - 00 (j)) (3.42)Equation 3.42 could be used for both reservoir routing and channel routing.For reservoir routing, KT - dS 1dQ etFor channel routing, the value ofKTis derived from an area-discharge relationshipas below,KT = L^idA STKdQ At , n which STK is 0.50 and L is channel's length.The inflow entered into the simple reservoir in eq. 3.42 (i.e., p i(j+1) ) is theupstream flow after it is translated to downstream (i.e., povi .1) in eq. 3.33 or3.34 ) with kinematic wave celerity. Therefore, eq. 3.42 is rewritten as,1 Qo (3.43)u+i) - Qo(J) + 1 + KT ( Qtyr(i+3.) - Qoul) ^45Chapter 4COMPARISON OF THE METHODS4.1. COMPARISON PROCEDUREFive flood routing models were compared to find the most reliable model.The accuracy of the results and the applicability (e.g., ease of use) of the modelswere examined when they were applied to channels with various bed slopes andfriction coefficients. The numerical models were written based on five differentmethods as mentioned in the previous chapter (i.e., the characteristic method, thedynamic wave, the kinematic wave, Muskingum-Cunge and UBC Flow methods).Appendix A presents flowcharts of the models. In order to get accurate results,t,x was chosen carefully so that the Courant number was close to one whenrunning each of the models. The error tolerances used in the iteration were0.0001 meter for y (flow depth) and 0.0001 m/sec for u (flow velocity).The above models were applied to channels with various bed slopes andfriction coefficients. Their results were compared for each particular combinationof bed slope and friction coefficient of channel. Ideally, the comparisons shouldbe done using data from natural rivers, so that observed data at both theupstream and downstream boundaries, bed slope and friction coefficient areneeded for each channel. However, because of the difficulties in obtainingsufficient data from natural rivers, most of the comparisons are done in an46Chapter 4. Comparison of the Methodsartificial prismatic channel.The comparison consists of two steps. The first one was a check on the fiveprograms by applying them to Fraser River flows, and comparing computed flowswith observed outflows. The second involved comparison of results of theseprograms for a uniform channel with various bed slopes and friction factors(Manning's n). For each combination of bed slope and Manning's n value, thecomputed outflows were compared with each other and with a referencedischarge. The latter was obtained from the dynamic wave method, i.e., themethod which gave the closest result to the observed outflow in the first step.In choosing the dynamic wave as the standard method, reports and analyses fromprevious studies (e.g., Price, 1974, Garbercht & Brunner, 1991) were also takeninto account.In the first procedure, inflow data was obtained from recorded flow datafrom the Fraser River at Marguerite hydrometric station (May - July 1988).Recorded flow data at Texas creek gauging station (270 km downstream) wasused as observed outflows (see figure 4.1.). The average bed slope betweenthese stations is 0.00079. A constant friction coefficient (Manning's n = 0.055)and an average channel width = 150 m were used for the computationthroughout the reach.Since the observed outflows included lateral flows entering along the channelreach, the lateral inflows were subtracted from the observed outflows.47ACTIVE STATIONSHYDROMETRIC HYDROMÉTRIQUESSTATIONS EN SERVICELEGEND LEGENDE^SMEwoiLOW ^•^DOOIT^SEO.A00' a STREAuflOw^•^SEDO•ENT a DEBITve■Teil LEV, •oVEAU OVAIJ^DRAINAGE INDEX^D^INDEX DE DRAINAGE.4.04 DrOSION .S1ON PIONC.0,4ESOO qv.. ^ SuODIVISIP-1SUB-SUB OnOVON — — — SOUS-SUBdwiSIONALBERTAancouvel okirArU. MOORS& "EOM" STATIONS KTMWFLED./CA•Tt - F WDOCE4.01^0•• 0OCOAORETexas creekIt■•Chapter 4. Comparison of the MethodsFig. 4.1. Location of Texas creek and Marguerite hydrometric sta.,Fraser River, British Columbia, Canada.484000 —‘.^• cr-o-0Fig. 4.2. FRASER RIVER HYDROGRAPHS (May — July 1988)^i ^■^■, ^I^\I ^1 I , `^* ^/^I^I —‘^ / \^ /^'id,^i l^‘/^\1 ^I,/ ^--\^/ •^i ^\ •^I ,\ ■ , ,‘^/^ il^• , / "/ ; *L .,^I\^•4 /^\^‘^. A ■^t i^\,/ Pf .^/ f^*--./‘\\i . \^i^li■\‘ .`k^ ii^'''*-'..^' 1^11 ‘ /05`ss)s I IIA^/0., 11 't^\/0 `o^ sl'.,,,5000P4 2000 —\ 5- 0-0-2'vg1000 —00^ Inflow Outflow• Outflow• Outflow1111111)111111111 11^II^11111111^III12^24^36 48^60TIME (days)96recorded at Marguerite)recorded at Texas Creek)version 1)version 2)I^1^III^1111^172^84Chapter 4. Comparison of the MethodsTwo different approximations were applied here :1. x1 1 — X0 ifX0i — fX.Ti( 3.^3. ^)N (4.1)t X0i - t Xii2. X2 i — X0i^( 1^1 ^) xi 1 ^ (4.2)f x_ri1in which : X1 1 = outflow version 1 at time iX2, = outflow version 2 at time iX0 = outflow data (recorded at Texas creek)XI = inflow data (recorded at Marguerite)N = total number of dataIn version 1, each outflow is reduced by an average lateral inflow. In version 2,the reduction for a particular outflow is proportional to the inflow. The biggerthe inflow, the more reduction is applied to the pertinent outflow. The results arepresented in figure 4.2.In the second procedure as mentioned before, the five numerical modelswere applied to a prismatic channel with various bed slopes and friction factors.The channel used was 80 km long, 20.0 m wide at the bottom and had a bankslope of 1 : 1 (horizontal vs. vertical). A discharge hydrograph from 100 cu.m/sto 250 cu.m/s with 5.00 cu.m/s increment and back to 100 cu.m/s with 3.73 cu.m/sdecrement was used as the inflow data. The time interval was 30 minutes for theboth processes. Initial flow was 100 cu.m/s. The results of the computation were50Chapter 4. Comparison of the Methodscompared in two series of comparisons. The first one used a constant bed slopewith various friction coefficients (Manning's n), and the second series used aconstant Manning's n with various bed slopes.The first series of comparisons used bed slope = 0.0010 and Manning's nvaried from 0.015 to 0.055. For each combination of variables, the computedoutflows were compared. Discrepancies between the computed outflows and thereference discharge were calculated for each time interval. Distribution of thediscrepancies through hydrograph curve was presented in two ways. First, thediscrepancies were divided by the reference flows for each time interval, andsecond, the discrepancies were divided by the reference peak flow.Di1. DRF ^ (4.3)Q R2. DPFDi (4.4)Q PRin which : D RFi = discrepancies were divided by the reference flowD pFi = discrepancies were divided by the reference peak flowDi = discrepancy at time iQRi = the reference flow at time iQPR = the reference peak flowA discrepancy indication is determined by formula as follow,0.5t( QC OR)2DI = [^in which : DI = discrepancy indication(4.5)51Chapter 4. Comparison of the MethodsQc = computed outflowsQR = reference dischargeThe second series of comparisons used Manning's n = 0.035 and the bed slopeswere varied from 0.0005 to 0.0050. The procedures for the comparison are thesame as those in the first series.4 2 RESULTS AND DISCUSSIONResults of all the computations are presented in figures 4.3 to 4.9.b. Beforediscussing the results, there are some notes about how the compared modelsobtain their results. The dynamic wave method uses an implicit finite differencescheme to approximate the complete unsteady flow equations. The results areobtained by solving simultaneously a system of non-linear equations throughoutthe reach. A generalized Newton Raphson iteration method is used here, and aset of initial trial values is needed to find the solutions for every time interval ofthe hydrograph. The trial value data should be close to the exact solution,otherwise the iteration can not converge. Regarding this condition, it wassupposed that when running the model it would take time to find the properinitial trial values. However, a subroutine can help model users cope with theproblem (e.g., an approximation procedure with two constants input for creatinga set of initial trial values). By choosing carefully the constants for the initial52Chapter 4. Comparison of the Methodsguess, everything proceeds smoothly. Convergence usually occurs in the third orfourth iteration cycle. From various combinations of bed slopes, Manning's n and AXused in this study, the dynamic wave method has shown its applicability andconvergence.The characteristic method transforms the complete unsteady equations intothe characteristic forms, then it approximates them by an explicit finite differencescheme. The method needs a small time step, so that the upstream discharge isobtained from the inflow data by linear interpolation. The results are convertedto become those at the original time step by Lagrangian 4-degree polynomialinterpolation. A Newton Raphson iteration method is used to obtain results thatfulfil conditions at both the upstream and downstream boundaries. In this study,the characteristic method with an explicit scheme failed to obtain the resultswhen it was applied to a steep bed slope channel (e.g., So z 0.0050). Themethod needed a bigger discharge or flow depth for its initial condition. Thisrequirement limits its application. The method also needs longer computing timethan the dynamic wave method does.The kinematic wave and UBC Flow methods simplify the unsteady flowequations. The kinematic wave method neglects both the inertia and the depthslope terms, while the UBC Flow model still considers the depth slope term.However, the latter model obtains the results by a hydrological approach. There53Chapter 4. Comparison of the Methodswas no difficulty in using both the kinematic wave method and UBC Flow model.In the Muskingum-Cunge method, it was somewhat difficult to fulfil therequirements set in equation 3.31 (section 3.4) for a steep bed slope. The timeinterval should be made small enough so that ex can fulfil the Courant-Friedrichs-Lewy condition and ensure positive outflows. The above requirementscaused the Muskingum-Cunge method failed to obtained its result when it isapplied to a channel with So = 0.0050 and Manning's n = 0.035.The computed Fraser River outflows from the five models and the observedoutflow (version 1) are presented in figure 4.3. Although none of the computedoutflows perfectly matched with the observed outflows, all hydrographs from thecomputed and the observed outflows had good agreement. It also has to beconsidered that the observed outflow (version 1) was obtained from anapproximation. From the data mentioned above, it is impossible to know thereal lateral inflow for each time interval. No reliable guide was available forestimating the lateral inflows. Among the computed outflows, it is clear that thehydrograph from the dynamic wave method is the closest one to that of theobserved outflows as can be seen in fig. 4.3. Except for the absolute value ofdischarge, the method usually reached the peak flow coincident with the observedoutflows. On the other hand, the peak flows of the other methods arrived earlierthan that of the observed outflows. This condition causes the discrepancies544500^ Inflow (recorded at Marguerite) Recorded Outflow (ver.1)• Kinematic Wave Method• Muskingum—Cunge Method• — Characteristic MethodFlood Dynamic Wave Method (0=0.65)+++++ UBC Flow Model3500/■U 2500 —Fig. 4.3. COMPUTED OUTFLOW HYDROGRAPHSFraser River at Texas Creek gauging sta. : (May — July 1988)I^I^I^I^I^1^I^I^I^I^I^1^I^I^1^1^I^1^11^1^I^I^1^I^I^I^I^I^1^I^F^[^I^I^1^1^1^i^1^1^1^iiii^1^12 24 36 48 60 72 84 96TIME (days)1500TLX = 80 kmB = 20.0 mZ = 1.0S. = 0.001- Inflow^ Outflow with n = 0.015.4.4 ^ Outflow with n = 0.025••-4-•• Outflow with n = 0.035+04+0 Outflow with n = 0.045*0000 Outflow with n = 0.055O^ 1^1^1^1^1^1^1^1^1^1^1^1^f^1^1^1^1^1^1^11 11111111^1^1^1^1^1^1^1^1^1^1^1^f^1^1^1^I^10 10 20 30 40 50TIME (hours)_d 200DYNAMIC WAVE METHOD WITH VARIOUS S oin artificial prismatic channel(a)- Inflow^ Outflow with S0 = 0.005 Outflow•-•-•-••with S. = 0.002with S0 = 0.001Outflow0040. Outflow with S. = 0.00050 1 1 1 1 1 1 1^1^1^1^1^1^1^1^1^1^1^1^1^1 1^1^1^1^1^1^1^1^1^10 10 20 30 40TIME (hours)50TLX = 80 km= 20.0 mZ^= 1.0n^= 0.035(b)300 ^Chapter 4. Comparison of the MethodsDYNAMIC WAVE METHOD WITH VARIOUS MANNING'S Nin artificial prismatic channel300Fig. 4.4. Computed outflow hydrographs using the dynamic wave method4.4.a. : So =0.0010 with various Manning's n values4.4.b. : Manning's n =0.035 with various bed slopes56Chapter 4. Comparison of the Methodsbetween the observed outflows and the computed outflow of the dynamic wavemethod to be less than those between the observed outflow and the valuescomputed with the other methods. Therefore, the dynamic wave method is mostsuitable for obtaining the reference discharge. This results agrees with reportsand analyses from the previous studies (e.g., Price, 1974, and Garbercht &Brunner, 1991); and in fact, it is what we would expect.Figures 4.4.a and 4.4.b present the results from the dynamic wave methodfor a channel with various bed slopes and friction factor coefficients. Figure 4.4.ashows the computed outflows for a constant bed slope with various Manning's nvalues, and figure 4.4.b is for a constant value of Manning's n with various bedslopes. From figure 4.4.a, it is clear that the time to peak increased and the peakflow decreased when Manning's n was raised. The same situation occurred whenbed slope became milder. These phenomena are consistent with the actualsituations in natural rivers. Increment of the time to peak and decrement of thepeak flow were not proportional to the increment both of bed slope andManning's n values. The rate of change of both the time to peak and the peakflow tends to be smaller for a bigger value of Manning's n and for a steeper bedslope.The results of the first series of comparisons for a prismatic channel arepresented in figures 4.5 and 4.6. Figures 4.5.a and 4.5.b show the computed57Chapter 4. Comparison of the Methodsoutflows for a constant bed slope (S0 =0.0010) with two values of Manning's n,i.e., n =0.015 and n =0.025 respectively. The discrepancies of the computedoutflows from the reference discharge are presented in two ways as shown infigures 4.6.a and 4.6.b. In the first case, the discrepancies were divided by thereference flows, while in the second case, the discrepancies were divided by thereference peak flow. Figures 4.7 and 4.8 display the results of the second seriesof comparisons, i.e., the computation using a constant value of Manning's n withvarious bed slopes. Complete figures of both comparisons are provided inappendix B.Based on figures 4.5, 4.7 and those in appendix B, it can be concluded thatresults of the dynamic wave method usually reaches its peak flow later than theother methods. This phenomenon is similar to the result of the test with theFraser River flows (see fig. 4.3), in which the dynamic wave method produced theclosest result to the observed outflows. Figure 4.3 showed that the peaks and thetoes of the hydrograph curve from the dynamic wave method coincide with thoseof the observed outflows. Considering these facts, one can draw a conclusionthat a bigger discrepancy from the results with the dynamic wave method wouldalso mean a bigger discrepancy with the observed outflows in a natural river. Forexample, with results of the dynamic wave as the reference discharge, results ofthe kinematic wave method in fig. 4.5.a and 4.5.b gave bigger discrepancies than58TLX = 80 kmB = 20.0 mZ = 1.0So = 0.001n = 0.015- Inflow• vvv. Kinematic Wave Method00000 Muskingum-Cunge MethodCharacteristic Method■-•-•-•-• Dynamic Wave Methodseeee UBC Flow ModelTLX = 80 kmB = 20.0 mZ = 1.0So = 0.001n = 0.025- Inflowv.v^ Kinematic Wave Method^ Muskingum-Cunge MethodCharacteristic Method•-•-•-•-• Dynamic Wave Method4444+ UBC Flow Model401^I^1^1^1^1^1^I^I^1^I^1^1^1^1^1^1^1^1^1^1^1^1^1^1^1^1^1^1^1^1^1^1^1^1^1^1^1^1^1^1^I0^10 20 30TIME (hours)3000^111111111111I)111111■1111111111111111^1^1^1^1^1^1^1^1^1^1^1^10 10^ 20^ 30^ 40 50TIME (hours)3000)_2000 f50200Chapter 4. Comparison of the MethodsINFLOW AND COMPUTED OUTFLOW HYDROGRAPHSin artificial prismatic channelFig. 4.5. Comparison with various Manning's n values4.5.a. : So =0.0010 and Manning's n =0.0154.5.b. : So =0.0010 and Manning's n =0.02559o°%;484ostgoo-000400000...ft+. AAAAAAA422Raea AAAAA^m 4,11,9.altIteo•4:4188484W1-^ +1"Chapter 4. Comparison of the Methods0.60DISCREPANCIES FROM THE REFERENCE DISCHARGE• 40r11.414':114IVegter.*44444.44uuttiiitutstuntiintt'.%4^% ,taaLat4A•aao......o '1111111111111,111111/111110^20^ 30^40^50TIME (hours)TLX = 80 kmB = 20.0 mZ^= 1.0S. = 0.001n 0.015^ Kinematic Wave Method••• ^ Muskingum—Cunge MethodCharacteristic Methodoo•seo UBC Flow ModelTK time when the reference dischargecornea to peak flowO• 0.4w 0.2Ix• 0.0—0.2Cl)rxFl —0.40 0.3.16 0.2• 0.1z• 0.0Cr.1Ix—0.1z'44ja, —0.2f:4Cr)ral,▪ —0.3—0.4011111- 11111111111110^20^ 30^40^50TIME (hours)Fig. 4.6. Comparison of discrepancies with So =0.0010 and n=0.0154.6.a. : divided by the reference flow4.6.b. : divided by the reference peak flow(a)(b)TLX = 80 kmB = 20.0 mZ = 1.0S. = 0.001n^= 0.015 Tr••••• Kinematic Wave Method• .••• Muskingum—Cunge MethodCharacteristic Method00000 UBC Flow ModelTw : time when the reference dischargecornea to peak flow60TLX = 80 kmB^20.0 rnZ^=^1.00.035So = 0.0005^ InflowKinematic Wave Method^Muskingum-Cunge Method Characteristic Method••-•-•-•-• Dynamic Wave Method00000 UBC Flow Model3000^111111111j111111111-11111111111111MIII1IIIr^ITII10 10^ 20^ 30^ 40^ 50TIME (hours)N" 200Chapter 4. Comparison of the MethodsINFLOW AND COMPUTED OUTFLOW HYDROGRAPHSin artificial prismatic channel300COEn 200f=X 100 -Co^-TLX = 80 kmB^= 20.0 mZ^= 1.0n^= 0.035So = 0.001^ InflowKinematic Wave Method...a. Muskingum-Conge Method Characteristic MethodDynamic Wave Methodeeeee UBC Flow Model0 r^r^i^i^i^i^r^r^I^r^r^r^i^r^i^r^r^r^I^r^i^i^r^r.^r^i^r^I0 10 20 30 40TIME (hours)Fig. 4.7. Comparison with various bed slopes4.7.a. : Manning's n =0.035 and So =0.00054.7.b. : Manning's n =0.035 and So =0.001050610.60(a)O• 0.4-aZ 0.21:4• 0.0C.)—0.2r.r.1V)fE1 —0.4Chapter 4. Comparison of the MethodsDISCREPANCIES FROM THE REFERENCE DISCHARGETLX = 80 kmB = 20.0 mZ^= 1.0n^= 0.035So = 0.0005te.^...no.......^ Kinematic Wave Methodeons ^ Muskingum—Cungs MethodCharacteristic Method4.4+44 UBC Flow ModelTT : time when the reference dischargecomes to peak flowI^i^i^i^i^i^i^l^i^i^i^i^I^I^ 11111111111110 20 30TIME (hours)..... ....^....0 °^.4.4.4++4......„...^,•^ a7,1104..°•0^ ...-•4 t st:•.°III '^t',,,,, . ............ ,I...22:wi....t.Tse40 500.30.20.0.1z4., 0.0F.1— 0.1al —0.2C/1Q —0.3(b)0. ^0...^....^.••• ^?:01:0 ^t• •:76444Stit • ...^; 1 WTLX = 80 kmB = 20.0 m2^= 1.0n^= 0.035So = 0.0005^ Kinematic Wave Method00000 Muskingum—Cunge MethodCharacteristic Method000e. UBC Flow ModelTrr : time when the reference dischargeTse^comes to peak flow500.40^ 10^20^30^40TIME (hours)Fig. 4.8. Comparison of discrepancies with So =0.0005 and n=0.0354.8.a. : divided by the reference flow4.8.b. : divided by the reference peak flow62Chapter 4. Comparison of the Methodsthose of the characteristic method. The same situation is also true for the caseof Fraser River flows as shown in fig. 4.3, that the characteristic wave methodgave closer results to the observed outflow than the kinematic wave did.Based on figures 4.5 - 4.8 and those in appendix B, the author would like togive some highlights as follow,For all the methods that were compared here, discrepancies to thereference discharge are bigger for a milder bed slope as well as for abigger value of Manning's n, and vice versa.The difference between the time to peak of the reference discharge andthose that are computed by the kinematic wave, UBC Flow andMuskingum-Cunge methods decreases as the Manning's n value becomessmaller or the bed slope becomes steeper.In general, the method of characteristics gave the closest results to thereference discharge, followed by the kinematic wave, UBC Flow andMuskingum-Cunge method. The absolute peak flow value which is closeto the reference peak flow was that from the UBC Flow model for achannel with steep bed slope or small value of Manning's n, and wasthat from the kinematic wave method for a channel with mild bed slopeor big value of Manning's n.Discrepancies of the characteristic method from the reference discharge63Chapter 4. Comparison of the Methodsachieved extreme values at the peak and toes of the hydrograph curve.For the kinematic wave method and the UBC Flow model, almost thesame values of discrepancies spread along the increasing and thedeclining slope of the hydrograph curves. Discrepancies with thekinematic wave method were usually less than those with the UBC Flowmodel.The maximum discrepancy from the reference flow at the same timestep of the characteristic method was 5.15%. Those values of thekinematic wave, UBC Flow, and Muskingum-Cunge methods were9.04%, 16.69%, and 25.05% respectively. All these maximumdiscrepancies occurred when the methods were applied to the mildestchannel's bed slope of this study (e.g., channel with Manning's n = 0.035and So=0.0005).- If those maximum discrepancies were compared to the reference peakflow, all those values mentioned above become smaller. Those values ofthe characteristic, the kinematic wave, UBC Flow model, andMuskingum-Cunge methods were 3.43%, 7.11%, 11.16%, and 13.36%.These maximum discrepancies were smaller for a steeper bed slope aswell as for a smaller value of Manning's n, and vice versa.Discrepancy indicator for all the methods for various bed slopes and64)-(-)a_CC 10U -(n-,o=ssAia Kinematic Wave MethodU13C Flow Model00080 Muskingum—Cunge Method***ww Characteristic MethodManning's n = 0.035z -O --< 20 -C_)^_E5^-z -Chapter 4. Comparison of the MethodsDISCREPANCIES FROM THE REFERENCE DISCHARGE30 ^z -O -e&o.so. Kinematic Wave MethodON* UBC Flow Model00800 Muskingum—Cunge Methodwww*-* Characteristic MethodSo = 0.0010< 20 -a0------^----------- 00^IIIIITIIIIIIIIIIIIIifIrf11111111111111111IIIIIIII0.010 0.020^0.030^0.040^0.050^0.060FRICTION FACTOR (Manning's n)300 tlrrlttttltt it t Itit irt tit^ ttitimtlit nlrr urtt t0.0^1.0^2.0^3.0^4.0^5.0^6.0BED SLOPE (X 1/1000)Fig. 4.9. Comparison of flood routing discrepancies4.9.a. : So =0.0010 with various Manning's n values4.9.b. : Manning's n = 0.035 with various bed slopes65Chapter 4. Comparison of the MethodsManning's n values are presented in figure 4.9.a and 4.9.b. The discrepancyindications from all methods decreased as the bed slope became steeper orManning's n value became smaller.From the results of this study and considering that the dynamic wave and thecharacteristic methods use the complete unsteady flow equations, it can beconcluded that the inertia and the depth slope terms of the equations influencedgreatly the results of the computations, especially for a channel with a mild bedslope or with a big value of friction coefficient (e.g., rough bed, covered by densevegetation). However, their effects were small for a steep or a smooth linedchannel. For example, the maximum discrepancy of the computed outflows ofthe kinematic wave for a channel with So =0.0005 and Manning's n =0.035 was7.11% compared with the reference peak flow. The discrepancy indication was9.370. On the other hand, those values for a channel with So = 0.0050 andManning's n =0.035 were 3.56% for discrepancy compared with the referencepeak flow and 4.557 for the discrepancy indication (see fig. 6.4.a2, 6.4.d2 and4.9.b).66Chapter 5SUMMARY AND CONCLUSIONSResults of five flood routing models have been compared to find the mostreliable model when they are applied to channels with various bed slopes and frictionfactor coefficients. The models represent five different methods. The first threemodels are the dynamic wave, the characteristic and the kinematic wave method.They are based on a hydraulic approach, however, they use different methods tosolve the non-linear partial different equations of unsteady flows. The fourth modelis Muskingum-Cunge method based on a hydrological approach, and the last one isa hybrid model that was developed by the University of British Columbia. Since theabove models use different ways to obtain their results, they produce different qualityof results.The dynamic wave method -- mathematically the most complicated one amongthem -- obtains the results by solving simultaneously the complete unsteady flowequations throughout the reach. The characteristic method simplifies the equationsthrough characteristic forms, while the kinematic wave method eliminates the inertiaand depth slope terms from the equations. The UBC Flow and Muskingum-Cungemethods use different equations that are based on the hydrological approach. Sincethe five models have different features, originally they need different inputparameters. They are also have different boundary conditions. With some67Chapter 5. Summary and Conclusionsmodifications, however, it is possible to compare them for the purpose of this study.All the above models were first applied to Fraser River flows, and then to aprismatic channel with various bed slopes and friction factor coefficients. The resultswere compared with each other and with the observed outflows. The computedoutflows of all the models had a good agreement to the Fraser River flows, but themost accurate result was that come from the dynamic wave method. In general, thecharacteristic, the kinematic wave, UBC Flow and Muskingum-Cunge methodsrespectively followed the dynamic wave method for their accuracy. The absolutepeak flow value which is close to the reference peak flow was that from the UBCFlow model for a channel with steep bed slope or small value of Manning's n, andwas that from the kinematic wave method for a channel with mild bed slope or bigManning's n value.In obtaining the results from various combinations of bed slopes and frictioncoefficients, the dynamic wave method is flexible, fast and easy to find convergence.The UBC Flow and the kinematic wave are also fast and flexible, but these methodsproduce less accurate results than the dynamic wave method does, especially for anearly horizontal channel. Results from the characteristic method are very close tothose of the dynamic wave method, except for the peak and toes of the hydrograph-curve. For a steep bed slope channel, however, the method requires a bigger valueof discharge or flow depth for initial condition than the other method do. Its68Chapter 5. Summar), and Conclusionscomputing time is also longer than those of the other method. The Muskingum-Cunge method produces the largest discrepancy among the compared methods.The maximum discrepancies of the characteristic method from the referencepeak flow for a nearly horizontal channel (e.g., S 0=0.0005 and Manning's n =0.035)is 3.43%. Those values for the kinematic wave, UBC Flow, and Muskingum-Cungemethods are 7.11%, 11.16% and 13.36% respectively. These maximum discrepanciesdecrease for a smaller value of the friction coefficient or a steeper bed slope channel.Although the UBC Flow model is not the most accurate one in flood routing, thismodel is simple and flexible, especially when the friction coefficient varies along withthe fluctuation of discharge in a channel. The input data can be obtained directlyfrom field measurements (e.g., velocity-discharge curve) that reduce error. Also, itis easy to adapt the method for solving problems in which the floods flow throughrivers and lakes (channel and reservoir routing).Results of the Muskingum-Cunge method are greatly influenced by its parameters(i.e., C1 , C2 and C3). Since the method is originally based on an inflow-outflowrelationship, the best way to obtain its parameters is from direct data (e.g., largehistorical inflow and outflow data). Unless using them, the method may giveinaccurate results. If there is not sufficient data available (e.g., remote area channel,man-made channel to improve the river flows), calculation should be done by usinga method which is based on the hydraulic approach (e.g., dynamic wave method).69Chapter 5. Summary and ConclusionsThe inertia and the depth slope terms of the unsteady flow equations influencednot only the way of solving the equations, but also the results. By neglecting them,to solve the equations becomes much easier. However, the results were greatlyaffected, especially for a channel with big value of friction coefficient (e.g., coveredby dense vegetation, rough bed) or with nearly horizontal bed slope. Their effectswere small for a smooth lined or steep slope channel.The flood routing models omitting the inertia and the depth slope terms aresimpler than those using the complete unsteady equations. However, the lattermodels gave more accurate results than the former ones. By using the completeunsteady flow equations, the dynamic wave model has a complex feature, but, thismodel is flexible, applicable, and reliable.It is recommended that further research be carried out on drainage networks inestuary areas affected by tides. These areas usually have greater probabilities ofreceiving floods that are caused not only by a large value of upstream discharge, butalso by the combination of the river flood and the high sea tide. However, locationsof villages and towns are often on low ground in estuary areas, so that the floodscause great damaging effects. The interaction between river flows and the tidesshould be considered. The additional stochastic nature of the river system as well asthe tides could be included in the simulations.70Bibliography1. Abbott, M.B. (1979). Computational Hydraulics : Elements of the Theory of FreeSurface Flows, Pitman Publsh.Ltd., London.2. Aldama, A.A. (1990). "Least-Square Parameter Estimation for Muskingum FloodRouting", Journal of Hydraulic Engineering, ASCE, 116(4), p.580-586.3. Amien, M. and Fang, C.S. (1970). "Implicit Flood Routing in Natural Channels",Journal of the Hydraulic Division, ASCE, 96(HY12), p.2481-2500.4. Amien, M., and Chu, H.L. (1975). "Implicit Numerical Modelling of UnsteadyFlows", Journal of the Hydraulic Division, ASCE, 101(HY6), p.717-'731.5. Anderson, D.A., Tannehill, J.C., and Pletcher, R.H. (1984). Computational FluidMechanics and Heat Transfer, Hemisphere Publsh.Co., New York.6. Ascher, U.M., Mattheij, R.M.M., and Russell, R.D. (1988). Numerical Solutionof Boundary Value Problems for ODEs, Prentice Hall, Englewood Cliffs, NewJersey.7. Atkinson, K.E. (1989). An Introduction to Numerical Analysis, 2nd Ed., JohnWiley & Sons Inc., New York.8. Ball, J.E. (1985). "An Algorithm for Routing Unsteady Flows in Urban DrainageNetworks", Journal of Hydraulic Research, 23(1), p.327-341.9. Blatzer.,R.A., Lai, C. (1968). "Computer Simulation of Unsteady Flows inWaterways", Journal of the Hydraulic Division, ASCE, 94(HY4), p.1083-1117.10. Boyd, M.J., and Bufill, M.C. (1989). "Determining Runoff Routing ModelParameters without Rainfall Data", Journal of Hydrology, 108, p.281-294.11. Carnahan, B., Luther, H.A., and Wilkes, J.O. (1969).Applied Numerical Methods,John Wiley & Sons Inc., New York.12. Chaudhry, M.H., and Contractor, D.N. (1973). "Application of Implicit Methodto Surge in Channels", Water Resources Research, 9(6), p.1605-1612.71Bibliography13. Chaudhry, M.H., and Schulte, A.M. (1986). "Computation of Steady-State,Gradually Varied Flow in Parallel Channel", Canadian Journal of CivilEngineering, 13, p.39-45.14. Chiu, C.L., and Isu, E.O. (1978). "Kalman filter in open channel flow estimation",Journal of the Hydraulic Division, ASCE, 104(HY8), p.1137-1152.15. Chow, V.T. (1959). Open Channel Hydraulics, McGrawHill-Kogakusha Ltd.,Tokyo.16. Chow, V.T. (1988). Hand Book of Applied Hydrology, McGraw Hill Co., NewYork, N.Y.17. Cooley, R.L., and MoM, S.A. (1976). "Finite Element Solution of Saint-Venantequations", Journal of the Hydraulic Division, ASCE, 102(HY6), p.759-775.18. Corps of Engineer, U.S.Army (1987). HEC-1: Flood Hydrograph Package UserManual, The Hydrologic Engineering Centre, Davis, C.A.19. Cunge, J.A., Holly Jr., F.M., and Verwey, A. (1980). Practical Aspects ofComputational River Hydraulics, Pitman Advanced Pblsh., Boston.20. Dawdy, D.R., Goldman,D., Merkel, W.H., Unkrich, C.L., Woolhiser, D.A.,Goodrich, D.C., Hromadka II, T.V., and De Vries, J.J. (1989). Discussion andClosure on:"Kinematic Wave Routing and Computational Error", Journal ofHydraulic Engineering, 115(2), p.278-289.21. Diskin, M.H. (1967). "On the Solution of the Muskingum Flood RoutingEquation", Journal of Hydrology, 5, p.286-289.22. Donnell, T.O., Pearson, C.P., and Woods, R.A. (1988). "Improved Fitting forThree-Parameter Muskingum Procedure", Journal of Hydraulic Engineering,114(5), p.516-528.23. Dooge, J.C.I. (1973). "Linier Theory of Hydrologic System",Agricultural ResearchService Technical Bulletin, no.1468.24. Dronkers, J.J. (1964). Tidal Computations in Rivers and Coastal Waters, NorthHoland Pblsh.Co., Amsterdam.72Bibliography25. Duthateau, P. and Zachman, D.W. (1986). Partial Defferential Equations,McGraw Hill Book Co., Toronto.26. Eagleson, P.S. (1970). Dynamic Hydrology, Mc.Graw Hill, New York.27. Field, W.G., and Williams, B.J. (1983). "A Generalised One-DimensionalKinematic Catchment Model", Journal of Hydrology, 60, p.25-45.28. Fleming, G. (1975). Computer Symulation Techniques in Hydrology, AmericanElsevier Publishing Co., New York.29. Fox, R.W., and McDonald, A.T. (1985). Introduction to Fluid Mechanics, 3rd Ed.,John Wiley & Sons Inc., Toronto.30. Fread, D.L. (1973). "Technique for Implicit Dynamic Routing in River withTributaries", Water Resources Research, 9(4), pp.919-925.31. Fread, D.L. (1980). "Capabilities of NSW Model to Forecast Flash FloodsCaused by Dam Failures", In: Proc. of 26th Annual Hydr.Div.Specialty Conf.,ASCE, College Park, Maryland, p.455-464.32. Garbrecht, J., and Brunner, G. (1991). "Hydrologic Channel-Flow Routing forCompound Sections", Journal of Hydraulic Engineering, ASCE, 117(5), p.629-642.33. Gill, M.A. (1977). "Routing Flood in River Channels", Nordic Hydrology, 8,p.163-170.34. Gill, M.A. (1978). "Flood Routing by the Muskingum Method", Journal ofHydrology, 36, p.353-363.35. Gill, M.A. (1979). "Critical Examination of the Muskingum Method", NordicHydrology, 10, p.261-270.36. Henderson, F.M. (1963). "Flood Waves in Prismatic Channels", Journal of theHydraulic Division, ASCE, 89(HY4), p.39-67.37. Henderson, F.M. (1966). Open Channel Flow, The MacMillan Co., New York.38. Hinwood, J.B., and Wallis, I.G. (1975). "Review of Models of Tidal Water",Journal of the Hydraulic Division, ASCE, 101(HY11), p.1405-1421.73Bibliography39. Hjelmfelt Jr., A.T. (1984). "Convolution and the Kinematic Wave Equations",Journal of Hydrology, 75, p.301-309.40. Hjelmfelt, Jr.A.T. (1985). "Negative Outflow from Muskingum Flood Routing",Journal of Hydraulic Engineering, 111(6), p.1010-1014.41. Hoggan, D.H. (1989). Computer-Assisted Floodplain Hydrology and Hydraulics,Featuring the US Army Corps of Engineers' Hec-1 and Hec-2 System Software,MacGraw Hill Co., New York.42. Hromadka II, T.V., and De Vries, J.J. (1988). "Kinematic Wave Routing andComputational Error", Journal of Hyraulic Engineering, 114(2), p.207-217.43. Hunt, B. (1987). "A Perturbation Solution of Flood Routing Problem", Journalof Hydraulic Research, 25(2), p.215-234.44. Jayawardena, A.W., White, J.K. (1977). "A Finite Element Distributed CatchmentModel, I : Analytical Basis", Journal of Hydrology, 34, p.269-286.45. Koussis, A. (1978). "Theoritical Estimation of Flood Routing Parameters",Journal of the Hydraulics Division, ASCE, 104(HY1), p.109-115.46. Koussis, A.D. (1980). "Comparison of Muskingum Method Difference Scheme",Journal of the Hydraulic Division, ASCE, 106(HY5), p.925-929.47. Kreyszig, E. (1983). Advanced Engineering Mathematics, 5th Ed., John Wiley &Sons Inc., Toronto.48. Kulandaiswany, V.C. (1966). "A Note on Muskingum Method of Flood Routing",Journal of Hydrology, 4, p.273-276.49. Liggett, J.E., and Woolhiser, D.A. (1967).^"Difference Solutions ofShallow-Water Equation", Journal of Mechanics Division, ASME, 93(EM3),p.39-71.50. Linsley, R.K., Kohler, M.A. and Paulhus, J.L. (1982). Hydrology for Engineers, 3ndEd., McGraw Hill Co., New York, N.Y.51. Lyn, D.A., and Goodwin, P. (1987), "Stability on General Preissmann Scheme",Journal of Hydraulic Engineering, ASCE, 113(HY1), p.16-28.74Bibliography52. Morris, E.M. (1979). "The Effect of Small Slope Aproximation and LowerBoundary Conditions on Solutions of the Saint Venant Equations", Journal ofHydrology, 40, p.31-47.53. Nash, J.E. (1959). "A note on the Muskingum Method of Flood Routing", Journalof Geophysical Res., 64, p.1053-1056.54. Parlange, J.Y., Rose, C.W., and Sander, G. (1981). "Kinematic FlowAppoximation of Runoff on a Plane : An Exact Analitycal Solution", Journal ofHydrology, 52, p.171-176.55. Ponce, V.H., and Simons, D.B. (1977). "Shallow Wave Propagation in OpenChannel Flow", Journal of the Hydraulic Div., ASCE, 103(HY12), p.1461-1476.56. Ponce, V.H., Li, R.M., and Simons, D.B. (1978a). "Applicability of Kinematic andDiffusion Models", Journal of the Hydraulic Div., ASCE, 104(HY3), p.353-360.57. Ponce, V.H., and Yevjevich, V. (1978b). "Muskingum-Cunge Method withVariable Parameters", Journal of the Hydraulic Div., ASCE, 104(HY12),p.1663-1667.58. Ponce, V.H., Indlekofer,H.,and Simons, D.B. (1979). "The Convergence ofImplicit Bed Transient Models", Journal of the Hydraulic Div., ASCE, 105(HY4),p.351-363.59. Price, R.K. (1974). "Comparison of Four Numerical Methods for Flood Routing",Journal of the Hydraulic Division, ASCE, 100(HY7), p.879-899.60. Overton, D.E. (1966). "Muskingum Flood Routing of Upland Stream Flow",Journal of Hydrology, 4, p.185-200.61. Quick, M.C., and Pipes, A. (1975). "Nonlinear Channel Routing by Computer",Journal of the Hydraulic Division, ASCE, 101(HY6), p.651-664.62. Quick, M.C., and Pipes, A. (1976). "A Combined Snowmelt and Rainfall RunoffModel", Canadian Journal of Civil Engineering, National Research CouncilCanada, 3(3), p.449-460.63. Roache, P.J. (1982).^Computational Fluid Dynamics, Hermosa Pblsh.,Albuquerque, USA.75Bibliography64. Samuels, P.G., and Skeels, C.P. (1990). "Stability Limits for Preisssman's Scheme",Journal of Hydraulic Engineering, ASCE, 116(8), p.997-1012.65. Schulte, A.M., and Chaudhry, M.H. (1987). "Gradually-Varied Flow in OpenChannel Networks", Journal of Hydraulic Research, 25(3), p.357-371.66. Singh, V.P., and McCann, R.C. (1980). "Some Note on Muskingum Method ofFlood Routing", Journal of Hydrology, 48, p.343-361.67. Singh, V.P., and Panagiotis, D.S. (1987). "Analysis of Non-Linear MuskingumFlood Routing", Journal of Hydraulic Engineering, 113(1), p.61-79.68. Stepien, Ireneusz (1987). "On the Numerical Solution of the Saint VenantEquations", Journal of Hydrology, 67, p.1-11.69. Strupczewski, W., and Kundzewicz, Z. (1980). "Translatory Characteristics of theMuskingum Method of Flood Routing", Journal of Hydrology, 48, p.363-368.70. Tung, Y.K. (1985). "River Flood Routing by Nonlinear Muskingum Method",Journal of Hydraulic Engineering, 111(12), p.1447-1460.71. Vallentine, H.R. (1969). Applied Hydrodynamics, SI ed., Butterworth & Co.Ltd.,London.72. Weinmann, P.E., and Laurenson, E.M., (1979). "Approximate Flood RoutingMethods : A Review", Journal of the Hydraulic Division, ASCE, 105(HY12),p.1521-1536.76Appendix A- Flowchart of the Dymanic Wave method- Flowchart of the Characteristic method- Flowchart of the Kinematic Wave method- Flowchart of Muskingum-Cunge method- Flowchart of UBC Flow model77Appendix AFlowchart of the Dynamic Wave Method :78la thereany stagnancy ofAY or All)yesEnhancethe Iterationyr^i-Ays;= LP; +Au;Print :`Change DrPrint: No Covergence`Appendix Ano0 ITMA)noilAY IN>YERR '\ yesorL4 I > UERR)Call FDW9flood wave speedCs DX suitable for noaccuracyYesPrint OUTPUT IChange etc age memoryfor the next time step^ > JMAX)noYesChange ee=e+ Incrementle>e„„,,)OryeaEnd79noDefine TQI and TLITO1 = sum of DTQITLI = TQItime t = JL = 0VC JH = 1i noConvert the Increment of DTrelated to DXCall CH8//INPUT:DX,DT, and channels parametersUpstream HydrographInitial ConditionsINT = 0InItal mode = 00Call CH1INT = 1Initial mode = YOCall CH2Appendix AFlowchart of the Characteristic Method :time = J +A tL = 1!STEP = 0 / C ISTEP = 0) yesnoSet for condition A :II =2IN = ND-1101 = 2ION= NCall CH3 (upstream)Call CH4Upstream dischargeIteration for upstreamSet for condition B :II =2IN = NDIGI = 2ION= N®l-080Appendix ACall CH5Calculation for Interior partof channel(STEP = 0)Call CH6Iteration for downstreamSet variable TC, YC, UCto convert time intervalat upstream & downstreamby Lagrang's 4-point interpolationISTEP = 0L = L+ 1!STEP = 1Change storage memoryfor the next time steprControl the last time stepL > LMAX orKIV,IGN) > TOG1MX-.1H+1yi yesCall CH7convert time intervalby Lagrang's interpolationPrint OUTPUT1C)81noi JHnor^time:t.)L = 0INT = 0Initial mode = 00Call KWI (Y0)INT = 1Initial mode = YO00 - I (YO)^iChannel's variableatt=lMOBOUP = OIN(L)OUP=OIN(L) +00OUP = OIN(L).10^/INPUT:DX,DT, and channel's parametersUpstream HydrographInitial ConditionConvert the Increment of DTrelated to DXCall KW3Appendix AFlowchart of the Kinematic Wave Method :®82Appendix A®Call KW1Call KW2CI and Y at time: J + Mnode: ICK at time: I + A tCalculate : a , gCs DX suitable foraccuracy ? Print :^—1.-'Change Dr^-----yes...< Next distance step >node = I + A xCalculate 0 and Yat node = I + A xnonode > NyesPrint OUTPUTChange storage memoryfor the next time stepC > JNIAXyes83Flowchart of Muskingum-Gunge Method :/INPUT :DX,DT, and channel's parametersUpstream HydrographInitial Conditions'i^H = 1 ) yesf noConvert the Increment of DTrelated to DXCall Mus3r^time :1 = JL = 0Ccp®time .1 +,6 t )L - 1, JMAXnL = 1t OUP = OIN(L)OUP=OIN(L)+00OUP = OIN(L)Cell Must0 and V at lime : J + Atnode: I84Appendix AnoC 00 ir 0i yseINT = 0 INT = 1Initial mode = 00 Initial mode = Yt)Cell Must (l0) CO = f (Y0)1Channel's variableelt=1CKOBAppendix AnoNext distance stepnode =I +AxCall MustCalculate :CK013,CIO/E,OBVEK, X, Cl, C2, C3Cs DX suitable foraccuracy ? Print :'Change DX'1, yesCalculate 0 and Yat node = I + A xCall Mustno Cnode > NyesPrint OUTplffChange storage memoryfor the next time step> MAX )isEnd85Appendix AFlowchart of UBC Flow Model :INPUT :Enter choice :Data, Run, or Quit I RunINPUT:Query 2 to 4 :Input data manuallyPrint file of dataUse data from fileINITIAL Filetotal reaches LRmaxname of constant filename of Inflow fileReach # LRLR 1 , LRmaxReadfileread Constant file # LAInflow :read Inflow file # LRPart 1 :calculate travel timecalculate translated flow QOTO86Print fileprint outflow 00as Inflow # (LR+1)------Appendix AT(time step J )J - 1 ,NPart2 :route 00T through a reservoircalculate outflow 00^C J > N )YesLR . LRmax--)yesC)nono87Appendix BComparison of Computed Outflow Hydrographsfor a channel with a constant Manning's nand various bed slopesComparison of Discrepancies from the Reference Dischargefor a channel with a constant Manning's nand various bed slopesComparison of Computed Outflow Hydrographsfor a channel with a constant bed slopeand various Manning's n valuesComparison of Discrepancies from the Reference Dischargefor a channel with a constant bed slope andvarious Manning's n valuesTLX = 80 km8 = 20.0 mZ = 1.0S. = 0.001n^= 0.015^ Inflow Kinematic Wave Method^ Muskingum-Cunge MethodCharacteristic Method•-•-v-v• Dynamic Wave Method.4.44 UBC Flow Model01^ 110^20 30^40TIME (hours)TLX = 80 kmB = 20.0 mZ = 1.0S. = 0.001n = 0.025- Inflow^ Kinematic Wave Method Muskingum-Cunge MethodCharacteristic Method•-•-v-v-tt Dynamic Wave Methodeti-ott• UBC Flow Model400 1^1^1^I^l^1^1^1^1^1^1^I^1^l^I^I^l^I^1^I^I^I^1^1^1^0 10 20 30TIME (hours)3003005)G 20005050200Appendix BINFLOW AND COMPUTED OUTFLOW HYDROGRAPHSin artificial prismatic channelFig. 6.1. Comparison with various Manning's n values6.1.a. : So =0.0010 and Manning's n =0.0156.1.b. : So =0.0010 and Manning's n =0.02589cn^-200 --•Appendix BINFLOW AND COMPUTED OUTFLOW HYDROGRAPHSin artificial prismatic channel300-200 -W0loo -UU)caTLX = 80 kmB = 20.0 mZ = 1.0= 0.001n^= 0.035- Inflow••••• Kinematic Wave Method^ Muskingum-Cunge MethodCharacteristic Method.-0•-d-n Dynamic Wave Method00000 UBC Flow Model0 r,[11111111111IT1111111111111[1111,1111f0 10^ 20^ 30^ 40^ 50TIME (hours)300TLX = 80 kmB = 20.0 mZ^= 1.0Se = 0.001n^= 0.045^ Inflow Kinematic Wave Methodon^ Muskingum-Congo Method---o Characteristic Methodo-o-.-o-• Dynamic Wave Methodeoeee UBC Flow Model0II^f I^I0 10^ 20^ 30^ 40^ 50TIME (hours)Fig. 6.1. Comparison with various Manning's n values6.1.c. : So =0.0010 and Manning's n =0.0356.1.d. : So =0.0010 and Manning's n =0.04590TLX = 80 kmB = 20.0 mZ = 1.0S. = 0.001n = 0.055^ Inflow••••• Kinematic Wave Method••••• Muskingum-Cunge MethodCharacteristic Method•-•-•-•-• Dynamic Wave Methodo•ooso UBC Flow Model50111111■1111111IIIII11111111III11111111111(1111110^10^20^30^40TIME (hours)300U1^-200 --U0Appendix BINFLOW AND COMPUTED OUTFLOW HYDROGRAPHSin artificial prismatic channelFig. 6.1. Comparison with various Manning's n values6.1.e. : So =0.0010 and Manning's n =0.05591,,A+° 6127,=="4-•eeP wag 11-6-seit ■-■-■'.1',131:1444;;1444444.4.44145/448 74.4 .opapp.010 'Appendix BDISCREPANCIES FROM THE REFERENCE DISCHARGE0 0.4Crl0.20.0z<4 —0.2C.C/D5 —0.4 #44t1Itsnrit"44." ...^......TLX = 80 kmB = 20.0 mZ = 1.0So = 0.001n^= 0.015^ Kinematic Wave Method...so Muskingum—Cunge MethodCharacteristic Method44441. UBC Flow ModelTw : time when the reference dischargecomes to peak flow0.60I^1^1^I^I^I^I^1^1i^I^I^I^I^I^I^I^11^I10 20 30TIME (hours)0.30.2Cm1a.0.1za. 0.0—0.1zra., —0.2Cs1IXC/)—0.340^50TLX = 80 km8 = 20.0 mZ = 1.0S. = 0.001n^= 0.015^ Kinematic Wave MethodMuskingum—Cunge Method--- Characteristic Method000eo UBC Flow ModelT,., : time when the reference dischargecomes to peak flow0.40I^I^t^I^1^t^I^I^1^I^1^1^I^1^I^I^I^I^1^1^1^I^i^1^I^t^I^1^I^I^I^I^I^I^I^1^1^i^1^1^I^I^I^111111 0 20 30 40 50TIME (hours)Fig. 6.2.a. Comparison of discrepancies with So =0.0010 and n=0.0156.2.a1.: divided by the reference flow6.2.a2. : divided by the reference peak flow(a2)92•50Appendix BDISCREPANCIES FROM THE REFERENCE DISCHARGEo 0.4 -.-aC...Fs1C.)^-Z 0.2 -fraX^/.::**tieo 44,__F.T.1F... 04* ....... :;.'::*oo............C=1^-^.11:w^° et ...... 4o,C4 0.0 - .4•444444.'..---#""flieg:"5'412''..--.evil'-- -''s`›.C.) **°.....°' -Z^-.•=4 -0.2 -0...^-rs1XC)Cil^--Kinematic Wave MethodA -0.4 - TLX = 80 km Muskingum-Cunge Method.....---- Characteristic MethodB^= 20.0 m^ 4.400 UBC Flow Model- Z^= 1.0- So = 0.001 Tit : time when the reference discharge- n^= 0.025^ 'Fp,^ comes to peak flow0.6^0^111,1111111^IIIIIIII'ITT111,111111111111111711111(bl) 10 20 30^ 40TIME (hours)O 0.3-a4 0.2 -GaC.^-C.) 0.1 -(246 .T.,W 0.0 - ^ :!:8*.44444;88444.oeeoeeeesitkoeirg:!z^-r:4<4a, -0.2 --TB" ==Z^= 1.08200 . 0km m ^ Kinematic Wave Method^ Muskingum-Cunge Method= Characteristic Method-0.3 eeeee UBC Flow ModelS. = 0.001^ Tw : time when the reference dischargen^= 0.025 Tr^ cornea to peak flow-0 . 4^i^T^T^T^T^T^i^T^T^1^i^r^i^T^T^i^T^i^T^l i ti^1111111 31IIIIIT111 4I1,11111I10 10 20 30^ 40^50(b2) TIME (hours)Fig. 6.2.b. Comparison of discrepancies with S0 =0.0010 and n =0.0256.2.bl. : divided by the reference flow6.2.b2. : divided by the reference peak flow93Appendix BDISCREPANCIES FROM THE REFERENCE DISCHARGEO• 0.4rr.3C.)Z 0.204col040.0›-C)z •444.3.8.... ::.totTLX = 80 kmB = 20.0 mZ^= 1.0n^= 0.035Sa = 0.0005111 11111111110^ Kinematic Wave Method...a ^ Muskingum—Gunge MethodCharacteristic Methodoe000 UBC Flow Model: time when the reference dischargecome. to peak flow1111111111111120^ 30^40TIME (hours).... . ..°^44+4.4,-. .... :::..°.. . .500.4(a2)^0Fig. 6.4.a. Comparison of discrepancies with So =0.0005 and n=0.0356.4.al. : divided by the reference flow6.4.a2. : divided by the reference peak flow99^visewww.svr•o...1444444._... ..^-Trr^ Kinematic Wave Method Muskingum-Cunge Method•-w.w^ Characteristic Method444^ UBC Flow ModelTrr time when the reference dischargecomes to peak flowTLX = 80 kmB = 20.0 mZ^= 1.0n^= 0.035Se = 0.001Appendix BDISCREPANCIES FROM THE REFERENCE DISCHARGE0.60..°°-.00.+"•44444g^.14 ........^.4400,,^A4 ^AAAAAA aRgaigh6:4..... .°4.■TLX = 80 kmB = 20.0 mZ^= 1.0n^= 0.035^Se = 0.001^Trr^ Kinematic Wave Methodea °° Muskingum-Cunge Method---•-•••■ Characteristic Methodoeeee UBC Flow ModelTw : time when the reference dischargecomes to peak flaw0.4Cs.1C.)0.2rs1Ls.Cs10.0zo -0.2l -0.4111111111111111111111111111111111111111111111111110 20^ 30^40^50TIME (hours)0.4 1^1111111111111111111^1111^1^1^1111^111^1^1^1^111111111^1^1^r . 50Fig. 6.4.b. Comparison of discrepancies with So =0.0010 and n=0.0356.4.bl. : divided by the reference flow6.4.b2. : divided by the reference peak flow0 10^ 20 30 40TIME (hours)(b1)(b2)10011414.441111121111111122tratia..04:..o ...... o ........20^ 30TIME (hours)40 50[ Tr.^1111111111111IIIIIIIIIII^ Kinematic Wave Method Muskingum—Cunge Method--- Characteristic Method0000e UBC Flow ModelTrr time when the reference dischargecomes to peak flowAppendix BDISCREPANCIES FROM THE REFERENCE DISCHARGE• 0.4rm.(-)Cr)• 0.2 —0404• 0.0>-Uz—0.2 -a..Cc)04to-0.4 -TmIIIIIIIII(1111111111111110^ 20^ 30^ 40^ 50TIME (hours)0.60TLX = 80 kmB = 20.0 mZ^= 1.0n^= 0.035So = 0.002^ Kinematic Wave Methodst.o ^ Muskingum—Cunge MethodCharacteristic Methodooeee UBC Flow Model: time when the reference dischargecomes to peak flow...44.44$2444ilitte,• 0.31=40.2Crl0'C..)^0.1z1=1• 0.0-0.1z0..• —0.2Cr)04(/)0-0.3TLX = 80 kmB^= 20.0 mZ^= 1.0n^= 0.035So = 0.0020.4(c2)^0 10Fig. 6.4.c. Comparison of discrepancies with So =0.0020 and n=0.0356.4.c1. : divided by the reference flow6.4.c2. : divided by the reference peak flow101..e4 TLX = 80 kmB = 20.0 mZ^= 1.0n^= 0.035So = 0.0050 0.4Cc.Z 0.2lx 0.0zU—0.2F1 —0.4414 ♦ ^ "'Ittgitews•••°4440444646$04$488=1:4:::•°.*+...,2444444004+4■44404044++1,64.*****t.-^ ° 4 r•^ef' AAAA$0.•••• 'e e•••112*4444044444444=U4=#$$344.40.Appendix BDISCREPANCIES FROM THE REFERENCE DISCHARGE(dl)^0.6 0 1111111111111111111111111111111111111111/ 11110^20^ 30^40^50TIME (hours)^ Kinematic Wave Methodoyee• UBC Flow ModelTw : time when the reference dischargecomes to peak flow• •O• 0.3-3F.vw0.2Cvl0..0.1z• o.o-0.1z—0.2Cc2txC.)Ci)• —0.3111111111111111111111111111111111111111111111111110^20^ 30^40TIME (hours)Fig. 6.4.d. Comparison of discrepancies with So =0.0050 and n=0.0356.4.d1. : divided by the reference flow6.4.d2. : divided by the reference peak flow0.4(d2)^0TLX = 80 kmB = 20.0 mZ^= 1.0n^= 0.0355„ = 0.005 T,,^ Kinematic Wave Method• 0000 U8C Flow ModelT,, : time when the reference dischargecomes to peak flow50102Appendix C- Jacobian matrix for the dynamic wave methodExample data for the UBC Flow model( area-discharge and velocity-discharge relationships )103Appendix C1. Jacobian matrix for the dynamic wave method0 E _ A ) +^(y1+3--y14-1) +^(1-0) (y11 -3715 )2^+1^2aF(u,y)^Axi^10(1_ A aT, .1+1—^e^+ „ayil+ 1^2 At 2^T2 0Y+ (1-0) [leo__ A aT,j+12^T2 ay) j ± 1aF(u,y) - 0[(-lS)^-1 0 (yil+? -34+1) +^(1-4)) (yil+1-yil)2 2aF(u,y) _ 1 AX1 4_ 0 ELT+ 10 (1_ „ n,r, _;+,.1-1 v.:. ) - - ( ili::1._ u.14-1 ) ]ay1:1-^2 A t^2^T2 ay- i+1p_ 0(1 _, A aT i+i+ (1-0) (u •^•T2 a,)^1-74. 1 -up]2 -1"^i+3.aauG(l+1^2 etiu,y)^•^ •^•Ax^0 [-LT+ —2 0 (L11.--u .1 4-1 )] + (1-0) [-2 0(uf7 - u:7 )]/-fi•S -1 + 1+ OgAxiuaG (u,y) ay1+1aG (u,y)aLd++11- 0g- [-1 + 2 Ax.S -1 + 1 ( aP/ay_ aA/a_y ) j+ 1 ]3^ A^i2 e ( uli++11 - up')^2+ ( 1 -0) [ 0 (^ul) ]exi + 0 [E +2 etS j+1+ OgAxi (u 1+1aG^2^.+1 ap/ay 44/ay r1+1 ]( u‘Y ) = eg[1 +^AxiS4 +1 (^p - 3 3 A^i+1ayl++3.1104Appendix Caz-, (u,y)aul+ 1aFi (u , y)ay1+1aFN(u,y) auji++11aFN(u,y) 8)4: 11- Al +1.up. ( aA ) j+1ay i12 t u \ j+ 1 , aA D ap,J+1_ ,—,--,3 PR i+i ay ay j+1105....-...0^ 1^I^1111111CU I^I^1^111111I^1^1^111111^1111111til 2.0 1^I^1^11111.......Appendix C(a)AREA—DISCHARGE CURVEfor UBC Flow Model0^200^ 400^ 600^800DISCHARGE (cu.m./sec.)VELOCITY—DISCHARGE CURVEfor UBC Flow Model3.0I^I^11111111^I^I^1^11111-^I^I 11111^1^11111111^11 11 1 111^I^111111111111111111111I 1 111111I 1 1111111 1111111 1111111111111111 I 1^11^1^11I I 1 11^1^[III 1 1111111,--. _ _L_ ____I __L. 4_1 4-4 i 4__^.._.^+__^_4^__4_ 4__ 1 _l_i 1^1^ 1^II I 1 t^1I I 1 iilI 1I 1^ 1111I 11 111,- - - 1- - 1 -1 -1 -1- 1-1 1 ^I 1^1111111111I 1111^111^111111I^1^[^t111 II^1I^t^t^1^1^:1i :^1^tI^11 1111I^It'^IIII 1- - 1-- 1- -1 -1-I^1111111^111111I^I^11111^111111->-i_ I^I^a^1^1^1^111E- 1^[^1^1^1[11E3 -^I^I^1^1^11111o.-^- 1^I^11111114] 1.0 --- - -- - 1- - -1 - -t -1 - 1 -1 -1 ti- - ->^-^I^1^111111111^I^11111^ 11111I^111111 1^I^111111111^1^1111111^I^I^1111111(b)1^I^1111111 I^I^I^IIIIIIIS0=10.001 I ^11I I^1111111 I1I^I^111111,n^.=-1, 0.035: I 11I11I I^I^111110.0^I^I^1^111111^I^I^11^11111 i^111111}1^I^61111^10 100DISCHARGE (cu.m./sec.)Fig. 6.7. Example data for the UBC Flow model(a) Area-discharge relationship(b) Velocity-Discharge relationship106